实验题目 2 Romberg 积分法
摘要
考虑积分
欲求其近似值，可以采用如下公式：
（复化）梯形公式
（复化）辛卜生公式
（复化）柯特斯公式
这里，梯形公式显得算法简单，具有如下递推关系
因此，很容易实现从低阶的计算结果推算出高阶的近似值，而只需要花费较少的附加函数
计算。但是，由于梯形公式收敛阶较低，收敛速度缓慢。所以，如何提高收敛速度，自然
是人们极为关心的课题。为此，记 为将区间 进行 等份的复化梯形积分结果，
为将区间 进行 等份的复化辛卜生积分结果， 为将区间 进行 等份的
复化柯特斯积分结果。根据李查逊（Richardson）外推加速方法，可得到
可以证明，如果 充分光滑，则有
（ 固定）
这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
该方法的计算可按下表进行
…
…
…
… …
很明显，龙贝格计算过程在计算机上实现时，只需开辟一个一维数组，即每次计算的结果
，可存放在 位置上，其最终结果 是存放在 位置上。
前言
利用龙贝格积分法计算积分
程序设计流程：
1．准备初值，计算
且 （ 为等份次数）
2．按梯形公式的递推关系，计算
3．按龙贝格公式计算加速值
4．精度控制。对给定的精度 ，若
则终止计算，并取 作为所求结果；否则 ，重复 2~4 步，直到满足精度
为止。
问题 1
（1） 程序运行如下：
I = RombergInterg(inline('x.^2.*exp(x)'),0,1,20,1e-6,1)
0.7200 0 0 0 0
0.7184 0.7179 0 0 0
0.7183 0.7183 0.7183 0 0
0.7183 0.7183 0.7183 0.7183 0
0.7183 0.7183 0.7183 0.7183 0.7183
I = 0.7183
（2） 程序运行如下：
I = Romberginterg(inline('exp(x).*sin(x)'),1,3,25,1e-6,1)
10.9391 0 0 0 0 0
10.9495 10.9529 0 0 0 0
10.9500 10.9502 10.9500 0 0 0
10.9501 10.9502 10.9502 10.9502 0 0
10.9502 10.9502 10.9502 10.9502 10.9502 0
10.9502 10.9502 10.9502 10.9502 10.9502 10.950
I = 10.9502
（3） 程序运行如下：
I = Romberginterg(inline('4./(1+x.^2)'),0,1,25,1e-6,1)
3.1413 0 0 0 0
3.1416 3.1417 0 0 0
3.1416 3.1416 3.1416 0 0
3.1416 3.1416 3.1416 3.1416 0
3.1416 3.1416 3.1416 3.1416 3.1416
I = 3.1416
（4） 程序运行如下：
I = Romberginterg(inline('1./(1+x)'),0,1,25,1e-6,1)
0.6932 0 0 0 0
0.6932 0.6931 0 0 0
0.6931 0.6931 0.6931 0 0
0.6931 0.6931 0.6931 0.6931 0
0.6931 0.6931 0.6931 0.6931 0.6931
I = 0.6931
问题 2
（1） 程序运行如下：
I = Romberginterg(inline('cos(x)./sqrt(1-x)'),0,1,25,1e-6)
I = NaN
因为函数在 x=0 处出现了 0/0 的情况，极限为 1，所以 Matlab 的结果显示为 NaN 非数，
解决方法是把下限 0 改为一个小数 eps。
I = Romberginterg(inline('sin(x)./x'),eps,1,25,1e-6)
I = 0.9461
（2） 程序运行如下：
I = Romberginterg(inline('cos(x)./sqrt(1-x)'),0,1,25,1e-6)
I = NaN
与（1）的原因相同，函数在 x=1 处出现了 0/0 的情况，结果为无穷，此时应选择变换
t=
√
x
，将积分变为
I=2
∫
0
1
cos t
2
√
1−t
2
t dt
，再进行变换
t=sinu
，将积分变为
I=2
∫
0
π / 2
sinu∗cos (si n
2
u)du
，变换后输入命令：
I = RombergInterg(inline('2*sin(x).*cos((sin(x)).^2)'),0,pi/2,25,1e-6)
I = 1.4996
（3） 程序运行如下：
I = Romberginterg(inline('cos(x)./sqrt(x)'),0,1,25,1e-6)
I = NaN
函数在 x=0 处出现了 1/0 的情况，结果为无穷。先将积分变为
I=2
∫
0
1
cosxd
(
√
x
)
，再做
变换
t=
√
x
，
I=2
∫
0
1
cos t
2
dt
I = 2*Romberginterg(inline('cos(x.^2)'),0,1,25,1e-6)
I = 1.8090
（4） 程序运行如下：
I = Romberginterg(inline('x.*sin(x)./(1-x.^2)'),-1,1,25,1e-6)
I = NaN
函数在 x=1,-1 处出现了 sin1/0 的情况，结果为无穷。被积函数为偶函数，做变换
x=sint
，积分变为
I=2
∫
0
π / 2
sint∗sin
(
sint
)
dt
I = 2*Romberginterg(inline('sin(x).*sin(sin(x))'),0,pi/2,25,1e-6)
I = 1.3825
