龙贝格积分实验报告

二、R o m b e r g 积分法
1.变步长Romberg 积分法的原理
复化求积方法对于提高精度是行之有效的方法，但复化公式的一个主要缺点在于要事先估计出部长。

若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。

而用复化公式的截断误差来估计步长，其结果是步长往往过小，而且''()f x 和(4)()f x 在区间[,]a b 上的上界M 的估计是较为困难的。

在实际计算中通常采用变步长的方法，即把步长逐次分半（也就是把步长二等分），直到达到某种精度为止，这种方法就是Romberg 积分法的思想。

在步长的逐步分半过程中，要解决两个问题：
1. 在计算出N T 后，如何计算2N T ,即导出2N T 和N T 之间的递推公式；
2. 在计算出N T 后，如何估计其误差，即算法的终止的准则是什么。

首先推导梯形值的递推公式，在计算N T 时，需要计算1N +个点处的函数值在计算出N T 后，在计算2N T 时，需将每个子区间再做二等分，共新增N 个节点。

为了避免重复计算，计算2N T 时，将已计算的1N +个点的数值保留下来，只计算新增N 个节点处的值。

为此，把2N T 表示成两部分之和，即
由此得到梯形值递推公式 因此
由复化梯形公式的截断误差有
若''()f x 变化不大时，即''''12()()f f ηη≈，则有
式（2）表明，用2N T 作为定积分I 的近似值，其误差大致为21
()3
N N T T -，因
此其终止条件为
其中ε是预先给定的精度。

积分公式
将上述方法不断推广下去，可以得到一个求积分的序列，而且这个序列很快收敛到所求的定积分。

记
(0)N N T T =，将区间N 等分的梯形值。

(1)N N T S =，将区间N 等分的Simpson (2)N N T C =，将区间N 等分的Cotes 。

(3)N N T R =，将区间N 等分的Romberg 。

由其可构造一个序列(){}k N T ，次序列称为Romberg 序列，并满足如下递推关系： 以上递推公式就是Romberg 积分递推公式。

积分程序
1. 置1N =，精度要求ε，1h b a =-；
2. 计算(0)1[()()]2
b a
T f a f b -=+； 3. 置22
N N
h h =，并计算(0)
(0)211((21))222N N N k b a b a T T f a k N N =--=++-∑； 4. 置,2,1;M N N N K ===
5. 计算(1)(1)
2441
k k k k M M
M
k T T T ---=-；
6. 若 1M =，则转（7）；否则置2
M
M =
，1k k =+转（5）； 7. 若()(1)11k k T T ε--≤，则停止计算（输出()1k T ），否则转（3）。

积分法的应用
function [T,n] = romb(f,a,b,eps) double R ;
if nargin<4,eps=1e-8;
end
h=b-a;R(1,1)=(h/2)*(feval(f,a)+feval(f,b));
n=1;J=0;err=1;
while (err>eps)
J=J+1;h=h/2;S=0;
for i=1:n
x=a+h*(2*i-1);
S=S+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*S;
for k=1:J
R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1);
end
err=abs(R(J+1,J+1)-R(J+1,J));
n=2*n;
end
R;
T=R(J+1,J+1)
End
其中输入项：f为被积函数，ab为积分区间的端点值，ep为积分精度；输出项：T 是逐次积分表值，n是迭代次数，R是最后积分值。

程序调用
可以将被积分函数编成函数文件，也可以直接使用内联函数来表示被积分函数，
示例
如下: >>f=inli
ne('1/(1+x.^2)','x'); >> [T,n,R]=romb(f,2,9,1e-9)
运行后得出其迭代次数，最终积分结果以及龙贝格积分矩阵如表2-1所示， 迭代次数N=64，最终的积分值R=.
表2-1 龙贝格积分矩阵
3.课本例题求解
1 当迭代精度ep=1e-9的条件下，迭代次数N=32，迭代结果R=
表2-2 式1对应的龙贝格积分矩阵
2 当迭代精度ep=1e-9的条件下，迭代次数N=32，迭代结果R=.
表2-3 式2对应的龙贝格积分矩阵
3对于积分
1
0ln(1)x
dx
x
+
⎰，由于积分下限0为其奇点，理论上无法进行数值积分，本题中近似取下限为1*10-9来进行计算。

当迭代精度ep=1e-9的条件下，迭代次数N=16，迭代结果R=.
表2-4 式3对应的龙贝格积分矩阵
4.对于积分2
0sin()x
dx x
π
⎰，同样积分下限0为积分函数的奇点，理论上无法进行数值积分运算，本题中仍取积分下限近似为1*10-9进行计算。

当迭代精度ep=1e-9的条件下，迭代次数N=16，迭代结果R=.
表2-5 式4对应的龙贝格积分矩阵。