数值积分之龙贝格积分

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

对于梯形公式，我们知道Ｔ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下运行的结果，看来龙贝格公式果然高精度啊！