数值积分之龙贝格积分

除了复化求积外，这里用龙贝格积分法进行近似求积，其原理与埃特金插值有些类似，进行线性整合后使结果具有高精度的求积效果。在实际过程中，由于对于评判合理步长的困难，我们常采取变步长的办法进行计算，使结果满足精度要求。龙贝格运算过程：梯形公式装载数据→辛甫生化→柯特斯化→龙贝格高精度化。

 

对于梯形公式，我们知道Ｔ1=[Ｆ(Ｘｋ)+Ｆ(Ｘｋ+1)]*ｈ/2，Ｔ2=[Ｆ(Ｘｋ)+2*Ｆ(Ｘｋ+1/2)+Ｆ(Ｘｋ+1)]*ｈ/4，显然得出Ｔ2=Ｔ1/2+∑Ｆ(Ｘｋ+１/２)*ｈ/2的关系。其中的ｈ=(ｂ-ａ)/Ｎ，在下一步计算时，要进行步长二分。当取得一定的梯形公式数据点后，就可以直接进行数据高精度化了。

辛甫生化：Ｓｎ=(4*Ｔ2ｎ-Ｔｎ)/3；

柯特斯化：Ｃｎ=(１６*Ｓ2ｎ-Ｓｎ)/１５；

龙贝格化：Ｒｎ=(６４*Ｃ2ｎ-Ｃｎ)/６3。

 

建立数组Ｔ[]、Ｓ[]、Ｃ[]、Ｒ[]分别存放各自的数据，用函数Ｆ(ｘ)=Ｓｉｎ(ｘ)/ｘ在区间[0,1]上进行检验，Ｆ(0)=1。

h=b-a;T2=T1=(F(a)+F(b))*h/2;//进行梯形公式数据装载for(int i=0;i<N+3;i++)//N为龙贝格数据个数{T[i]=T2;sum=0;x=a+h/2;while(x<b){sum=F(x)+sum;x+=h;}T2=(T1+h*sum)/2;h=h/2;T1=T2;}//进行梯形公式到辛甫生的转化for(int i=0;i<N+2;i++)S[i]=(4*T[i+1]-T[i])/3;//进行辛甫生到柯特斯的转化for(int i=0;i<N+1;i++)C[i]=(16*S[i+1]-S[i])/15;//进行柯特斯到龙贝格的转化for(int i=0;i<N;i++)R[i]=(64*C[i+1]-C[i])/63;

最后，当然是验证输出结果了：

在VS2010下运行的结果，看来龙贝格公式果然高精度啊！