龙贝格积分
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
k
开始计算时，取
龙贝格算法是由递推算法得来的。由梯形公式得出辛普森公式得出柯特斯公式最后得到龙贝格公式。
根据梯形法的误差公式，积分值将减至四分之一，即有
T n 的截断误差大致与h 2成正比，因此步长减半后误差
1-T 2n 1
≈1-T n 4
将上式移项整理，知
1
1-T 2n ≈(T 2n -T n )
3
由此可见，只要二分前后两个积分值
T n 和T 2n 相当接近，就可以保证计算保证结果计算
结果
T 2n 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。
1
(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 。
41
S n =T 2n -T n
33
4
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
2
⎪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. 计算结果：