龙贝格求积法

T2
4

1
1 3
T1

3.133333
C1

16 15
S2
1 15
S1

3.142118
S2 3 T4 3 T2 3.141569
T8

1 2 T4

1 8

f
(
1 8
)

f
(
3 8
)

f
(
5 8
)

f
(
7 8
)

3.138989
S4

4 3
T8

1 3
T4

3.141593
C2
， 类似推导可得：
Rn

64 63
C2n

1 63
Cn

43C2n Cn 43 1
（4.4.8）
即在Cotes序列 C2k 的基础上，产生了一个称为Romberg 序列的新序列 R2k : R1 、R2 、R4 ...
我们在变步长的过程中运用公式(4.4.4-4.4.8)(也称它们为
(由上节梯形公式)

h 6

f
n1
a 2
k 1
f
xk
n1
4
k 0
f


x
k

1 2

f
b

Sn
（由上节Simpson公式）
故（4.4.6）式成立.
同理，由上节近似式
I

S2n

1 15
S
2n
Sn
类似推导可得：
C2n
C k2 2
R k3 2
O O.9207355 1 O.9397933 2 O.9445135 3 0.9456909
O.9461459 O.9460869 O.9460833
O.9460830 O.9460831
O.9460831
我们看到，这里利用二分3次的数据(它们的精度都很差，只有二三位
是有效数字)，通过三次加速求得 R1 =0.9460831，这个结果的每一位
根据梯形法的误差公式，可知积分值 Tn的截断误差大致与 h2
成正比，因此当步长二分后，截断误差将减至原有误差的1/4，
即有 I T2n 1 将上式移项整理，可得
I Tn 4
I
T2n

1 3
(T2
n
Tn )
（4.4.3）
由此可见，只要二分前后的两个积分值Tn 与T2n 相当接近，就可以保
f
(
1 2
)

1 3 2
1 16 25

3.1
S1

4 3
T2

1 3
T1

3.133333
T4

1 2
T2

1 4

f
(
1 4
)

f
(
3 4
)

1 3.1 2
1 (3.764706 2.56) 4

3.131177
S2

4 3 T4

1 3 T2

3.141569
S1

4 3
数字都是有效数字，可见加速的效果是十分显著的。三次外推后达 到6位有效数字。
三、Romberg 算法计算公式的简化
注意 教材中介绍的Richardson外推法，为便于上机计算，引用记
号 Tmk 来表示各近似值，其中k仍代表积分区间的二分次数，而
下标m则指出了近似值 所在序列的性质。如m=1在梯形序列中，
定理4.5 假设F(h)逼近F(0)的余项为
F (h) F (0) a1h p1 a2h p2 a3h p3 L ,
其中，p1 p2 p3 L , ak (k 0,1, 2,L )是与h无关的非零常数。
则由
F0 (h)

F (h),
Fk 1 (h)

q pk
Fk
(
S8
C4
R2
5 25=32
T32
S16
C8
R4
…
…
…
…
…
例1: 利用Romberg算法计算

14 0 1 x2 dx
计算到R1.
解:由题意
a

0, b

1,
f
(x)

4 1 x2
T1

1f
2
(0)
f
(1)
1 (4 2
2)

3
T1,T2 ,T4 ,T8.
T2

1 2
T1

1 2
n 1, 2, 4,8,...
的递推化公式求得
龙贝格求积算法也可用下表来表示：
区间等分数 梯形序列 辛普森序列 柯特斯序列 龙贝格序列
k
2k
T2k
S 2 k 1
C 2k2
R2 k 3
0 20=1
T1
1 21=2
T2
S1
2 22=4
T4
S2
C1
3 23=8
T8
S4
C2
R1
4 24=16
T16
1, 2,...
T (k 1) m1

4mTm(k ) Tm(k 1) 4m 1
,m
1, 2,..., k
1, 2,...,i
由此可逐行构造 出一个三角形数表---称为T数表
外推次数
分
k
半i
T1k
次 数
0 T1 = T10 1 T2 = T11
2 3
T4 = T12 T8 = T13
Romberg序列为止。
2、Tmk 可用二维数组来存放并参加运算，也可用一维数组。
3、 对于积分限为无穷的积分 ，可利用变量代换化成有限区
间的积分然后再进行计算。例如： dx
1 1 x2 x3
令
x1 t
则
dx


1 t2
dt
代入得：
dx
1 1 x2 x3
0
(42.4,.33-4,.4.8) 可以将粗糙的近似值 逐步地“加工”成越来
越精确的近似值 Sn ,Cn , Rn , 也就是将收敛速度缓慢
的梯形序列 T2k 逐步地“加工”成收敛速度越来越快的新
新序列 S2k , C2k , R2k ,
这种加速的方法称为 Romberg算法。其“加工”过程如下图， 其中圆圈中号码表示计算顺序。
值小于事先给定的精度要求。取最后一次的数值积分值作为积分的
近似值。
四、几点说明：
4m
1
1、在上面“加工”过程中的系数4m 1 和4m 1 ，当
m
4
时， 1
4m 1

1 255
而另一个 系数则接近于1，也就是
新公式与原公式差别不大，Tmk Tmk1 但工作量却大增。
因此，在实际计算中常规定 m 3，即计算到出现
可将上述外推思想推广到一般情况。设F(h)是计算F(0)的一 种近似算式，带截断误差的表达式为
F(h) F(0) aphp O(hs ), s p,
其中，a p 与h无关。如果我们用h和h/q (q>1)两种步长分别计算
F(h)和F(h /q)，则有
F
(
h q
)

F
(0)

ap
(
h q

dx

1
t 2 1
.
可能是wk.baidu.com好的结果。
（4.4.4）
1 Tn T2n 3 (T2n Tn )
这就是说，用梯形法二分前后两个积分值 Tn 和T2n作线性组合
即
41 Tn 3 T2n 3 Tn
有可能比
T2n
更好地接近于积分
b
a
f
xdx
的真值
I
(4.4.5)
这种新近似值 Tn实质上又是什么呢？可以验证：

16 15
S4
1 15
S2

3.1415946
64 1 R1 63 C2 63 C1 3.141586292
1 sin x
算例结2 果用见R表om4-b5(ekr代g算表法二计分算次I数)。0 计x算值dx.的得误到差的不梯超形过值，计
0.510-6.
表4-5
k
T2k
S k 1 2

120
1 4
h4
L
可见计算的近似值的算法F(h)的截断误差是O (h2), 而算法F1(h)的截断误差是O (h4)。外推一次，精度就提 高了。这就是外推法的基本思想。若重复以上过程，不
断外推，即不断折半步长h, 得到计算的算法序列 {Fk(h)}。随着k的增加，算法的截断误差阶越来越高，计 算精度越来越好。
h p
)

Fk
q pk 1
(h)
,k

0,1, 2,L
定义的序列{Fk(h)}有
Fk (h) F (0) an(n)1h pn1
a h (n) pn2 n2
a h (n) pn3 n3
L
,
其中an(n)k (k 1, 2, 3,L )与h 无关，q>1.
Richardson外推法应用非常广泛且有效，下面介绍应用 于数值积分的情形。
① T1 ② T2 ④ T4
③ S1 ⑤ S2
⑥ C1
⑨ C2
⑩ R1
⑦ T8
⑾ T16

⑧S4 ⑿ S8
⒀ C4

⒁ R2

Sn

43T2n

1 3 Tn , Cn

16 15
S2n
1 15
Sn , Rn

64 63 C2n

1 63 Cn ,
T1,T2、T4、T6 ……由梯 形公式、复化梯形公式

5
5!n4
L
若记 h , F (h) 6sin ，则有
6
6
F(h) h2 h4 L , F(h) 1 h2 1 h4 L ,
6 120
2
6 4 120 16
由此构造新的表达式
F1(h)

4F(h) 2 3
F (h)



§4.4 外推原理与Romberg求积方法
4.4.1 外推原理
在科学与工程计算中，很多算法与步长h有关，特别是数值 积分、数值微分和微分方程数值解的问题。对于这些算法，我 们可以通过外推技巧提高计算精度。
例1 计算的近似值。
由函数 sin x 的Taylor展开式有
n sin
n


3
3!n2
4T2n Tn
3

Tn
2ba n
n k 1
f
a 2k

1
1 2

b
n
a

1 3

1 h
3

2

f
a
n1
2
k 1
f
xk

f
b
n
2h
k 1
f

a



k

1 2

h

m=2在Simpson序列中，m=3在Cotes序列中，Tmk 引入上面记号
后，Romberg算法所用到的各个计算公式可统一化为：
T (0) 1

ba
2
f
(a)
f
(b)
T1i

1 2
T1i
1

b 2i
a
2i1 j1
f
(a

2
j 1 2i (b

a) , i

16 15
S2n
1 15
Sn
42
S2n Sn 42 1
（4.4.7）
即将Simpson序列 S2k 按(4.4.7) 作线性组合就可产生收敛
速度更快的新序列 C2k : C1,C2,C4,K
------------Cotes序列
由近似式
I

C2n

1 63
(C2
n
Cn )
Tn Sn 即：
Sn

4 3
T2
n

1 3
Tn

4T2n Tn 4 1
(4.4.6)
这表明在收敛缓慢的梯形数列 T2k 的基础上，若将T2n 与 Tn
按(4.4.6)作线性组合就可产生收敛速度较快的Simpson序列：
S2k ：S1 、S2 、S4...
关于(4.4.3.6)的证明：由(4.4.1)可知：
加速公式) ，就能将粗糙的梯形值Tn逐步加工成精度较高的辛 普森值Sn 、柯特斯值Cn和龙贝格值Rn .
可以证明：当f（x）满足一定条件时，Romberg序列 R2k
比Cotes序列 C2k
能更快地收敛到积分 b f xdx 的真值I。 a
综上可知：在积分区间逐次分半过程中利用公式
证 T2n 计算结果的误差很小。这样直接用计算结果来估计误差的方法
通常称作误差的事后估计法。
按式(4.4.3),积分近似值与 T2n
的误差大致等于
1 3
T2
n
Tn
因此如果
用这个误差值作为 T2n 的一种补偿，可以期望所得到的
T

T2 n

1 3
(T2n
Tn )

4 3
T2
n

1 3
Tn
二、Romberg算法
变步长梯形求积法算法简单，但精度较差，收敛速度较慢， 但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙 贝格求积公式。龙贝格公式又称逐次分半加速法。
Romberg算法是在积分区间逐次分半的过程中，对用复化梯 形法产生的近似值进行加权平均，以获得准确度较高的近似值的 一种方法，具有公式简练，使用方便，结果较可靠的优点。
T2k
< ?
S1 = T20 S2 = T21 S4 = T22
T3k
< ?
C1 = T30 C2 = T31
T4k
< ?
R1 = T40
… …
… … …
Simpson序列
Cotes序列
Romberg序列
注：实际计算中常常只计算到第4列,只使用3次理查逊外推法。
Romberg算法中止准则，一般取同列或同行相邻两数值的误差绝对
)
p

O(hs
),
削去截断误差的主项，得新的算法
q pF( h) F(h)
F1(h)
q qp 1
F (0) O(hs ),
我们称这个计算过程为Richardson外推法。 F1(h)逼近F(0)的截 断误差是 O(hs ).
只要知道F(h)的更加完整的关于h幂的展开式，而无需知 道展开式中各个系数的具体数值，就能重复使用Richardson外 推法，直到截断误差达到容许误差。用归纳法可以证明下面更 一般的定理。