2004-11-1 1
第四章数值积分与数值微分
§1 数值积分概述
§2 Newton Cotes 公式
§3 Romberg求积法
§4 Gauss型求积公式
§5 数值微分
2004-11-1 2
§4.1 数值积分概述由积分学基本定理Newton-Leibenize公式有：
)()()( aFbFdxxf
b
a
=
∫
但在应用中常会碰到如下情况：
①f(x)的原函数无法用初等函数给出
②f(x)用表格形式给出
③虽然f(x)的原函数能用初等函数表示,
但表达式过于复杂。
2004-11-1 3
∑
∫
=
=≈=
n
k
nkk
b
a
IxfAdxxfI
0
)()(
∫
∑
=
+=+==
b
a
n
k
nkk
fRIfRxfAdxxfI
0
][][)()(
一.求积公式其中R[f]称为求积公式的余项，
),2,1,0( nkx
k
L=
称为求积节点，
),2,1,0( nkA
k
L=
称为求积系数。
k
x
k
A
仅与求积节点的选取有关，而不依赖与被积函数f(x)的具体形式。
2004-11-1 4
二.求积公式的代数精确度衡量一个求积公式好坏的标准。
∫
∑
=
≈
b
a
n
k
kk
xfAdxxf
0
)()(
定义：如果求积公式对于一切不高于m次的代数多项式准确成立，
而对于某个m+1次多项式并不准确成立，
则称上述求积公式具有m次代数精确度，或称为具有m次代数精度。
2004-11-1 5
求积公式的代数精确度（续）
如果要构造具有m次代数精度的求积公式，
m
xxxxf,,,,1)(
2
L=
只要令它对于都能准确成立即可。
∑
∫
=
+++=+++
n
k
m
kmkk
b
a
m
m
xaxaaAdxxaxaa
0
1010
)()( LL
∫∫∫
+++?
b
a
m
m
b
a
b
a
dxxbaxdxadxa L
10
∑∑∑
===
+++=
n
k
m
kkm
n
k
n
k
kkk
xAaxAaAa
000
10
L
2004-11-1 6
求积公式的代数精确度（续）
∫
∑
=
=
b
a
n
k
k
Adx
0
∫
∑
=
=
b
a
n
k
kk
xAxdx
0
∫
∑
=
=
b
a
n
k
m
kk
m
xAdxx
0
M
即若求积公式具有m次代数精度，
k
A应满足上述方程组。反之亦然。
2004-11-1 7
特别地,若取n = m，即取m+1个节点,
m
xx L
0
求积公式的代数精确度（续）
则可通过给定的m+1个节点得到上述含m+1个未知数、m+1个方程的方程组。
若求积节点互异，则
0
111
det
10
10
≠=
m
m
mm
m
xxx
xxx
A
L
MMM
L
L
从而可得如下定理：
2004-11-1 8
求积公式的代数精确度（续）
定理：
使求积公式至少有m次代数精度。
在区间[a,b]上,对于给定m+1个互异节点，
bxxa
m
≤<<≤ L
0 m
AAA L,,
10
，总存在求积系数,
Remark：定出，则求积公式至少具有m次代数精度，但并不一定它具有m次代数精度,要将代入求积公式，如果等式不准确成立，即求积公式具有m 次代数精度，
否则代数精度将大于m。
),,1,0( mkA
k
L=
1+m
x
2004-11-1 9
求积公式的代数精确度（续）
Remark1：代数精度越高，求积公式的适应性越强。
Remark2：凡至少具有零次代数精度的求积公式
∫
∑
=
≈
b
a
n
k
kk
xfAdxxf
0
)()(
一定满足
∫
∑
=
=?
b
a
n
k
k
Adx
0
11
abA
n
k
k
=
∑
=0
从而有即求积系数之和等于积分区间长度，这是求积系数的一个基本特性。
2004-11-1 10
求积公式的代数精确度（续）
例：确定求积公式
)]()0([
2
)]()0([
)(
2
0
hffh
hffh
dxxf
h
′
′
+
+
≈
∫
α
中的待定参数，使其代数精度尽量高，并指明所构造的求积公式所具有的代数精度。
解：求积公式中含一个待定参数当f(x)=1,x时,有
∫
+≡
h
h
dx
0
]11[
2
2004-11-1 11
求积公式的代数精确度（续）
∫
++=
h
hh
h
xdx
0
2
]11[]0[
2
α
故令求积公式对f(x)=x
2
成立,即
∫
×++=
h
hhh
h
dxx
0
222
)202(]0[
2
α
12
1
=α
得
2004-11-1 12
求积公式的代数精确度（续）
3
)( xxf =
将代入已求得的求积公式，显然
]40[
12
]0[
2
3
2
4
0
4
h
h
h
h
dxx
h
++≠
∫
∫
′
′
++≈
h
hff
h
hff
h
dxxf
0
2
)]()0([
12
)]()0([
2
)(
故具有三次代数精度。
]30[
12
]0[
2
2
2
3
0
3
h
h
h
h
dxx
h
++=
∫
令
4
)( xxf =
2004-11-1 13
三.收敛性与稳定性
),0][lim()()(lim
0
0
==
∞→
→
=
∞→
∑
∫
fRdxxfxfA
n
h
n
k
b
a
kk
n
如果如果求积公式对舍入误差不敏感（误差能够控制），则称该求积公式是稳定的。
其中，则称该求积公式是收敛的。
)(max
1
1
≤≤
=
ii
ni
xxh
一个求积公式首先应该是收敛的，其次应该是稳定的。
2004-11-1 14
收敛性与稳定性（续）
则
∑∑
==
=
n
k
n
k
kkkkn
xfAxfAIe
00
*
)()()(
∑∑
==
≤=
n
k
kk
n
k
kk
AA
00
εε
∑
=
=≤
n
k
kn
abAIe
0
)()( εε
则设计算时有绝对误差，即
k
ε)(
k
xf
.)()(
*
kkk
xfxf ε=?
n
nk
εε
≤≤
=
0
max
),,1,0(0 nkA
k
L=>
取若
k
A故当系数全为正时，求积公式是稳定的。
2004-11-1 15
§4.2 Newton Cotes 公式一.插值型求积公式找到一个足够精度的简单函数p(x)代替f(x),于是
∫∫
≈
b
a
b
a
dxxpdxxf )()(
有,把p(x)取成插值多项式，
则可得到插值型求积公式。
设给定节点
bxxxxa
n
≤<<<<≤ L
210
),,1,0()( nkxf
k
L=
并已知这些节点上的函数值
2004-11-1 16
插值型求积公式（续）
由
)(
)!1(
)(
)()()(
1
1
0
x
n
f
xfxlxf
n
n
n
k
kk +
+
=
+
+=
∑
ω
ξ
dxx
n
f
dxxfxldxxf
b
a
b
a
n
n
b
a
n
k
kk
∫∫∫
∑ +
+
=
+
+= )(
)!1(
)(
)()()(
1
)1(
0
ω
ξ
∫
∑ +
+
=
+
+=
b
a
n
n
n
k
kk
dxx
n
f
xfA )(
)!1(
)(
)(
1
)1(
0
ω
ξ
∫
∏
∫
≠
=
==
b
a
n
ki
i
ik
i
b
a
kk
dx
xx
xx
dxxlA
0
)(
)(
)(其中，系数
2004-11-1 17
插值型求积公式（续）
∫
=
b
a
kk
dxxlA )(
当求积系数由截断误差
∫
∏
=
+
+
=
b
a
n
i
i
n
dxxxf
n
fR
0
)1(
)()(
)!1(
1
][ ξ
所唯一确定时，所得的求积公式称为插值型求积公式。
Remark：由截断误差可知，插值型求积公式至少具有n次代数精度。
2004-11-1 18
二,Newton-cotes公式将[a，b]分为n等份，
取节点（k＝0,1,…,n）
khax
k
+=
nabh )(?=
由Lagrange插值公式，可得
∫∫
+
+
′
==
b
a
b
a
knk
n
kk
dx
xxx
x
dxxlA
)()(
)(
)(
1
1
ω
ω
系数可以进一步表示：
k
A
2004-11-1 19
Newton-cotes公式（续）
),,1,0( nkkhax
k
LQ =+=
令x=a+th即有dx=hdt
],0[],[ nba?
hktxx
k
)(?=?
)()1()()()(
1
01
nttthxxxxx
n
nn
==
+
+
LLω
故
)())(()()(
1101 nkkkkkkkn
xxxxxxxxx=
′
+?+
LLω
))(()2)(1(1)1( knkkh
n
= LL
)!()1(! knkh
knn
=
2004-11-1 20
Newton-cotes公式（续）
故
∫
=
+
n
knn
n
k
hdt
knkhhkt
nttth
A
0
1
)!()1(!)(
)()1( L
∫
+
=
n
kn
dtntktkttt
knk
h
0
)()]1()][1([)1(
)!(!
)1(
LL
∫
+
=
n
kn
dtntktkttt
knnk
ab
0
)()1)(1()1(
)!(!
)1(
)( LL
)(
)(
n
k
cab?=
2004-11-1 21
Newton-cotes公式（续）
故求积公式可写为,其
∫
∑
=
+?≈
b
a
n
k
n
k
khxfcabdxxf
0
0
)(
)()()(
)(n
k
c
中称为柯特斯系数，上式称为Newton----Cotes公式。
可以证明。
)()( n
kn
n
k
cc
=
2004-11-1 22
Newton-cotes公式（续）
当n=1时，
2
1
2
)1(
)1()1(
!0!11
)1(
1
0
1
0
201
)1(
0
=
×?=?
=
∫
t
dttc
∫
==
=
1
0
1
0
2
11
)1(
1
2
1
2
1
!0!11
)1(
ttdtc
∫
+
≈
b
a
bfaf
ab
dxxf )]()([
2
)(
)(
上式称为梯形公式。
2004-11-1 23
Newton-cotes公式（续）
n=2可计算
6
1
6
4
6
1
)2(
2
)2(
1
)2(
0
=== ccc
)](
6
1
)
2
(
6
4
)(
6
1
)[()( bf
ba
fafabdxxf
b
a
+
+
+?≈
∫
故有它称为辛浦生（Simpson）公式或抛物线公式
n=4 Newton—Cotes公式为
)](
90
7
)(
90
32
)(
90
12
)(
90
32
)(
90
7
[
432
10
xfxfxf
xfxfabdxxf
b
a
+++
+?≈
∫
）（）（
其中，
)4,,1,0(
0
L=+= kkhax
k
这个公式特别称为柯特斯公式。
2004-11-1 24
Newton-cotes公式（续）
Remark：对n＋1个节点的Newton—Cotes公式，当n
为奇数时，Newton—Cotes公式的代数精度至少为n次；
当n为偶数时，Newton—Cotes公式的代数精度至少为
n＋1次。
只要验证，当n为偶数时，Newton—Cotes公式对
∑
+
=
+
=
1
0
1
)(
n
k
k
kn
xaxP
准确成立即可。
1
)1(
1
)!1()(
+
+
+
+=
n
n
n
anxP
由
∫
+
+
+
+
=?=
b
a
n
n
n
nn
dxx
n
p
IIPR )(
)!1(
)(
][
1
)1(
1
1
ω
ξ
＋
2004-11-1 25
Newton-cotes公式（续）
∫
∏
=
++
=
b
a
n
j
jnn
dxxxaPR
0
11
][）（即
∫
=
+
++
n
n
nn
dsnssshaPR
0
2
11
)(1][ LL）（
令x= a+sh，并注意
jhax
j
+=则，
当n为偶数，即n=2k，k为正整数，再令s=u+k（u=s-k）
∫∫
+=
nn
dsksksksksssdsnsss
00
)2)(12()1)(()1()(1 LLL）（
∫
+++
k
k
dukukuuukuku ))(1()1()1)(( LL＝
2004-11-1 26
Newton-cotes公式（续）
))(1()1()1)(()( kukuuukukuu?+++= LL?记
))(1()1)(()1)(()( kukuuukukuu++?+?=? LL?
))(1()1()1)((1
12
kukuuukuku
k
+?+++=
+
）（LL
)()(1
12
uu
k
=?=
+
）（
从而，为奇函数。
)(u?
故，即当n为偶数时，Newton—Cotes
公式至少有n+1次代数精度。
0][
1
=
+n
PR
2004-11-1 27
三.几种低阶求积公式的截断误差
1.若f(x)在[a,b]上有二阶连续导数，梯形求积公式有下列误差估计：
baf
ab
bfaf
ab
dxxffR
b
a
≤≤
=
+
=
∫
ηη）（
）（
）（）（）（
''
12
][
2
][
3
证：
,))((
!2
)(''
][ dxbxax
f
fR
b
a
=
∫
ξ
由
x依赖于ξ
2004-11-1 28
几种低阶求积公式的截断误差（续）
)(ξf
′′
注意到是倚赖于x的函数，且在[a,b]上连续，故运用积分中值定理，
在[a,b]上存在一点使得：
,0))(( ≤ bxax
,η
)())((''
6
1
))(()(''))()((''
3
baabf
dxbxaxfdxbxaxf
b
a
b
a
≤≤=
=
∫∫
ηη
ηξ
3
)(
12
)(''
][ ab
f
fR=∴
η
2004-11-1 29
几种低阶求积公式的截断误差（续）
2.若f(x)在[a，b]上有四阶连续导数，
辛浦生公式有下列误差估计：
）（
）（
）（）（
）（）（）（
）（
）（
）（
η
η
4
5
)4(4
2880
2180
]
2
4[
6
][
f
ab
f
abab
bf
ba
faf
ab
dxxffR
b
a
=
=
+
+
+
=
∫
ba ≤≤η
2004-11-1 30
几种低阶求积公式的截断误差（续）
）（）（bfbHafaH == )()(
33
)
2
()
2
(
22
33
ba
f
ba
H
ba
f
ba
H
+
′
=
+
′
+
=
+
）（）（
)]()
2
(4)([
6
)(
3333
bH
ba
HaH
ab
dxxH
b
a
+
+
+
=
∫
由于
)()
2
)((
!4
)(
)()(
2
)4(
3
bx
ba
xax
f
xHxf?
+
+=
ξ
3≤
由于辛浦生公式的代数精度为3，为此构造次数的多项式，使满足：)(
3
xH
证：
2004-11-1 31
几种低阶求积公式的截断误差（续）
)]()
2
(4)([
6
)(
)]()
2
(4)([
6
)(][
333
bH
ba
HaH
ab
dxxf
bf
ba
faf
ab
dxxffR
b
a
b
a
+
+
+
=
+
+
+
=
∫
∫
∫
=
b
a
dxxHxf )]()([
3
dxbx
ba
xax
f
b
a
)()
2
)((
!4
)(
2
)4(
+
=
∫
ξ
2004-11-1 32
几种低阶求积公式的截断误差（续）
dxbxcxaxf
dxcxbxaxf
b
a
b
a
)())(()(
))()()((
2)4(
2)4(
=
∫
∫
η
ξ
4)4(
)
2
)()((
180
1
][
ab
abffR
= η
)())((
2880
1
5)4(
baabf ≤≤= ηη
η
由于是依赖于x的函数，在[a,b]上连续，
故可运用积分中值定理，
在[a,b]上存在一点，使
,0)())((
2
≤ bxcxax
)(
)4(
ξf
2004-11-1 33
几种低阶求积公式的截断误差（续）
3.若f(x)在[a，b]上有六阶连续导数，
Cotes公式有下列误差估计：
）（
）（
）（
η
66
)
4
(
945
2
)(][ f
abab
nCIfR
=?=
ba ≤≤η
2004-11-1 34
四.复化求积公式当n 7时，Cotes系数均为正，但从n=8开始，Cotes
系数有正有负，这会使计算误差得不到控制，稳定性得不到保证。另外，节点增多的高次插值多项式舍入误差很大，导致插值型求积公式的舍入误差增大。
≤
因此，实际计算时，一般不采用n较大的Newton-
Cotes公式，而是将区间[a,b]等分为n个小区间，其长度为h=(b-a)/n,在每个小区间上应用低阶的公式，
然后对所有小区间上的计算结果求和，这样得出的
Newton-Cotes公式称为复化求积公式。
2004-11-1 35
1.复化梯形公式
∑∑
=
=
+
′′
+=
1
0
1
0
3
1
)(
12
)]()([
2
n
k
n
k
kkk
f
h
xfxf
h
η
∑
∫∫
=
+
==
1
0
1
)()(
n
k
x
x
b
a
k
k
dxxfdxxfI
由将[a,b]等分为n个子区间
.1,,1,0],[
1
=
+
nkxx
kk
L
][)(
)(
fRnT
nT
+=
),(
1+
∈
kkk
xxη
2004-11-1 36
复化梯形公式（续）
∑
=
+
+=
1
0
1
))()((
2
)(
n
k
kk
xfxf
h
nT
∑
=
++=
1
1
)]()(2)([
2
n
k
k
bfxfaf
h
∑
=
′′
=?=
1
0
3
)(
)(
12
)(][
n
k
knT
f
h
nTIfR η
其中余项若f″(x)在[a,b]上连续,设m为f ″(x)的最小值,M
为f ″(x)的最大值,则
2004-11-1 37
复化梯形公式（续）
故由介值定理,一定存在(a,b)上一点使
η
.)("
1
1
0
Mf
n
m
n
k
k
≤≤
∑
=
η
∑
=
=
1
0
)("
1
)("
n
k
k
f
n
f ηη
)(
12
)("
12
][
2
1
0
3
)(
ηη fh
nh
f
h
fR
n
k
knT
′′
=?=
∑
=
故
),( ba∈η)("
12
)(
2
ηfh
ab?
=
2004-11-1 38
2.复化Simpson公式
n
ab
h
=
(将[a,b]n等分(偶数份),子区间长度).
∑
∫∫
=
+
==
1
2
0
)1(2
2
)()(
n
k
x
x
b
a
k
k
dxxfdxxfI由
∑∑
=
+
=
+
++=
1
2
0
)4(
5
22
1
2
0
122
)(
2880
)2(
)]()(4)([
6
2
n
k
kk
n
k
kk
f
h
xfxfxf
h
η
L++++++= )()(4)()()(4)([
6
2
43221
xfxfxfxfxfaf
h
∑
=
+++
1
2
0
)4(
5
12
)(
2880
)2(
)]()(4)(
n
k
knn
f
h
bfxfxf η
2004-11-1 39
复化Simpson公式（续）
∑∑∑
=
==
+++=
1
2
0
)4(
5
1
2
1
2
2
1
12
)(
2880
)2(
)]()(2)(4)([
6
2
n
k
k
n
k
k
n
k
k
f
h
bfxfxfaf
h
η
][)(
)(
fRnS
nS
+=
其中称为
)]()(2)(4)([
3
)(
1
2
1
2
2
1
12
bfxfxfaf
h
nS
n
k
k
n
k
k
+++=
∑∑
==
复化辛浦生（Simpson）公式。
∑
=
=?=
1
2
0
)4(
5
)(
)(
90
)(][
n
k
knS
f
h
nSIfR η
余项
2004-11-1 40
复化Simpson公式（续）
M
n
m
n
k
k
f
≤≤
∑
=
1
2
0
)4(
)(
2
η
故由介值定理,一定在(a,b)中有一使
η
∑
=
=
1
2
0
)4()4(
)(
2
)(
n
k
k
ff
n
ηη
设在在[a,b]上连续,设m为的最小值,M
为的最大值,则
)(
)4(
xf )(
)4(
xf
)(
)4(
xf
)(
180
)(
)(
290
)(
90
][
)4(4
)4(
5
1
2
0
)4(
5
)(
η
ηη
fh
ab
f
nh
f
h
fR
n
k
knS
=
×?=?=
∑
=
故：
),( ba∈η
2004-11-1 41
3.复化Cotes公式
n
ab
h
=(将[a,b]n等分(n为4的倍数),子区间长度).
)(
945
)(2
)](7)(14)(32
)(12)(32)(7[
90
4
)(
)6(6
1
4
1
4
1
4
0
34
1
4
0
24
1
4
0
14
ηfh
ab
bfxfxf
xfxfaf
h
dxxf
n
k
k
n
k
k
n
k
k
n
k
k
b
a
+++
++≈
∑∑
∑∑
∫
=
=
+
=
+
=
+
Remark:复化梯形公式，复化辛浦生(Simpson)
公式和复化柯特斯(Cotes)公式当时均收敛到所求积分值I。
0→h
2004-11-1 42
例1
若取9个节点，用复化梯形公式、复化辛浦生公式和复化柯特斯公式计算积分，其步长以及与9
个节点所对应的求积系数分别是多少？
解：复化梯形公式：n=8，h=(b-a)/8，对应的求积系数为1、2、2、2、2、2、2、2、1。
复化辛浦生公式：n=8，h=(b-a)/8，对应的求积系数为1、4、2、4、2、4、2、4、1。
复化柯特斯公式：n=8，h=(b-a)/8，对应的求积系数为7、32、12、32、14、32、12、32、7。
2004-11-1 43
用积分计算ln2,要使所得结果的近似值具有5位有效数字。问用复化梯形公式，复化Simpson公式时，至少要取多少个节点?
,2ln2
1
8
2
=
∫
dx
x
例2
解：由且
∫
=
8
2
1
2
1
2ln dx
x
375.0
8
3
8
1
2
11
2
1
8
2
8
2
==>
∫∫
dxdx
x
故计算ln2时，要使误差不超过也即计算2ln2，其误差不超过。
,10
2
1
5?
×
5
10
eln2ln375.0 <<
即
2004-11-1 44
例2（续）
其中
)(
max
2
xfM
bxa
′′
=
≤≤
2
1
)('
x
xf =
3
2
)("
x
xf =
4
1
2
22
max
33
82
2
===
≤≤
x
M
x
2
2
3
2
2
)(
12
)(
12
)(
][1 M
n
ab
Mh
ab
fR
nT
=
≤）由
5
2
3
2
2
3
)(
10
4
1
12
6
12
)(
][
≤×=
≤
n
M
n
ab
fR
nT
令
（
....8203933.67010
412
6
5
3
=×
×
≥n
故区间应取671个，节点至少应取672。
2004-11-1 45
例2（续）
4
4
5
4
4)4(4
)(
180
)(
180
)(
180
][)2( M
n
ab
Mh
ab
fh
ab
fR
nS
=
≤
= η由
)(max
)4(
4
xfM
bxa ≤≤
=
其中由
5
)4(
4
)3(
32
24
)(,
6
)(,
2
)(",
1
)('
x
xf
x
xf
x
xf
x
xf =
==
=
4
3
2
2424
max
55
82
4
=≤=
≤≤
x
M
x
2004-11-1 46
例2（续）
5
4
5
4
4
5
)(
10
4
3
180
6
180
)(
][
≤×=×
=
nn
Mab
fR
nS
令
42640687.4260707106781.0
10625.010
4180
1036
}10
4180
36
{
4
4
5
4
1
5
5
=×=
××=×
×
××
=×
×
×
≥n
故区间n应取44（偶数份），即45个节点。
2004-11-1 47
例3
下面我们来看一个数值算例：
1
4
1
0
2
=
+
∫
dx
x
求解
2004-11-1 48
五.区间逐次分半求积法
1.误差的事后估计法复化求积公式是提高精度的一种有效方法，但在使用复化求积公式之前，必须根据复化求积公式的余项进行先验估计，以确定节点数目，从而确定合适的等分步长。
因为余项表达式中包含了被积函数的导数，而估计各阶导数的最大值往往是很困难的，且估计的误差上界往往偏大。所以实际中，常常使用“事后估计误差”的方法，
通过区间的逐次分半，在步长逐次分半的过程中，反复利用复化求积公式进行计算，查看相继两次计算结果的差值是否达到要求，直到所求得的积分值满足精度要求。
2004-11-1 49
误差的事后估计法（续）
对梯形公式，假定区间分为n等份时，由公式
∫
∑
=++≈=
=
b
a
n
k
k
nTbfxfaf
h
xfI )()]()(2)([
2
)(
1
1
其中算出的积分近似值为T(n),因而有
.khax
k
+=
)("
12
)(
1
2
ηfh
ab
nTI
=
再把各个小区间分为二等份，得积分的近似值为T(2n)，则积分值为
)(")
2
(
12
)2(
2
2
ηf
hab
nTI
=
),(
1
ba∈η
),(
2
ba∈η
2004-11-1 50
误差的事后估计法（续）
22
)
2
(
12
)(
)2(
12
)(
)(
hab
nTI
h
ab
nTI
≈
假定f′′(x)有[a,b]上变化不大，即有，则
)(")("
21
ηη ff ≈
4
)2(
)(
≈
nTI
nTI
即上式可改写为
))()2((
14
1
)2(
))()2((
3
1
)2(
nTnTnT
nTnTnTI
+=
+≈
2004-11-1 51
误差的事后估计法（续）
上式说明用T(2n)做为积分Ｉ的近似值时，其误差大约为
)]()2([
3
1
nTnT?
计算时只需检验是否满足？若不满足，则再把区间分半进行计算，直到满足要求为止。这样就实现了步长得自动选取。
ε<? )()2( nTnT
2004-11-1 52
误差的事后估计法（续）
类似的，还可以得到下面的结论：
对于Simpson公式，假定在在[a,b]上变化不大，则有
)(
)4(
xf
)(
180
)(
)(
1
)4(4
ηfh
ab
nSI
=由),(
1
ba∈η
)()
2
(
180
)(
)2(
2
)4(4
ηf
hab
nSI
= ),(
2
ba∈η
))()2((
14
1
)2())()2((
15
1
)2(
2
nSnSnSnSnSnSI?
+=?+≈
24
42
)2(
)(
=≈
nSI
nSI
变形即得
2004-11-1 53
误差的事后估计法（续）
对于Cotes公式，假定在[a,b]上变化不大，则由
)(
)6(
x
f
))()2((
14
1
)2(
))()2((
63
1
)2(
3
nCnCnC
nCnCnCI
+=
+≈
)(
945
)(2
)(
1
)6(6
ηfh
ab
nCI
=
)()
2
(
945
)(2
)2(
2
)6(6
ηf
hab
nCI
=
36
42
)2(
)(
=≈
nCI
nCI
),(
1
ba∈η
),(
2
ba∈η
得变形即得
2004-11-1 54
2.区间逐次分半的梯形公式
)]()(2)([
2
)(
1
1
bfxfaf
h
nT
n
k
++=
∑
=
由于
)]()(2)([
2
'
)2(
12
1
bfxfaf
h
nT
n
k
k
+
′
+=
∑
=
2
h
h =
′
∑∑
=
=
+++=
n
k
k
n
k
k
xfxfbfaf
h
nT
1
2
1
1
1
)](2)(2)()([
4
)2(故
∑
=
+=
n
k
k
xf
h
nT
1
2
1
)(
2
)(
2
1
∑
=
+
+=
n
k
n
ab
kaf
n
ab
nTnT
1
)
2
)12((
2
)(
2
1
)2(即
2004-11-1 55
区间逐次分半的梯形公式（续）
这就是复化梯形公式当区间逐次分半时的计算公式。由T(n)计算T(2n)时只需计算新增加节点处的函数值。
==
+
+=
+
=
=
∑
L,3,2,1,2
))12(
2
(
2
)(
2
1
)2(
)]()([
2
)1(
1
1
kn
i
n
ab
af
n
ab
nTnT
bfaf
ab
T
k
n
i
算例求解故有递推公式
2004-11-1 56
§4.3 Romberg求积算法一.对低精度公式经过组合构造高精度公式用T(2n)做为积分Ｉ的近似值时，其误差大约为
)]()2([
3
1
nTnT?
若我们不是把误差舍掉，而是把误差补偿到
T(2n)后，我们会得到精度更好的求积公式－
Simpson公式：
)(
14
1
)2(
14
4
)(
3
1
)2(
3
4
)2( nTnTnTnTnS
=?=
2004-11-1 57
对低精度公式经过组合构造高精度公式（续）
)]()(4)(2
)(2)(4)(2)(4)([
6
2
2222
43210
nnn
xfxfxf
xfxfxfxfxf
h
+++
+++++=
L
事实上，
∑
=
++
++=
1
0
22122
)]()(4)([
6
2
)2(
n
k
kkk
xfxfxf
h
nS
)](2)(4)(4)(4)(4)(2[
6
2
2123210 nn
xfxfxfxfxfxf
h
++++++=
L
)]()(2)(2)(2)([
6
2
222420 nn
xfxfxfxfxf
h
+++++?
L
2004-11-1 58
对低精度公式经过组合构造高精度公式（续）
])(2)()([
6
2
])(2)()([
6
4
1
1
2
12
1
∑∑
=
=
++?++=
n
k
k
n
k
k
xfbfaf
h
xfbfaf
h
)(
3
1
)2(
3
4
nTnT?=
14
)()2(4
)(
3
1
)2(
3
4
)2(
=?=
nTnT
nTnTnS即
（1）
Reamrk：两个二阶精度的梯形公式组合起来就得到一个四阶精度的Simpson公式，截断误差提高了二阶，
起到了加速收敛的作用。
2004-11-1 59
对低精度公式经过组合构造高精度公式（续）
)(
14
1
)2(
14
4
))()2((
15
1
)2()2(
22
2
nSnS
nSnSnSnC
=
+=
类似地，可以证明：
（2）
如果再用柯特斯公式作线性组合，
)(
14
1
)2(
14
4
)]()2([
63
1
)2()2(
33
3
nCnC
nCnCnCnR
=
+=
（3）
这个公式称为Romberg公式。
2004-11-1 60
对低精度公式经过组合构造高精度公式（续）
由（1），（2），（3）组成的方法称为Romberg
算法。
序列{T(n)},{S(n)},{C(n)}和{R(n)}分别称为梯形序列，
Simpson序列，Cotes序列和Romberg序列。
上述用若干个积分近似值推算出更为精确的积分近似值的方法，称为外推算法。得到Romberg序列后还可以继续外推，得到新的求积序列，称为Richardson外推算法。但由于在新的求积序列中，其线性组合的系数分别为：
2004-11-1 61
对低精度公式经过组合构造高精度公式（续）
0
14
1
1
14
4
≈
≈
mm
m
因此，新的求积序列与前一个序列结果相差不大，故通常外推到Romberg序列为止。
可以证明，梯形序列，Simpson序列，Cotes序列和
Romberg序列均收敛到积分值，且每次外推可使误差阶提高二阶。
2004-11-1 62
二.Romberg算法的实现
T数表：
R（32）C（32）S（32）T（32）32
R（16）C（16）S（16）T（16）16
R（8）C（8）S（8）T（8）8
C（4）S（4）T（4）4
S（2）T（2）2
T（1）1
R（n）C（n）S（n）T（n）
区间等分数n=2
k
MMMMM
2004-11-1 63
Romberg算法的实现（续）
对上面的T数表作计算，一直到Romberg序列中前后两项之差的绝对值不超过给定的误差限为止。
Remark:Romberg求积算法从最简单的梯形序列开始逐步进行线性加速，具有占用内存少，
精确度高的优点，是实际中最常用的算法之一。
算例求解
2004-11-1 64
§4.4 Gauss型求积公式问题：
若求积公式
∑
∫
=
≈=
n
k
kk
b
a
xfAdxxfI
0
)()(
),,2,1,0(
,
nkAx
kk
L=
中含有2n+2个待定参数我们能否通过节点的选择将求积公式的代数精度从n或者n+1提高到2n+1？
2004-11-1 65
一．带权积分的插值型求积公式
0)( ≥xρ
dxxxfI
b
a
))(（
∫
= ρ
考虑更一般的带权积分,其中为[a,b]上的给定的权函数。
类似于前面插值型求积公式的建立过程，有
)(
)!1(
)(
)()()(
1
)1(
0
x
n
f
xfxlxf
n
n
n
k
kk +
+
=
+
+=
∑
ω
ξ
)
~
,
~
( Mm∈ξ
由
)()(
)!1(
)(
)()}()({)()(
1
)1(
0
xx
n
f
xxfxlxfx
n
n
n
k
kk
ρω
ξ
ρρ
+
+
=
+
+=
∑
得
2004-11-1 66
带权积分的插值型求积公式（续）
积分
)()()()()(
0
k
n
k
b
a
k
b
a
xfdxxxldxxfx
∑
∫∫
=
= ρρ
∫
+
+
+
+
b
a
n
n
dxx
n
f
x )(
)!1(
)(
)(
1
)1(
ω
ξ
ρ
[]
∑
=
+=
n
k
kk
fRxfA
0
)(
[]
+
=
=
∫
∫
+
+
b
a
n
n
b
a
kk
dxx
n
f
xfR
dxxlxA
)(
)!1(
)(
)(
)()(
1
)1(
ω
ξ
ρ
ρ
其中上述公式称为计算的插值型求积公式。
∫
=
b
a
dxxfxI )()(ρ
2004-11-1 67
二．高斯公式和高斯点我们知道，n＋1个插值节点的插值型求积公式的代数精确度不可能低于n。由于在带权插值求积公式中未知的节点及系数共2n＋2个参数。因此，当节点数n＋1确定时，只要适当选取节点位置x
k
及系数A
k
，就可以使上面的插值型求积公式对f(x)=1,x,…,x
2n+1
精确成立，即使其代数精度提高到2n+1。若某个求积公式（n＋1个插值节点）的代数精度达到2n+1，则称该公式为Gauss
型求积公式。
2004-11-1 68
1.定义若是上的一组互异节点，且求积公式
[ ]ba，
n
xxx L,
1,0
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
具有2n+1次代数精度，则称该求积公式为guass
型求积公式，其求积节点（k=0，1，……n）
称为高斯点，系数称为高斯系数。
k
x
k
A
2004-11-1 69
2.引理求积公式至少具有
n次代数精度的充要条件是它是插值型的。
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
证明：充分性：如果求积公式为插值型的，
由
[]
∫
+
+
+
=
b
a
n
n
dxxx
n
f
fR )()(
)!1(
)(
1
)1(
ωρ
ξ
可知，对于任意次数≤n的多项式f(x),R[f]=0。
这时求积公式至少具有n次精度。
必要性：
设求积公式具有n
次代数精度。
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
2004-11-1 70
续则对于不高于n次的多项式，等式精确成立。
),,1,0()( njxl
j
L=对n次插值函数仍有，
特别地，由于为常数，
k
A
∫
∑
=
=
b
a
n
k
kjkj
xlAdxxlx
0
)()()(ρ
),,1,0()()( njdxxlxA
b
a
jj
L==
∫
ρ
根据插值型求积公式定义知，其求积公式为插值型求积公式。
2004-11-1 71
两条结论：
①．高斯型求积公式一定是插值型求积公式，其系数由高斯点唯一确定。
②．高斯型求积公式是精度最高的求积公式。
2004-11-1 72
续事实上，在中，取
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
,()()(
0
22
1 ∏
=
+
==
n
i
in
xxxxf）ω
，则0)()(
2
1
>
∫
+
b
a
n
dxxx ωρ
.0)()(
00
2
1∑∑
==
+
==
n
k
n
k
knkkk
xAxfA ω
，即对)()(
2
1
xxf
n+
=ω
∑
∫
=
++
≠
n
k
knk
b
a
n
xAdxxx
0
2
1
2
1
)()()( ωωρ
也即插值型求积公式的精度不可能达到2n+2。
2004-11-1 73
3.高斯点的选取定理：插值型求积公式中的节点是高斯点的充要条件是，在[a,b]上，以这些点为零点的n+1次多项式与任意次数不超过n的多项式P(x)带权正交，即
),,1,0( nkx
k
L=
∏
=
+
=
n
j
jn
xxx
0
1
)()(ω
)(xρ
∫
=
+
b
a
n
dxxPxx 0)()()(
1
ωρ
2004-11-1 74
定理证明必要性：
设是高斯点，于是对任意次数不超过n
的多项式P(x)，
),,1,0( nkx
k
L=
的次数不超过2n+1。
)()()(
1
xxPxf
n+
= ω
故有
∫
∑
=
++
==
b
a
n
k
kknkn
xPxAdxxPxx
0
11
0)()()()()( ωωρ
充分性：
设对于任意次数不超过
2n+1的多项式设除f(x)的商为P(x),余项为q(x)。
∫
=
+
b
a
n
dxxPxx 0)()()(
1
ωρ
),(xf
)(
1
x
n+
ω
2004-11-1 75
定理证明（续）
)()()()(
1
xqxxPxf
n
+=
+
ω即
nxqxP ≤的次数其中)(),(
∫∫ ∫
+=
+
b
a
b
a
b
a
n
dxxqxdxxxPxdxxxf )()()()()()()(
1
ρωρρ
∫
=
+
b
a
n
dxxxPx，由条件0)()()(
1
ωρ
所给的求积公式是插值型的，其代数精度至少为n。
∫
∑
=
=
b
a
n
n
kk
xqAdxxqx
0
)()()(ρ故
2004-11-1 76
定理证明（续）
Remark：
)()()()()(
1 kkknkk
xqxqxxPxf =+=
+
ω
∑∑
∫∫
==
==
=
n
k
n
k
kkkk
b
a
b
a
xfAxqA
dxxqxdxxxf
00
)()(
)()()()( ρρ
即求积公式具有2n+1次代数精度，
从而为Gauss点。
),,1,0( nkx
k
L=
2004-11-1 77
Reamrk：
当为正交多项式系中的n+1次多项式时，
)(
1
xf
n+
()()0)(),(
1
)(,
1
)1(
1
1
==
+
+
+
xPxf
A
xP
n
n
n
ω
)(
1
x
n+
ω
)1(
1
1
1
)(
)(
+
+
+
=
n
n
n
A
xf
xω
取，则有n+1个互异的零点，且对任意次数不超过n的多项式有
①[a,b]上带权正交的n+1次多项式的零点就是高斯型求积公式的高斯点。
)(xρ
2004-11-1 78
Reamrk（续）
②当高斯点确定以后，高斯系数
),,1,0( nkA
k
L=
也可以由插值型求积公式中的系数公式确定。
∫
=
b
a
kk
dxxlxA )()(ρ
=
=
=
∫
∑
∫
∑
∫
∑
=
=
=
b
a
n
k
n
kk
n
b
a
n
k
kk
b
a
n
k
k
xAdxxx
xAxdxx
Adxx
0
0
0
)(
)(
)(
ρ
ρ
ρ
M
确定。即可由线性方程组
2004-11-1 79
4.高斯型求积公式的收敛性和稳定性定理：高斯型求积公式总是稳定的。
证明：
只需证明Gauss系数全为正即可。事实上，由于高斯型求积公式对次数不超过2n+1的多项式准确成立，故可取（j＝0,1,…,n），
其中l
j
(x)是n次Lagrange插值基函数，故有
)()(
2
xlxf
j
=
∫
∑
=
>==
b
a
n
k
jkjkj
AxlAdxxlx
0
22
0)()()(ρ
（j=0,1,…n）
故高斯求积系数全为正。因此高斯求积公式是稳定的。
k
A
2004-11-1 80
高斯型求积公式的收敛性和稳定性（续）
Remark：当求积节点数目增加时，Newton-
Cotes求积公式的系数变得有正有负（n≥8），
从而不能保证计算的稳定性。但对于Gauss型求积公式，不论节点数n＋1有多大，Gauss型求积公式总是稳定的。
定理：设f(x)∈C[a,b]，则Gauss型求积公式是收敛的。
2004-11-1 81
5.高斯型求积公式的余项定理：
设在内有直到2n+2阶导数，则高斯型求积公式的余项为：
)(xf
[ ]ba,
[]
∫
+
+
+
=
b
a
n
n
dxxx
n
f
fR )()(
)!22(
)(
2
1
)22(
ωρ
ξ
( )ba，∈ξ
证明：
设为满足
)(
12
xH
n+
=
=
+
+
)()(
)()(
''
12
12
kkn
kkn
xfxH
xfxH
),,1,0( nk L=
的Hermite插值多项式，则的次数。)(
12
xH
n+
12 +≤ n
2004-11-1 82
高斯型求积公式的余项（续）
)(
)!22(
)(
)()(
2
1
22
12
x
n
f
xHxf
n
n
n +
+
+
+
=? ω
η
且( )ba，∈η
由于高斯型求积公式的代数精度为2n+1，故
[]
∫
∑
=
=
b
a
n
k
kk
xfAdxxfxfR
0
)()()(ρ
∫
∑
=
+
=
b
a
n
k
knk
xHAdxxfx
0
12
)()()(ρ
∫∫
+
=
b
a
b
a
n
dxxHxdxxfx )()()()(
12
ρρ
2004-11-1 83
高斯型求积公式的余项（续）
∫
+
+
+
=
b
a
n
n
dxx
n
f
x )(
)!22(
)(
)(
2
1
22
ω
η
ρ
∫
+
+
+
=
b
a
n
n
dxxx
n
f
)()(
)!22(
)(
2
1
22
ωρ
ξ
( )ba,∈ξ
证毕
Remark：高斯型求积公式具有代数精度高、
且总是收敛、稳定的优点。在实用中也可构造复化高斯求积公式。
2004-11-1 84
三.几种特殊的高斯型求积公式
1.高斯—勒让德（Gauss—Legendre）求积公式
]1,1[],[,1)(?== baxρ
公式：
∫
∑
=
≈
1
1
0
)()(
n
k
kk
xfAdxxf
为的零点。
),,1,0( nkx
k
L=
12
1
1
1
1
)1(
2)!1(
1
)(
+
+
+
+
+
+
=
n
n
n
n
n
x
dx
d
n
xP
[]
2
'
1
2
)()1(
2
knk
k
xPx
A
+
= ),,1,0( nk L=
2004-11-1 85
Gauss—Legendre求积公式（续）
[]
[ ]
[]
)(
)!22()32(
)!1(2
22
3
4
32
ξ
+
+
++
+
=
n
n
f
nn
n
fR）（1,1?∈ξ
Remark：当积分区间为时，可通过线性变换
[ ]ba，
22
ab
t
ab
x
+
+
=
[]ba，[ ]11，?
将区间变换为标准化区间，从而有
∫∫
+
+
=
1
1
)
22
(
2
)( dt
ba
t
ab
f
ab
dxxf
b
a
从而可以使用Gauss-Legendre求积公式来计算。
2004-11-1 86
2.高斯—切比雪夫(Gauss-Chebyshev)求积公式
2
1
1
)(
x
x
=ρ
∫
∑
=
≈
1
1
0
2
)()(
1
1
n
k
kk
xfAdxxf
x
高斯点为n+1次切比雪夫多项式的零点.
),1,0( nkx
k
L=
π
22
12
cos
+
+
=
n
k
x
k
1+
=
n
A
k
π
系数
[] )(
)!22(2
)22(
12
ξ
π
+
+
+
=
n
n
f
n
fR
2004-11-1 87
3.高斯—拉盖尔(Gauss-Laugerre)求积公式
[ ] ),0[,,)( +∞==
baex
x
ρ
∫
∑
∞+
=
≈
0
0
)()(
n
k
kk
x
xfAdxxfe
的零点。为)()(
1
1
1
1
xn
n
n
x
nk
ex
dx
d
exLx
+
+
+
+
=
[]
[ ]
)(
)!22(
)!1(
22
2
ξ
+
+
+
=
n
f
n
n
fR )0[ ∞+∈，ξ
[ ]
[]
2
'
1
2
)(
)!1(
knk
k
xLx
n
A
+
+
= ),,1,0( nk L=
2004-11-1 88
4.高斯—埃尔米特(Gauss-Hermite)求积公式
[ ] ),(,,)(
2
+∞?∞==
baex
x
ρ
∑
∫
=
∞
∞+
≈
n
k
kk
x
xfAdxxfe
0
)()(
2
的零点。为)()1(
22
1
1
1
1
x
n
n
xn
nk
e
dx
d
eHx
+
+
+
+
=
[]
2
'
1
2
)(
)!1(2
kn
n
k
xH
n
A
+
+
+
=
π
[] )(
)!22(2
)!1(
)22(
1
ξ
π
+
+
+
+
=
n
n
f
n
b
fR
( )∞+∞?∈，ξ
关于高斯点和求积系数有详细的表格可供参阅。
2004-11-1 89
四．一般情况下的高斯型求积公式的构造例：求高斯型求积公式
)()()(
2211
1
0
xfAxfAdxxfx +≈
∫
的系数及节点。
21
,AA
21
,xx
解：上面的求积公式对函数f(x)=1，
应该是精确成立的。
32
,,xxx
于是有下述公式成立：
2004-11-1 90
续
5
2
1
0
2211
==+
∫
dxxxAxAx
3
2
1
0
21
==+
∫
dxxAA
∫
==+
1
0
2
2
2
21
2
1
7
2
dxxxAxAx
9
2
1
0
3
2
3
21
3
1
==+
∫
dxxxAxAx
2004-11-1 91
续
=
9
2
7
2
1
3
2
3
1
2
2
2
1
2
1
xx
xx
A
A
写成矩阵有：
=
5
2
3
2
11
1
212
1
xxA
A
上面两式右端写成等式，并化简有：
2004-11-1 92
续
+=
+?++=
9
2
)(
7
2
5
2
)(
9
2
)(
7
2
)(
3
2
2121
21
2
221
2
1
2
21
xxxx
xxxxxxxx
=
=
0
9
1
5
9
7
1
0
3
1
7
1
9
1
7
1
22
uv
uuvv
得
vxx =+
21
uxx =
21
令
21
5
=u
9
10
=v
2004-11-1 93
续
389111.0
1
=A
277556.0
2
=A
821162.0
1
=x 289949.0
2
=x
从而有
()
()289949.0277556.0
821162.0389111.0)(
1
0
f
fdxxfx
+
≈
∫
Remark：正交多项式构造的困难，是
Gauss型求积公式的不足之处。
2004-11-1 94
§4.5 数值微分以离散数据近似表达在节点处的微分,通常称为数值微分。
))(,(
kk
xfx
),......2,1,0( nk =
)(xfy =
k
x
一,插值型求导公式给出列表函数，可建立插值多项式，取作为的近似函数，
则称为插值型求导公式。
)(xfy =
)(xp
n
),()( xpxf
n
′
≈
′
)(xf
′
)(' xp
n
2004-11-1 95
插值型求导公式（续）
)(
)!1(
)(
)()(
1
)1(
x
n
f
xpxf
n
n
n +
+
+
+= ω
ξ
由][
n
xx,
0
∈ξ
)(
)!1(
)(
)(
)!1(
)(
)()(
)1(
1
1
)1(
ξ
ω
ω
ξ
+
+
+
+
+
+
′
+
=
′
′
n
n
n
n
n
f
dx
d
n
x
x
n
f
xpxf得若要求节点上的导数值，有余项k
x
)(
)!1(
)(
)()(
1
)1(
kn
n
knk
x
n
f
xpxf
+
+
′
+
=
′
′
ω
ξ
2004-11-1 96
插值型求导公式（续）
若f(x)的各阶导数有界，即
),2,1,0()(
)1(
L=≤
+
nMxf
n
),,2,1,0(],,[ nkbax
k
L=∈
节点
)(0)(
)!1(
)()( ∞→→?
+
≤
′
′
nab
n
M
xpxf
n
kk则误差
),,1,0()()( nkxpxf
knk
L=
′
≈
′
从而有为讨论方便，假定所给节点是等距的，我们限定求节点上的各阶导数值。
2004-11-1 97
1.两点公式
[])()(
1
)(
101
xfxf
h
xp +?=
′
得
01
xxh?=其中由余项
)()()(
1
01
0
0
10
1
1
xf
xx
xx
xf
xx
xx
xp
+
=由
)(
!2
)(
)()(
21 kkk
x
f
xpxf ω
ξ
′
′′
=
′
′
),(
10
xx∈ξ
[]
∈
′′
+?=
′
∈
′′
=
′
),()(
2
)]()([
1
)(
),(
!2
)(
)()(
1
)(
1022011
101
1
010
xxf
h
xfxf
h
xf
xx
fh
xfxf
h
xf
ξξ
ξ
ξ
得
2004-11-1 98
2.三点公式由
)(
))((
))((
)(
))((
))((
)(
))((
))((
)(
2
1202
10
1
2101
20
0
2010
21
2
xf
xxxx
xxxx
xf
xxxx
xxxx
xf
xxxx
xxxx
xp
+
+
=
)(
2
)1(
)(
1
)2(
)(
2
)2)(1(
)(
2
1002
xf
tt
xf
tt
xf
tt
thxp
+
+
=+
hxx 2
02
+= thxx +=
0
hxx +=
01
令
2004-11-1 99
三点公式（续）
由两边对t求导
[]
[][])1(
2
)(
)2()(
)2()1(
2
)(
)(
2
1
0
02
++?+?
+?=+
′
tt
xf
ttxf
tt
xf
thxph
)}()12()()44()()32{(
2
1
)(
21002
xftxftxft
h
thxp?+=+
′
),())((
!3
)(
))((
!3
)(
))((
!3
)(
)(')(
201020
212
xxxxxx
f
xxxx
f
xxxx
f
xpxf
kk
∈
′′′
+
′′′
+
′′′
=?
′
ξ
ξξ
ξ
2004-11-1 100
三点公式（续）
注意在中，当时得
thxx +=
0 210
,,xxx
2,1,0=t
{})2)((
!3
)(
)()(4)(3
2
1
)(
1
2100
hh
f
xfxfxf
h
xf
′′′
+?+?=
′
ξ
)(
3
)()(4)(3(
2
1
1
2
200
ξf
h
xfxfxf
h
′′′
+?+?=
{} ))((
!3
)(
)()(
2
1
)('
2
2010
hh
f
xfxf
h
xf?
′′′
++?=
ξ
{}
2
2
20
6
)(
)()(
2
1
h
f
xfxf
h
ξ
′′′
+?=
2004-11-1 101
三点公式（续）
{}
2
3
2102
2
!3
)(
)(3)(4)(
2
1
)( h
f
xfxfxf
h
xf
ξ
′′′
++?=
′
)(
3
)(3)(4)(
2
1
3
2
210
ξf
h
xfxfxf
h
′′′
++?=
),(,,
20321
xx∈ξξξ
其中类似的，我们还可以建立高阶数值微分公式：
L,2,1)()(
)()(
=≈ kxpxf
k
n
k
2004-11-1 102
三点公式（续）
如可以得到二阶三点公式
[])(
6
)()()(2)(
1
)(
1
)4(
2
1210
2
0
ηξ f
h
fhxfxfxf
h
xf +
′′′
+?=
′′
)(
12
)()(2)(
1
)(
2
)4(
2
210
2
1
ηf
h
xfxfxf
h
xf?+?=
′′
)](
6
)()()(2)([
1
)(
3
)4(
2
2211
2
2
ηξ f
h
fhxfxfxf
h
xf +
′′′
++?=
′′
),(,,,,
2032121
xx∈ηηηξξ
其中
2004-11-1 103
二.Taylor展开法根据Taylor展开式可得
)(
2
)()(
)(
1
ξf
h
h
xfhxf
xf
kk
k
′′
+
=
′
)(
2
)()(
)(
2
ξf
h
h
hxfxf
xf
kk
k
′′
+
=
′
)(
62
)()(
)(
3
2
ξf
h
h
hxfhxf
xf
kk
k
′′′
+
=
′
)(
12
)()(2)(
)(
4
)4(
2
2
ξf
h
h
hxfxfhxf
xf
kkk
k
+?+
=
′′
),(
1
hxx
kk
+∈ξ
),(
2 kk
xhx?∈ξ
),(
3
hxhx
kk
+?=ξ
),(
4
hxhx
kk
+?=ξ
（1）
（2）
（3）
（4）
2004-11-1 104
Taylor展开法（续）
以(3)为例
)
~
(
!3
)(
2
)()()(
1
32
ξf
h
xf
h
xfhxfhxf
kkkk
′′′
+
′′
+
′
+=+
)
~
(
!3
)(
2
)()()(
2
32
ξf
h
xf
h
xfhxfhxf
kkkk
′′′
′′
+
′
=?
),(
~
,
~
2021
xx∈ξξ
)(a
)(b
其中两式相减得
))
~
()
~
((
6
)(2)()(
21
3
ξξ ff
h
xfhhxfhxf
kkk
′′′
+
′′′
+
′
=+
2004-11-1 105
Taylor展开法（续）
[ ],)
~
()
~
(
2
1
)(
213
ξξξ fff
′′′
+
′′′
=
′′′
从而有：
[])(
6
)()(
2
1
)(
3
2
ξf
h
hxfhxf
h
xf
kkk
′′′
+=
′
),(
3
hxhx
kk
+?∈ξ
若连续，则
)(xf
′′′
Mffm ≤
′′′
+
′′′
≤ )
~
()
~
((
2
1
21
ξξ
（为的最小值与最大值）故由连续函数的介值定理，必存在使得
Mm,
)(xf
′′′
,
3
ξ
由(a)可以得到(1)；由(b)可以得到(2)。
2004-11-1 106
Taylor展开法（续）
若将(a)(b)写为
)
~
(
!4
)(
!3
)(
2
)()()(
1
)4(
43
2
ξf
h
xf
h
xf
h
xfhxfhxf
k
kkkk
+
′′′
+
′′
+
′
+=+
)
~
(
!4
)(
!3
)(
2
)()()(
2
)4(
43
2
ξf
h
xf
h
xf
h
xfhxfhxf
k
kkkk
+
′′′
′′
+
′
=?
两式相加可以类似的得到（4）
2004-11-1 107
Remark
Remark1：在数值微分计算中，并非步长越小精度越高。这是因为数值微分对舍入误差非常敏感，
它随步长h的缩小而增大，导致计算不稳定。
Remark2：在数值微分计算中，当插值多项式收敛到函数f(x)时，P′
n
（x）不一定收敛到f ′(x)。
Remark3：为了避免上述问题，可以用样条插值函数的导函数来代替f(x)的导函数。
第四章数值积分与数值微分
§1 数值积分概述
§2 Newton Cotes 公式
§3 Romberg求积法
§4 Gauss型求积公式
§5 数值微分
2004-11-1 2
§4.1 数值积分概述由积分学基本定理Newton-Leibenize公式有：
)()()( aFbFdxxf
b
a
=
∫
但在应用中常会碰到如下情况：
①f(x)的原函数无法用初等函数给出
②f(x)用表格形式给出
③虽然f(x)的原函数能用初等函数表示,
但表达式过于复杂。
2004-11-1 3
∑
∫
=
=≈=
n
k
nkk
b
a
IxfAdxxfI
0
)()(
∫
∑
=
+=+==
b
a
n
k
nkk
fRIfRxfAdxxfI
0
][][)()(
一.求积公式其中R[f]称为求积公式的余项，
),2,1,0( nkx
k
L=
称为求积节点，
),2,1,0( nkA
k
L=
称为求积系数。
k
x
k
A
仅与求积节点的选取有关，而不依赖与被积函数f(x)的具体形式。
2004-11-1 4
二.求积公式的代数精确度衡量一个求积公式好坏的标准。
∫
∑
=
≈
b
a
n
k
kk
xfAdxxf
0
)()(
定义：如果求积公式对于一切不高于m次的代数多项式准确成立，
而对于某个m+1次多项式并不准确成立，
则称上述求积公式具有m次代数精确度，或称为具有m次代数精度。
2004-11-1 5
求积公式的代数精确度（续）
如果要构造具有m次代数精度的求积公式，
m
xxxxf,,,,1)(
2
L=
只要令它对于都能准确成立即可。
∑
∫
=
+++=+++
n
k
m
kmkk
b
a
m
m
xaxaaAdxxaxaa
0
1010
)()( LL
∫∫∫
+++?
b
a
m
m
b
a
b
a
dxxbaxdxadxa L
10
∑∑∑
===
+++=
n
k
m
kkm
n
k
n
k
kkk
xAaxAaAa
000
10
L
2004-11-1 6
求积公式的代数精确度（续）
∫
∑
=
=
b
a
n
k
k
Adx
0
∫
∑
=
=
b
a
n
k
kk
xAxdx
0
∫
∑
=
=
b
a
n
k
m
kk
m
xAdxx
0
M
即若求积公式具有m次代数精度，
k
A应满足上述方程组。反之亦然。
2004-11-1 7
特别地,若取n = m，即取m+1个节点,
m
xx L
0
求积公式的代数精确度（续）
则可通过给定的m+1个节点得到上述含m+1个未知数、m+1个方程的方程组。
若求积节点互异，则
0
111
det
10
10
≠=
m
m
mm
m
xxx
xxx
A
L
MMM
L
L
从而可得如下定理：
2004-11-1 8
求积公式的代数精确度（续）
定理：
使求积公式至少有m次代数精度。
在区间[a,b]上,对于给定m+1个互异节点，
bxxa
m
≤<<≤ L
0 m
AAA L,,
10
，总存在求积系数,
Remark：定出，则求积公式至少具有m次代数精度，但并不一定它具有m次代数精度,要将代入求积公式，如果等式不准确成立，即求积公式具有m 次代数精度，
否则代数精度将大于m。
),,1,0( mkA
k
L=
1+m
x
2004-11-1 9
求积公式的代数精确度（续）
Remark1：代数精度越高，求积公式的适应性越强。
Remark2：凡至少具有零次代数精度的求积公式
∫
∑
=
≈
b
a
n
k
kk
xfAdxxf
0
)()(
一定满足
∫
∑
=
=?
b
a
n
k
k
Adx
0
11
abA
n
k
k
=
∑
=0
从而有即求积系数之和等于积分区间长度，这是求积系数的一个基本特性。
2004-11-1 10
求积公式的代数精确度（续）
例：确定求积公式
)]()0([
2
)]()0([
)(
2
0
hffh
hffh
dxxf
h
′
′
+
+
≈
∫
α
中的待定参数，使其代数精度尽量高，并指明所构造的求积公式所具有的代数精度。
解：求积公式中含一个待定参数当f(x)=1,x时,有
∫
+≡
h
h
dx
0
]11[
2
2004-11-1 11
求积公式的代数精确度（续）
∫
++=
h
hh
h
xdx
0
2
]11[]0[
2
α
故令求积公式对f(x)=x
2
成立,即
∫
×++=
h
hhh
h
dxx
0
222
)202(]0[
2
α
12
1
=α
得
2004-11-1 12
求积公式的代数精确度（续）
3
)( xxf =
将代入已求得的求积公式，显然
]40[
12
]0[
2
3
2
4
0
4
h
h
h
h
dxx
h
++≠
∫
∫
′
′
++≈
h
hff
h
hff
h
dxxf
0
2
)]()0([
12
)]()0([
2
)(
故具有三次代数精度。
]30[
12
]0[
2
2
2
3
0
3
h
h
h
h
dxx
h
++=
∫
令
4
)( xxf =
2004-11-1 13
三.收敛性与稳定性
),0][lim()()(lim
0
0
==
∞→
→
=
∞→
∑
∫
fRdxxfxfA
n
h
n
k
b
a
kk
n
如果如果求积公式对舍入误差不敏感（误差能够控制），则称该求积公式是稳定的。
其中，则称该求积公式是收敛的。
)(max
1
1
≤≤
=
ii
ni
xxh
一个求积公式首先应该是收敛的，其次应该是稳定的。
2004-11-1 14
收敛性与稳定性（续）
则
∑∑
==
=
n
k
n
k
kkkkn
xfAxfAIe
00
*
)()()(
∑∑
==
≤=
n
k
kk
n
k
kk
AA
00
εε
∑
=
=≤
n
k
kn
abAIe
0
)()( εε
则设计算时有绝对误差，即
k
ε)(
k
xf
.)()(
*
kkk
xfxf ε=?
n
nk
εε
≤≤
=
0
max
),,1,0(0 nkA
k
L=>
取若
k
A故当系数全为正时，求积公式是稳定的。
2004-11-1 15
§4.2 Newton Cotes 公式一.插值型求积公式找到一个足够精度的简单函数p(x)代替f(x),于是
∫∫
≈
b
a
b
a
dxxpdxxf )()(
有,把p(x)取成插值多项式，
则可得到插值型求积公式。
设给定节点
bxxxxa
n
≤<<<<≤ L
210
),,1,0()( nkxf
k
L=
并已知这些节点上的函数值
2004-11-1 16
插值型求积公式（续）
由
)(
)!1(
)(
)()()(
1
1
0
x
n
f
xfxlxf
n
n
n
k
kk +
+
=
+
+=
∑
ω
ξ
dxx
n
f
dxxfxldxxf
b
a
b
a
n
n
b
a
n
k
kk
∫∫∫
∑ +
+
=
+
+= )(
)!1(
)(
)()()(
1
)1(
0
ω
ξ
∫
∑ +
+
=
+
+=
b
a
n
n
n
k
kk
dxx
n
f
xfA )(
)!1(
)(
)(
1
)1(
0
ω
ξ
∫
∏
∫
≠
=
==
b
a
n
ki
i
ik
i
b
a
kk
dx
xx
xx
dxxlA
0
)(
)(
)(其中，系数
2004-11-1 17
插值型求积公式（续）
∫
=
b
a
kk
dxxlA )(
当求积系数由截断误差
∫
∏
=
+
+
=
b
a
n
i
i
n
dxxxf
n
fR
0
)1(
)()(
)!1(
1
][ ξ
所唯一确定时，所得的求积公式称为插值型求积公式。
Remark：由截断误差可知，插值型求积公式至少具有n次代数精度。
2004-11-1 18
二,Newton-cotes公式将[a，b]分为n等份，
取节点（k＝0,1,…,n）
khax
k
+=
nabh )(?=
由Lagrange插值公式，可得
∫∫
+
+
′
==
b
a
b
a
knk
n
kk
dx
xxx
x
dxxlA
)()(
)(
)(
1
1
ω
ω
系数可以进一步表示：
k
A
2004-11-1 19
Newton-cotes公式（续）
),,1,0( nkkhax
k
LQ =+=
令x=a+th即有dx=hdt
],0[],[ nba?
hktxx
k
)(?=?
)()1()()()(
1
01
nttthxxxxx
n
nn
==
+
+
LLω
故
)())(()()(
1101 nkkkkkkkn
xxxxxxxxx=
′
+?+
LLω
))(()2)(1(1)1( knkkh
n
= LL
)!()1(! knkh
knn
=
2004-11-1 20
Newton-cotes公式（续）
故
∫
=
+
n
knn
n
k
hdt
knkhhkt
nttth
A
0
1
)!()1(!)(
)()1( L
∫
+
=
n
kn
dtntktkttt
knk
h
0
)()]1()][1([)1(
)!(!
)1(
LL
∫
+
=
n
kn
dtntktkttt
knnk
ab
0
)()1)(1()1(
)!(!
)1(
)( LL
)(
)(
n
k
cab?=
2004-11-1 21
Newton-cotes公式（续）
故求积公式可写为,其
∫
∑
=
+?≈
b
a
n
k
n
k
khxfcabdxxf
0
0
)(
)()()(
)(n
k
c
中称为柯特斯系数，上式称为Newton----Cotes公式。
可以证明。
)()( n
kn
n
k
cc
=
2004-11-1 22
Newton-cotes公式（续）
当n=1时，
2
1
2
)1(
)1()1(
!0!11
)1(
1
0
1
0
201
)1(
0
=
×?=?
=
∫
t
dttc
∫
==
=
1
0
1
0
2
11
)1(
1
2
1
2
1
!0!11
)1(
ttdtc
∫
+
≈
b
a
bfaf
ab
dxxf )]()([
2
)(
)(
上式称为梯形公式。
2004-11-1 23
Newton-cotes公式（续）
n=2可计算
6
1
6
4
6
1
)2(
2
)2(
1
)2(
0
=== ccc
)](
6
1
)
2
(
6
4
)(
6
1
)[()( bf
ba
fafabdxxf
b
a
+
+
+?≈
∫
故有它称为辛浦生（Simpson）公式或抛物线公式
n=4 Newton—Cotes公式为
)](
90
7
)(
90
32
)(
90
12
)(
90
32
)(
90
7
[
432
10
xfxfxf
xfxfabdxxf
b
a
+++
+?≈
∫
）（）（
其中，
)4,,1,0(
0
L=+= kkhax
k
这个公式特别称为柯特斯公式。
2004-11-1 24
Newton-cotes公式（续）
Remark：对n＋1个节点的Newton—Cotes公式，当n
为奇数时，Newton—Cotes公式的代数精度至少为n次；
当n为偶数时，Newton—Cotes公式的代数精度至少为
n＋1次。
只要验证，当n为偶数时，Newton—Cotes公式对
∑
+
=
+
=
1
0
1
)(
n
k
k
kn
xaxP
准确成立即可。
1
)1(
1
)!1()(
+
+
+
+=
n
n
n
anxP
由
∫
+
+
+
+
=?=
b
a
n
n
n
nn
dxx
n
p
IIPR )(
)!1(
)(
][
1
)1(
1
1
ω
ξ
＋
2004-11-1 25
Newton-cotes公式（续）
∫
∏
=
++
=
b
a
n
j
jnn
dxxxaPR
0
11
][）（即
∫
=
+
++
n
n
nn
dsnssshaPR
0
2
11
)(1][ LL）（
令x= a+sh，并注意
jhax
j
+=则，
当n为偶数，即n=2k，k为正整数，再令s=u+k（u=s-k）
∫∫
+=
nn
dsksksksksssdsnsss
00
)2)(12()1)(()1()(1 LLL）（
∫
+++
k
k
dukukuuukuku ))(1()1()1)(( LL＝
2004-11-1 26
Newton-cotes公式（续）
))(1()1()1)(()( kukuuukukuu?+++= LL?记
))(1()1)(()1)(()( kukuuukukuu++?+?=? LL?
))(1()1()1)((1
12
kukuuukuku
k
+?+++=
+
）（LL
)()(1
12
uu
k
=?=
+
）（
从而，为奇函数。
)(u?
故，即当n为偶数时，Newton—Cotes
公式至少有n+1次代数精度。
0][
1
=
+n
PR
2004-11-1 27
三.几种低阶求积公式的截断误差
1.若f(x)在[a,b]上有二阶连续导数，梯形求积公式有下列误差估计：
baf
ab
bfaf
ab
dxxffR
b
a
≤≤
=
+
=
∫
ηη）（
）（
）（）（）（
''
12
][
2
][
3
证：
,))((
!2
)(''
][ dxbxax
f
fR
b
a
=
∫
ξ
由
x依赖于ξ
2004-11-1 28
几种低阶求积公式的截断误差（续）
)(ξf
′′
注意到是倚赖于x的函数，且在[a,b]上连续，故运用积分中值定理，
在[a,b]上存在一点使得：
,0))(( ≤ bxax
,η
)())((''
6
1
))(()(''))()((''
3
baabf
dxbxaxfdxbxaxf
b
a
b
a
≤≤=
=
∫∫
ηη
ηξ
3
)(
12
)(''
][ ab
f
fR=∴
η
2004-11-1 29
几种低阶求积公式的截断误差（续）
2.若f(x)在[a，b]上有四阶连续导数，
辛浦生公式有下列误差估计：
）（
）（
）（）（
）（）（）（
）（
）（
）（
η
η
4
5
)4(4
2880
2180
]
2
4[
6
][
f
ab
f
abab
bf
ba
faf
ab
dxxffR
b
a
=
=
+
+
+
=
∫
ba ≤≤η
2004-11-1 30
几种低阶求积公式的截断误差（续）
）（）（bfbHafaH == )()(
33
)
2
()
2
(
22
33
ba
f
ba
H
ba
f
ba
H
+
′
=
+
′
+
=
+
）（）（
)]()
2
(4)([
6
)(
3333
bH
ba
HaH
ab
dxxH
b
a
+
+
+
=
∫
由于
)()
2
)((
!4
)(
)()(
2
)4(
3
bx
ba
xax
f
xHxf?
+
+=
ξ
3≤
由于辛浦生公式的代数精度为3，为此构造次数的多项式，使满足：)(
3
xH
证：
2004-11-1 31
几种低阶求积公式的截断误差（续）
)]()
2
(4)([
6
)(
)]()
2
(4)([
6
)(][
333
bH
ba
HaH
ab
dxxf
bf
ba
faf
ab
dxxffR
b
a
b
a
+
+
+
=
+
+
+
=
∫
∫
∫
=
b
a
dxxHxf )]()([
3
dxbx
ba
xax
f
b
a
)()
2
)((
!4
)(
2
)4(
+
=
∫
ξ
2004-11-1 32
几种低阶求积公式的截断误差（续）
dxbxcxaxf
dxcxbxaxf
b
a
b
a
)())(()(
))()()((
2)4(
2)4(
=
∫
∫
η
ξ
4)4(
)
2
)()((
180
1
][
ab
abffR
= η
)())((
2880
1
5)4(
baabf ≤≤= ηη
η
由于是依赖于x的函数，在[a,b]上连续，
故可运用积分中值定理，
在[a,b]上存在一点，使
,0)())((
2
≤ bxcxax
)(
)4(
ξf
2004-11-1 33
几种低阶求积公式的截断误差（续）
3.若f(x)在[a，b]上有六阶连续导数，
Cotes公式有下列误差估计：
）（
）（
）（
η
66
)
4
(
945
2
)(][ f
abab
nCIfR
=?=
ba ≤≤η
2004-11-1 34
四.复化求积公式当n 7时，Cotes系数均为正，但从n=8开始，Cotes
系数有正有负，这会使计算误差得不到控制，稳定性得不到保证。另外，节点增多的高次插值多项式舍入误差很大，导致插值型求积公式的舍入误差增大。
≤
因此，实际计算时，一般不采用n较大的Newton-
Cotes公式，而是将区间[a,b]等分为n个小区间，其长度为h=(b-a)/n,在每个小区间上应用低阶的公式，
然后对所有小区间上的计算结果求和，这样得出的
Newton-Cotes公式称为复化求积公式。
2004-11-1 35
1.复化梯形公式
∑∑
=
=
+
′′
+=
1
0
1
0
3
1
)(
12
)]()([
2
n
k
n
k
kkk
f
h
xfxf
h
η
∑
∫∫
=
+
==
1
0
1
)()(
n
k
x
x
b
a
k
k
dxxfdxxfI
由将[a,b]等分为n个子区间
.1,,1,0],[
1
=
+
nkxx
kk
L
][)(
)(
fRnT
nT
+=
),(
1+
∈
kkk
xxη
2004-11-1 36
复化梯形公式（续）
∑
=
+
+=
1
0
1
))()((
2
)(
n
k
kk
xfxf
h
nT
∑
=
++=
1
1
)]()(2)([
2
n
k
k
bfxfaf
h
∑
=
′′
=?=
1
0
3
)(
)(
12
)(][
n
k
knT
f
h
nTIfR η
其中余项若f″(x)在[a,b]上连续,设m为f ″(x)的最小值,M
为f ″(x)的最大值,则
2004-11-1 37
复化梯形公式（续）
故由介值定理,一定存在(a,b)上一点使
η
.)("
1
1
0
Mf
n
m
n
k
k
≤≤
∑
=
η
∑
=
=
1
0
)("
1
)("
n
k
k
f
n
f ηη
)(
12
)("
12
][
2
1
0
3
)(
ηη fh
nh
f
h
fR
n
k
knT
′′
=?=
∑
=
故
),( ba∈η)("
12
)(
2
ηfh
ab?
=
2004-11-1 38
2.复化Simpson公式
n
ab
h
=
(将[a,b]n等分(偶数份),子区间长度).
∑
∫∫
=
+
==
1
2
0
)1(2
2
)()(
n
k
x
x
b
a
k
k
dxxfdxxfI由
∑∑
=
+
=
+
++=
1
2
0
)4(
5
22
1
2
0
122
)(
2880
)2(
)]()(4)([
6
2
n
k
kk
n
k
kk
f
h
xfxfxf
h
η
L++++++= )()(4)()()(4)([
6
2
43221
xfxfxfxfxfaf
h
∑
=
+++
1
2
0
)4(
5
12
)(
2880
)2(
)]()(4)(
n
k
knn
f
h
bfxfxf η
2004-11-1 39
复化Simpson公式（续）
∑∑∑
=
==
+++=
1
2
0
)4(
5
1
2
1
2
2
1
12
)(
2880
)2(
)]()(2)(4)([
6
2
n
k
k
n
k
k
n
k
k
f
h
bfxfxfaf
h
η
][)(
)(
fRnS
nS
+=
其中称为
)]()(2)(4)([
3
)(
1
2
1
2
2
1
12
bfxfxfaf
h
nS
n
k
k
n
k
k
+++=
∑∑
==
复化辛浦生（Simpson）公式。
∑
=
=?=
1
2
0
)4(
5
)(
)(
90
)(][
n
k
knS
f
h
nSIfR η
余项
2004-11-1 40
复化Simpson公式（续）
M
n
m
n
k
k
f
≤≤
∑
=
1
2
0
)4(
)(
2
η
故由介值定理,一定在(a,b)中有一使
η
∑
=
=
1
2
0
)4()4(
)(
2
)(
n
k
k
ff
n
ηη
设在在[a,b]上连续,设m为的最小值,M
为的最大值,则
)(
)4(
xf )(
)4(
xf
)(
)4(
xf
)(
180
)(
)(
290
)(
90
][
)4(4
)4(
5
1
2
0
)4(
5
)(
η
ηη
fh
ab
f
nh
f
h
fR
n
k
knS
=
×?=?=
∑
=
故：
),( ba∈η
2004-11-1 41
3.复化Cotes公式
n
ab
h
=(将[a,b]n等分(n为4的倍数),子区间长度).
)(
945
)(2
)](7)(14)(32
)(12)(32)(7[
90
4
)(
)6(6
1
4
1
4
1
4
0
34
1
4
0
24
1
4
0
14
ηfh
ab
bfxfxf
xfxfaf
h
dxxf
n
k
k
n
k
k
n
k
k
n
k
k
b
a
+++
++≈
∑∑
∑∑
∫
=
=
+
=
+
=
+
Remark:复化梯形公式，复化辛浦生(Simpson)
公式和复化柯特斯(Cotes)公式当时均收敛到所求积分值I。
0→h
2004-11-1 42
例1
若取9个节点，用复化梯形公式、复化辛浦生公式和复化柯特斯公式计算积分，其步长以及与9
个节点所对应的求积系数分别是多少？
解：复化梯形公式：n=8，h=(b-a)/8，对应的求积系数为1、2、2、2、2、2、2、2、1。
复化辛浦生公式：n=8，h=(b-a)/8，对应的求积系数为1、4、2、4、2、4、2、4、1。
复化柯特斯公式：n=8，h=(b-a)/8，对应的求积系数为7、32、12、32、14、32、12、32、7。
2004-11-1 43
用积分计算ln2,要使所得结果的近似值具有5位有效数字。问用复化梯形公式，复化Simpson公式时，至少要取多少个节点?
,2ln2
1
8
2
=
∫
dx
x
例2
解：由且
∫
=
8
2
1
2
1
2ln dx
x
375.0
8
3
8
1
2
11
2
1
8
2
8
2
==>
∫∫
dxdx
x
故计算ln2时，要使误差不超过也即计算2ln2，其误差不超过。
,10
2
1
5?
×
5
10
eln2ln375.0 <<
即
2004-11-1 44
例2（续）
其中
)(
max
2
xfM
bxa
′′
=
≤≤
2
1
)('
x
xf =
3
2
)("
x
xf =
4
1
2
22
max
33
82
2
===
≤≤
x
M
x
2
2
3
2
2
)(
12
)(
12
)(
][1 M
n
ab
Mh
ab
fR
nT
=
≤）由
5
2
3
2
2
3
)(
10
4
1
12
6
12
)(
][
≤×=
≤
n
M
n
ab
fR
nT
令
（
....8203933.67010
412
6
5
3
=×
×
≥n
故区间应取671个，节点至少应取672。
2004-11-1 45
例2（续）
4
4
5
4
4)4(4
)(
180
)(
180
)(
180
][)2( M
n
ab
Mh
ab
fh
ab
fR
nS
=
≤
= η由
)(max
)4(
4
xfM
bxa ≤≤
=
其中由
5
)4(
4
)3(
32
24
)(,
6
)(,
2
)(",
1
)('
x
xf
x
xf
x
xf
x
xf =
==
=
4
3
2
2424
max
55
82
4
=≤=
≤≤
x
M
x
2004-11-1 46
例2（续）
5
4
5
4
4
5
)(
10
4
3
180
6
180
)(
][
≤×=×
=
nn
Mab
fR
nS
令
42640687.4260707106781.0
10625.010
4180
1036
}10
4180
36
{
4
4
5
4
1
5
5
=×=
××=×
×
××
=×
×
×
≥n
故区间n应取44（偶数份），即45个节点。
2004-11-1 47
例3
下面我们来看一个数值算例：
1
4
1
0
2
=
+
∫
dx
x
求解
2004-11-1 48
五.区间逐次分半求积法
1.误差的事后估计法复化求积公式是提高精度的一种有效方法，但在使用复化求积公式之前，必须根据复化求积公式的余项进行先验估计，以确定节点数目，从而确定合适的等分步长。
因为余项表达式中包含了被积函数的导数，而估计各阶导数的最大值往往是很困难的，且估计的误差上界往往偏大。所以实际中，常常使用“事后估计误差”的方法，
通过区间的逐次分半，在步长逐次分半的过程中，反复利用复化求积公式进行计算，查看相继两次计算结果的差值是否达到要求，直到所求得的积分值满足精度要求。
2004-11-1 49
误差的事后估计法（续）
对梯形公式，假定区间分为n等份时，由公式
∫
∑
=++≈=
=
b
a
n
k
k
nTbfxfaf
h
xfI )()]()(2)([
2
)(
1
1
其中算出的积分近似值为T(n),因而有
.khax
k
+=
)("
12
)(
1
2
ηfh
ab
nTI
=
再把各个小区间分为二等份，得积分的近似值为T(2n)，则积分值为
)(")
2
(
12
)2(
2
2
ηf
hab
nTI
=
),(
1
ba∈η
),(
2
ba∈η
2004-11-1 50
误差的事后估计法（续）
22
)
2
(
12
)(
)2(
12
)(
)(
hab
nTI
h
ab
nTI
≈
假定f′′(x)有[a,b]上变化不大，即有，则
)(")("
21
ηη ff ≈
4
)2(
)(
≈
nTI
nTI
即上式可改写为
))()2((
14
1
)2(
))()2((
3
1
)2(
nTnTnT
nTnTnTI
+=
+≈
2004-11-1 51
误差的事后估计法（续）
上式说明用T(2n)做为积分Ｉ的近似值时，其误差大约为
)]()2([
3
1
nTnT?
计算时只需检验是否满足？若不满足，则再把区间分半进行计算，直到满足要求为止。这样就实现了步长得自动选取。
ε<? )()2( nTnT
2004-11-1 52
误差的事后估计法（续）
类似的，还可以得到下面的结论：
对于Simpson公式，假定在在[a,b]上变化不大，则有
)(
)4(
xf
)(
180
)(
)(
1
)4(4
ηfh
ab
nSI
=由),(
1
ba∈η
)()
2
(
180
)(
)2(
2
)4(4
ηf
hab
nSI
= ),(
2
ba∈η
))()2((
14
1
)2())()2((
15
1
)2(
2
nSnSnSnSnSnSI?
+=?+≈
24
42
)2(
)(
=≈
nSI
nSI
变形即得
2004-11-1 53
误差的事后估计法（续）
对于Cotes公式，假定在[a,b]上变化不大，则由
)(
)6(
x
f
))()2((
14
1
)2(
))()2((
63
1
)2(
3
nCnCnC
nCnCnCI
+=
+≈
)(
945
)(2
)(
1
)6(6
ηfh
ab
nCI
=
)()
2
(
945
)(2
)2(
2
)6(6
ηf
hab
nCI
=
36
42
)2(
)(
=≈
nCI
nCI
),(
1
ba∈η
),(
2
ba∈η
得变形即得
2004-11-1 54
2.区间逐次分半的梯形公式
)]()(2)([
2
)(
1
1
bfxfaf
h
nT
n
k
++=
∑
=
由于
)]()(2)([
2
'
)2(
12
1
bfxfaf
h
nT
n
k
k
+
′
+=
∑
=
2
h
h =
′
∑∑
=
=
+++=
n
k
k
n
k
k
xfxfbfaf
h
nT
1
2
1
1
1
)](2)(2)()([
4
)2(故
∑
=
+=
n
k
k
xf
h
nT
1
2
1
)(
2
)(
2
1
∑
=
+
+=
n
k
n
ab
kaf
n
ab
nTnT
1
)
2
)12((
2
)(
2
1
)2(即
2004-11-1 55
区间逐次分半的梯形公式（续）
这就是复化梯形公式当区间逐次分半时的计算公式。由T(n)计算T(2n)时只需计算新增加节点处的函数值。
==
+
+=
+
=
=
∑
L,3,2,1,2
))12(
2
(
2
)(
2
1
)2(
)]()([
2
)1(
1
1
kn
i
n
ab
af
n
ab
nTnT
bfaf
ab
T
k
n
i
算例求解故有递推公式
2004-11-1 56
§4.3 Romberg求积算法一.对低精度公式经过组合构造高精度公式用T(2n)做为积分Ｉ的近似值时，其误差大约为
)]()2([
3
1
nTnT?
若我们不是把误差舍掉，而是把误差补偿到
T(2n)后，我们会得到精度更好的求积公式－
Simpson公式：
)(
14
1
)2(
14
4
)(
3
1
)2(
3
4
)2( nTnTnTnTnS
=?=
2004-11-1 57
对低精度公式经过组合构造高精度公式（续）
)]()(4)(2
)(2)(4)(2)(4)([
6
2
2222
43210
nnn
xfxfxf
xfxfxfxfxf
h
+++
+++++=
L
事实上，
∑
=
++
++=
1
0
22122
)]()(4)([
6
2
)2(
n
k
kkk
xfxfxf
h
nS
)](2)(4)(4)(4)(4)(2[
6
2
2123210 nn
xfxfxfxfxfxf
h
++++++=
L
)]()(2)(2)(2)([
6
2
222420 nn
xfxfxfxfxf
h
+++++?
L
2004-11-1 58
对低精度公式经过组合构造高精度公式（续）
])(2)()([
6
2
])(2)()([
6
4
1
1
2
12
1
∑∑
=
=
++?++=
n
k
k
n
k
k
xfbfaf
h
xfbfaf
h
)(
3
1
)2(
3
4
nTnT?=
14
)()2(4
)(
3
1
)2(
3
4
)2(
=?=
nTnT
nTnTnS即
（1）
Reamrk：两个二阶精度的梯形公式组合起来就得到一个四阶精度的Simpson公式，截断误差提高了二阶，
起到了加速收敛的作用。
2004-11-1 59
对低精度公式经过组合构造高精度公式（续）
)(
14
1
)2(
14
4
))()2((
15
1
)2()2(
22
2
nSnS
nSnSnSnC
=
+=
类似地，可以证明：
（2）
如果再用柯特斯公式作线性组合，
)(
14
1
)2(
14
4
)]()2([
63
1
)2()2(
33
3
nCnC
nCnCnCnR
=
+=
（3）
这个公式称为Romberg公式。
2004-11-1 60
对低精度公式经过组合构造高精度公式（续）
由（1），（2），（3）组成的方法称为Romberg
算法。
序列{T(n)},{S(n)},{C(n)}和{R(n)}分别称为梯形序列，
Simpson序列，Cotes序列和Romberg序列。
上述用若干个积分近似值推算出更为精确的积分近似值的方法，称为外推算法。得到Romberg序列后还可以继续外推，得到新的求积序列，称为Richardson外推算法。但由于在新的求积序列中，其线性组合的系数分别为：
2004-11-1 61
对低精度公式经过组合构造高精度公式（续）
0
14
1
1
14
4
≈
≈
mm
m
因此，新的求积序列与前一个序列结果相差不大，故通常外推到Romberg序列为止。
可以证明，梯形序列，Simpson序列，Cotes序列和
Romberg序列均收敛到积分值，且每次外推可使误差阶提高二阶。
2004-11-1 62
二.Romberg算法的实现
T数表：
R（32）C（32）S（32）T（32）32
R（16）C（16）S（16）T（16）16
R（8）C（8）S（8）T（8）8
C（4）S（4）T（4）4
S（2）T（2）2
T（1）1
R（n）C（n）S（n）T（n）
区间等分数n=2
k
MMMMM
2004-11-1 63
Romberg算法的实现（续）
对上面的T数表作计算，一直到Romberg序列中前后两项之差的绝对值不超过给定的误差限为止。
Remark:Romberg求积算法从最简单的梯形序列开始逐步进行线性加速，具有占用内存少，
精确度高的优点，是实际中最常用的算法之一。
算例求解
2004-11-1 64
§4.4 Gauss型求积公式问题：
若求积公式
∑
∫
=
≈=
n
k
kk
b
a
xfAdxxfI
0
)()(
),,2,1,0(
,
nkAx
kk
L=
中含有2n+2个待定参数我们能否通过节点的选择将求积公式的代数精度从n或者n+1提高到2n+1？
2004-11-1 65
一．带权积分的插值型求积公式
0)( ≥xρ
dxxxfI
b
a
))(（
∫
= ρ
考虑更一般的带权积分,其中为[a,b]上的给定的权函数。
类似于前面插值型求积公式的建立过程，有
)(
)!1(
)(
)()()(
1
)1(
0
x
n
f
xfxlxf
n
n
n
k
kk +
+
=
+
+=
∑
ω
ξ
)
~
,
~
( Mm∈ξ
由
)()(
)!1(
)(
)()}()({)()(
1
)1(
0
xx
n
f
xxfxlxfx
n
n
n
k
kk
ρω
ξ
ρρ
+
+
=
+
+=
∑
得
2004-11-1 66
带权积分的插值型求积公式（续）
积分
)()()()()(
0
k
n
k
b
a
k
b
a
xfdxxxldxxfx
∑
∫∫
=
= ρρ
∫
+
+
+
+
b
a
n
n
dxx
n
f
x )(
)!1(
)(
)(
1
)1(
ω
ξ
ρ
[]
∑
=
+=
n
k
kk
fRxfA
0
)(
[]
+
=
=
∫
∫
+
+
b
a
n
n
b
a
kk
dxx
n
f
xfR
dxxlxA
)(
)!1(
)(
)(
)()(
1
)1(
ω
ξ
ρ
ρ
其中上述公式称为计算的插值型求积公式。
∫
=
b
a
dxxfxI )()(ρ
2004-11-1 67
二．高斯公式和高斯点我们知道，n＋1个插值节点的插值型求积公式的代数精确度不可能低于n。由于在带权插值求积公式中未知的节点及系数共2n＋2个参数。因此，当节点数n＋1确定时，只要适当选取节点位置x
k
及系数A
k
，就可以使上面的插值型求积公式对f(x)=1,x,…,x
2n+1
精确成立，即使其代数精度提高到2n+1。若某个求积公式（n＋1个插值节点）的代数精度达到2n+1，则称该公式为Gauss
型求积公式。
2004-11-1 68
1.定义若是上的一组互异节点，且求积公式
[ ]ba，
n
xxx L,
1,0
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
具有2n+1次代数精度，则称该求积公式为guass
型求积公式，其求积节点（k=0，1，……n）
称为高斯点，系数称为高斯系数。
k
x
k
A
2004-11-1 69
2.引理求积公式至少具有
n次代数精度的充要条件是它是插值型的。
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
证明：充分性：如果求积公式为插值型的，
由
[]
∫
+
+
+
=
b
a
n
n
dxxx
n
f
fR )()(
)!1(
)(
1
)1(
ωρ
ξ
可知，对于任意次数≤n的多项式f(x),R[f]=0。
这时求积公式至少具有n次精度。
必要性：
设求积公式具有n
次代数精度。
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
2004-11-1 70
续则对于不高于n次的多项式，等式精确成立。
),,1,0()( njxl
j
L=对n次插值函数仍有，
特别地，由于为常数，
k
A
∫
∑
=
=
b
a
n
k
kjkj
xlAdxxlx
0
)()()(ρ
),,1,0()()( njdxxlxA
b
a
jj
L==
∫
ρ
根据插值型求积公式定义知，其求积公式为插值型求积公式。
2004-11-1 71
两条结论：
①．高斯型求积公式一定是插值型求积公式，其系数由高斯点唯一确定。
②．高斯型求积公式是精度最高的求积公式。
2004-11-1 72
续事实上，在中，取
∫
∑
=
≈
b
a
n
k
kk
xfAdxxfx
0
)()()(ρ
,()()(
0
22
1 ∏
=
+
==
n
i
in
xxxxf）ω
，则0)()(
2
1
>
∫
+
b
a
n
dxxx ωρ
.0)()(
00
2
1∑∑
==
+
==
n
k
n
k
knkkk
xAxfA ω
，即对)()(
2
1
xxf
n+
=ω
∑
∫
=
++
≠
n
k
knk
b
a
n
xAdxxx
0
2
1
2
1
)()()( ωωρ
也即插值型求积公式的精度不可能达到2n+2。
2004-11-1 73
3.高斯点的选取定理：插值型求积公式中的节点是高斯点的充要条件是，在[a,b]上，以这些点为零点的n+1次多项式与任意次数不超过n的多项式P(x)带权正交，即
),,1,0( nkx
k
L=
∏
=
+
=
n
j
jn
xxx
0
1
)()(ω
)(xρ
∫
=
+
b
a
n
dxxPxx 0)()()(
1
ωρ
2004-11-1 74
定理证明必要性：
设是高斯点，于是对任意次数不超过n
的多项式P(x)，
),,1,0( nkx
k
L=
的次数不超过2n+1。
)()()(
1
xxPxf
n+
= ω
故有
∫
∑
=
++
==
b
a
n
k
kknkn
xPxAdxxPxx
0
11
0)()()()()( ωωρ
充分性：
设对于任意次数不超过
2n+1的多项式设除f(x)的商为P(x),余项为q(x)。
∫
=
+
b
a
n
dxxPxx 0)()()(
1
ωρ
),(xf
)(
1
x
n+
ω
2004-11-1 75
定理证明（续）
)()()()(
1
xqxxPxf
n
+=
+
ω即
nxqxP ≤的次数其中)(),(
∫∫ ∫
+=
+
b
a
b
a
b
a
n
dxxqxdxxxPxdxxxf )()()()()()()(
1
ρωρρ
∫
=
+
b
a
n
dxxxPx，由条件0)()()(
1
ωρ
所给的求积公式是插值型的，其代数精度至少为n。
∫
∑
=
=
b
a
n
n
kk
xqAdxxqx
0
)()()(ρ故
2004-11-1 76
定理证明（续）
Remark：
)()()()()(
1 kkknkk
xqxqxxPxf =+=
+
ω
∑∑
∫∫
==
==
=
n
k
n
k
kkkk
b
a
b
a
xfAxqA
dxxqxdxxxf
00
)()(
)()()()( ρρ
即求积公式具有2n+1次代数精度，
从而为Gauss点。
),,1,0( nkx
k
L=
2004-11-1 77
Reamrk：
当为正交多项式系中的n+1次多项式时，
)(
1
xf
n+
()()0)(),(
1
)(,
1
)1(
1
1
==
+
+
+
xPxf
A
xP
n
n
n
ω
)(
1
x
n+
ω
)1(
1
1
1
)(
)(
+
+
+
=
n
n
n
A
xf
xω
取，则有n+1个互异的零点，且对任意次数不超过n的多项式有
①[a,b]上带权正交的n+1次多项式的零点就是高斯型求积公式的高斯点。
)(xρ
2004-11-1 78
Reamrk（续）
②当高斯点确定以后，高斯系数
),,1,0( nkA
k
L=
也可以由插值型求积公式中的系数公式确定。
∫
=
b
a
kk
dxxlxA )()(ρ
=
=
=
∫
∑
∫
∑
∫
∑
=
=
=
b
a
n
k
n
kk
n
b
a
n
k
kk
b
a
n
k
k
xAdxxx
xAxdxx
Adxx
0
0
0
)(
)(
)(
ρ
ρ
ρ
M
确定。即可由线性方程组
2004-11-1 79
4.高斯型求积公式的收敛性和稳定性定理：高斯型求积公式总是稳定的。
证明：
只需证明Gauss系数全为正即可。事实上，由于高斯型求积公式对次数不超过2n+1的多项式准确成立，故可取（j＝0,1,…,n），
其中l
j
(x)是n次Lagrange插值基函数，故有
)()(
2
xlxf
j
=
∫
∑
=
>==
b
a
n
k
jkjkj
AxlAdxxlx
0
22
0)()()(ρ
（j=0,1,…n）
故高斯求积系数全为正。因此高斯求积公式是稳定的。
k
A
2004-11-1 80
高斯型求积公式的收敛性和稳定性（续）
Remark：当求积节点数目增加时，Newton-
Cotes求积公式的系数变得有正有负（n≥8），
从而不能保证计算的稳定性。但对于Gauss型求积公式，不论节点数n＋1有多大，Gauss型求积公式总是稳定的。
定理：设f(x)∈C[a,b]，则Gauss型求积公式是收敛的。
2004-11-1 81
5.高斯型求积公式的余项定理：
设在内有直到2n+2阶导数，则高斯型求积公式的余项为：
)(xf
[ ]ba,
[]
∫
+
+
+
=
b
a
n
n
dxxx
n
f
fR )()(
)!22(
)(
2
1
)22(
ωρ
ξ
( )ba，∈ξ
证明：
设为满足
)(
12
xH
n+
=
=
+
+
)()(
)()(
''
12
12
kkn
kkn
xfxH
xfxH
),,1,0( nk L=
的Hermite插值多项式，则的次数。)(
12
xH
n+
12 +≤ n
2004-11-1 82
高斯型求积公式的余项（续）
)(
)!22(
)(
)()(
2
1
22
12
x
n
f
xHxf
n
n
n +
+
+
+
=? ω
η
且( )ba，∈η
由于高斯型求积公式的代数精度为2n+1，故
[]
∫
∑
=
=
b
a
n
k
kk
xfAdxxfxfR
0
)()()(ρ
∫
∑
=
+
=
b
a
n
k
knk
xHAdxxfx
0
12
)()()(ρ
∫∫
+
=
b
a
b
a
n
dxxHxdxxfx )()()()(
12
ρρ
2004-11-1 83
高斯型求积公式的余项（续）
∫
+
+
+
=
b
a
n
n
dxx
n
f
x )(
)!22(
)(
)(
2
1
22
ω
η
ρ
∫
+
+
+
=
b
a
n
n
dxxx
n
f
)()(
)!22(
)(
2
1
22
ωρ
ξ
( )ba,∈ξ
证毕
Remark：高斯型求积公式具有代数精度高、
且总是收敛、稳定的优点。在实用中也可构造复化高斯求积公式。
2004-11-1 84
三.几种特殊的高斯型求积公式
1.高斯—勒让德（Gauss—Legendre）求积公式
]1,1[],[,1)(?== baxρ
公式：
∫
∑
=
≈
1
1
0
)()(
n
k
kk
xfAdxxf
为的零点。
),,1,0( nkx
k
L=
12
1
1
1
1
)1(
2)!1(
1
)(
+
+
+
+
+
+
=
n
n
n
n
n
x
dx
d
n
xP
[]
2
'
1
2
)()1(
2
knk
k
xPx
A
+
= ),,1,0( nk L=
2004-11-1 85
Gauss—Legendre求积公式（续）
[]
[ ]
[]
)(
)!22()32(
)!1(2
22
3
4
32
ξ
+
+
++
+
=
n
n
f
nn
n
fR）（1,1?∈ξ
Remark：当积分区间为时，可通过线性变换
[ ]ba，
22
ab
t
ab
x
+
+
=
[]ba，[ ]11，?
将区间变换为标准化区间，从而有
∫∫
+
+
=
1
1
)
22
(
2
)( dt
ba
t
ab
f
ab
dxxf
b
a
从而可以使用Gauss-Legendre求积公式来计算。
2004-11-1 86
2.高斯—切比雪夫(Gauss-Chebyshev)求积公式
2
1
1
)(
x
x
=ρ
∫
∑
=
≈
1
1
0
2
)()(
1
1
n
k
kk
xfAdxxf
x
高斯点为n+1次切比雪夫多项式的零点.
),1,0( nkx
k
L=
π
22
12
cos
+
+
=
n
k
x
k
1+
=
n
A
k
π
系数
[] )(
)!22(2
)22(
12
ξ
π
+
+
+
=
n
n
f
n
fR
2004-11-1 87
3.高斯—拉盖尔(Gauss-Laugerre)求积公式
[ ] ),0[,,)( +∞==
baex
x
ρ
∫
∑
∞+
=
≈
0
0
)()(
n
k
kk
x
xfAdxxfe
的零点。为)()(
1
1
1
1
xn
n
n
x
nk
ex
dx
d
exLx
+
+
+
+
=
[]
[ ]
)(
)!22(
)!1(
22
2
ξ
+
+
+
=
n
f
n
n
fR )0[ ∞+∈，ξ
[ ]
[]
2
'
1
2
)(
)!1(
knk
k
xLx
n
A
+
+
= ),,1,0( nk L=
2004-11-1 88
4.高斯—埃尔米特(Gauss-Hermite)求积公式
[ ] ),(,,)(
2
+∞?∞==
baex
x
ρ
∑
∫
=
∞
∞+
≈
n
k
kk
x
xfAdxxfe
0
)()(
2
的零点。为)()1(
22
1
1
1
1
x
n
n
xn
nk
e
dx
d
eHx
+
+
+
+
=
[]
2
'
1
2
)(
)!1(2
kn
n
k
xH
n
A
+
+
+
=
π
[] )(
)!22(2
)!1(
)22(
1
ξ
π
+
+
+
+
=
n
n
f
n
b
fR
( )∞+∞?∈，ξ
关于高斯点和求积系数有详细的表格可供参阅。
2004-11-1 89
四．一般情况下的高斯型求积公式的构造例：求高斯型求积公式
)()()(
2211
1
0
xfAxfAdxxfx +≈
∫
的系数及节点。
21
,AA
21
,xx
解：上面的求积公式对函数f(x)=1，
应该是精确成立的。
32
,,xxx
于是有下述公式成立：
2004-11-1 90
续
5
2
1
0
2211
==+
∫
dxxxAxAx
3
2
1
0
21
==+
∫
dxxAA
∫
==+
1
0
2
2
2
21
2
1
7
2
dxxxAxAx
9
2
1
0
3
2
3
21
3
1
==+
∫
dxxxAxAx
2004-11-1 91
续
=
9
2
7
2
1
3
2
3
1
2
2
2
1
2
1
xx
xx
A
A
写成矩阵有：
=
5
2
3
2
11
1
212
1
xxA
A
上面两式右端写成等式，并化简有：
2004-11-1 92
续
+=
+?++=
9
2
)(
7
2
5
2
)(
9
2
)(
7
2
)(
3
2
2121
21
2
221
2
1
2
21
xxxx
xxxxxxxx
=
=
0
9
1
5
9
7
1
0
3
1
7
1
9
1
7
1
22
uv
uuvv
得
vxx =+
21
uxx =
21
令
21
5
=u
9
10
=v
2004-11-1 93
续
389111.0
1
=A
277556.0
2
=A
821162.0
1
=x 289949.0
2
=x
从而有
()
()289949.0277556.0
821162.0389111.0)(
1
0
f
fdxxfx
+
≈
∫
Remark：正交多项式构造的困难，是
Gauss型求积公式的不足之处。
2004-11-1 94
§4.5 数值微分以离散数据近似表达在节点处的微分,通常称为数值微分。
))(,(
kk
xfx
),......2,1,0( nk =
)(xfy =
k
x
一,插值型求导公式给出列表函数，可建立插值多项式，取作为的近似函数，
则称为插值型求导公式。
)(xfy =
)(xp
n
),()( xpxf
n
′
≈
′
)(xf
′
)(' xp
n
2004-11-1 95
插值型求导公式（续）
)(
)!1(
)(
)()(
1
)1(
x
n
f
xpxf
n
n
n +
+
+
+= ω
ξ
由][
n
xx,
0
∈ξ
)(
)!1(
)(
)(
)!1(
)(
)()(
)1(
1
1
)1(
ξ
ω
ω
ξ
+
+
+
+
+
+
′
+
=
′
′
n
n
n
n
n
f
dx
d
n
x
x
n
f
xpxf得若要求节点上的导数值，有余项k
x
)(
)!1(
)(
)()(
1
)1(
kn
n
knk
x
n
f
xpxf
+
+
′
+
=
′
′
ω
ξ
2004-11-1 96
插值型求导公式（续）
若f(x)的各阶导数有界，即
),2,1,0()(
)1(
L=≤
+
nMxf
n
),,2,1,0(],,[ nkbax
k
L=∈
节点
)(0)(
)!1(
)()( ∞→→?
+
≤
′
′
nab
n
M
xpxf
n
kk则误差
),,1,0()()( nkxpxf
knk
L=
′
≈
′
从而有为讨论方便，假定所给节点是等距的，我们限定求节点上的各阶导数值。
2004-11-1 97
1.两点公式
[])()(
1
)(
101
xfxf
h
xp +?=
′
得
01
xxh?=其中由余项
)()()(
1
01
0
0
10
1
1
xf
xx
xx
xf
xx
xx
xp
+
=由
)(
!2
)(
)()(
21 kkk
x
f
xpxf ω
ξ
′
′′
=
′
′
),(
10
xx∈ξ
[]
∈
′′
+?=
′
∈
′′
=
′
),()(
2
)]()([
1
)(
),(
!2
)(
)()(
1
)(
1022011
101
1
010
xxf
h
xfxf
h
xf
xx
fh
xfxf
h
xf
ξξ
ξ
ξ
得
2004-11-1 98
2.三点公式由
)(
))((
))((
)(
))((
))((
)(
))((
))((
)(
2
1202
10
1
2101
20
0
2010
21
2
xf
xxxx
xxxx
xf
xxxx
xxxx
xf
xxxx
xxxx
xp
+
+
=
)(
2
)1(
)(
1
)2(
)(
2
)2)(1(
)(
2
1002
xf
tt
xf
tt
xf
tt
thxp
+
+
=+
hxx 2
02
+= thxx +=
0
hxx +=
01
令
2004-11-1 99
三点公式（续）
由两边对t求导
[]
[][])1(
2
)(
)2()(
)2()1(
2
)(
)(
2
1
0
02
++?+?
+?=+
′
tt
xf
ttxf
tt
xf
thxph
)}()12()()44()()32{(
2
1
)(
21002
xftxftxft
h
thxp?+=+
′
),())((
!3
)(
))((
!3
)(
))((
!3
)(
)(')(
201020
212
xxxxxx
f
xxxx
f
xxxx
f
xpxf
kk
∈
′′′
+
′′′
+
′′′
=?
′
ξ
ξξ
ξ
2004-11-1 100
三点公式（续）
注意在中，当时得
thxx +=
0 210
,,xxx
2,1,0=t
{})2)((
!3
)(
)()(4)(3
2
1
)(
1
2100
hh
f
xfxfxf
h
xf
′′′
+?+?=
′
ξ
)(
3
)()(4)(3(
2
1
1
2
200
ξf
h
xfxfxf
h
′′′
+?+?=
{} ))((
!3
)(
)()(
2
1
)('
2
2010
hh
f
xfxf
h
xf?
′′′
++?=
ξ
{}
2
2
20
6
)(
)()(
2
1
h
f
xfxf
h
ξ
′′′
+?=
2004-11-1 101
三点公式（续）
{}
2
3
2102
2
!3
)(
)(3)(4)(
2
1
)( h
f
xfxfxf
h
xf
ξ
′′′
++?=
′
)(
3
)(3)(4)(
2
1
3
2
210
ξf
h
xfxfxf
h
′′′
++?=
),(,,
20321
xx∈ξξξ
其中类似的，我们还可以建立高阶数值微分公式：
L,2,1)()(
)()(
=≈ kxpxf
k
n
k
2004-11-1 102
三点公式（续）
如可以得到二阶三点公式
[])(
6
)()()(2)(
1
)(
1
)4(
2
1210
2
0
ηξ f
h
fhxfxfxf
h
xf +
′′′
+?=
′′
)(
12
)()(2)(
1
)(
2
)4(
2
210
2
1
ηf
h
xfxfxf
h
xf?+?=
′′
)](
6
)()()(2)([
1
)(
3
)4(
2
2211
2
2
ηξ f
h
fhxfxfxf
h
xf +
′′′
++?=
′′
),(,,,,
2032121
xx∈ηηηξξ
其中
2004-11-1 103
二.Taylor展开法根据Taylor展开式可得
)(
2
)()(
)(
1
ξf
h
h
xfhxf
xf
kk
k
′′
+
=
′
)(
2
)()(
)(
2
ξf
h
h
hxfxf
xf
kk
k
′′
+
=
′
)(
62
)()(
)(
3
2
ξf
h
h
hxfhxf
xf
kk
k
′′′
+
=
′
)(
12
)()(2)(
)(
4
)4(
2
2
ξf
h
h
hxfxfhxf
xf
kkk
k
+?+
=
′′
),(
1
hxx
kk
+∈ξ
),(
2 kk
xhx?∈ξ
),(
3
hxhx
kk
+?=ξ
),(
4
hxhx
kk
+?=ξ
（1）
（2）
（3）
（4）
2004-11-1 104
Taylor展开法（续）
以(3)为例
)
~
(
!3
)(
2
)()()(
1
32
ξf
h
xf
h
xfhxfhxf
kkkk
′′′
+
′′
+
′
+=+
)
~
(
!3
)(
2
)()()(
2
32
ξf
h
xf
h
xfhxfhxf
kkkk
′′′
′′
+
′
=?
),(
~
,
~
2021
xx∈ξξ
)(a
)(b
其中两式相减得
))
~
()
~
((
6
)(2)()(
21
3
ξξ ff
h
xfhhxfhxf
kkk
′′′
+
′′′
+
′
=+
2004-11-1 105
Taylor展开法（续）
[ ],)
~
()
~
(
2
1
)(
213
ξξξ fff
′′′
+
′′′
=
′′′
从而有：
[])(
6
)()(
2
1
)(
3
2
ξf
h
hxfhxf
h
xf
kkk
′′′
+=
′
),(
3
hxhx
kk
+?∈ξ
若连续，则
)(xf
′′′
Mffm ≤
′′′
+
′′′
≤ )
~
()
~
((
2
1
21
ξξ
（为的最小值与最大值）故由连续函数的介值定理，必存在使得
Mm,
)(xf
′′′
,
3
ξ
由(a)可以得到(1)；由(b)可以得到(2)。
2004-11-1 106
Taylor展开法（续）
若将(a)(b)写为
)
~
(
!4
)(
!3
)(
2
)()()(
1
)4(
43
2
ξf
h
xf
h
xf
h
xfhxfhxf
k
kkkk
+
′′′
+
′′
+
′
+=+
)
~
(
!4
)(
!3
)(
2
)()()(
2
)4(
43
2
ξf
h
xf
h
xf
h
xfhxfhxf
k
kkkk
+
′′′
′′
+
′
=?
两式相加可以类似的得到（4）
2004-11-1 107
Remark
Remark1：在数值微分计算中，并非步长越小精度越高。这是因为数值微分对舍入误差非常敏感，
它随步长h的缩小而增大，导致计算不稳定。
Remark2：在数值微分计算中，当插值多项式收敛到函数f(x)时，P′
n
（x）不一定收敛到f ′(x)。
Remark3：为了避免上述问题，可以用样条插值函数的导函数来代替f(x)的导函数。
