第 5章 数值积分
引言
在数学分析中，我们学习过微积分基
本定理 N e w ton - L e ibniz 公式,
( 5.0.1 )
其中，
( )Fx
是被积函数
( )fx
的原函数。
随着学习的不断深化，发现 N e w ton -
L e ibniz 公式有很大的局限性。
? ???
b
a
b
a aFbFxFdxxf )()()()(
引言
首先，遇到的是一类被积函数
( )fx
没有初
等函数有限形式的原函数，如
等 。 dxe
daL
x
?
?
?
??
1
0
2
0
2
2
si n14
正态分布函数；椭圆周长
?
??
引言
其次，被积函数 ( )fx 由表格形式给出，没有解析形式，也无
法使用 N e w ton - L e ibni z 公式；
第三，常常 ( )fx 本身形式并不复杂，而原函数 ( )Fx 推
导十分冗长，且表达式复杂，给计算结果带来十分不便。
引言
为克服上述许多缺点，定积分计算的数值求解能弥
补上述不足，并可带来满意的结果。
积分数值算法的思想是，首先求被积函数 ( )fx 的一个逼近函数
( )px,即 ( ) ( ) ( )f x p x r x=+,这里 ( )rx 为误差函数，于是
引言
? 由定积分定义
i
n
i
i
n
i
in
iiiii
in
b
a
i
n
i
i
x
xfsS
xxxxfs
bxxxa
xfdxxf
??
? ?
??
?
?
??
????
??????
?????
??
00
1
10
0
0
)()3(
)()2(
...)1(
)(l i m)(
?
?
?
求和
近似
分割
引言
?? ???
???
?
????
??
b
a
n
i
ii
x
n
x
i
ni
dxxfxfS
xx
)()(l i ml i m
}{m a x)4(
0
00
1
?
求极限
的近似值。也是
加权和，是权系数，其中
由此想到机械求积公式
?
?
?
?
?
?
??
????
b
a
i
n
i
iii
n
i
ii
nn
b
a
dxxf
xfxfAA
fRxfA
fRxfAxfAxfAdxxf
)(
)()(
][)(
][)(.,,)()()(
0
0
1100
5.1 Newton-Cotes求积公式
5.1.1 C ot e s 系数
首先，我们考察一种简单情况。设
()y f x=
用节点
(,( ) ),(,( ( ) )a f a b f b
的一次插值多项式代替，即
( ) ( ) ( )11f x L x r x=+
( 5, 1, 1 )
( ) ( ) ( ) ( ) ( )"
1
2
x b x a
f a f b f x a x b
a b a b
x
--
= + + - -
--
),( bax ?
),()(
12
][
][)]()([
2
))()((
2
1
)]()([)(
''
3
''
baf
ab
fR
fRbfaf
ab
dxbxaxf
dxbf
ab
ax
af
ba
bx
dxxf
T
T
b
a
b
a
b
a
?
?
??
??
?
?
???
?
?
?
?
?
?
?
? ?
??
?
）（
其中
所以
? 由 Lagrange插值,任何一的函数 都
可以近似的表示成
其中
)( xfy ?
)()()( xRxLxf nn ??
插值多项式。的是 L a g r a g exfyxlxL
n
j
jjn )()()(
0
?
?
?
? 为简便起见,取节点为等分
? 现在关键是求
njjhaxn abh j,.,,,2,1,0,?????
dxxRdxxLdxxf ba nba ba n ?? ? ?? )()()(
? ??
?
? ?
? ?? ??
?
?
?
?
??
?
?
?
?
?
?
??
?
??
??
b
a
n
ji
i ij
i
b
a
j
n
j
j
n
j
n
j
j
n
j
b
a
j
j
n
j
b
a
j
b
a
n
j
jj
b
a
n
dx
xx
xx
ab
dxxl
ab
C
xfCab
xfdxxl
ab
ab
ydxxldxyxldxxL
0
)(
0
)(
0
00
1
)(
1
)()(
)(])(
1
[)(
])([))(()(
因此就归结为求权
? ??
? ?? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
??????
?????
???
?
?
n
n
ji
i
n
ji
i
n
n
ji
i
b
a
n
ji
i ij
in
j
iji
i
j
dtit
ijn
hdt
ij
it
nh
dx
xx
xx
ab
C
ntbxtax
hijxxhitxx
hdtdxthaxihax
njjhax
n
ab
h
0
00
0
00
)(
)(
11
11;0
,)(,)(
,,
,.,,,2,1,0,
因此
。时时
，知
由
2
1
21
1
)0(
)!11(!11
)1(
2
1
2
)1(
1
1
)1(
)!01(!01
)1(
1
)(
)!(!
)1(
1
0
2
1
0
11
)1(
1
1
0
2
1
0
01
)1(
0
0
0
???
???
?
?
?
??
??
???
?
?
?
?
?
?
?
?
?
? ?
?
?
?
?
?
t
dttC
t
dttC
n
dtit
jnnj
n
n
ji
i
jn
时，仅有两个节点：当
6
1
6
4
6
1
])2(
2
1
)2(
3
1
[
4
1
)]2()2[(
4
1
)2)(1(
)!02(!02
)1(
2
)2(
2
)2(
1
2
0
23
2
0
2
2
0
02
)2(
0
??
?????
????
??
???
?
?
?
?
?
?
CC
tt
dttt
dtttC
n
，同理可得
时当
? 以此类推得 Cotes系数表,
n ( )n
kC
1
{ }1 1,12
2 { }1 1,4,1
6
3
{ }1 1,3,3,18
4
{ }1 7,32,12,32,790
5
{ }1 19,75,50,50,75,19288
6
{ }1 41,216,27,272,27,216,41840
7
{ }1 751,357 7,132 3,298 9,298 9,132 3,357 7,75117280
8
{ }1 989,588 8,928,104 96,454 0,104 96,928,588 8,9 8928350 - - -
Newton Cotes积分公式
][)(
)()()(
)()()(
)(
}{,)(),.,,,2,1,0(
,],[
],[)(
0
)(
0
fRfCab
dxxRdxxLdxxf
xRxLxf
l a g r a g exf
xfxfnj
khax
n
ab
hnba
baxf
n
j
j
n
j
b
a
b
a
b
a
nn
nn
n
jjj
j
???
??
??
??
??
?
?
?
? ? ?
?
则称
插值多项式，即的
为节点作以记
等分，取区间等分
上的连续函数，将是设定义
的截断误差。为
，（其中
称
系数；是其中
积分公式。
?
?
?
? ?
????
?
?
?
?
?
?
?
?
?
?
?
?
b
a
n
b
a
n
b
a
n
n
n
n
ji
i
jn
n
j
dxxf
xxxxxxx
dxxxxxxf
dxxf
n
fR
C o t e sdtit
jnnj
C
C o t e sN e w t o n
)(
)),,, ()()(
)(],.,,,,[
)()(
)!1(
1
][
)(
)!(!
)1(
10
10
)1(
0
0
)(
?
?
??
常用的几个积分公式
? 梯形公式 (n=1)
),()(
12
)(
][
))()((
2
)](
2
1
)(
2
1
)[(][
][][)(,
2
1
''
3
1
1
1
0
baf
ab
fR
bfaf
ab
bfafabfT
fRfTdxxfCC
T
b
a
T
?
?
??
?
?
?
???
????
?
??
且
。则因为
? Simpson公式 (n=2)
),()(
2 8 8 0
)(
][
))()
2
(4)((
6
][
][][)(
6
1
,
6
4
,
6
1
)4(
5
)2(
2
)2(
1
)2(
0
baf
ab
fR
bf
ba
faf
ab
fS
fRfSdxxf
CCC
S
b
a
S
?
?
??
?
?
?
?
?
??
???
?
??
且
所以
因为
? Newton公式 (n=3)
。其中
且
所以
，因为
3
))()2(3)(3)((
8
][
][][)(
8
1
8
3
,
8
3
,
8
1
)3(
3
)3(
2
)3(
1
)3(
0
ab
h
bfhafhafaf
ab
fN
fRfNdxxf
CCCC
b
a
N
?
?
?????
?
?
??
????
?
? Cotes公式 (n=4)
。其中
且
所以
，因为
4
))(7)2(32
)2(12)(32)(7(
90
][
][][)(
90
7
,
90
32
90
12
,
90
32
,
90
7
)4(
4
)4(
3
)4(
2
)4(
1
)4(
0
ab
h
bfhaf
hafhafaf
ab
fC
fRfCdxxf
CCCCC
b
a
C
?
?
???
????
?
?
??
?????
?
例题
756 9 3 1.0
756 9 3.0)
2
1
3
5
1
3
3
4
1
31(
8
1;4469.0)
2
1
2
3
1
4
1
1
(
6
1;57.0
1
1
2
1
2
1
6 9 3 1 4 7 1 8.02ln
1
1
2
1
2
1
?
?????
????
???
???
?
?
?
?
IC o t e s
IN e w t o n
IS i m p s o n
I
dx
x
I
L e i b n i zN e w t o n
dx
x
I
公式得由
公式由
公式由
）（由梯形公式
公式得解：由
。计算
5.1.2 Newton-Cotes公式截断误差及代数精度
],[,
()),,, (2)(1(
)!1(
)(
()),,, (2)(1)(
2
(
)!2(
)(
][
,
],[)(
2 8 8 0
)(
][
],,[)(1.5, 1
0
)1(2
0
)2(3
)4(
5
4
ba
n
ab
h
ndtnttt
n
fh
ndtnttt
n
t
n
fh
fR
C o t e sN e w t o nn
baf
ab
fR
S i m p s o nbaCxf
n
nn
n
nn
n
S
?
?
?
?
?
?
?
?
?
?
???
?
????
?
?
?
?
?
??
?
?
?
??
??
?
?
?
??
其中，
为奇即数）
为偶数）
公式的截断误差为一般地对任意
的截断误差为
公式则设定理
],[)(
4 8 3 8 4 0
)(][
4
)6(
7
bafabfR
C o t e sn
C ?
???
?
??
求积公式的截断误差得令
定义 5.1.1 如果定积分的求积公式对于所有不高于 n 次代数多项
式 f ( x ) 精度 成立，即截断误差 [ ]
0Rf =
，但对于至少 1 个 n + 1 次代
数多项式不能精确成立，则称该求积多项式具有 n 次代数精度。
几个常用的求积公式的代数精度
1.T 公式的代数精度
度公式具有一次的代数精所以
时当
时当
T
dxxfba
ab
bfaf
ab
fT
abxdxxdxxf
xxf
dxxfba
ab
bfaf
ab
fT
abxx d xdxxf
xxf
b
a
b
a
b
a
b
a
b
a
b
a
b
a
b
a
?
? ?
?
? ?
??
?
??
?
?
????
?
??
?
??
?
?
????
?
)()(
2
))()((
2
][
)(
3
1
3
1
)(
)(
)()(
2
))()((
2
][
)(
2
1
2
1
)(
)(
22
3332
2
222
2,S–公式的代数精度
成立所以
时当
][)(
)(
2
1
)
2
4(
6
))()
2
(4)((
6
][
)(
2
1
)(
)(
22
22
fSdxxf
abb
ba
a
ab
bf
ab
faf
ab
fS
abx d xdxxf
xxf
b
a
b
a
b
a
?
???
?
?
?
?
?
?
?
?
?
???
?
?
??
精确成立即
时当
][)(
)(
3
1
)242(
6
))
2
(4(
6
))()
2
(4)((
6
][
)(
3
1
)(
)(
3322
222
332
2
fSdxxf
abbaba
ab
b
ba
a
ab
bf
ab
faf
ab
fS
abdxxdxxf
xxf
b
a
b
a
b
a
?
????
?
?
?
?
?
?
?
?
?
?
?
?
???
?
?
??
精确成立即
时当
][)(
)(
4
1
)(
2
3
6
))33(
2
1
(
6
))
2
(4(
6
))()
2
(4)((
6
][
)(
4
1
)(
)(
443223
332233
333
442
3
fSdxxf
abbabbaa
ab
bbabbaaa
ab
b
ba
a
ab
bf
ab
faf
ab
fS
abdxxdxxf
xxf
b
a
b
a
b
a
?
?????
?
?
?????
?
?
?
?
?
?
?
?
?
?
?
?
???
?
?
??
? 因此 S-公式具有三次代数精度。
? 同理可得 N-公式具有三次代数精度,C-
公式具有五次代数精度。
定理 5.1.2 机械求积公式 (5.1.1 1) 至少具有 n 次代
数精度的充分必要条件是求积公式为插值型。
由定理 5.1.2 知,Ne wto n - C otes 公式至少具有 n 次代
数精度。由 S im pson 公式具有 3 次代数精度,C otes 公式
具有 5 次代数精度启发，对偶阶 Ne wton - C ote s 公式的代
数精度有如下结论。
定理 5.1.3 当 n 为偶数时,N e w ton - Cotes 公式具有 n + 1 次
代数精度。
证明 只需证明当 ( )
1nf x x +=
时，余式 [ ]
0Rf =
。由式 ( 5.1.6 )
并注意 到此时 ( ) ( )
1 1!nfn x+ =+
则
0)),,, (1(][
,2
)),,, (1()(][
22222
0
2
????
???
????
?
??
?
?
?
dumuuuhfR
mtumnn
dtnttthdxxfR
m
m
m
nb
a
n
则
为偶数，所以令因为
?
5.2 复化求积公式
? 将区间 [a,b]适当分割成若干个字区间，
对每个子区间使用求积公式，构成所谓
的复化求积公式，这是提高积分精度的
一个常用的方法。
5.2.1定步长复化求积公式
? 1.复化梯形求积公式
bbaa 2? ))()
2
(2)((
4
))()
2
((
2
2
))
2
()((
2
2
)
2
(
))()((
2
)(
bf
ba
faf
ab
bf
ba
f
ba
ba
faf
ba
h
T
bfaf
ba
hT
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? 一般地将 [a,b]区间 n等分，则
))()((
2
),.,,2,1(],[
),.,,2,1,0(,
2
1
1
jjj
jj
j
xfxf
h
S
T
njxx
njjhax
ba
h
??
?
???
?
?
?
?
公式有使用
对每个子区间
? 所以
)())(2)()((
2
))(2)()((
2
))()()()(
...)()()()((
2
))()((
2
],[)(
1
1
1
211
1
1
1
1
hTjhafbfaf
h
xfbfaf
h
bfbfbfxf
xfxfxfaf
h
xfxf
h
S
hfRSdxxf
n
n
j
n
j
j
n
n
j
jj
n
j
j
T
b
a
n
j
j
?????
???
????
?????
??
??
?
?
??
? ?
?
?
?
?
?
?
?
而
? 而
)(
12
)(
)(
12
],[
)(
1
)(
,
],[)(
),()(
12
],[
''
2
''
3
1
''''
1
1
''
3
??
??
?
??
f
hab
nf
h
hfR
f
n
f
ba
baxf
xxf
h
hfR
T
n
j
j
jjj
n
j
jT
?
????
?
?
???
?
?
?
?
?
故
）使得（则必存在一点
区间上的二阶连续函数是若
定步长复化梯形求积公式算法
S
hSS
njjhafSS
bfaf
S
n
ab
h
hT
nba
n
输出
计算
输入
.3
)3(;,.,,2,1)()2(;
2
)()(;)1(
)(.2;,,.1
?
????
?
?
?
?
? 2.复化 Simpson公式
类似于梯形公式,
],[)(
))()(4)((
6
,.,,2,1],[
1
2
11
1
hfRSdxxf
xfxfxf
h
S
njxx
S
n
j
j
b
a
j
j
jj
jj
??
???
?
??
?
?
?
?
则
上有在每个子区间
? ?
? ?
??
? ?
?
? ?
?
?
?
?
?
?
?
??
?
?
????
?????
???
???
???
n
j
n
j
j
j
n
j
n
j
j
j
nnn
n
n
n
j
j
j
j
n
j
j
xfxf
bfafh
xfxfbfaf
h
fffff
fff
fff
h
xfxfxf
h
S
1 1
2
1
1 1
2
1
2
11
2
2
31
1
2
10
1
2
11
1
))()(2
2
)()(
(
3
))(2)(4)()((
6
4
..,,,,,,,,,
4
4
6
))()(4)((
6
（
而
),()(
2880
)(
),()(
2880
],[
)(
)))
2
1
((2)((
2
)()(
(
3
))())
2
1
((2
2
)()(
(
3
)4(
4
1
1
)4(
5
1
1 1
baf
hab
xxf
h
hfR
hS
hjafjhaf
bfafh
jhafhjaf
bfafh
jji
n
j
iS
n
n
j
n
j
n
j
?
?
??
???
?
?????
?
?
?????
?
?
?
?
?
? ?
?
?
? ?
??
??
定步长复化 Simpson求积公式算法
S
S
h
S
njjhafhjafSS
bfaf
S
n
ab
h
hS
nba
n
输出
计算
输入
.3;
3
)3(
);,.,,2,1()())
2
1
((2)2(;
2
)()(;)1(
)(.2;,,.1
?
???????
?
?
?
?
例题
至少取多少？试问划分数
要使得截断误差不超过公式计算
用复化如果用复化梯形公式和
对于定积分例
n
x d xI
5
2
0
10
2
1
,
S i m p s o n
s i n
1.2.5
?
?
?
?
?
254
10
2
1
48
si n
224
],[
)( si n)
2
(
12
0
2
)(
12
)(
],[
5
2
''2''
2
?
????
?
??
?
??
?
n
nn
hfR
n
f
hab
hfR
T
T
解得
即
解：由截断误差有
?
?
??
?
??
?
。取故解得
即
由
6,1.5
31.0
2
,31.0
10
2
1
s i n
22 8 8 0
],[
)( s i n)(
2 8 8 0
0
2
)(
2 8 8 0
)(
],[
54
)4(4)4(
4
????
??
?
?
?
??
?
??
?
nnh
hhfR
hf
hab
hfR
S
S
?
?
?
?
?
?
5.2.2变步长求积公式
定步长复化求积公式的一个明显缺点是：事先
很难估计分划数 n 使结果达到预期精度。由于适当
加密分点，精度会有所改善，为此采用自动加密分
点的方法，并利用事后估计来控制加密次数，以判
断是否达到预期精度，从而停止计算。
首先我们讨论变步长梯形求积公式。
变步长梯形求积公式
设区间
[ ],ab
划分为 n 等分，即步长
ba
h
n
-
=
，
计算
( )nTh
：然后将区间
[ ],ab
分点加密一倍，即步长
缩小一半为
2
h
，再 计算出
2
2
()n hT
。如果
??? )()
2
(2 hT
h
T nn
则 取 S=
2
2
()n hT
作为定积分的近似值 。已知
( )nTh
，如何计算
2
2
()n hT
且计算量小？
j
j
j xxx
2
11
?
?
因为
)2(
2
2
2
11 j
j
jj
fff
h
S ???
?
?
所以
}2)({
4 1
2
1
1
1
1
2 ???
?
?
?
?
?
????
n
j
j
n
j
jj
n
j
jn
fff
h
ST
))()((
2
1
})(
2
{
2
1
1 2
1
1
1
hHhT
fhff
h
nn
n
j
j
n
j
jj
??
??? ??
?
?
?
?
其中
))
2
1
(()(
1
hjafhhH
n
j
n ?
?
???
算法 5.2.1变步长梯形求积算法
);(
2
1
)2(
);)
2
1
(()1(
)(.3;
2
))()((;;1
)(.2;,,,.1
12
1
2
1
n
n
j
n
n
n
HTT
hjafhH
hT
bfafh
Tabhn
hT
Mba
??
???
?
????
?
?
计算
计算
输入 ?
3;;
2;2.6
t h e nif.5;t h e nif.4
21
212
返回
计算；输出划分数过大，停止
输出
TT
h
hnn
Mn
TTT
???
?
?? ?
变步长 Simpson求积公式
))
2
(2)
2
((
3
1
)
2
(
))(2)((
3
1
)})
2
1
((2)](
2
)()(
[{
3
1
)))
2
1
((2)((
2
)()(
(
3
)(
222
11
1
h
H
h
T
h
S
hHhT
hjafhjhaf
bfaf
h
hjafjhaf
bfafh
hS
S i m p so n
nnn
nn
n
j
n
j
n
j
n
??
??
?????
?
?
?????
?
?
??
?
??
?
所以
公式由复化
重复上述过程。
将分点加密一倍，定积分的近似值，否则
为时，取当
程序实现的基本思想：
)
2
()()
2
(
)
2
()
2
()
2
(
)()()(
22
222
h
ShS
h
S
h
S
h
H
h
T
hShHhT
nnn
nnn
nnn
???
???
??
变步长 Sinmpson求积算法
);(
3
1
);)
2
1
(()2(;
2;2);(
2
1
)1(
)
2
(),
2
(),
2
(.3
);2(
3
1
);
2
()2(;
2
))()((;;1)1(
)(),(),(.2;,,,.1
22
1
12
222
11
1
n
n
j
n
n
nnn
nn
nnn
HTShjafhH
h
hnnHTT
h
S
h
H
h
T
HTS
ba
hfH
bfafh
Tabhn
hShHhT
Mba
?????
????
??
?
?
?
????
?
?
计算
计算
输入 ?
。返回
计算；输出分划数过大，停止
停止计算；输出
3;;.6
t h e nif.5
,t h e nif.4
2121
212
SSTT
Mn
SSS
??
?
?? ?
14
)()
2
(4
)(
14
)()
2
(4
)(
2
2
2
2
?
?
?
?
?
?
hS
h
S
hC
C o t e s
hT
h
T
hS
S i m p s o n
nn
n
nn
n
公式同理也可以推出复化
有如下关系
公式与复化梯形公式复化
例题
。于公式计算，使其误差小用变步长
计算，并估计误差；的复化用
公式；用
计算
5
1
0
10)3(
5)2(
)1(
1
?
?
?
? ?
S i m p so n
S i m p so nn
S i m p so n
x
dx
I
所以有：
则节点
）（
)5,4,3,2,1,0(0
2.0
5
1
)2(
6 9 4 4 4.0
36
25
)
2
1
2
3
1
41(
6
1
))1()
2
1
(4)0((
6
1
1
???
??
?
?
?????
???
iihx
n
ab
h
fffI
i
?
6 9 3 1 5.0
))
2
1
8.1
1
6.1
1
4.1
1
2.1
1
(2
2
1
(
30
1
)]
11
1
8.01
1
6.01
1
4.01
1
2.01
1
(2)
11
1
01
1
[(
5
1
6
1
)(
5
?
???????
?
?
?
?
?
?
?
?
?
??
?
?
?
?? hSI
5
44
5
)4(
103333.1
120
)2.0(
24
2880
],[
]1,0[24
1
24
)(
?
??
???
??
?
?
h
hfR
x
x
xf
S
所以
因为
5.3 Romberg求积公式
? 5.3.1外推法基本思想
以较小的计算量为代价,达到提高数值结
果的精度是外推法的中心思想,
展式起截断误差为由
。逼近
用中心差商设
T a y l o r
xf
h
h
xf
h
xf
hG
h
x
h
xCxf
k
)(
)
2
()
2
(
)(
],
2
,
2
[)(
'
22
???
?
???
?
)
2
()
2
(...)
2
()
2
()
2
()(
2
)()()(),()(
))
2
((
),.,,,2,1(
)!12(2
)(
)(...)()(
2
2
4
4
2
2
'
2'22
)12(
2
2
2
4
4
2
2
'
h
E
hhhh
Gxf
h
h
hOhGxfhOhE
T a y l a r
h
xfh
kj
j
xf
hEhhhhGxf
k
k
k
j
j
j
k
k
??????
???
?
?
?
?
??????
?
?
???
?
???
代得用上式中
。且
展式得的证明可以从无关的常数。是与
其中
)()()(
),.,,3,2()14(
3
1
,
14
)()(4
)(,
14
)()
2
(4
)(
)()()(...
))
2
()(())()((4
4
1
'
2
)1(
2
11
11
1
'
1
2
2
6
6
4
4
''
hOhGxf
kjr
hEhE
hE
hG
h
G
hG
hGxfhEhrhrhr
h
GxfhGxf
j
j
j
k
k
??
???
?
?
?
?
?
?
???????
???
??
此时
。
其中
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
,.,2,1,0
14
)
2
()(4
)(
)()(
),()()(
14
)
2
()(4
)(
)()
2
()(
1
1
1
0
6
2
'
2
11
2
2
211
m
h
GhG
hG
hGhG
hOhGxf
h
GhG
hG
hG
h
GhG
m
mm
m
m
依此类推有：这时
构造和同理可以有
外推法。
上述过程称为而
R i ch a r d s o n
hOhGxf mm ),()()( 22' ???
例 5.3.1 试用 Ri c ha r dson 外推算法式 ( 5.3.5 ),
计算
( )f x x=
在 1x = 的导数值。
解 由式 ( 5.3.1 ) 有
( )
22
hh
xx
Gh
h
+ - -
=
取 h = 0.4,并由式 ( 5.3.5 ) 可算得结果如 下 表所示,
与真值
( )' 1 0,5f =
相比，仅两次外推，效果
极为显著。
K
2
()
k
h
G
1
2
()
k
h
G
2
2
()
k
h
G
0 0.50254481 0.49998873 0.50000015
1 0.50062775 0.49999931
2 0.50015642
5.3.2 Romberg求积算法
首先，我们介绍复化梯形求积公式的截断误差的 Eu le r - Ma c l a ur in 公
式。
定理 5.3.1 设函数,则复化梯形求积公式余式
为
=
)()(],[ hTdxxfhfR nb
aT
?? ?
],[)( 22 baCxf k ??
Ehhh kk ??? 224422,.,???
其中,
n
ab
h
)( ?
?
?
?
??
?
?
n
j
n
jhaf
bfaf
hhT
1
)](
2
)()(
[)(；
)]()([ )12()12(22 afbfb jjjj ?? ????
为与 h 无关的常数 ；
2 jb
( )1,2,,,1j k k=+ L
为与 h,
( )fx
无关的常数
],[)()( )22(22 bafabbE kk ???? ?? ??
。
证明过于复杂，这里从略。
由 于复化梯形求积公式的截断误差有 Eule r - Ma c lar uin 级数形式,
且 j
h
2 项的系数是与 h 无关的常数，因此满足 R icha rdson 外推算法 的条
件,类似式 (5.3.5 ) 有
?
?
?
?
?
?
?
?
?
?
?
?
?
????
?
?
?
?????
?
?
??
?
?
?
?
?
?
)1,2,.,,,1;,.,,2,1(
14
)
2
()
2
(4
)
2
(
.,,3,2,1)]
2
1
(()
2
([
2
1
)
2
(
)]()([
2
)()(
1
1
1
1
1
11
1
mkjkm
h
T
h
T
h
T
kihafh
h
T
h
T
bfaf
ab
hThT
m
j
m
j
m
m
j
m
n
i
kk
(5.3.7 )
其中，
ba
h
n
-
=
,1
2 kn -=
，式 (5.3.7 ) 称为 R omber g 求积公式。
表 5, 3,2 式 ( 5, 3,7 ) 的表格形式
k 1
2
()
k
h
T
2
2
()
k
h
T
3
2
()
k
h
T
4
2
()
k
h
T
0
( 1)
( )1Th
( 3)
( )2Th
( 6)
( )3Th
( 10)
( )4Th
1
( 2)
1
2
()
h
T
( 5)
2
2
()
h
T
( 9)
3
2
()
h
T
( 14)
4
2
()
h
T
2
( 4)
1
4
()
h
T
( 8)
2
4
()
h
T
( 13)
3
4
()
h
T
M
3
( 7)
1
8
()
h
T
( 12)
2
8
()
h
T
M
4
( 1 1)
1
16
()
h
T
M
M
M
说明 表 5.3.2 中 ( i ) ( )1,2,i = L 表示计算次序。当 ( ) ( )
1mmT h T h e+ -<
时,取 ( )
1mTh+
作为满足精度的结果。
例 5,3.2 试用 R omber g 算法计算定积分
)102
1
(s i n
51
0
?
??? ?，dx
x
x
解 由 R omber g 算法 计算结果如 下 表所示。
由于
5
34 102
1)0()0( ???? TT,故取 ( )
4 0 0,94 60 83 06 9T =
作为该定积分的满足精度的结果。
表 5.3.3 例 5.3.2 计算结果
k
1 2() k
hT
2 2() k
hT
3 2() k
hT
4 2() k
hT
0 0.920735492 0.946145881 0.946083003 0.946083069
1 0.939793284 0.946086933 0.946083068
2 0.944513521 0.94608331
3 0.945690863
我们指出,Rombe r g 算法外推次数也不宜过多，一般以 m = 4 或 5
为宜。如果 m 过大，则，
0
14
1
,1
14
4
?
?
?
? mm
m
，
也就是说，外推所得新序列与前一次序列值无多大差别，即精度不
再有明显改善。
Romberg求积算法
1, 输入：
,,ab x；
2,计算 T ( h )
（ 1 ） 0k = ； 1n = ； h b a=- ；
（ 2 ）
)]()([
2
)0( bfaf
h
T ??；
( )0cT =；
3, 计算
)
2
(1
k
h
T
（ 1 ） 1kk=+ ；
（ 2 ）
)]
2
1
(()1([
2
1
)(
1
?
?
?????
n
i
ihafhkTkT；
3 ) 1d = ；
4, 对于
,1,,2,1j k k=- L
做
（ 1 ） 4dd= ；
（ 2 ）
( ) ( ) ( ) ( )( ) ( )1 1 1T j T j T j T j d- = + - - -；
5, if
( )0cT e-<
,输出
( )0T
，停止计算；
6, 2nn= ;
2hh=;
( )0cT=; 返回 3 。
引言
在数学分析中，我们学习过微积分基
本定理 N e w ton - L e ibniz 公式,
( 5.0.1 )
其中，
( )Fx
是被积函数
( )fx
的原函数。
随着学习的不断深化，发现 N e w ton -
L e ibniz 公式有很大的局限性。
? ???
b
a
b
a aFbFxFdxxf )()()()(
引言
首先，遇到的是一类被积函数
( )fx
没有初
等函数有限形式的原函数，如
等 。 dxe
daL
x
?
?
?
??
1
0
2
0
2
2
si n14
正态分布函数；椭圆周长
?
??
引言
其次，被积函数 ( )fx 由表格形式给出，没有解析形式，也无
法使用 N e w ton - L e ibni z 公式；
第三，常常 ( )fx 本身形式并不复杂，而原函数 ( )Fx 推
导十分冗长，且表达式复杂，给计算结果带来十分不便。
引言
为克服上述许多缺点，定积分计算的数值求解能弥
补上述不足，并可带来满意的结果。
积分数值算法的思想是，首先求被积函数 ( )fx 的一个逼近函数
( )px,即 ( ) ( ) ( )f x p x r x=+,这里 ( )rx 为误差函数，于是
引言
? 由定积分定义
i
n
i
i
n
i
in
iiiii
in
b
a
i
n
i
i
x
xfsS
xxxxfs
bxxxa
xfdxxf
??
? ?
??
?
?
??
????
??????
?????
??
00
1
10
0
0
)()3(
)()2(
...)1(
)(l i m)(
?
?
?
求和
近似
分割
引言
?? ???
???
?
????
??
b
a
n
i
ii
x
n
x
i
ni
dxxfxfS
xx
)()(l i ml i m
}{m a x)4(
0
00
1
?
求极限
的近似值。也是
加权和，是权系数，其中
由此想到机械求积公式
?
?
?
?
?
?
??
????
b
a
i
n
i
iii
n
i
ii
nn
b
a
dxxf
xfxfAA
fRxfA
fRxfAxfAxfAdxxf
)(
)()(
][)(
][)(.,,)()()(
0
0
1100
5.1 Newton-Cotes求积公式
5.1.1 C ot e s 系数
首先，我们考察一种简单情况。设
()y f x=
用节点
(,( ) ),(,( ( ) )a f a b f b
的一次插值多项式代替，即
( ) ( ) ( )11f x L x r x=+
( 5, 1, 1 )
( ) ( ) ( ) ( ) ( )"
1
2
x b x a
f a f b f x a x b
a b a b
x
--
= + + - -
--
),( bax ?
),()(
12
][
][)]()([
2
))()((
2
1
)]()([)(
''
3
''
baf
ab
fR
fRbfaf
ab
dxbxaxf
dxbf
ab
ax
af
ba
bx
dxxf
T
T
b
a
b
a
b
a
?
?
??
??
?
?
???
?
?
?
?
?
?
?
? ?
??
?
）（
其中
所以
? 由 Lagrange插值,任何一的函数 都
可以近似的表示成
其中
)( xfy ?
)()()( xRxLxf nn ??
插值多项式。的是 L a g r a g exfyxlxL
n
j
jjn )()()(
0
?
?
?
? 为简便起见,取节点为等分
? 现在关键是求
njjhaxn abh j,.,,,2,1,0,?????
dxxRdxxLdxxf ba nba ba n ?? ? ?? )()()(
? ??
?
? ?
? ?? ??
?
?
?
?
??
?
?
?
?
?
?
??
?
??
??
b
a
n
ji
i ij
i
b
a
j
n
j
j
n
j
n
j
j
n
j
b
a
j
j
n
j
b
a
j
b
a
n
j
jj
b
a
n
dx
xx
xx
ab
dxxl
ab
C
xfCab
xfdxxl
ab
ab
ydxxldxyxldxxL
0
)(
0
)(
0
00
1
)(
1
)()(
)(])(
1
[)(
])([))(()(
因此就归结为求权
? ??
? ?? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
??????
?????
???
?
?
n
n
ji
i
n
ji
i
n
n
ji
i
b
a
n
ji
i ij
in
j
iji
i
j
dtit
ijn
hdt
ij
it
nh
dx
xx
xx
ab
C
ntbxtax
hijxxhitxx
hdtdxthaxihax
njjhax
n
ab
h
0
00
0
00
)(
)(
11
11;0
,)(,)(
,,
,.,,,2,1,0,
因此
。时时
，知
由
2
1
21
1
)0(
)!11(!11
)1(
2
1
2
)1(
1
1
)1(
)!01(!01
)1(
1
)(
)!(!
)1(
1
0
2
1
0
11
)1(
1
1
0
2
1
0
01
)1(
0
0
0
???
???
?
?
?
??
??
???
?
?
?
?
?
?
?
?
?
? ?
?
?
?
?
?
t
dttC
t
dttC
n
dtit
jnnj
n
n
ji
i
jn
时，仅有两个节点：当
6
1
6
4
6
1
])2(
2
1
)2(
3
1
[
4
1
)]2()2[(
4
1
)2)(1(
)!02(!02
)1(
2
)2(
2
)2(
1
2
0
23
2
0
2
2
0
02
)2(
0
??
?????
????
??
???
?
?
?
?
?
?
CC
tt
dttt
dtttC
n
，同理可得
时当
? 以此类推得 Cotes系数表,
n ( )n
kC
1
{ }1 1,12
2 { }1 1,4,1
6
3
{ }1 1,3,3,18
4
{ }1 7,32,12,32,790
5
{ }1 19,75,50,50,75,19288
6
{ }1 41,216,27,272,27,216,41840
7
{ }1 751,357 7,132 3,298 9,298 9,132 3,357 7,75117280
8
{ }1 989,588 8,928,104 96,454 0,104 96,928,588 8,9 8928350 - - -
Newton Cotes积分公式
][)(
)()()(
)()()(
)(
}{,)(),.,,,2,1,0(
,],[
],[)(
0
)(
0
fRfCab
dxxRdxxLdxxf
xRxLxf
l a g r a g exf
xfxfnj
khax
n
ab
hnba
baxf
n
j
j
n
j
b
a
b
a
b
a
nn
nn
n
jjj
j
???
??
??
??
??
?
?
?
? ? ?
?
则称
插值多项式，即的
为节点作以记
等分，取区间等分
上的连续函数，将是设定义
的截断误差。为
，（其中
称
系数；是其中
积分公式。
?
?
?
? ?
????
?
?
?
?
?
?
?
?
?
?
?
?
b
a
n
b
a
n
b
a
n
n
n
n
ji
i
jn
n
j
dxxf
xxxxxxx
dxxxxxxf
dxxf
n
fR
C o t e sdtit
jnnj
C
C o t e sN e w t o n
)(
)),,, ()()(
)(],.,,,,[
)()(
)!1(
1
][
)(
)!(!
)1(
10
10
)1(
0
0
)(
?
?
??
常用的几个积分公式
? 梯形公式 (n=1)
),()(
12
)(
][
))()((
2
)](
2
1
)(
2
1
)[(][
][][)(,
2
1
''
3
1
1
1
0
baf
ab
fR
bfaf
ab
bfafabfT
fRfTdxxfCC
T
b
a
T
?
?
??
?
?
?
???
????
?
??
且
。则因为
? Simpson公式 (n=2)
),()(
2 8 8 0
)(
][
))()
2
(4)((
6
][
][][)(
6
1
,
6
4
,
6
1
)4(
5
)2(
2
)2(
1
)2(
0
baf
ab
fR
bf
ba
faf
ab
fS
fRfSdxxf
CCC
S
b
a
S
?
?
??
?
?
?
?
?
??
???
?
??
且
所以
因为
? Newton公式 (n=3)
。其中
且
所以
，因为
3
))()2(3)(3)((
8
][
][][)(
8
1
8
3
,
8
3
,
8
1
)3(
3
)3(
2
)3(
1
)3(
0
ab
h
bfhafhafaf
ab
fN
fRfNdxxf
CCCC
b
a
N
?
?
?????
?
?
??
????
?
? Cotes公式 (n=4)
。其中
且
所以
，因为
4
))(7)2(32
)2(12)(32)(7(
90
][
][][)(
90
7
,
90
32
90
12
,
90
32
,
90
7
)4(
4
)4(
3
)4(
2
)4(
1
)4(
0
ab
h
bfhaf
hafhafaf
ab
fC
fRfCdxxf
CCCCC
b
a
C
?
?
???
????
?
?
??
?????
?
例题
756 9 3 1.0
756 9 3.0)
2
1
3
5
1
3
3
4
1
31(
8
1;4469.0)
2
1
2
3
1
4
1
1
(
6
1;57.0
1
1
2
1
2
1
6 9 3 1 4 7 1 8.02ln
1
1
2
1
2
1
?
?????
????
???
???
?
?
?
?
IC o t e s
IN e w t o n
IS i m p s o n
I
dx
x
I
L e i b n i zN e w t o n
dx
x
I
公式得由
公式由
公式由
）（由梯形公式
公式得解：由
。计算
5.1.2 Newton-Cotes公式截断误差及代数精度
],[,
()),,, (2)(1(
)!1(
)(
()),,, (2)(1)(
2
(
)!2(
)(
][
,
],[)(
2 8 8 0
)(
][
],,[)(1.5, 1
0
)1(2
0
)2(3
)4(
5
4
ba
n
ab
h
ndtnttt
n
fh
ndtnttt
n
t
n
fh
fR
C o t e sN e w t o nn
baf
ab
fR
S i m p s o nbaCxf
n
nn
n
nn
n
S
?
?
?
?
?
?
?
?
?
?
???
?
????
?
?
?
?
?
??
?
?
?
??
??
?
?
?
??
其中，
为奇即数）
为偶数）
公式的截断误差为一般地对任意
的截断误差为
公式则设定理
],[)(
4 8 3 8 4 0
)(][
4
)6(
7
bafabfR
C o t e sn
C ?
???
?
??
求积公式的截断误差得令
定义 5.1.1 如果定积分的求积公式对于所有不高于 n 次代数多项
式 f ( x ) 精度 成立，即截断误差 [ ]
0Rf =
，但对于至少 1 个 n + 1 次代
数多项式不能精确成立，则称该求积多项式具有 n 次代数精度。
几个常用的求积公式的代数精度
1.T 公式的代数精度
度公式具有一次的代数精所以
时当
时当
T
dxxfba
ab
bfaf
ab
fT
abxdxxdxxf
xxf
dxxfba
ab
bfaf
ab
fT
abxx d xdxxf
xxf
b
a
b
a
b
a
b
a
b
a
b
a
b
a
b
a
?
? ?
?
? ?
??
?
??
?
?
????
?
??
?
??
?
?
????
?
)()(
2
))()((
2
][
)(
3
1
3
1
)(
)(
)()(
2
))()((
2
][
)(
2
1
2
1
)(
)(
22
3332
2
222
2,S–公式的代数精度
成立所以
时当
][)(
)(
2
1
)
2
4(
6
))()
2
(4)((
6
][
)(
2
1
)(
)(
22
22
fSdxxf
abb
ba
a
ab
bf
ab
faf
ab
fS
abx d xdxxf
xxf
b
a
b
a
b
a
?
???
?
?
?
?
?
?
?
?
?
???
?
?
??
精确成立即
时当
][)(
)(
3
1
)242(
6
))
2
(4(
6
))()
2
(4)((
6
][
)(
3
1
)(
)(
3322
222
332
2
fSdxxf
abbaba
ab
b
ba
a
ab
bf
ab
faf
ab
fS
abdxxdxxf
xxf
b
a
b
a
b
a
?
????
?
?
?
?
?
?
?
?
?
?
?
?
???
?
?
??
精确成立即
时当
][)(
)(
4
1
)(
2
3
6
))33(
2
1
(
6
))
2
(4(
6
))()
2
(4)((
6
][
)(
4
1
)(
)(
443223
332233
333
442
3
fSdxxf
abbabbaa
ab
bbabbaaa
ab
b
ba
a
ab
bf
ab
faf
ab
fS
abdxxdxxf
xxf
b
a
b
a
b
a
?
?????
?
?
?????
?
?
?
?
?
?
?
?
?
?
?
?
???
?
?
??
? 因此 S-公式具有三次代数精度。
? 同理可得 N-公式具有三次代数精度,C-
公式具有五次代数精度。
定理 5.1.2 机械求积公式 (5.1.1 1) 至少具有 n 次代
数精度的充分必要条件是求积公式为插值型。
由定理 5.1.2 知,Ne wto n - C otes 公式至少具有 n 次代
数精度。由 S im pson 公式具有 3 次代数精度,C otes 公式
具有 5 次代数精度启发，对偶阶 Ne wton - C ote s 公式的代
数精度有如下结论。
定理 5.1.3 当 n 为偶数时,N e w ton - Cotes 公式具有 n + 1 次
代数精度。
证明 只需证明当 ( )
1nf x x +=
时，余式 [ ]
0Rf =
。由式 ( 5.1.6 )
并注意 到此时 ( ) ( )
1 1!nfn x+ =+
则
0)),,, (1(][
,2
)),,, (1()(][
22222
0
2
????
???
????
?
??
?
?
?
dumuuuhfR
mtumnn
dtnttthdxxfR
m
m
m
nb
a
n
则
为偶数，所以令因为
?
5.2 复化求积公式
? 将区间 [a,b]适当分割成若干个字区间，
对每个子区间使用求积公式，构成所谓
的复化求积公式，这是提高积分精度的
一个常用的方法。
5.2.1定步长复化求积公式
? 1.复化梯形求积公式
bbaa 2? ))()
2
(2)((
4
))()
2
((
2
2
))
2
()((
2
2
)
2
(
))()((
2
)(
bf
ba
faf
ab
bf
ba
f
ba
ba
faf
ba
h
T
bfaf
ba
hT
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? 一般地将 [a,b]区间 n等分，则
))()((
2
),.,,2,1(],[
),.,,2,1,0(,
2
1
1
jjj
jj
j
xfxf
h
S
T
njxx
njjhax
ba
h
??
?
???
?
?
?
?
公式有使用
对每个子区间
? 所以
)())(2)()((
2
))(2)()((
2
))()()()(
...)()()()((
2
))()((
2
],[)(
1
1
1
211
1
1
1
1
hTjhafbfaf
h
xfbfaf
h
bfbfbfxf
xfxfxfaf
h
xfxf
h
S
hfRSdxxf
n
n
j
n
j
j
n
n
j
jj
n
j
j
T
b
a
n
j
j
?????
???
????
?????
??
??
?
?
??
? ?
?
?
?
?
?
?
?
而
? 而
)(
12
)(
)(
12
],[
)(
1
)(
,
],[)(
),()(
12
],[
''
2
''
3
1
''''
1
1
''
3
??
??
?
??
f
hab
nf
h
hfR
f
n
f
ba
baxf
xxf
h
hfR
T
n
j
j
jjj
n
j
jT
?
????
?
?
???
?
?
?
?
?
故
）使得（则必存在一点
区间上的二阶连续函数是若
定步长复化梯形求积公式算法
S
hSS
njjhafSS
bfaf
S
n
ab
h
hT
nba
n
输出
计算
输入
.3
)3(;,.,,2,1)()2(;
2
)()(;)1(
)(.2;,,.1
?
????
?
?
?
?
? 2.复化 Simpson公式
类似于梯形公式,
],[)(
))()(4)((
6
,.,,2,1],[
1
2
11
1
hfRSdxxf
xfxfxf
h
S
njxx
S
n
j
j
b
a
j
j
jj
jj
??
???
?
??
?
?
?
?
则
上有在每个子区间
? ?
? ?
??
? ?
?
? ?
?
?
?
?
?
?
?
??
?
?
????
?????
???
???
???
n
j
n
j
j
j
n
j
n
j
j
j
nnn
n
n
n
j
j
j
j
n
j
j
xfxf
bfafh
xfxfbfaf
h
fffff
fff
fff
h
xfxfxf
h
S
1 1
2
1
1 1
2
1
2
11
2
2
31
1
2
10
1
2
11
1
))()(2
2
)()(
(
3
))(2)(4)()((
6
4
..,,,,,,,,,
4
4
6
))()(4)((
6
（
而
),()(
2880
)(
),()(
2880
],[
)(
)))
2
1
((2)((
2
)()(
(
3
))())
2
1
((2
2
)()(
(
3
)4(
4
1
1
)4(
5
1
1 1
baf
hab
xxf
h
hfR
hS
hjafjhaf
bfafh
jhafhjaf
bfafh
jji
n
j
iS
n
n
j
n
j
n
j
?
?
??
???
?
?????
?
?
?????
?
?
?
?
?
? ?
?
?
? ?
??
??
定步长复化 Simpson求积公式算法
S
S
h
S
njjhafhjafSS
bfaf
S
n
ab
h
hS
nba
n
输出
计算
输入
.3;
3
)3(
);,.,,2,1()())
2
1
((2)2(;
2
)()(;)1(
)(.2;,,.1
?
???????
?
?
?
?
例题
至少取多少？试问划分数
要使得截断误差不超过公式计算
用复化如果用复化梯形公式和
对于定积分例
n
x d xI
5
2
0
10
2
1
,
S i m p s o n
s i n
1.2.5
?
?
?
?
?
254
10
2
1
48
si n
224
],[
)( si n)
2
(
12
0
2
)(
12
)(
],[
5
2
''2''
2
?
????
?
??
?
??
?
n
nn
hfR
n
f
hab
hfR
T
T
解得
即
解：由截断误差有
?
?
??
?
??
?
。取故解得
即
由
6,1.5
31.0
2
,31.0
10
2
1
s i n
22 8 8 0
],[
)( s i n)(
2 8 8 0
0
2
)(
2 8 8 0
)(
],[
54
)4(4)4(
4
????
??
?
?
?
??
?
??
?
nnh
hhfR
hf
hab
hfR
S
S
?
?
?
?
?
?
5.2.2变步长求积公式
定步长复化求积公式的一个明显缺点是：事先
很难估计分划数 n 使结果达到预期精度。由于适当
加密分点，精度会有所改善，为此采用自动加密分
点的方法，并利用事后估计来控制加密次数，以判
断是否达到预期精度，从而停止计算。
首先我们讨论变步长梯形求积公式。
变步长梯形求积公式
设区间
[ ],ab
划分为 n 等分，即步长
ba
h
n
-
=
，
计算
( )nTh
：然后将区间
[ ],ab
分点加密一倍，即步长
缩小一半为
2
h
，再 计算出
2
2
()n hT
。如果
??? )()
2
(2 hT
h
T nn
则 取 S=
2
2
()n hT
作为定积分的近似值 。已知
( )nTh
，如何计算
2
2
()n hT
且计算量小？
j
j
j xxx
2
11
?
?
因为
)2(
2
2
2
11 j
j
jj
fff
h
S ???
?
?
所以
}2)({
4 1
2
1
1
1
1
2 ???
?
?
?
?
?
????
n
j
j
n
j
jj
n
j
jn
fff
h
ST
))()((
2
1
})(
2
{
2
1
1 2
1
1
1
hHhT
fhff
h
nn
n
j
j
n
j
jj
??
??? ??
?
?
?
?
其中
))
2
1
(()(
1
hjafhhH
n
j
n ?
?
???
算法 5.2.1变步长梯形求积算法
);(
2
1
)2(
);)
2
1
(()1(
)(.3;
2
))()((;;1
)(.2;,,,.1
12
1
2
1
n
n
j
n
n
n
HTT
hjafhH
hT
bfafh
Tabhn
hT
Mba
??
???
?
????
?
?
计算
计算
输入 ?
3;;
2;2.6
t h e nif.5;t h e nif.4
21
212
返回
计算；输出划分数过大，停止
输出
TT
h
hnn
Mn
TTT
???
?
?? ?
变步长 Simpson求积公式
))
2
(2)
2
((
3
1
)
2
(
))(2)((
3
1
)})
2
1
((2)](
2
)()(
[{
3
1
)))
2
1
((2)((
2
)()(
(
3
)(
222
11
1
h
H
h
T
h
S
hHhT
hjafhjhaf
bfaf
h
hjafjhaf
bfafh
hS
S i m p so n
nnn
nn
n
j
n
j
n
j
n
??
??
?????
?
?
?????
?
?
??
?
??
?
所以
公式由复化
重复上述过程。
将分点加密一倍，定积分的近似值，否则
为时，取当
程序实现的基本思想：
)
2
()()
2
(
)
2
()
2
()
2
(
)()()(
22
222
h
ShS
h
S
h
S
h
H
h
T
hShHhT
nnn
nnn
nnn
???
???
??
变步长 Sinmpson求积算法
);(
3
1
);)
2
1
(()2(;
2;2);(
2
1
)1(
)
2
(),
2
(),
2
(.3
);2(
3
1
);
2
()2(;
2
))()((;;1)1(
)(),(),(.2;,,,.1
22
1
12
222
11
1
n
n
j
n
n
nnn
nn
nnn
HTShjafhH
h
hnnHTT
h
S
h
H
h
T
HTS
ba
hfH
bfafh
Tabhn
hShHhT
Mba
?????
????
??
?
?
?
????
?
?
计算
计算
输入 ?
。返回
计算；输出分划数过大，停止
停止计算；输出
3;;.6
t h e nif.5
,t h e nif.4
2121
212
SSTT
Mn
SSS
??
?
?? ?
14
)()
2
(4
)(
14
)()
2
(4
)(
2
2
2
2
?
?
?
?
?
?
hS
h
S
hC
C o t e s
hT
h
T
hS
S i m p s o n
nn
n
nn
n
公式同理也可以推出复化
有如下关系
公式与复化梯形公式复化
例题
。于公式计算，使其误差小用变步长
计算，并估计误差；的复化用
公式；用
计算
5
1
0
10)3(
5)2(
)1(
1
?
?
?
? ?
S i m p so n
S i m p so nn
S i m p so n
x
dx
I
所以有：
则节点
）（
)5,4,3,2,1,0(0
2.0
5
1
)2(
6 9 4 4 4.0
36
25
)
2
1
2
3
1
41(
6
1
))1()
2
1
(4)0((
6
1
1
???
??
?
?
?????
???
iihx
n
ab
h
fffI
i
?
6 9 3 1 5.0
))
2
1
8.1
1
6.1
1
4.1
1
2.1
1
(2
2
1
(
30
1
)]
11
1
8.01
1
6.01
1
4.01
1
2.01
1
(2)
11
1
01
1
[(
5
1
6
1
)(
5
?
???????
?
?
?
?
?
?
?
?
?
??
?
?
?
?? hSI
5
44
5
)4(
103333.1
120
)2.0(
24
2880
],[
]1,0[24
1
24
)(
?
??
???
??
?
?
h
hfR
x
x
xf
S
所以
因为
5.3 Romberg求积公式
? 5.3.1外推法基本思想
以较小的计算量为代价,达到提高数值结
果的精度是外推法的中心思想,
展式起截断误差为由
。逼近
用中心差商设
T a y l o r
xf
h
h
xf
h
xf
hG
h
x
h
xCxf
k
)(
)
2
()
2
(
)(
],
2
,
2
[)(
'
22
???
?
???
?
)
2
()
2
(...)
2
()
2
()
2
()(
2
)()()(),()(
))
2
((
),.,,,2,1(
)!12(2
)(
)(...)()(
2
2
4
4
2
2
'
2'22
)12(
2
2
2
4
4
2
2
'
h
E
hhhh
Gxf
h
h
hOhGxfhOhE
T a y l a r
h
xfh
kj
j
xf
hEhhhhGxf
k
k
k
j
j
j
k
k
??????
???
?
?
?
?
??????
?
?
???
?
???
代得用上式中
。且
展式得的证明可以从无关的常数。是与
其中
)()()(
),.,,3,2()14(
3
1
,
14
)()(4
)(,
14
)()
2
(4
)(
)()()(...
))
2
()(())()((4
4
1
'
2
)1(
2
11
11
1
'
1
2
2
6
6
4
4
''
hOhGxf
kjr
hEhE
hE
hG
h
G
hG
hGxfhEhrhrhr
h
GxfhGxf
j
j
j
k
k
??
???
?
?
?
?
?
?
???????
???
??
此时
。
其中
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
,.,2,1,0
14
)
2
()(4
)(
)()(
),()()(
14
)
2
()(4
)(
)()
2
()(
1
1
1
0
6
2
'
2
11
2
2
211
m
h
GhG
hG
hGhG
hOhGxf
h
GhG
hG
hG
h
GhG
m
mm
m
m
依此类推有：这时
构造和同理可以有
外推法。
上述过程称为而
R i ch a r d s o n
hOhGxf mm ),()()( 22' ???
例 5.3.1 试用 Ri c ha r dson 外推算法式 ( 5.3.5 ),
计算
( )f x x=
在 1x = 的导数值。
解 由式 ( 5.3.1 ) 有
( )
22
hh
xx
Gh
h
+ - -
=
取 h = 0.4,并由式 ( 5.3.5 ) 可算得结果如 下 表所示,
与真值
( )' 1 0,5f =
相比，仅两次外推，效果
极为显著。
K
2
()
k
h
G
1
2
()
k
h
G
2
2
()
k
h
G
0 0.50254481 0.49998873 0.50000015
1 0.50062775 0.49999931
2 0.50015642
5.3.2 Romberg求积算法
首先，我们介绍复化梯形求积公式的截断误差的 Eu le r - Ma c l a ur in 公
式。
定理 5.3.1 设函数,则复化梯形求积公式余式
为
=
)()(],[ hTdxxfhfR nb
aT
?? ?
],[)( 22 baCxf k ??
Ehhh kk ??? 224422,.,???
其中,
n
ab
h
)( ?
?
?
?
??
?
?
n
j
n
jhaf
bfaf
hhT
1
)](
2
)()(
[)(；
)]()([ )12()12(22 afbfb jjjj ?? ????
为与 h 无关的常数 ；
2 jb
( )1,2,,,1j k k=+ L
为与 h,
( )fx
无关的常数
],[)()( )22(22 bafabbE kk ???? ?? ??
。
证明过于复杂，这里从略。
由 于复化梯形求积公式的截断误差有 Eule r - Ma c lar uin 级数形式,
且 j
h
2 项的系数是与 h 无关的常数，因此满足 R icha rdson 外推算法 的条
件,类似式 (5.3.5 ) 有
?
?
?
?
?
?
?
?
?
?
?
?
?
????
?
?
?
?????
?
?
??
?
?
?
?
?
?
)1,2,.,,,1;,.,,2,1(
14
)
2
()
2
(4
)
2
(
.,,3,2,1)]
2
1
(()
2
([
2
1
)
2
(
)]()([
2
)()(
1
1
1
1
1
11
1
mkjkm
h
T
h
T
h
T
kihafh
h
T
h
T
bfaf
ab
hThT
m
j
m
j
m
m
j
m
n
i
kk
(5.3.7 )
其中，
ba
h
n
-
=
,1
2 kn -=
，式 (5.3.7 ) 称为 R omber g 求积公式。
表 5, 3,2 式 ( 5, 3,7 ) 的表格形式
k 1
2
()
k
h
T
2
2
()
k
h
T
3
2
()
k
h
T
4
2
()
k
h
T
0
( 1)
( )1Th
( 3)
( )2Th
( 6)
( )3Th
( 10)
( )4Th
1
( 2)
1
2
()
h
T
( 5)
2
2
()
h
T
( 9)
3
2
()
h
T
( 14)
4
2
()
h
T
2
( 4)
1
4
()
h
T
( 8)
2
4
()
h
T
( 13)
3
4
()
h
T
M
3
( 7)
1
8
()
h
T
( 12)
2
8
()
h
T
M
4
( 1 1)
1
16
()
h
T
M
M
M
说明 表 5.3.2 中 ( i ) ( )1,2,i = L 表示计算次序。当 ( ) ( )
1mmT h T h e+ -<
时,取 ( )
1mTh+
作为满足精度的结果。
例 5,3.2 试用 R omber g 算法计算定积分
)102
1
(s i n
51
0
?
??? ?，dx
x
x
解 由 R omber g 算法 计算结果如 下 表所示。
由于
5
34 102
1)0()0( ???? TT,故取 ( )
4 0 0,94 60 83 06 9T =
作为该定积分的满足精度的结果。
表 5.3.3 例 5.3.2 计算结果
k
1 2() k
hT
2 2() k
hT
3 2() k
hT
4 2() k
hT
0 0.920735492 0.946145881 0.946083003 0.946083069
1 0.939793284 0.946086933 0.946083068
2 0.944513521 0.94608331
3 0.945690863
我们指出,Rombe r g 算法外推次数也不宜过多，一般以 m = 4 或 5
为宜。如果 m 过大，则，
0
14
1
,1
14
4
?
?
?
? mm
m
，
也就是说，外推所得新序列与前一次序列值无多大差别，即精度不
再有明显改善。
Romberg求积算法
1, 输入：
,,ab x；
2,计算 T ( h )
（ 1 ） 0k = ； 1n = ； h b a=- ；
（ 2 ）
)]()([
2
)0( bfaf
h
T ??；
( )0cT =；
3, 计算
)
2
(1
k
h
T
（ 1 ） 1kk=+ ；
（ 2 ）
)]
2
1
(()1([
2
1
)(
1
?
?
?????
n
i
ihafhkTkT；
3 ) 1d = ；
4, 对于
,1,,2,1j k k=- L
做
（ 1 ） 4dd= ；
（ 2 ）
( ) ( ) ( ) ( )( ) ( )1 1 1T j T j T j T j d- = + - - -；
5, if
( )0cT e-<
,输出
( )0T
，停止计算；
6, 2nn= ;
2hh=;
( )0cT=; 返回 3 。
