龙贝格算法

龙贝格算法 一、问题分析

1.1龙贝格积分题目

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

二、方法原理

2.1龙贝格积分原理

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

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

b -a

, k =0,1, ，n 。这里用T n 表示复化分，则一共有n+1个等分点x k =a +kh , h =n

梯形法求得的积分值，其下标n 表示等分数。

先考察下一个字段[x k , x k +1]，其中点x

k +1

2

=

1

(x k +x k +1)，在该子段上二分前后2

两个积分值

T 1=

h

f (x k )+f (x k +1)⎤⎡ ⎣⎦2

T 2=

h ⎡

⎢f (x k )+4⎢⎣⎤⎛⎫

f x 1⎪+f (x k +1)⎥

⎥⎝k +2⎭⎦

显然有下列关系

1h

T 2=T 1+

22

⎛⎫

f x 1⎪ ⎝k +2⎭

将这一关系式关于k 从0到n-1累加求和，即可导出下列递推公式

⎫1h n -1⎛

T 2n =T n +∑f x 1⎪

22k =0⎝k +2⎭

需要强调指出的是，上式中的h =

b -a

代表二分前的步长，而 n

12

x

k +

1⎫⎛

=a + k +⎪h

2⎭⎝

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

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

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

1-T 2n 1

≈ 1-T n 4

将上式移项整理，知

1

1-T 2n ≈(T 2n -T n )

3

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

法。

按上式，积分值T 2n 的误差大致等于(T 2n -T n ) ，如果用这个误差值作为T 2n 的一种

13

补偿，可以期望，所得的

T =T 2n +

141(T 2n -T n )=T 2n -T n 333

应当是更好的结果。

按上式，组合得到的近似值T 直接验证，用梯形二分前后的两个积分值T n 和T 2n 按

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

41S n =T 2n -T n

33

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

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

I -S 2n 1

≈ I -S n 16

由此得

I ≈

161S 2n -S n 1515

不难验证，上式右端的值其实就等于C n ，就是说，用辛普生法二分前后的两个积分值

S n 和S 2n ，在按上式再做线性组合，结果得到柯特斯法的积分值C n ，既有

161

S 2n -S n 1515

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

641R n ≈C 2n -C n

6363

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

C n ≈

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

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

三、算法设计

3.1龙贝格积分算法

就是求出T 1, 再走一遍求出T 2，根据T 1T 2求出S 1, 再走一遍求出T 4，根据T 2T 4求出S 2，

根据S 1S 2求出C 1, 再走一遍程序求出T 8，根据T 4T 8得出S 4，根据S 2S 4得出C 2，再根据得出T 16，根据T 8T 16得出S 8，根据S 4S 8得出C 4，再由C 2C 4C 1C 2得出R 1, 再走一边程序，

得出R 2。再根据R 1R 2相减的绝对值小于其精度。那其中R 2为求出的值。

四、案例分析

4.1龙贝格积分分析

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

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

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

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

五、总结

5.1龙贝格积分总结

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

附录

龙贝格积分：

#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",&a);

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

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

printf("请输入L ："); scanf("%f",&l);

printf("请输入T ："); scanf("%f",&t); 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",i,T1,S1,C1,R1); do {

s=0; x=a+h/2; while(x

s+=f(x); x+=h; }

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",i+1,T2,S2,C2,R2); 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",i,T2,S2,C2,R2); }while (1);

getchar();

return 0; }

float f(float x) {

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

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

}