《龙贝格算法.doc》由会员分享，可在线阅读，更多相关《龙贝格算法.doc（7页珍藏版）》请在三一文库上搜索。
1、龙贝格算法一、问题分析1.1龙贝格积分题目要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y(x)=lsin(tx),则给定雨篷的长度后，求所需要平板材料的长度）。二、方法原理2.1龙贝格积分原理龙贝格算法是由递推算法得来的。由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。在变步长的过程中探讨梯形法的计算规律。设将求积区间a，b分为n个等分，则一共有n+1个等分点，n。这里用表示复化梯形法求得的积分值，其下标n表示等分数。先考察下一个字段，其中点，在该子段上二分前后两个积分值显然有下列关系将这一关系式关于k从0到n-1累加求和，即可导出下列递推公式需要强调指出的是，上式中的代
2、表二分前的步长，而梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计算量，自然式人们极为关心的。根据梯形法的误差公式，积分值的截断误差大致与成正比，因此步长减半后误差将减至四分之一，既有将上式移项整理，知由此可见，只要二分前后两个积分值和相当接近，就可以保证计算保证结果计算结果的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。按上式，积分值的误差大致等于，如果用这个误差值作为的一种补偿，可以期望，所得的应当是更好的结果。按上式，组合得到的近似值直接验证，用梯形二分前后的两个积分值和按式组合，结果得到辛普生法的积分值。再考察辛普生法。其截断误差与成正比。因此，若
3、将步长折半，则误差相应的减至十六分之一。既有由此得不难验证，上式右端的值其实就等于，就是说，用辛普生法二分前后的两个积分值和，在按上式再做线性组合，结果得到柯特斯法的积分值，既有重复同样的手续，依据斯科特法的误差公式可进一步导出龙贝格公式应当注意龙贝格公式已经不属于牛顿柯特斯公式的范畴。在步长二分的过程中运用公式加工三次，就能将粗糙的积分值逐步加工成精度较高的龙贝格，或者说，将收敛缓慢的梯形值序列加工成熟练迅速的龙贝格值序列，这种加速方法称龙贝格算法。三、算法设计3.1龙贝格积分算法就是求出,再走一遍求出，根据求出,再走一遍求出，根据求出，根据求出,再走一遍程序求出，根据得出，根据得出，再根据
4、得出,再走一边程序，得出，根据得出，根据得出，再由得出。再根据相减的绝对值小于其精度。那其中为求出的值。四、案例分析4.1龙贝格积分分析a积分下限b积分上限n区间个数e积分值要求达到的精度s用以存放除积分区间两端点以外的其他各节点函数值的累加和p积分区间两端点函数值之和h步长值T1 、T2分别存放二分区间前后梯形积分值S1 、S2分别存放二分区间前后辛普生积分值C1 、C2分别存放二分区间前后斯科特积分值R1 、R2分别存放二分区间前后龙贝格积分值五、总结5.1龙贝格积分总结通过本次试验，了解了龙贝格算法的计算过程，了解了龙贝格公式的计算收敛过程，用变步长的方法，逐步减小步长，反复积分，逐步得
5、到所求积分值满足精度要求。一步步从梯形法的递推到辛普森到柯特斯法，最后到龙贝格，让精度逐步升高。附录龙贝格积分：#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(*nn); printf(请输入积分的下限：);scanf(%f,&a);printf(n请输入积分的上限：);scan
6、f(%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 R2n); printf(%d %10.7f %10.7f %10.7f %10.7fn,i,T1,S1,C1,R1); do s=0; x=a+h/2; while(xb) s+=f(x); x+=h; T2=T1/2+s*h/2; S2=T2+(T2-T1)/3; if (k=1)
7、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.7fn,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.7fn,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;