龙贝格算法

Romberg算法算法思路1）梯形法的递推化为了提高求积精度，实际计算时若精度不够可以将步长逐次分半，以此求积分abf(x)的近似值。首先将[a,b]分为n等份，共有n+1个分点，注意到每个子区间[xk,xk+1]经过二分只增加了一个分点：xk+12=12(xk,+xk+1)。然后利用复合梯形求积公式可以推导出二分前后的积分值递推关系（h代表二分前的步长）：T2n=12Tn+h2k=0n-1f(xk+12)。（1.1）2）理查森外推加速从梯形公式出发，将区间[a,b]逐次二分可以提高求积公式精度，当[a,b]分为n等份时，若记Tn=Th，当区间[a,b]分为2n等份时，则有T2n=Th2。再有泰勒公式展开为：Th=I+α1h2+α2h4+…+αlh2l+…，然后用h/2代替h有Th2=I+α1h24+α2h416+…+αlh22l+…。再记Sh=4Th2-Th3，这将复合梯形公式的的误差阶O(h2)提高到了O(h4)，并且易知Sh=Sn，即将[a,b]分为n等份得到的复合辛普森公式。与上述做法类似，从Sh出发，当n在增加一倍，即h减少一半时得到Sh2，然后记Ch=16Sh2-Sh3，易知Ch=Cn，即将[a,b]分为n等份得到的复合柯特斯公式……如此继续下去就可以得到龙贝格公式。我们重新引入记号T0h=Th,T1h=Sh等，从而可以将上述公式写成统一的形式：Tmh=I+δ1h2(m+1)+δ2h2(m+2)+⋯.上述处理方法就称为理查森外推加速法。3）龙贝格求积算法设以T0(k)表示二分k次后求得的梯形值，且以Tm(k)表示序列T0(k)的m次加速值，则依理查森外推加速公式，可以得到：Tm(k)=4m4m-1Tm-1(k-1)-14m-1Tm-1k，k=1，2⋯.（1.2）这就是龙贝格求积算法。算法的程序设计（matlab实现）1）我们在matlab中定义函数function[T]=Romberg(f,a,b,e)。令T=zeros(10,10)，将其初始化为一个10阶的零元方阵，以此保存Tm(k)。f代表被积函数，a,b分别是积分限，e是我们要求的绝对误差；2）取k=1,h=b-a,求T1,1=T0(0)=h2fa+f(b).表示初始时的梯形值，k→10逐渐往后递推求梯形值；3）求梯形值T1b-a2(k-1)，即按照递推公式（1.1）计算Tk,1=T0k；4）求加速值，按公式（1.2）逐个求出T矩阵的第k行其余各元素Tk,j(j=2,3⋯k.)；5）若Tk,k-Tk-1,k-15%从第6行开始判断下面的精度要求，如果没有这个判断，计算可能发生错误if(abs(T(k,k)-T(k-1,k-1))=11)%超过矩阵大小则自动提示disp('溢出');end1）计算积分I=01x3/2dx注：利用函数f=inline(‘’)输入函数，再用Romberg（f,a,b,e）实现函数调用来计算积分值。2）计算积分I=02πxsinxdx四、Romberg公式讨论龙贝格公式是在梯形公式，辛普森公式以及柯特斯公式关系的基础上，构造出的一中加速计算积分的方法。相对于利用递推的复合梯形公式计算积分时，需要多次二分才能达到预定的精度。龙贝格积分利用外推加速可以极大的提高收敛的速度。也就是说在不增加计算量的前提下极大地提高了误差的精度。而且龙贝格公式在计算积分时，区间不断二分，这样前一次分半以后的函数值仍然能够保存下来，可以继续利用，这一点非常利于编程。当然对于f(x)不是特别光滑的函数也可以用龙贝格算法计算，不过收敛速度相对会慢一点。这种情况下直接使用复合辛普森公式也可以计算。比如说上述计算的积分I=01x3/2dx。