1. 復化梯形法公式以及遞推化
復化梯形法是一種有效改善求積公式精度的方法。將[a,b]區間n等分，步長h = (b-a)/n，分點xk = a + kh。復化求積公式就是將這n等分的每一個小區間進行常規的梯形法求積，再將這n的小區間累加求和。 公式如下：
使用復化梯形法積分時，可以將此過程遞推化，以更方便的使用計算機實現。設積分區間[a,b]，將此區間n等分，則等分點共有n+1個，使用復化梯形積分求得Tn。進行二分，二分結果記為T2n，則有：
2. 龍貝格積分公式
龍貝格積分實際上是提高收斂速度的一種算法。由於復化梯形法步長減半后誤差減少至 ，即有：
整理得：
根據此思路，將收斂緩慢的梯形值序列Tn加工成收斂迅速的龍貝格值序列Rn，這就是龍貝格算法，加工算法流程如下：
實現：
1 #include<stdio.h> 2 #include<math.h> 3 #include<iostream> 4 #include<cstdio> 5 using namespace std; 6 int Rk=0; 7 int Tk=0; 8 double fx(double x) //被積函數 9 { 10 //if(x==0.0)return 1.0; 11 return 3*x*x*x+2*x*x+1 + sin(x); 12 } 13 double getReal(double a,double b){ 14 double r1 = 3.0/4.0 * b*b*b*b + 2.0/3.0*b*b*b + b - cos(b); 15 double r2 = 3.0/4.0 * a*a*a*a + 2.0/3.0*a*a*a + a - cos(a); 16 return r1 - r2; 17 } 18 double getS(double a,double b,double h) 19 { 20 double res=0.0; 21 for(double i=a+h/2.0; i<b; i+=h){ 22 res+=fx(i); 23 } 24 25 return res; 26 } 27 double Romberg(double a,double b,double e) 28 { 29 int k=1; 30 double T1,T2,S1,S2,C1,C2,R1,R2; 31 double h=b-a; 32 double s; 33 T1=(fx(a)+fx(b))*h/2.0; 34 int counter=0; 35 while(1) 36 { 37 Rk++; 38 counter++; 39 s=getS(a,b,h); 40 T2=(T1+h*s)/2.0; 41 S2=(4.0*T2-T1)/3.0; 42 h/=2.0; 43 T1=T2; 44 S1=S2; 45 C1=C2; 46 R1=R2; 47 if(k==1) 48 { 49 k++; 50 continue; 51 } 52 C2=(16.0*S2-S1)/15.0; 53 if(k==2) 54 { 55 k++; 56 continue; 57 } 58 R2=(64.0*C2-C1)/63.0; 59 if(k==3) 60 { 61 k++; 62 continue; 63 } 64 if(fabs(R1-R2)<e||counter>=100)break; 65 } 66 return R2; 67 } 68 double Tn(double a,double b,double e) 69 { 70 double T1,T2; 71 double h=b-a; 72 T1=(fx(a)+fx(b))*h/2.0; 73 while(1) 74 { 75 Tk++; 76 double s=getS(a,b,h); 77 T2=(T1+h*s)/2.0; 78 if(fabs(T2-T1)<e)break; 79 h/=2.0; 80 T1=T2; 81 } 82 return T2; 83 } 84 int main() 85 { 86 double a,b,e; 87 printf("輸入積分限和精度： a b e:"); 88 //輸入區間[a,b],和精度e 89 scanf("%lf%lf%lf",&a,&b,&e); 90 double t=Romberg(a,b,e); 91 //分別輸出龍貝格算法和梯形法的計算結果和相應二分次數 92 printf("\nRomberg:積分值：%.7lf -- 二分次數：%d\n",t,Rk); 93 t=Tn(a,b,e); 94 printf(" Tn:積分值：%.7lf -- 二分次數：%d\n",t,Tk); 95 double tf = getReal(a,b); 96 printf(" Real:%.7lf",tf); 97 return 0; 98 }