龙贝格算法

求积分的龙贝格算法

计算f(x)=1/x在[1,3]上的积分：

#include <iostream>
#include <cstring>
#include <cmath>
using namespace std;

double f(double x){
    return 1.0/x;
}

int main(){
    double a,b;//积分区间的上界和下届 
    a=1;
    b=3;
    
    double t[100];
    memset(t,-1,sizeof(t));
    t[1]=(b-a)/2.0*( f(a)+f(b) ) ;
    for(int n=1;n<=8;n=2*n )
    {
        double temp=0;
        for(int i=0;i<=n-1;i++){
            temp+=f( a+(2*i+1)*(b-a)/(2*n)  );
        }
        temp=temp*(b-a)/(2*n);
        t[2*n]=t[n]/2.0 + temp;
    }
    
    for(int i=1;i<=16;i=i*2){
        printf("t(%d)=%2.6f
",i,t[i]);
    }
    printf("
");
    
    double s[100];
    memset(s,-1,sizeof(s));
    for(int n=1;n<=8;n*=2)
    {
        s[n]=4.0/3.0*t[2*n]-t[n]/3.0;
    }
    for(int i=1;i<=8;i=i*2){
        printf("s(%d)=%2.6f
",i,s[i]);
    }
    printf("
");
    
    double c[100];
    memset(c,-1,sizeof(c));
    for(int n=1;n<=4;n*=2)
    {
        c[n]=16.0/15.0*s[2*n]-s[n]/15.0;  
    } 
    for(int i=1;i<=4;i=i*2){
        printf("c(%d)=%2.6f
",i,c[i]);
    }
    printf("
");
    
    double r[100];
    memset(r,-1,sizeof(r));
    for(int n=1;n<=2;n*=2)
    {
        r[n]=64.0/63.0*c[2*n]-c[n]/63.0;
    }
    for(int i=1;i<=2;i=i*2){
        printf("r(%d)=%2.8f
",i,r[i]);
    }
    
    
}

输出结果为：