龙贝格算法

龙贝格算法 一、问题分析

1.1龙贝格积分题目

要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y(x)=lsin(tx),则给定雨篷的长度后，求所需要平板材料的长度）。

二、方法原理

2.1龙贝格积分原理

龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。

在变步长的过程中探讨梯形法的计算规律。设将求积区间[a，b]分为n个等

ba,k0,1,，n。这里用Tn表示复化分，则一共有n+1个等分点xkakh,hn梯形法求得的积分值，其下标n表示等分数。

先考察下一个字段[xk,xk1]，其中点xk121xkxk1，在该子段上二分前后2两个积分值

T1hfxkfxk1 2hT2fxkfx1fxk1

4k2显然有下列关系

1hT2T122fx1 k2将这一关系式关于k从0到n-1累加求和，即可导出下列递推公式

1 / 7

真诚为您提供优质参考资料，若有不当之处，请指正。

1hn1T2nTnf22k0需要强调指出的是，上式中的hxk1 2ba代表二分前的步长，而 n12xk1akh

2梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计

算量，自然式人们极为关心的。

根据梯形法的误差公式，积分值Tn的截断误差大致与h2成正比，因此步长减半后误差将减至四分之一，既有

1T2n1 1Tn4将上式移项整理，知

11T2n(T2nTn)

3由此可见，只要二分前后两个积分值Tn和T2n相当接近，就可以保证计算保证结果计算结果T2n的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计

法。

按上式，积分值T2n的误差大致等于(T2nTn)，如果用这个误差值作为T2n的一种

13补偿，可以期望，所得的

TT2n应当是更好的结果。

141T2nTnT2nTn 333按上式，组合得到的近似值T直接验证，用梯形二分前后的两个积分值Tn和T2n按

式组合，结果得到辛普生法的积分值Sn。

41SnT2nTn

33

再考察辛普生法。其截断误差与h4成正比。因此，若将步长折半，则误差相

应的减至十六分之一。既有

IS2n1 ISn16由此得

2 / 7

真诚为您提供优质参考资料，若有不当之处，请指正。

I161S2nSn 1515不难验证，上式右端的值其实就等于Cn，就是说，用辛普生法二分前后的两个积分值

Sn和S2n，在按上式再做线性组合，结果得到柯特斯法的积分值Cn，既有

161S2nSn 1515重复同样的手续，依据斯科特法的误差公式可进一步导出龙贝格公式

1RnC2nCn

6363应当注意龙贝格公式已经不属于牛顿—柯特斯公式的范畴。

Cn

在步长二分的过程中运用公式加工三次，就能将粗糙的积分值Tn逐步加工成

精度较高的龙贝格Rn，或者说，将收敛缓慢的梯形值序列Tn加工成熟练迅速的龙贝格值序列Rn，这种加速方法称龙贝格算法。

三、算法设计

3.1龙贝格积分算法

就是求出T1,再走一遍求出T2，根据T1T2求出S1,再走一遍求出T4，根据T2T4求出S2，

根据S1S2求出C1,再走一遍程序求出T8，根据T4T8得出S4，根据S2S4得出C2，再根据

C1C2得出R1,再走一边程序，得出T16，根据T8T16得出S8，根据S4S8得出C4，再由C2C4得出R2。再根据R1R2相减的绝对值小于其精度。那其中R2为求出的值。

3 / 7

真诚为您提供优质参考资料，若有不当之处，请指正。

四、案例分析

4.1龙贝格积分分析

a—积分下限 b—积分上限 n—区间个数

e—积分值要求达到的精度

s—用以存放除积分区间两端点以外的其他各节点函数值的累加和 p—积分区间两端点函数值之和 h—步长值

T1 、T2分别存放二分区间前后梯形积分值 S1 、S2分别存放二分区间前后辛普生积分值 C1 、C2分别存放二分区间前后斯科特积分值 R1 、R2分别存放二分区间前后龙贝格积分值

五、总结

5.1龙贝格积分总结

通过本次试验，了解了龙贝格算法的计算过程，了解了龙贝格公式的计算收敛过程，用变步长的方法，逐步减小步长，反复积分，逐步得到所求积分值满足精度要求。一步步从梯形法的递推到辛普森到柯特斯法，最后到龙贝格，让精度逐步升高。

4 / 7

真诚为您提供优质参考资料，若有不当之处，请指正。

附录

龙贝格积分：

#include \"stdio.h\" #include \"math.h\"

float l; float t;

int main(void){ float f(float);

float a,b,e,h,T1=0,T2=0,S1=0,S2=0,C1=0,C2=0,R1=0,R2=0,k,s,x; int i=0;

printf(\"\\n****************************************\\n\"); printf(\"****************龙贝格算法**************\\n\"); printf(\"****************************************\\n\\n\"); printf(\"请输入积分的下限：\"); scanf(\"%f\

printf(\"\\n请输入积分的上限：\"); scanf(\"%f\

printf(\"\\n请输入允许误差：\"); scanf(\"%f\

printf(\"请输入L：\"); scanf(\"%f\

printf(\"请输入T：\"); scanf(\"%f\ k=1; h=b-a;

T1=h*(f(a)+f(b))/2;

printf(\"----------------------\\n\");

printf(\"k T2 S2 C2 R2\\n\"); printf(\"%d %10.7f %10.7f %10.7f %10.7f\\n\ do {

s=0; x=a+h/2; while(xs+=f(x); x+=h; }

5 / 7

真诚为您提供优质参考资料，若有不当之处，请指正。

T2=T1/2+s*h/2; S2=T2+(T2-T1)/3; if (k==1) { k=k+1; h=h/2; T1=T2; S1=S2; }

else if (k==2) {

C2=S2+(S2-S1)/15;C1=C2; k=k+1; h=h/2; T1=T2; S1=S2; }

else if (k==3) {

R2=C2+(C2-C1)/63;C2=S2+(S2-S1)/15;C1=C2;k=k+1;h=h/2;T1=T2;S1=S2; } else {

C2=S2+(S2-S1)/15; R2=C2+(C2-C1)/63; if (fabs(R2-R1)printf(\"%d %10.7f %10.7f %10.7f %10.7f\\n\ break; } else {

R1=R2; C1=C2; k=k+1; h=h/2; T1=T2; S1=S2; } }

i++;

printf(\"%d %10.7f %10.7f %10.7f %10.7f\\n\}while (1);

getchar();

return 0; }

float f(float x) {

6 / 7

真诚为您提供优质参考资料，若有不当之处，请指正。

float y=0; if(x==0.0) return 1;

y=(float)l*sin(t*x)/x; return y; }

7 / 7