龙贝格算法

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

二、方法原理
2.1龙贝格积分原理
龙贝格算法是由递推算法得来的。

由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。

在变步长的过程中探讨梯形法的计算规律。

设将求积区间[a ，b]分为n 个等
分，则一共有n+1个等分点,k x a kh =+,0,1,b a
h k n
-==，n 。

这里用n T 表示复化梯形法求得的积分值，其下标n 表示等分数。

先考察下一个字段[1,k k x x +]，其中点()1
12
1
2
k k k x
x x ++=
+，在该子段上二分前后两个积分值
()()112
k k h
T f x f x +=
+⎡⎤⎣⎦ ()()21124k k k h T f x f x f x ++⎡
⎤⎛⎫=++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦
显然有下列关系
2112122
k h T T f x +⎛⎫
=+
⎪⎝⎭
将这一关系式关于k 从0到n-1累加求和，即可导出下列递推公式
1
210
2122n n n k k h T T f
x -+=⎛⎫
=+ ⎪⎝⎭
∑ 需要强调指出的是，上式中的b a
h n
-=
代表二分前的步长，而 12
12k x
a k h +
⎛
⎫=++ ⎪⎝
⎭
梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计
算量，自然式人们极为关心的。

根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，既有
211
14
n n T T -≈- 将上式移项整理，知
221
1()3
n n n T T T -≈-
由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计
法。

按上式，积分值2n T 的误差大致等于21()3
n n T T -，如果用这个误差值作为2n T 的一种
补偿，可以期望，所得的
()222141333
n n n n n T T T T T T =+
-=- 应当是更好的结果。

按上式，组合得到的近似值T 直接验证，用梯形二分前后的两个积分值n T 和2n T 按
式组合，结果得到辛普生法的积分值n S 。

24133
n n n S T T =-
再考察辛普生法。

其截断误差与4h 成正比。

因此，若将步长折半，则误差相
应的减至十六分之一。

既有
21
16
n n I S I S -≈- 由此得
21611515
n n I S S ≈
- 不难验证，上式右端的值其实就等于n C ，就是说，用辛普生法二分前后的两个积分值
n S 和2n S ，在按上式再做线性组合，结果得到柯特斯法的积分值n C ，既有
21611515
n n n C S S ≈
- 重复同样的手续，依据斯科特法的误差公式可进一步导出龙贝格公式
26416363
n n n R C C ≈-
应当注意龙贝格公式已经不属于牛顿—柯特斯公式的范畴。

在步长二分的过程中运用公式加工三次，就能将粗糙的积分值n T 逐步加工成
精度较高的龙贝格n R ，或者说，将收敛缓慢的梯形值序列n T 加工成熟练迅速的龙贝格值序列n R ，这种加速方法称龙贝格算法。

三、算法设计
3.1龙贝格积分算法
就是求出1T ,再走一遍求出2T ，根据12T T 求出1S ,再走一遍求出4T ，根据2T 4T 求出2S ，
根据1S 2S 求出1C ,再走一遍程序求出8T ，根据4T 8T 得出4S ，根据2S 4S 得出2C ，再根据
1C 2C 得出1R ,再走一边程序，得出16T ，根据8T 16T 得出8S ，根据4S 8S 得出4C ，再由2C 4
C 得出2R 。

再根据1R 2R 相减的绝对值小于其精度。

那其中2R 为求出的值。

四、案例分析
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<b)
{
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)<e)
{
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;
}。