1、6.4 龙贝格求积公式 龙贝格积分(Romberg Integration)以逐次对分区间和复化求积为出发点，逐次修正近似值的方法构造出一种新的求积算法，称为龙贝格算法. 6.4.1 梯形法的递推化 在区间 , ba 上取步长knnabh 2, ，k可取 ,2,1,0 正整数.用复化梯形公式求积，可得 12012)()(2kkjjjxfxfhT . 将每个子区间对分，即 , ba 区间122kn 等分，112kabh ，每个子区间的新增分点为)(21jxf ，则 )()()()(212112021121jjjjjxfxfxfxfhTkk )(22)()(221120211201 kkjjjjj
2、xfhxfxfh 1202112120212)(22)(22kkkkjjkjjxfabTxfhT， 即对 ,2,1,0k ，求得梯形法的递推计算，称为变步长梯形公式 )()(202bfafabT ， 12021122)(221kkkjjkxfabTT . （7.18） 由递推公式（7.18）可知，在已经算出kT2的基础上再计算12kT 时，只要计算k2 个区间上新增分点的函数值.这与直接利用复化梯形公式求12kT 相比较，计算量几乎节省了一半. 例5 利用递推公式（7.18）重新计算积分10214dxx . 解 1） 首先对区间0，1使用梯形公式，得 3)24(21)1()0(2011 ffT
3、 . 2） 将区间二等分，新增分点 2/1x ，由递推公式（7.18）得 1.3)21(212112 fTT . 3） 再将各小区间二等分，新增分点 4/3,4/1x ，由递推公式（7.18）得 13117647.3)43()41(412124 ffTT . 4） 将区间八等分，新增分点 8/7,8/5,8/3,8/1x ，由递推公式（7.18）得 13898849.3)87()85()83()81(812148 ffffTT . 这样不断二分下去，计算结果见表7-4，其中k代表二分的次数，区间等分数kn 2 . 表7-4 梯形法递推计算值 k 0 1 2 3 4 5 6 nT 3 3.1 3
4、.13117647 3.13898849 3.14094161 3.14142989 3.14155196 k 7 8 9 nT 3.14158248, 3.14159011, 3.14159202 表7-4说明用复化梯形公式计算积分，将区间二分9次，即有分点513个，达到7位有效数字，计算量很大. 变步长梯形算法：由 )()(2)0( bfafhT 开始，由如下递推公式可生成一个梯形公式 )( jT 的序列： ,2,1,)(2)1()(112jxfhjTjTnkk，其中jabh 2/)( ，12jn ， khaxk . 变步长梯形公式的MATLAB程序：Rctrap.m function T
5、=Rctrap(f,a,b,n) %求积分的变步长梯形算法，其中， % f为被积函数； % a和b为积分区间的下限、上限； % n为复化区间个数； % 输出项T为变步长梯形公式计算值. m=1;h=b-a;T=zeros(1,n+1); T(1)=h*(feval(f,a)+feval(f,b)/2; for j=1:n m=2*m;h=h/2;s=0; for k=1:m/2 x=a+h*(2*k-1);s=s+feval(f,x); end T(j+1)=T(j)/2+h*s; end 对于例5，利用符号积分计算：exact=int( 4/(1+x2) ,0,1),exact=vpa(ex
6、act,10) 得到符号积分为exact = pi，10位有效数字的近似值为exact = 3.141592654. 利用变步长梯形求积函数 Rctrap 计算，首先编写被积函数f(x) f=inline( 4/(1+x2) ); 将积分区间对分9次，利用函数 Rctrap a=0;b=1;n=9;T=Rctrap(f,a,b,n);T=vpa(T,9) 计算得T = 3., 3.10000000, 3.13117647, 3.13898849, 3.14094161, 3.14142989, 3.14155196, 3.14158248, 3.14159011, 3.14159202. 6.
7、4.2 龙贝格求积公式 梯形法的递推计算简单但收敛慢，本节将讨论利用龙贝格算法提高收敛速度以节省计算量.龙贝格算法是在积分区间逐次分半的过程中，对用变步长梯形法产生的近似值进行加权平均，以获得准确程度较高的近似值的一种方法，具有公式简练、使用方便、结果较可靠等优点. 1、第一次外推计算（将梯形公式用余项修正为辛普森公式） 对于badxxfI )( ，由复化梯形公式的余项，得到 )(12222kkfhabTI ， ),(2bak ； )()2(1211222kkfhabTI ， ),(12bak . 当 )(xf 连续时， )()(122kkff ，因此 41221kkTITI. 移项整理可得
8、)(3122211 kkkTTTI . 【注】上式说明，用12kT 近似积分，其余项大约是对分区间前后两次结果之差的31倍. 用余项来修正梯形求积公式，则有 kkkkkTTTTTI222223134)(31111. 通过验证可知， )(231)4(31120212221kkkkjjxfhTTTI )(2)()(231120211120kkjjjjjxfhxfxfh )(4)(2)()(612021120kkjjjjxfxfbfafh， 这恰好是复化辛普森公式kS2，即 kkkTTS22231341. （7.19） 【注】 将12kT 用余项修正后得到辛普森公式kS2，加速了收敛，提高了计算精度
9、. 2、第二次外推计算（将12kS 用余项修正为柯特斯公式） 因为 ),(),()2(18022)4(42bafhabSIkkk ， ),(),()4(180)()2(180111122)4(42)4(412bafhabfhabSIkkkk， 所以161)/()(221 kkSISI ，因此 )(15122211 kkkSSSI ，从而 kkkkkSSSSSI222221511516)(151111 . （7.20） 可验证此式就是柯特斯公式kC2，则 kkkSSC22215115161. （7.21） 3、第三次外推计算（将柯特斯公式修正为龙贝格公式） 类似上述的加速公式，可得 kkkCCR
10、22263163641. （7.22） 【注】 （1） 在积分区间逐次分半的过程中利用公式（7.19）、（7.21）与（7.22），将粗糙的近似值kT2逐步地“加工”成越来越精确的近似值kS2、kC2、kR2，也就是说，将收敛速度缓慢的梯形序列逐步地“加工”成收敛速度越来越快的新序列.公式（7.19）、（7.21）与（7.22）结合就构成了龙贝格算法. （2） 若由龙贝格公式求出的积分精度不够时，可以将求积区间对分，求出12kR ，对给定的允许误差，当 kkRR221时，12kR 就是所求的解，否则再对分区间（k增加1），直至满足要求为止. 例6 用龙贝格算法计算积分10214dxx . 解
11、用龙贝格加速算法公式（7.19）、（7.20）与（7.21）加工例5得到的梯形值，计算结果见表7-5. 表7-5 龙贝格计算值 k kT2 12kS 22kC 32kR 0 1 2 3 4 3 3.1 3.131176 3.138988 3.140941 3.133333 3.141569 3.141593 3.141592 3.142118 3.141594 3.141592 3.141586 3.141593 由表7-5可知，这里利用二分4次的数据（它们的精度都很差，只有两三位有效数字），通过三次加速求得 141593.31R ，这个结果的每一位都是有效数字，加速的效果十分显著.而对于加速
12、算法的计算量，因只需做少量的四则运算，没有涉及到求函数值，可以忽略不计. 龙贝格算法：龙贝格积分 ),( kjR 保存在下三角矩阵中，第0列的元素 )0,( jR 用基于j2个 , ba 子区间的变步长梯形公式计算，然后利用龙贝格公式计算 ),( kjR . 第 j行的元素为： jkkjRkjRkjRkjRk 1,14)1,1()1,()1,(),( . 通过生成 kj 的龙贝格积分表 ),( kjR ，并以 )1,1( jjR 为最终解来逼近积分：),()( jjRdxxfba. 龙贝格公式的MATLAB程序：Romberg.m function R,quad,err,h=Romberg(f
13、,a,b,n,tol) %求积分的龙贝格算法，其中， % f为被积函数； % a和b为积分区间的下限、上限； % n为复化区间个数； % tol为容差； % 输出项R为龙贝格计算表； % quad为积分值； % err为误差估计； % h为最小步长. m=1; h=b-a; err=1; j=0; R=zeros(4,4); R(1,1)=h*(feval(f,a)+feval(f,b)/2; while (errtol)&(jn)|(j4) j=j+1; h=h/2; s=0; for k=1:m x=a+h*(2*k-1); s=s+feval(f,x); end R(j+1,1)=R(j
14、,1)/2+h*s; m=2*m; for k=1:j R(j+1,k+1)=R(j+1,k)+(R(j+1,k)-R(j,k)/(4k-1); end err=abs(R(j,j)-R(j+1,k+1); end quad=R(j+1,j+1); 对于例6，利用龙贝格函数Romberg计算 f=inline(4/(1+x2);a=0;b=1;n=4;tol=0.0001; R,quad,err,h=Romberg(f,a,b,n,tol) 计算得 R = 3.0000 0 0 0 0 3.1000 3.1333 0 0 0 3.1312 3.1416 3.1421 0 0 3.1390 3.1416 3.1416 3.1416 0 3.1409 3.1416 3.1416 3.1416 3.1416 quad = 3.1416 err = 6.8815e-006 h = 0.0625
