龙贝格求解积分

龙贝格公式就是逐次对分积分区间的方法，可以把前面计算的结果作为一个整体带入对分后的计算公式中，只需要增加新的分点的函数值。

以)(0k T表示],[b a 二分k 次后的梯形值，)(k mT 表示)(0k T 的m 次加速值，则有),2,1(,141144)(14)1(1)( =---=-+-k TTTk m k m mmk m龙贝格算法：步1 初始化：计算)]()([2)0(0b f a f ab T+-=，置1:=k （k记录区间],[b a 的二分次数），a b h -=；步2（二分） 计算积分值： ∑-=+---+=122/11)1(0)(01)(221k j j k k k x f h TT；步3（加速） 求加速值：计算)(1)1(1)(141144j k j jj k j jjj k jTTT--+------=，),,2,1(k j =；步4 精度检验：对指定的精度ε，若ε<--)0(1)0(k kTT，则终止计算，并取)0(kT作为所求的结果；否则置1:+=k k ，hh 21:=，转步2。

计算次序：(0)0T第1次循环 二分：(1)0T加速：(0)1T 第2次循环二分：(2)1T加速：(1)(0)12,T T第3次循环 二分：(3)2T加速：(2)(1)(0)123,,T T T第4次循环 二分：(4)3T加速：(3)(2)(1)(0)1234,,,T T TT……在命令窗口输入： a = 0; b = 1;epsilon = 5e-6;f = @(x)sin(x);%@(X)申请变量空间,计算sin 积分y = romberg(f,a,b,epsilon) %函数调用安回车，出结果。

一、实验目的及题目1. 实验目的：（1） 学习用复化辛普森公式及龙贝格方法求解积分并掌握这种方法。

（2）了解这些辛普森公式及龙贝格方法的概念，参考课本写出用复化辛普森算法以及龙贝格方法计算目标题目的程序，在matlab 中实现，并用matlab 内置的函数计算出结果，并提出存在的问题。

2. 题目：利用复化辛普森公式和龙贝格方法计算下列积分：（1）dx e x ⎰-5.002（2）dx x x ⎰202sin )2sin(cos π二、实验用仪器设备、器材或软件环境计算机、matlab 软件。

三、实验原理、程序框图、程序代码1.实验原理：根据微积分学基本定理，若被积函数f(x)在区间[a,b]上连续，只要能找到f(x)的一个原函数F(x)，便可利用牛顿-莱布尼茨公式求得积分值。

但会经常遇到如下问题：找不到用初等函数，找到了原函数，但因表达式过于复杂而不便计算等等。

此时则不能用牛顿-莱布尼茨公式，因此有必要研究如下公式。

1）复化求积公式及原理由于高阶插值的不稳定性，为了提高计算积分的精度，可把积分区间分为若干个小区间，将()I f 写成这些小区间上的积分之和，然后对每一个小区间上的积分应用到辛普森公式，或柯特斯公式，并把每个小区间上的结果累加，所得到的求积公式就称为复化求积公式。

辛普森公式的数值积分公式为：⎰+++-≈ba b f b a f a f a b dx x f )]()2(4)([6)(它的集合意义为用通过三点))(,()),2(,2()),(,(b f b b a f b a a f a =+的抛物线围城的曲边形面积来代替给定函数的积分。

同梯形公式一样，也有复化辛普森公式：)()(),()()]()(4)([6)(010121b f x f a f x f x f x x f h dx x f n n k k k k ba ==++≈∑⎰-=++ 其中 n ab h x x xk k k -=+=++,2121。

龙贝格(Romberg ）求积法1.算法理论Romberg 求积方法是以复化梯形公式为基础，应用Richardson 外推法导出的数值求积方法。

由复化梯形公式 )]()(2)([2222b f h a f a f h T +++=可以化为)]()]()([2[212112h a f h b f a f hT +++=＝)]([21211h a f h T ++一般地，把区间［a,b ］逐次分半k －1次，（k ＝1，2,……,n)区间长度（步长)为kk m a b h -=，其中mk ＝2k －1。

记k T ＝)1(k T由)1(k T ＝]))12(([21211)1(1∑=---++km j k k k h j a f h T 从而⎰badxx f )(＝)1(kT－)(''122k f h a b ξ- (1)按Richardson 外推思想，可将（1）看成关于k h ，误差为)(2k h O 的一个近似公式，因而，复化梯形公式的误差公式为⎰badxx f )(－)1(k T ＝......4221++kkh K h K ＝∑∞=12i i k i h K （2）取1+k h ＝k h 21有 ⎰ba dx x f )(－)1(1+k T ＝∑∞=+121221i ik ii hK （3）误差为)(2jh O 的误差公式 )(j kT＝)1(-j kT＋141)1(1)1(------j j k j k T T2。

误差及收敛性分析（1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差,然后将n 个子区间上的误差相加作为整个积分区间上的误差。

(2）收敛性,记h x i =∆，由于∑=++=ni i i n x f x f h f T 01))]()([2)(＝))()((21101∑∑-==∆+∆n i ni i i i i x x f x x f上面两个累加式都是积分和，由于)(x f 在区间],[b a 上可积可知，只要],[b a 的分划的最大子区间的长度0→λ时,也即∞→n 时,它们的极限都等于积分值)(f I 。

龙贝格积分1. 算法原理采用复化求积公式计算时，为使截断误差不超过ε，需要估计被积函数高阶导数的最大值，从而确定把积分区间[]b a ,分成等长子区间的个数n 。

首先在整个区间[]b a ,上应用梯形公式，算出积分近似值T1；然后将[]b a ,分半，对 应用复化梯形公式算出T2；再将每个小区间分半，一般地，每次总是在前一次的基础上再将小区间分半，然后利用递推公式进行计算，直至相邻两个值之差小于允许误差为止。

实际计算中，常用ε≤-n n T T 2作为判别计算终止的条件。

若满足，则取n T f I 2][≈；否则将区间再分半进行计算，知道满足精度要求为止。

又经过推导可知，∑=-++=ni i i n n x x f h T T 112)2(221，在实际计算中，取kn 2=，则k a b h 2-=，112)1*2(2++--+=+k i i ab i a x x 。

所以，上式可以写为∑=++--+-+=+kk i k k ab i a f a b T T 211122)2)12((2211k开始计算时，取())()(21b f a f ab T +-=龙贝格算法是由递推算法得来的。

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

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

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

数值分析实验报告（七）一.实验名称：利用龙贝格法求解积分二.实验目的：学会使用龙贝格法求解积分，体会其基本思想三．题目：用龙贝格法计算积分：311dyy⎰，要求()()500110k kT T---〈四．程序：#include<stdio.h>#include<math.h>double Romb(double a,double b,double eps) {double Rombf();int m,n,i,k;double y[10],h,ep,p,x,s,q;h=b-a;y[0]=h*(Rombf(a)+Rombf(b))/2.0;m=1;n=1;ep=eps+1.0;while((ep>=eps)&&(m<=9)){p=0.0;for(i=0;i<=n-1;i++){x=a+(i+0.5)*h;p=p+Rombf(x);}p=(y[0]+h*p)/2.0;s=1.0;for(k=1;k<=m;k++){s=4.0*s;q=(s*p-y[k-1])/(s-1.0);y[k-1]=p;p=q;}ep=fabs(q-y[m-1]);m=m+1;y[m-1]=q;n=n+n;h=h/2.0;}return q;}double Rombf(double x){double y;y=1/x;return y;}main(){double a=1;double b=3;double eps=1e-5;double t=Romb(a,b,eps);printf("The result is:%.3f\n",t);}五．运行结果：The result is:1.099六．体会：通过这次实验上机，使我更进一步了解了龙贝格法的计算思想，其在精度上很有保证，收敛较快，是解积分问题的有效方法。

信计050213 孟松见。