龙贝格积分实验报告

二、Romberg 积分法
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 表示成两部分之和，即
21
2221
122211
122211
1
[()()2()]21
[()()2()2((21))]21[()()2()]((21))]22N N
N N k N N N N N k k N N
N N N N k k T h f a f b f a kh h f a f b f a kh f a k h h f a f b f a kh h f a k h -=-==-===+++=+++++-=+++++-∑∑∑∑∑ 由此得到梯形值递推公式
2221
1
((21))12N
N
N N N k T T h f a k h ==++-∑
因此
1
11212122,[()()],211
,()
22h h b a T f a f b h h T T h f a h =-=
+==++
由复化梯形公式的截断误差有
2''
112''2222(),12(),12
N N N
N b a I T h f a b b a I T h f a b
ηηηη--=-≤≤--=-≤≤ 若''()f x 变化不大时，即''''12()()f f ηη≈，则有
22241
()
2413
N N N N N T T I T T T -≈
=+--
式（2）表明，用2N T 作为定积分I 的近似值，其误差大致为21
()3
N N T T -，
因此其终止条件为
2N N T T ε-≤
其中ε是预先给定的精度。

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

记
(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 序列，并满足如下递推关系：
(0)
(0)(0)1
211[()()],((21)),2222N N N k b a b a b a T
f a f b T T f a k N N
=---=+=++-∑ (1)(1)
24,1,2,41k k k k N N
N
k
T T T k ---==
-
以上递推公式就是Romberg 积分递推公式。

3.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）。

4.Romberg 积分法的应用
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 是最后积分值。

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

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

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