龙贝格算法

龙贝格积分

1. 算法原理

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

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

实际计算中，常用

T 2n -T n ≤ε

作为判别计算终止的条件。若满足，则取

I [f ]≈T 2n ；

否则将区间再分半进行计算，知道满足精度要求为止。

x +x 1h n

T 2n =T n +∑f (i -1i ) k

222n =2i =1又经过推导可知，，在实际计算中，取，则h =

b -a x i +x i +1b -a

=a +(2*i -1)

2k ，22k +1。所以，上式可以写为

T 2k +1

1b -a 2b -a

=T 2k +k +1∑f (a +(2i -1) k +1) 222i =1

T 1=

b -a

(f (a ) +f (b ) )2

开始计算时，取

龙贝格算法是由递推算法得来的。由梯形公式得出辛普森公式得出柯特斯公式最后得到龙贝格公式。

根据梯形法的误差公式，积分值将减至四分之一，即有

T n 的截断误差大致与h 2成正比，因此步长减半后误差

1-T 2n 1

≈1-T n 4

将上式移项整理，知

1-T 2n ≈(T 2n -T n )

由此可见，只要二分前后两个积分值

T n 和T 2n 相当接近，就可以保证计算保证结果计算

结果

T 2n 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。

(T 2n -T n )

T T 32n 按上式，积分值的误差大致等于，如果用这个误差值作为2n 的一

种补偿，可以期望，所得的

T =T 2n +

141

(T 2n -T n )=T 2n -T n 333

应当是更好的结果。

按上式，组合得到的近似值T 直接验证，用梯形二分前后的两个积分值

T n 和T 2n 按

式组合，结果得到辛普生法的积分值

S n 。

S n =T 2n -T n

h 再考察辛普森法。其截断误差与成正比。因此，若将步长折半，则误差相应的

减至十六分之一。即有

I -S 2n 1

≈I -S n 16

由此得

I ≈

161S 2n -S n 1515

不难验证，上式右端的值其实就等于和

C n ，S 就是说，用辛普生法二分前后的两个积分值n

S 2n ，在按上式再做线性组合，结果得到柯特斯法的积分值C n ，即有

C n ≈

161

S 2n -S n 1515 641C 2n -C n 6363

重复同样的手续，依据斯科特法的误差公式可进一步导出龙贝格公式

R n ≈

龙贝格积分法的计算公式如下：

b -a ⎧T =(f (a )+f (b )) ⎪1

⎪k

1b -a 2b -a ⎪T T f (a +(2i -1) ) , k =0, 1, 2, ⋯⋯k +1=k +∑⎪2

222k +1i =12k +1

⎪

1⎪

S =T +(T 2k +1-T 2k ) ，k =0, 1, 2，⋯⋯k k +1⎨22

4-1⎪

1⎪

C (S k +1-S 2k ) ，k =0, 1, 2, ⋯⋯k =S k +1+⎪22

42-12

⎪

⎪R =C +1(C -C )，k =0, 1, 2, ⋯⋯k +1 ⎪2k 2k +12k

43-12

⎩

2. 程序框图

3. 源代码

%龙贝格积分公式

function [T,n]=romberg(f,a,b,fa,fb,eps)

%f为被积函数，a 和b 为区间左右端点，fa 为函数在左端点的值，fb 为函数在右端点的值 %eps为精度要求

%返回近似积分结果T 和区间等分数n h=b-a;

R(1,1)=(h/2)*(fa+fb); n=1;err=1;J=0;k=0; while (err>eps) k=k+1;J=J+1;s=0; for i=1:n

s=s+feval(f,a+h*(2*i-1)/(2^k)); end

R(J+1,1)=R(J,1)/2+s*(h/(2^k)); for m=1:J

R(J+1,m+1)=R(J+1,m)+(R(J+1,m)-R(J,m))/(4^m-1); end if J

err=abs(R(J+1,1)-R(J,1)); else

err=abs(R(J+1,4)-R(J,4)); end n=2*n; end R

T=R(J+1,4); end

%复化辛普森求积

function s=simp(f,a,b,fa,fb,n)

%f为被积函数，a 、b 是积分区间的左右端点，fa 为函数在左端点的取值，fb 为函数在右端点的取值，n 为区间的等分数，s 为返回的积分近似值 h=(b-a)/n;

x=linspace(a,b,2*n+1); for i=2:2*n y=feval(f,x); end

s=(h/6)*(fa+fb+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))); end

%复化梯形求积

function s=trap(f,a,b,fa,fb,n)

%f为被积函数，a 、b 是积分区间的左右端点，fa 为函数在左端点的取值，fb 为函数在右端点的取值，n 为区间的等分数，s 为返回的积分近似值 h=(b-a)/n;

x=linspace(a,b,n+1); for i=2:n y=feval(f,x); end

s=0.5*h*(fa+fb+2*sum(y(2:n))); end

%习题6.2 f1=@(x)1./(1+x);

[T1,n1]=romberg(f1,0,1,1,0.5,1.e-7) T1=vpa(T1,7)%有效数字为七位 f2=@(x)log(1+x)./(1+x^2);

[T2,n2]=romberg(f2,0,1,0,log(2)/2,1.e-7) T2=vpa(T2,7) f3=@(x)log(1+x)./x;

[T3,n3]=romberg(f3,0,1,0,log(2),1.e-7) T3=vpa(T3,7) f4=@(x)sin(x)./x;

[T4,n4]=romberg(f4,0,1,1,2/pi,1.e-7) T4=vpa(T4,7) f1=@(x)1./(1+x);

s11=vpa(trap(f1,0,1,1,0.5,8),7) s12=vpa(simp(f1,0,1,1,0.5,8),7) f2=@(x)log(1+x)./(1+x.^2);

s21=vpa(trap(f2,0,1,0,log(2)/2,8),7) s22=vpa(simp(f2,0,1,0,log(2)/2,8),7) f3=@(x)log(1+x)./x;

s31=vpa(trap(f3,0,1,0,log(2),8),7) s32=vpa(simp(f3,0,1,0,log(2),8),7) f4=@(x)sin(x)./x;

s41=vpa(trap(f4,0,1,1,2/pi,8),7) s42=vpa(simp(f4,0,1,1,2/pi,8),7)

4. 计算结果：