龙贝格算法

求积分的龙贝格算法

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

 1 #include <iostream>
 2 #include <cstring>
 3 #include <cmath>
 4 using namespace std;
 5
 6 double f(double x){
 7     return 1.0/x;
 8 }
 9
10 int main(){
11     double a,b;//积分区间的上界和下届
12     a=1;
13     b=3;
14
15     double t[100];
16     memset(t,-1,sizeof(t));
17     t[1]=(b-a)/2.0*( f(a)+f(b) ) ;
18     for(int n=1;n<=8;n=2*n )
19     {
20         double temp=0;
21         for(int i=0;i<=n-1;i++){
22             temp+=f( a+(2*i+1)*(b-a)/(2*n)  );
23         }
24         temp=temp*(b-a)/(2*n);
25         t[2*n]=t[n]/2.0 + temp;
26     }
27
28     for(int i=1;i<=16;i=i*2){
29         printf("t(%d)=%2.6f\n",i,t[i]);
30     }
31     printf("\n");
32
33     double s[100];
34     memset(s,-1,sizeof(s));
35     for(int n=1;n<=8;n*=2)
36     {
37         s[n]=4.0/3.0*t[2*n]-t[n]/3.0;
38     }
39     for(int i=1;i<=8;i=i*2){
40         printf("s(%d)=%2.6f\n",i,s[i]);
41     }
42     printf("\n");
43
44     double c[100];
45     memset(c,-1,sizeof(c));
46     for(int n=1;n<=4;n*=2)
47     {
48         c[n]=16.0/15.0*s[2*n]-s[n]/15.0;
49     }
50     for(int i=1;i<=4;i=i*2){
51         printf("c(%d)=%2.6f\n",i,c[i]);
52     }
53     printf("\n");
54
55     double r[100];
56     memset(r,-1,sizeof(r));
57     for(int n=1;n<=2;n*=2)
58     {
59         r[n]=64.0/63.0*c[2*n]-c[n]/63.0;
60     }
61     for(int i=1;i<=2;i=i*2){
62         printf("r(%d)=%2.8f\n",i,r[i]);
63     }
64
65
66 }

输出结果为：