计算方法第五讲2

16 0.9459850 0.0002941 512 0.9460830 0.0000003
§4.3 龙贝格 （Romberg）算法
一、问题 复化梯形公式简单，但计算精度较差 有没有办法改善复化梯形公式的精度呢?
由上节关于复化梯形公式的误差公式
I
T2n
1 3
(T2
n
Tn )
可得
I
T2n
1 3
T2n
4I 4T2n I Tn
I
T2n
1 3
(T2
n
Tn )
I
T2n
1 3
(T2n
Tn )
T2n
4
1
1
(T2
n
Tn )
实际计算时，给定误差，若
1 3 (T2n
Tn )
,
即
(T2n
Tn )
3，则停止计算，取I
=T2
；
n
否则，区间继续减半计算T4n，并检验误差。
类似地可以推得
S2n作为I的近似值的截断误差约 为
x 导数，并估计其最大值M。 为此将函数变为如下形式：
f (x)= sin x
1
cos txdt
x
0
L f (x) 1t2 cos txdt 0
1
f (x) 0 tsintxdt
f (k) (x)
1t k cos
0
tx
k 2
dt
于是有
f (k) (x)
1 tk cos
0
比较三种复化公式的的余项
R(Tn
)
I
Tn
(b a) 12
h2
f
( )
o(h2 )
R(Sn ) I Sn b a h4 f (4) ()
2880
o(h4 )
R(Cn )
I
Cn
2(b a) 945
h 4
6
f
( 6 ) ( )
o(h6 )
分别是h的2，4，6阶无穷小量
即当n 时，Tn , Sn ,Cn都收敛于I
§5.3 复化求积公式
一、问题 当积分区间[a, b]的长度较大时，
(1) 如果使用低阶的Newton-Cotes公式，误差将会较大； (2) 如果使用高阶的Newton-Cotes公式，则系数复杂，
且公式的稳定性和收敛性也不能保证。
为了提高公式的精度,又使算法简单易行,往往使用复化求积法.
所谓复化求积法，就是
（1）计算f (a)，f (b)，按梯形公式计算T1。
（2）将区间[a,
b]对半分，计算f
(
ab 2
)和T2，进而
计算S1。
（3）将区间再对半分，计算T4，进而计算S2，C1。
（4）将区间再对半分，计算T8，进而计算S4，C2，R1。 （5）将区间再对半分，重复第4步过程，计算T16，
进而计算S8，C4，R2。反复进行此过程，直到 满足精度 R2k1 R2k 10m为止。此时，取R2k1 为积分近似值，精度达到10m.
精度最低 精度次高 精度最高
三种方法的工作量基本相同。收敛速度一个比一 个快。实际应用时，较多地应用辛普森公式。
四、自动选取积分步长
事前确定步长的问题 (1) 高阶导数的估计往往是很困难的； (2) 这种估计往往是很保守的，得到的n往往偏大。
为了改正上述缺点，实际常采用“事后估计法”
“事后估计法”的基本思想是 (1) 求数值积分时，将区间逐次分半； (2) 利用前后两次的计算结果来判断误差是否满足精度要求，
在子区间[xk , xk1](k 0,1,L , n 1)上使用Newton-Cotes公式
二、复化求积公式
复化梯形公式
b
I f (x)dx
n1
xk1 f ( x)dx
a
k 0 xk
n1 xk1 xk
k 0
2
f (xk ) f (xk1)
n1 h k0 2
f (xk ) f (xk1)
C2
1 [7 f 180
(0)
1
[32 f
k 0
(
x
k
1
)
12
f
4
(
x
k
2
)
32
f
4
(
x
k
3
)]
4
1
14 f (xk ) 7 f (1)] k 1
0.9460830
三个公式的结果比较 原积分的精确值为
I 1sin x dx 0.946083070367183 0x
T8 0.9456908 S4 0.9460832 C2 0.9460830
1 0.9397933 0.9460869 0.9460831 0.9460830
2 0.9445135 0.9460834 0.9460830
3 0.9456909 0.9460830
4 0.9459850
这里只用了T1，T2，T4，T8，T16 就得到了精度较高 的 R2 ，效果显著。
二、龙贝格求积公式的计算步骤
945 4
h3 12
f (k )
h5 2880
f
(4) (k )
8 945
h 4
7
f (6) (k )
复化梯形公式的截断误差
1. 设函数f (x)在[a,b]连续,则
I
Tn
n1
[
h3
k 0 12
f (k )]
h3 12
n1 k 0
f
(k
)
由于
min { f (x)} n1 f (k ) max f (x)
(1) 首先把整个积分区间分成若干个子区间（通常采用等分) ； (2) 然后在每个小区间上使用低阶Newton-Cotes公式;
(3) 最后将每个小区间上的积分的近似值相加。
Байду номын сангаас
将定积分 b f (x)dx的积分区间[a,b]分割为n等份 a
各节点为 xk a kh , k 0,1, , n
h ba n
0.9088516 0.8771925 0.8414709
分别由复化梯形、复化辛普森。复化柯特斯公式有
T8
1[ 16
f
(0)
7
2
k 1
f
( xk
)
f (1)]
0.9456908
S4
1 [ f (0) 4 3
24
k 0
3
f
(
x
k
1 2
)
2
k
1
f (xk )
f (1)]
0.9460832
),
4
3
T2n
1 3
Tn
2 3
Tn
2 h n1 3 k0
f
1
(
x k
1 2
)
3
Tn
Sn
h
n1
2 n1
[ 6
f
(a) 2
k 1
f
(xk )
f
(b)]
3
k 0
f
(
xk
1 2
)
h
n1
n1
[ 6
f
(a) 2
k 1
f
(xk ) 4
k 0
f
(
xk
1 2
)
f
(b)]
为了便于记忆，将上式写成
4
x
0
1
2
x3
x 1
1
2
x
0
3
4
x4
x2
x1
x5
x 2
1
2
x 1
1
4
x6
x3
x 1
1
2
x7 x8
x 3
1
2
x4
x 1
3
4
x2
xi
0 0.125
0.25 0.375
0.5 0.625
0.75 0.875
1
f ( xi )
1 0.9973978
0.9896158 0.9767267
0.9588510 0.9361556
b f (x)dx,而且
a
收敛速度一个比一个快。
例1 要用复化梯形公式计算
2
3ln xdx
的值，要求误差
1
不超过102，则要将区间分成多少份？
解: f x 3ln x f 'x 3 1
x
f
''
x
3
1 x2
M2
max
1x2
f '' x
=3
R
Tn
(b a) 12
h2
f
(
)
RTn
I Sn
n1 k 0
h5 2880
f
(4) (k
)
ba 2880
h4
f
(4) ()
R(Sn )
复化cotes公式的截断误差
若函数f (6) (x)在[a,b]上连续，则
I
Cn
n1 k 0
8 945
h 4
7
f
(6)
(k
)
2(b a) 945
h 4
6
f (6)()
R(Cn )
I
S2n
1 42
1
(S2
n
Sn )
C2n作为I的近似值的截断误差约 为
I
C2n
1 43
1
(C2n
Cn )
例3 使用复化梯形公式和步长逐次减半法计算
I 1 sin x dx，要求误差不超过 1 106。
0x
2
解: 计算结果见下表
2n T2n
|T2n-Tn|
1 0.9207355
2n T2n
f
(
xk
2 4
)
32
f
(
x k
3 4
)]
14
k 1
f
(xk ) 7
f
(b)]
Cn
三、复化求积公式的截断误差
我们上一节知道，三个求积公式的截断误差分别为
单纯的求积公式
复化求积公式在每个小区间
R(T ) (b a)3 f ()
12
R(S) b a5 f (4) ()
2880
R(C ) 8 (b a )7 f (6) ()
从而确定n. 下面以复化梯形公式为例来介绍这种步长逐次减半求积法
将区间[a,b]
n等分，步长
h
ba n
，复化梯形公式为
h
n1
Tn
[ 2
f
(a)
2
k 1
f
(xk )
f
(b)],
再将步长减半，增加新分点
xk
1 2
1 2
xk
xk1 ，
公式变为
h
n1
n1
T2n
[ 4
f
(a) 2
k 1
f
(xk ) 2
Tn
4 3
T2n
1 3
Tn
此式可视作用误差对 T2n 作一种补偿，可以提高精度。
对上式进行计算
Tn
h[ 2
f
n1
(a) 2
k 1
f
(xk )
f
(b)],
T2n
h[ f 4
n1
(a) 2
k 1
f
n1
(xk ) 2
k 0
f
(
xk
1 2
)
f
(b)],
1
h n1
2 Tn
2
k 0
f
(
xk
1 2
1 12
1 n
2
3 102
n2 25
n取5，即将区间分成5份即可。
例2 使用各种复化求积公式计算定积分 I 1 sin x dx， 0x 要求按复化辛普森公式计算时误差不超过 1 106。 2
解: 首先，由于lim sin x 1，所以该积分不是广义积分。
x0 x 其次，为了满足精度要求，需要求函数f (x)= sin x 的4阶
k 0
f
(
xk
1 2
)
f
(b)],
1
h n1
2 Tn
2
k 0
f
(
xk
1 2
),
如何根据Tn和T2n来确定误差是否满足要求？
I Tn
(b a) h2 f ()
12
I T2n
b a ( h)2 f ( )
12 2
如果二阶导数在区 间[a,b]上变化不大
则有 即
I T2n 1 I Tn 4
Sn
复化柯特斯公式
n1
I
xk1 f (x)dx
k 0 xk
记 xk
i 4
i xk 4 h
（i 1，2，3）,
h 90
n1
[7
k 0
f
( xk
)
32
f
(
x
k
1
4
)
12
f
(
x
k
2
4
)
32
f
(
x
k
3
)
4
7
f
( xk 1 )]
h
n1
n1
[7 90
f
(a)
[32
k 0
f
(
xk
1 4
)
12
1
Sn 4 1T2n 4 1Tn
类似地可以推得
42
1
Cn 42 1 S2n 42 1 Sn
Rn
43 43 1 C2n
1 43 1 Cn
梯形加速公式 抛物线加速公式 龙贝格求积公式
说明
（1）根据上述公式，可以由序列 {T2k }求得序列{S2k } ， {C2k } ， {R2k } 。
（2）上述将变步长的梯形法得到的积分值加工成精度 较高的积分值的求积方法称为龙贝格求积算法。
32 0.9460586
|T2n-Tn|
0.0000736
2 0.9397933 0.0190578 64 0.9460769 0.0000183
4 0.9445135 0.0047202 128 0.9460815 0.0000046
8 0.9456909 0.0011774 256 0.9460827 0.0000012
（3）利用两个系数
4m 4m
1
与
1 4m
1
可以构造出新的求积公式，
但当m>=4时，计算结果差别不大，但增加了计算量， 故一般不用。
例1 下面使用Romberg算法计算 I 1 sin x dx， 0x 要求误差不超过 1 106。 2
解: 计算结果见下表
k T2k
S2k
C2k
R2k
0 0.9207355 0.9461459 0.9460830 0.9460831
h 2
f
n1
a 2
k 1
f
(xk )
f
b
Tn
复化辛普森公式
记
xk
1 2
xk
1 h, 2
n1
I
k 0
xk 1 xk
f (x)dx
n1 k 0
h 6
f
( xk
)
4
f
( xk
1 2
)
f
(
xk
1
)
h
n1
n1
6 f
a 4
f
(
xk
1 2
)
2
f (xk ) f
k 0
k 1
b
例2
用Romberg算法计算积分 I
1 0
4 x2
1
dx，
解: 按照上述步骤计算：
（1）f
(x)
x
4 2
，a 1
0, b
1,
f
(0)
4，f
(1)
2，
T1
1
2
f
(0)
f
(1)
3
（2）f
(
1 2
)
16 5
tx
k 2
dt 1 k 1
取 f（4）(x) 的一个上界M 1。 5
R(Sn
)
ba 2880
h4
f
(
4)
(
)
则只要取 n4 2106 1 138.9, 即n=4即可。 2880 5
各节点的值如下面表格
各节点的值如下表
Trapz Simp. Cotes
x0
x0
x0
x1
x
0
1
2
x
0
1
4
x2
x1
a xb
k0 n
a xb
由介值定理 , [a,b], 使得 n1 f (k ) f ()
k0 n
即有
I
Tn
nh3 12
n1 k 0
f (k )
n
nh3 f () (b a) h2 f ()
12
12
R(Tn )
复化simpson公式的截断误差
若函数f（4）(x)在[a, b]上连续，则