计 算 方 法 学 习 指 导
温师院数学与信息科学学院
2003.7
目 录
第一章 引论...........................................................................................................................1
1.1 内容提要...................................................................................................................1
1.2 典型例题................................................................................................................3
第二章 解线性方程组的直接法..............................................................................7
2.1 内容提要...................................................................................................................7
2.2 典型例题..............................................................................................................14
第三章 插值法与最小二乘法..................................................................................21
3.1 内容提要.................................................................................................................21
3.2 典型例题.................................................................................................................27
第四章 数值积分与微分.............................................................................................47
4.1 内容提要.................................................................................................................47
4.2 典型例题..............................................................................................................50
第五章 常微分方程数值解法..................................................................................63
5.1 内容提要.................................................................................................................63
5.2 典型例题..............................................................................................................68
第六章 逐次逼近法........................................................................................................79
6.1 内容提要.................................................................................................................79
6.2 典型例题..............................................................................................................84
第一章 引论
1.1 内容提要
1,绝对误差与绝对误差限
1
*
设为准确值，x x是的一个近似值，称e为近似值
*
x xx?=
*
x的绝对误差，简称误差。如果有常数ε使得εe ≤，则称ε为近似值x的绝对误差限，简称误差限。
2,相对误差和相对误差限
设是准确值，
*
x x是的一个近似值，称为近似值
*
x
**
/)( xxx? x的相对误差，记作。
在实际计算中，是未知的，常以
r
e
*
x xx /)
*
xe
r
(=作为相对误差。如果有常数使得
r
ε≤
r
e (或
r
ε
r
e ≤ )，则称
r
ε为相对误差限。
r
ε
3,有效数字
如果近似数x的绝对误差不超过其上某一位数字的半个单位，且该位数字到x的第一位非零数字共有位，则称n x近似时具有位有效数字，简称
*
x n x有位有效数字。 n
定理1.1 设x是的近似值，它的表达式为，则
*
x
n
k
aaax L
21
*
.010 ×±= x的有效数字位数与x的相对误差之间有如下关系,
⑴若x具有位有效数字时，n x的相对误差
r
ε满足
1
1
10
2
1
+?
×≤
n
r
a
ε
⑵若x的相对误差
r
ε满足
1
1
10
)1(2
1
+?
×
+
≤
n
r
a
ε时，x至少具有位有效数字。 n
4,数据误差的影响
给定函数。设分别为的近似值，则 ),(
21
xxfy =
21
,xx
*
2
*
1
,xx
2
2
21
1
1
21
2
*
2
2
21
1
*
1
1
21
21
*
2
*
1
),(),(
)(
),(
)(
),(
),(),()(
e
x
xxf
e
x
xxf
xx
x
xxf
xx
x
xxf
xxfxxfye
+
=
+?
≈
=
rrr
e
y
x
x
xxf
e
y
x
x
xxf
y
ye
ye
2
2
2
21
1
1
1
21
),(),()(
)(
+
≈=
由以上两式得
() )()(
2121
xexexxe +≈+,)()()(
2
21
2
1
21
1
21
xe
xx
x
xe
xx
x
xxe
rrr
+
+
+
≈+
() )()(
2121
xexexxe?≈?,)()()(
2
21
2
1
21
1
21
xe
xx
x
xe
xx
x
xxe
rrr
≈?
() )()(
211221
xexxexxxe +≈,)()()(
2121
xexexxe
rrr
+≈
)(
)(
2
2
2
1
2
1
2
1
xe
x
x
x
xe
x
x
e?≈
,)()()(
21
2
1
xexe
x
x
e
rrr
≈，0
2
≠x
5,设计算法时要注意的几个问题
⑴应用数值稳定的递推公式
设是由递推公式 LL,,,,
21 n
xxx
（1.1）
==
给定
,2,1),(
0
1
x
nxFx
nn
L
得到的。若有误差e，则实际上只能得到
0
x
0
=
==
000
1
~
,2,1),
~
(
~
exx
nxFx
nn
L
记
nnn
xxe
~
=，L,2,1=n。如果存在不依赖于的常数C使得 n
0
eCe
n
≤，L,2,1=n
则称递推公式（1.1）是数值稳定的。否则称为数值不稳定的。
⑵注意简化运算步骤，减少运算次数
⑶要避免相近数相减
⑷多个数相加，应先将绝对值较小的数相加之后，再依次与绝对值较大的数相加
⑸避免小数作除数和大数作乘数
2
1.2 典型例题
例1.1 问
7
22
,141.3,142.3分别作为π的近似值各具有几位有效数字？
解 L14159265.3=π，记142.3
1
=x，141.3
2
=x，
7
22
3
=x。
由L00040.0
1
=? xπ知
3
1
4
10
2
1
10
2
1
×≤?<× xπ，因而具有4位有效数字。
1
x
由L00059.0
2
=? xπ知
2
2
3
10
2
1
10
2
1
×≤?<× xπ，因而具有3位有效数字。
2
x
由L00126.0
3
=? xπ知
23
10
2
1
7
22
10
2
1
×≤?<× π，因而具有3位有效数字。
3
x
例1.2 已知近似数x有两位有效数字，试求其相对误差限。
解 有效数字，所以相对误差限为2=n %510
2
1
10
2
1
121
1
=×≤×≤
+?+?n
r
a
ε。
例1.3 已知近似数x的相对误差限为0.3%，问x至少有几位有效数字？
解 设x具有位有效数字，则n
1
1
2
10
)1(2
1
102
1
1000
3
%3.0
×
+×
≤
×
<=
a
所以1，即。 1?=? n 2=n
例1.4 为使70的近似数的相对误差小于0.1%，问查开方表时，要取几位有效数字？
解 设查表时取位有效数字，由于n 9708 ≤≤，所以8
1
=a,
由公式%1.010
2
1
1
1
≤×≤
+?n
r
a
ε知，只需取3=n。
例1.5 求方程的两个根，使它们至少具有4位有效数字，取0156
2
=+? xx
982.27783 =？
解 方程的两个根为
98.5578328
1
≈+=x，0.01786
982.55
1
78328
1
78328
2
≈≈
+
=?=x
注：此处若直接利用018.0982.272878328
2
=?≈?=x，则只能得到两位有效数字。
3
一般地，
xx
xx
++
=?+
1
1
1
2
2
当1<<x时，
2
)
2
(sin2cos
x
x =?1,
当1>>x时，
)1(
1
1
11
+
=
+
xxxx
例1.6 计算
6
)12(?=a，取4.12 ≈，利用下列等式计算,
⑴
6
)12(
!
+；⑵27099?；⑶
3
)223(?；⑷
3
)223(
1
+
问哪一个得到的结果最好？
解 显然以上四个算式是恒等的，但当取4.12 ≈计算时，四个算式的误差是不同的。
记2 4
*
=x，，.1=x xxx?=?
*
⑴设
6
11
)1(
1
)(
+
==
t
tfy，则按⑴式计算产生的误差为
xxfxfxfye?′≈?= )()()()(
11
*
11
同理,
⑵设，则按⑵式计算产生的误差为 ttfy 7099)(
22
==
xxfxfxfye?′≈?= )()()()(
22
*
22
⑶设，则按⑶式计算产生的误差为
3
33
)23()( ttfy?==
xxfxfxfye?′≈?= )()()()(
33
*
33
⑷设
3
44
)23(
1
)(
t
tfy
+
==，则按⑷式计算产生的误差为
xxfxfxfye?′≈?= )()()()(
44
*
44
估计以上四个式子中的)(,)(,)(,)(
4321
xfxfxfxf ′′′′，即可得出⑷式误差最小的结论。
具体计算结果如下,
4
⑴
3
6
102.5
)12(
!
×≈
+；
⑵0.127099 ≈?；
⑶
33
100.8)223(
×≈?；
⑷
3
3
101.5
)223(
1
×≈
+
例1.7 数列{满足递推公式}
n
x ),2,1(,110
1
L=?=
nxx
nn
，若取41.12
0
≈=x，
问按上述递推公式，从计算到时误差有多大？这个计算过程稳定吗？
0
x
10
x
解 记作为41.1
*
0
=x 2
0
=x的近似值，则
2*
000
105.0
×<?= xxε，由
110,110
*
0
*
101
=?= xxxx
得
0
*
00
*
00
*
11
1010)110()110( ε<?==? xxxxxx
0
2*
11
*
11
*
22
1010)110()110( ε<?==? xxxxxx
M
8
0
10*
1010
105.010 ×<<? εxx
显然，误差积累很大，按递推公式110
1
=
nn
xx求时，会把初始误差
10
x
0
ε扩大10倍，
使计算精度受到严重影响，因此，这个计算过程不稳定。
10
5
第二章 解线性方程组的直接法
2.1 内容提要
1,直接法概述
目前，计算机上解线性代数方程组的数值方法尽管很多，但归纳起来，大致可以分为两大类：一类是直接法；另一类是迭代法。一般，对中、低阶以及高阶带形线性代数方程组，
用直接法有效，而对于高阶线性代数方程组，特别是对于系数矩阵为大型稀疏矩阵的线性代数方程组，用迭代法有效。
2,高斯(Gauss)消去法
高斯消去法是一个古老的求解线性代数方程组的直接法，但由它的改进、变形得到的如选主元消去法、三角分解法等，仍然是目前计算机上解中、低阶稠密矩阵方程组的常用有效方法。
7
n设有阶代数线性方程组
=+++
=+++
=+++
nnnnnn
nn
nn
bxaxaxa
bxaxaxa
bxaxaxa
L
L
L
L
2211
22222121
11212111
写为矩阵形式
bx =A,
其中
=
nnnn
n
n
aaa
aaa
aaa
A
L
MMMM
L
L
21
22221
11211
，,
=
n
x
x
x
M
2
1
x
=
n
b
b
b
M
2
1
b
为保证方程组有唯一解，设系数矩阵为非奇异阵。将A bx =A记为，高斯消去法描述如下,
)1()1(
bx =A
⑴消元过程
第次消元(k )11?≤≤ nk：假定已完成1?k步消元，得到,
)()( kk
A bx =
=
)()()(
)()()(
)2(
2
)2(
2
)2(
22
)1(
1
)1(
1
)1(
12
)1(
11
)()(
),(
k
n
k
nn
k
nk
k
k
k
kn
k
kk
n
n
kk
baa
baa
baa
baaa
A
L
MM
L
MMO
LLL
LLL
b
8
)k
k n
k
x
此时，的右下角阶矩阵必非奇异。
(
A 1+?kn
不妨设（称为主元素），消去方程组的第至个方程中的，计算公式为,
0
)(
≠
k
kk
a
)(k
kk
a
)()( kk
A bx = 1+
(a)计算行乘数
)(
)(
k
kk
k
ik
ik
a
a
m =，),,1( nki L+=
(b)第i行元素减去第行对应元素乘以m，即 k
ik
)()()1( k
kjik
k
ij
k
ij
amaa?=
+
,
)()()1( k
kik
k
i
k
i
bmbb?=
+
,
),,1,,,1( nkjnki LL +=+=
得到，此时的右下角
)1()1( ++
=
kk
A bx
)1( +k
A kn?阶矩阵必非奇异。当完成第次消元后，得到，其中
1?n
)()( nn
A bx =
=
)()(
)()()(
)2(
2
)2(
2
)2(
22
)1(
1
)1(
1
)1(
12
)1(
11
)()(
),(
n
n
n
nn
k
k
k
kn
k
kk
n
n
nn
ba
baa
baa
baaa
A
M
L
MM
LL
LL
b，。 0
)(
≠
n
nn
a
原方程组变形为
=
)(
)2(
2
)1(
1
2
1
)(
)2(
2
)2(
22
)1(
1
)1(
12
)1(
11
n
nn
n
nn
n
n
b
b
b
x
x
x
a
aa
aaa
MMMO
L
L
以上过程称为消元过程。
⑵回代过程
求解上述三角形方程组，得到
=?=
=
∑
+=
)1,2,,2,1(,/)(
/
)(
1
)()(
)()(
Lnniaxabx
abx
i
ii
n
ij
j
i
ij
i
ii
n
nn
n
nn
此过程称为回代过程。
Gauss消去法算法简述,
①输入：系数矩阵( ),bA
②消去：对于 1,,2,1?= nk L
a.如果，找非零元，交换两行，若找不到非零元，则算法停止。 0
)(
=
k
kk
a
b.消元
)(
)(
k
kk
k
ik
ik
a
a
m =，( ),,1 nki L+=
,
)()()1( k
kjik
k
ij
k
ij
amaa?=
+
),,1,,,1( nkjnki LL +=+=
b，(
)()()1( k
kik
k
i
k
i
bmb?=
+
),,1,,,1 nkjnki LL +=+=
③回代,
a.如果，算法停止。 0
)(
=
n
nn
a
b.按x的顺序回代求解。
11
,,,xx
nn
L
=?=
=
∑
+=
)1,2,,2,1(,/)(
/
)(
1
)()(
)()(
Lnniaxabx
abx
i
ii
n
ij
j
i
ij
i
ii
n
nn
n
nn
④输出：结果。
n
xxx,,,
21
L
Gauss消去法的运算量
第步消元时，共运算乘法(k )1)( + knkn次，除次kn?次，故完成步消元乘除法总数为
1?n
6
5
23
)2)((
231
1
nnn
knkn
n
k
+=+
∑
=
回代过程乘除法总数为
22
)1(
2
1
nn
in
n
i
+=+?
∑
=
总数为)(
333
2
3
2
3
nO
nn
n
n
+=?+
3,Gauss列主元素消去法
9
10
)(
)(
k
nk
k
kk
a
a
M
显然，Gauss消去法能够进行的条件是各步消元的主元素，如果在消元过程中发现某个主元素为零，即使
0
)(
≠
k
kk
a
0)det( ≠A，则消元过程将无法进行；其次，即使主元素不为零，
但如果主元素的绝对值很小，用它们作除数，将使该步消元的乘数绝对值很大，势必造成舍入误差的严重扩散，以致于方程组解的精度受到严重的影响。
)
k
)
)k (
i
a
) a
(k
kk
a
因此，为避免小主元作除数，在Gauss消去法中加入选主元过程，即在第步消元时 1,,1(?= nk L
=
)()(
)()(
)2(
2
)2(
2
)2(
22
)1(
1
)1(
1
)1(
12
)1(
11
)()(
),(
k
n
k
nn
k
k
k
kn
n
n
kk
ba
ba
baa
baaa
A
L
M
L
MMO
LLL
LLL
b
首先在的第列主对角元以下元素中挑选绝对值最大者
，并通过交换的第行与第行对应元素，使位于主对角线上，
仍记为，然后再进行消元计算。称每一步都按列选主元的Gauss消去法为Gauss列主元素消去法。
(
A
ni
k
≤
)(k
kk
a
k )(
)(
nika
k
ik
≤≤
k
i
)k
k
k
(k ≤ ),(
)()( kk
A b k
)(k
ki
k
4,系数矩阵的三角分解
从矩阵变换的角度来看，每一次消元都相当于矩阵A左乘一个初等行变换阵，最终将原方程组系数矩阵化为简单的上三角阵。而所有初等行变换阵之积的逆是一个单位下三角阵，所以，不带行交换的Gauss消元过程产生了一个单位下三角阵和一个上三角阵U，
它们的乘积等于，即
L
A
LUA =
称之为矩阵的分解，又称为Doolittel分解。 A LU
当然也可以将A分解为一个下三角阵L与一个单位上三角阵U的乘积，此时的LU分解又称为Crout分解。
当A为n阶对称正定阵，则必有分解
T
LLA =，其中L是对角元全为正的下三角阵，
称为正定阵的Cholesky分解。
定理2.1 设阶方阵的前n A 1?n个顺序主子式都不为零，则的分解存在且唯一。
A LU
定理2.2 设A为非奇异矩阵，则必存在排列矩阵P以及单位下三角阵和非奇异上三角阵U，使。
L
A
LUPA =
5,直接三角分解法
直接三角分解法的基本想法是，一旦实现了矩阵A的分解，那么求解方程组的问题就等价于求解两个三角形方程组
LU bx =
by =L和yx =U。
直接三角分解法的步骤,
⑴实现Doolittle分解，其中为单位下三角阵，U为上三角阵。其分解公式LUA = L
为：对，计算 nr,,3,2,1 L=
by =L
yx =
b=
r,2,1 L=
a
ir
←
),,1,(
1
1
nrriulau
r
k
kirkriri
L+=?=
∑
=
),,2,1(,
1
1
nrri
u
ula
l
rr
r
k
krikir
ir
L++=
=
∑
=
⑵求解，其计算公式为,
=?=
=
∑
=
),,3,2(
1
1
11
nrylby
by
r
i
irirr
L
⑶求解U，其计算公式为,
=
=
=
∑
+=
)1,2,,2,1(
/
1
Lnnr
u
xuy
x
uyx
rr
n
ri
irir
r
nnnn
事实上，计算U和解y的公式是一致的，因此系数矩阵的Doolittle分解与解方程组可以同时进行，都归结为对数组yL ),( bAA =的分解计算，然后再解U，称这种求解方式为紧凑格式。
yx =
6,部分选主元的Doolittle分解
部分选主元的Doolittle分解法与列主元消去法在理论上是等价的。具体步骤如下,
⑴对于进行（部分）选主元的分解计算。 n,
(1.1)计算中间量S作为可能的主元，并存入 ),1,( nrri
i
L+=
ir
a
),,1,(
1
1
1
1
nrriaaaulaS
r
k
krikir
r
k
krikiri
L+=?=?=
∑∑
=
=
(1.2)挑选绝对值最大的，即确定行号i，使满足
i
S
r i
nir
i
SS
r
≤≤
= max。
(1.3)交换两行：当时，交换A的第ri
r
≠ r行与第行的对应元素，交换后的元素以它在A中所处的新位置记。
r
i
11
(1.4)分解计算（以下两步可调换次序，因已计算）
rr
a
)1,,1(
1
1
1
1
++=?=?=←
∑∑
=
=
nriaaaulaua
r
k
kirkri
r
k
kirkririri
L
12
1此处,,1 ++= nr Li意味着对增广矩阵(进行分解计算，即紧凑格式 ),bA
),,2,1(,
1
1
nrri
a
a
S
S
u
ula
la
rr
ir
r
i
rr
r
k
krikir
irir
L++===
=←
∑
=
按上述方法完成分解，便实现了系数矩阵分解LUPA =，P是排列阵，即是将分解过程中所有的行交换依次作用在单位阵上的结果。
⑵求解U yx =
=
=
=
==
∑∑
+=
++
+
+=
++
)1,2,,2,1(
//
1
1,1,
1,
1
1,1,
Lnnr
a
aaa
a
u
xuy
x
aaauyx
rr
n
ri
nirinr
nr
rr
n
ri
irir
r
nnnnnnnnnn
7,平方根法
平方根法是用于解系数矩阵为对称正定阵的线性代数方程组，其计算步骤如下,
⑴对矩阵A进行Cholesky分解
T
LLA =，按列确定下三角阵的元素，即对于
，计算
L
nj,,2,1 L=
∑
=
=
1
1
2/12
)(
j
k
jkjjjj
lal
),,2,1(,
1
1
njji
l
lla
l
jj
j
k
jkikij
ij
L++=
=
∑
=
⑵求解下三角形方程组,by =L
=
=
=
∑
=
),,3,2(
/
1
1
1111
ni
l
ylb
y
lby
ii
i
j
jiji
i
L
⑶求解上三角形方程组,yx =
T
L
=
=
=
∑
+=
)1,2,,2,1(
/
1
Lnni
l
xly
x
lyx
ii
n
ij
jjii
i
nnnn
平方根法是目前在计算机上解对称正定方程组的一个有效方法，应用较广泛，其计算量和存贮量都比用消去法约节省近一半，且不需要选主元，便能求得较精确的数值解。
8,解三对角方程组的追赶法
在一些实际问题中，例如用差分法解二阶线性常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等，都会要求解系数矩阵为对角占优的三对角线方程组
=
+
+
n
n
n
i
i
i
n
n
n
i
i
i
nn
nnn
iii
f
f
f
f
f
f
f
f
f
x
x
x
x
x
x
x
x
x
ba
cba
cba
cba
cb
1
2
1
1
3
2
1
1
2
1
1
3
2
1
111
222
11
M
M
M
M
OOO
OOO
OOO
OOO
OOO
OOO
记为。其中A满足下列条件：（称为对角占优的三对角线矩阵） fx =A
（1）0
11
>> cb
（2）)1,,3,2,0(,?=≠+≥ nicacab
iiiii
L
（3）0>>
nn
ab
追赶法的计算步骤如下,
⑴对A进行Crout分解。设将A分解为如下两个三角形矩阵的乘积
L＝
nn
ii
αγ
αγ
αγ
α
OO
OO
OO
OO
22
1
13
U＝
1
1
1
1
1
1
1
2
2
1
n
i
i
i
β
β
β
β
β
β
O
OO
O
OO
则其中元素的计算公式为
==
==
=?==
),,3,2(
)1,2,1(/
),,3,2(,
111
nia
nic
niabb
ii
iii
iiii
L
L
L
γ
αβ
βαα
14
i
正好已处于的位置上，故只需计算
i
a
ii
βα,，其次序为 γ
nnn
αβαβαβα →→→→→→→
112211
L
⑵求解，其计算公式为 fy =L
=
=
=
),,3,2(
/
1
111
ni
yaf
y
fy
i
iii
i
L
α
α
⑶求解U，其计算公式为 yx =
=?=
=
+
)1,2,,2,1(
1
Lnnixyx
yx
iiii
nn
β
上述解三对角方程组的方法称为追赶法。追赶法不但方法简练，运算量少，而且是数值稳定的。
2.2 典型例题
例2.1 用Gauss消去法解方程组
=
2
6
15
15
1113
1111
11318
43312
4
3
2
1
x
x
x
x
解 消元过程相当于对增广矩阵进行初等变换
21113
61111
1511318
1543312
→?
14
13
12
4
1
12
1
)
2
3
(
rr
rr
rr
4
7
0
4
7
4
7
0
4
19
3
2
4
3
4
5
0
2
15
5
2
7
2
3
0
1543312
→?
24
23
)
6
7
(
)
6
5
(
rr
rr
7
11
2
15
15
6
35
3
7
00
6
29
3
11
00
5
2
7
2
3
0
43312
→?
3
4
11
7
rr
0
33
91
000
11
6
29
3
11
00
2
15
5
2
7
2
3
0
1543312
得同解三角方程组
=
=+
=++?
=++?
0
33
91
11
6
29
3
11
2
15
5
2
7
2
3
1543312
4
43
432
4321
x
xx
xxx
xxxx
回代得,，，0
4
=x 3
3
=x 2
2
=x 1
1
=x。
例2.2 用Gauss列主元消去法解如下方程组。
=
6
7
4
515
0710
623
3
2
1
x
x
x
解 对增广矩阵进行选主元及消元，加框的元素是主元素
→?
21
6515
70710
4623
rr
6515
4623
70710
－
－
→?
13
12
2
1
)
10
3
(
rr
rr
2
5
5
10
61
6
10
1
70710－
→?
32
rr
10
61
6
2
5
5
70
10
1
2
5
710－
3
)
25
1
( rr
→
2
5
31
5
31
2
5
5
2
5
70710－
2
5
15
得同解三角方程组
=
=+
=?
5
31
5
31
2
5
5
2
5
7710
3
32
21
x
xx
xx
回代求解得 1
3
=x，，1
2
=x 0
1
=x
数方程组：
回代求解得对增广矩阵进行选主元及消元，加框的元素是主元素例2.3 采用四位有效数字，分别用Gauss消去法和列主元素Gauss消去法求解线性代
=
000.1
8338.0
6867.0
100.1210.1331.1
000.1000.1000.1
9000.08100.07290.0
3
2
1
x
x
x
解 Gauss消元过程相当于对增广矩阵进行初等变换
000.1100.1210.1331.1
8338.0000.1000.1000.1
6867.09000.08100.07290.0
→?
13
12
826.1
372.1
rr
rr
2540.05430.02690.0
1084.02350.01110.0
6867.09000.08100.07290.0
→?
23
423.2 rr
008700.002640.0
1084.02350.01110.0
6867.09000.08100.07290.0
3295.0
3
=x，2790.0
2
=x，2252.0
1
=x
000.1100.1210.1331.1
8338.0000.1000.1000.1
6867.09000.08100.07290.0
→?
13
12
5477.0
7513.0
rr
rr
.0
.0
.1331.1
→?
31
rr
2975.01473
1736.00909
100.1210
6867.09000.08100.0729.0
8338.0000.1000.1000.1
000.1100.1210.1331.1
1390.0
08250.0
000.1
→?
32
rr
16
08250.01736.00909.0
1390.02975.01473.0
000.1100.1210.1331.1
,001000.0
.02975.01473.0
.1100.1210.1331.1
3280.0
3
→?
23
6171.0 rr
003280
1390
000
=x 2812.0
2
=x 2246.0
1
=x
3279.0= 2814.0
2
= 2245.0
1
=
2011612
6384
10278
5124
→
201163
6381
10272
5124
→
201103
6321
0032
5124
1403
121
12
1
3
24
14321 x
→
03
21
32
24
1
12
00
51
回代求解得,,
与精确解，x，
3
x x比较，显然列主元素消去法比一般的Gauss消去法得到的解误差要小。
例2.4 用Doolittle法对矩阵A作LU分解。
A＝
解 从上至下，从左至右，分解过程如下,
2011612
6384
10278
5124
1403
1221
0032
5124
204
12
00
51
→
每一次只变化加框部分的元素，具体计算公式参见Doolittle分解算法。
所以A＝LU，其中
L＝，U＝
例2.5 利用Doolittle分解的紧凑格式解下列方程组
=
13
8
142
210
3
2
1
x
x
17
解 所谓紧凑格式就是对增广矩阵进行Doolittle分解，过程如下,
→
13142
8210
14321
13142
8210
14321
→
13102
8210
14321
→
15502
8210
14321
即系数矩阵可分解为
=
5
21
321
102
10
1
142
210
321
得同解三角方程组为
=?
=+
=++
155
82
1432
3
32
321
x
xx
xxx
回代求解得 3
3
=x，，2
2
=x 1
1
=x
例2.6 用部分选主元的Doolittle法紧格式解矩阵方程组BAX =，其中
=
25681161
642781
16941
4321
A,
=
24082
5628
1210
24
B
解 对增广矩阵进行分解计算，每一次只变化加框部分的元素，过程如下,
2408225681161
5628642781
121016941
244321
（第一步选主元之后，矩阵各元素不变）→?
列的主元计算第2
2408225681141
5628642761
121016921
244321
12101697/11
562864277/31
2387825278141
244321
→?
42
rr
计算第3
121016921
5628642761
2408225681141
244321
→
主元
167/367/11
647/667/31
25278141
4321
→?
分解计算
1210
5628
23878
24
列的
18
→?
分解计算
7/11
7/31
141
21
7/11
7/31
2387825278141
244321
11/24011/2411/6
487/66487/66
2387825278
2443
=
111/67/11
17/31
11
1
L
=
14
21
=
0010
0100
1000
0001
P
=
=
=++
=+++
011/24
7/66487/66
782527814
4432
4
43
432
4321
x
xx
xxx
xxxx
+
+
66
7814
2
2
21
x
xx
=
=
=
=
0
1
0
1
4
3
2
1
x
x
x
x
=
=
=
=
1
0
1
0
4
3
2
1
x
x
x
x
=
10
01
10
01
X
→?
分解计算
24
7/66
78
3
+
++
11/24
487/
252
43
4
43
43
43
x
xx
xx
xx
12101611/6
487/66487/66
所以PA＝LU，其中
，U
11/
48
252
4
得以下两个三角方程组
,
=
=
=
=
11/24
48
238
2
回代得
,
所以
例2.7 用平方根法解对称正定方程组
19
=
7
3
5
1203
022
323
3
2
1
x
x
x
解 可以验证系数矩阵对称正定，故可用平方根法来解。
→
1203
22
3
1203
23/2
3
→
1263
3/23/2
3
→
363
3/23/2
3
得三角方程组
=
7
3
5
363
3/23/2
3
3
2
1
y
y
y
，
=
3
2
1
3
2
1
3
63/2
33/23
y
y
y
x
x
x
解得
3/5
1
=y，6/1
2
=y，3/1
3
=y
3/1
3
=x，2/1
2
=x，1
1
=x
例2.8 举例说明一个非奇异矩阵不一定存在LU分解。
解 取，显然A是非奇异矩阵。若A存在LU分解，则有
=
01
10
A
+
=
=
=
dacab
cb
d
cb
a
A
01
01
01
10
,
从而得，则，这与0=b 0=ab 1=ab矛盾。
例2.9 设是实对称阵，且
nnij
aA
×
= )( 0
11
≠a，经过Gauss消去法一步后，A约化为
，其中，是
2
11
A
a
0
α
),,
1n
aL(
12
a=α
2
A 1?n阶方阵，证明是对称阵。
2
A
解 设，显然也是对称阵。则经过Gauss消去法一步后A化为
=
1
11
A
a
A
T
α
α
1
A
T
αα0
α
11
1
11
1
a
A
a
显然
T
αα
11
12
1
a
AA?=是对称阵。
20
第三章 插值法与最小二乘法
3.1 内容提要
21
)
1,多项式插值的概念
设函数在区间[上有定义，且已知在一组互异点上的函数值，若一个不超过n次的多项式函数
，满足
(xfy =
x
n
≤<<L
],ba
y,
0
bxxa <≤
10
)(xP
n
n
yy,,
1
L
),,2,1,0()( niyxP
iin
L== （3.1）
称为的次插值多项式；相应的插值问题称为n次代数插值问题；称为插值节点；称[为插值区间；称式（3.1）为插值条件。
) )
)
n
n
(xP
n
(ix
i
=
(xf
,,2,L
n
1,0 n ],ba
定理3.1 次代数插值问题的解是存在且唯一的。
2,拉格朗日(Lagrange)插值多项式
次代数插值问题的解可表示为如下Lagrange插值多项式,
∑ ∏∑
=
≠
==
==
n
j
n
ji
i ij
i
j
n
j
jjn
xx
xx
yxlyxL
0 00
)()( (3.2)
其中
∏
≠
=
=
=
n
ji
i ij
i
j
nj
xx
xx
xl
0
),,1,0()( L
称为拉格朗日插值基函数，仅与插值节点有关。
3,牛顿(Newton)插值多项式
(1)差商的定义
定义3.1 设在互异的节点上的函数值)(xf
n
xxx,,,
10
L ),,1,0(),( nixff
ii
L==，
称
)(],[ ik
xx
ff
xxf
ik
ik
ki
≠
=
为关于的一阶均差（差商），称 )
k
(xf
i
xx,
)(
],[],[
],,[ kji
xx
xxfxxf
xxxf
jk
jiki
kji
≠≠
=
为关于的二阶均差（差商）。一般地，称 )
k
(xf
ji
xxx,,
1
1210210
10
],,,,[],,,,[
],,,[
=
kk
kkkk
k
xx
xxxxfxxxxf
xxxf
LL
L
为关于的阶均差（差商）。 )
)
(xf
k
xxx,,,
10
L k
定理3.2 (i)关于的阶均差是在这些点上的函数值的线性组合，即
(xf
k
xxx,,,
10
L k )(xf
∑ ∏
=
≠
=
=
k
j
k
ji
i ij
jk
xx
xfxxxf
0 0
10
1
)(],,,[ L
（ii）阶均差关于是对称的，即与节点的排列次序无关。 k ],,,[
10 k
xxxf L
k
xxx,,,
10
L
),2,1,0,(],,,,,,[],,,,,,[
00
kjixxxxfxxxxf
kijkji
LLLLLLL ==且ji ≠。从而可知在定义阶均差时，两个k 1?k均差的选取可以是任意的,
(2)差商与导数之间有如下关系
若在含有的区间上具有k阶导数，则在这一区间内至少有点)(xf
k
xxxx,,,,
10
L ξ，
使得
!
)(
],,[
)(
0
k
f
xxf
k
k
ξ
=L，ξ在之间
k
xxxx,,,,
10
L
(3)差商表
均差的计算过程如表
k
(
k
xfx )
]
,[
1+kk
xxf
一阶差商
],,[
21 ++ kkk
xxxf
二阶差商
],,,[
321 +++ kkkk
xxxxf
三阶差商
],,,,[
4321 ++++ kkkkk
xxxxxf
四阶差商
],,[
321
xxxf ],,,,[
43210
xxxxxf
],[
32
xxf ],,,[
4321
xxxxf
],,[
432
xxxf
],[
43
xxf
0
x
0
f
],[
10
xxf
1
f ],,[
210
xxxf
],[
21
xxf ],,,[
3210
xxxxf
1
x
22
2
x
2
f
3
x
3
f
4
x
4
f
(4) Newton插值多项式
n次代数插值问题的解还可表示为如下Newton插值多项式,
∑ ∏
=
=
+=
++?+=
n
k
k
j
jk
nnn
xxxxxfxf
xxxxxxxxxfxxxxfxfxN
1
1
0
100
110100100
)(],,,[)(
)())(](,,,[)](,[)()(
L
LLL
对于同一个插值条件的次代数插值问题，Newton插值公式与Lagrange插值公式所给出的是同一个多项式，只是形式上不同而已。但Newton插值公式较Lagrange插值公式的一个明显的优点是具有递推性,
n
)())(](,,,[)()(
101101 nnnn
xxxxxxxxxfxNxN+=
++
LL
(5) 差分与等距节点插值公式
定义 设在等距节点)(xf ),,1,0(
0
nkkhxx
k
L=+=上的函数值为，称
k
f
kkk
fff?=?
+1
，
1?
=?
kkk
fff
分别为在处的一阶向前差分和一阶向后差分。符号)
k
(xf xx =,分别称为向前差分算子和向后差分算子。类似地，可定义在)(xf
k
xx =处的阶向前差分和阶向后差分 m
m
,(
0
x
m
1?k
m
f
)(
)
ξ
m 1
1
k
fh=
,1=]
0
,,,
1
L
[= xxxf
k
k
k
f
0
!
+= f
0
)+x
0
(
+=
n
f
)(xf
+
n
thx(
k
k
f
!
),n
1+n
),( ba
)
)!1
)(ξ
(x
],
)(
x
f
x
nn
n
=
ωL
)(
0
xxx
,[
)(
0
xxf
xf?
() xx?=
)
1
M
n+
=
(
)
n
M
≤
,
i
y
∑
=
=
n
j
x)(
k
m
k
m
fff
1
+
=?，?
11
=
k
m
kk
ff
差分与导数的关系为,
)
(
0 k
kk
xf ∈? ξ
差分与差商的关系为,
),2(
!
,
!
],,,[
1
0
10
LL
=
= k
hk
f
x
hk
f
xxxf
k
k
k
k
k
k
k
称
∑ ∏
=
=
=
n
k
k
j
nn
jtthNxN
1
1
0
)()(
为等距节点的Newton向前（差分）插值公式。
称
∑ ∏
=
=
+
=
n
k
k
j
n
nn
jtNxN
1
1
0
)())(
为等距节点的Newton向后（差分）插值公式。
4,插值多项式的余项
定理3.3 （插值多项式余项公式）设在含节点1,0(ix
i
L=的区间[上次可微，是关于给定的
],ba
1 )+n (xP
n
)(xf个节点的次插值多项式，则对于
，存在与
n
],[ ba∈x? x有关的∈ξ，使得
(,
)(
(
)(
1
1
)1(
x
x
n
PxR
n
n
n
+
+
+
=
+
= ω
))(
11 nn
x
+
Lω
若在[上有界,) ](
)1(
xf
n+
,ba
(max
)1(
xf
n
bxa
+
≤≤
,
则
)(
)!1
(
1
1
xxR
n
n
n +
+
+
ω (3.3)
5,Hermite插值多项式
(1)Hermite插值问题：已知函数f(x)在节点bxxxa
n
≤≤≤≤ L
10
)(xH )(
i
xH =
处有，
，。要求插值多项式，使得，
。
ii
yx =)(
'
)(
ii
yx =
f
'
i
y
n
)(
i
xf =′
i,1,0 L=
ni,,1,0 L= H′
,
定理3.4 上述Hermite插值问题的解是存在且唯一的，其表达式为 )(
12
xH
n+
+
′+
jjjjn
xyxyH
0
12
))()(( βα
23
其中
)()21()(
2
0
xl
xx
x
x
j
n
ji
i ij
j
j ∑
≠
=
=α
)()()(
2
xlxxx
jjj
=β
)x
H
i
1
(l
j
为相应于节点的Lagrange插值基函数。
n
xxx,,,
10
L
(2) 两点三次Hermite插值问题,求一个三次Hermite插值多项式，满足
，，。
)(
3
x
i
yxH =)(
3 ii
yxH ′=′ )(
3
)1,0( =i
当时，两点三次Hermite插值问题的解为：,0
10
== xx
)()()()()(
110011003
xyxyxyxyxH ββαα ′+′++=
2
0
)1)(21()(?+= xxxα，，，。
2
1
)23()( xxx?=α
2
0
)1()(?= xxxβ
2
1
)1()( xxx?=β
一般情形下两点三次Hermite插值多项式为
)()()()()(
110011003
xByxByxAyxAyxH ′+′++=
)()(
01
0
00
xx
xx
xA
=α，)()(
01
0
11
xx
xx
x
A
=α，)()(
01
0
00
xx
xx
x
B
= β，)()(
01
0
11
xx
xx
x
= βB
(3) Hermite插值多项式的余项公式
定理3.5 设在上存在且连续，则对于任意的都存在)(
)4(
xf ],[
10
xx ],[
10
xx∈x
),(
10
xx∈ξ，使得
2
1
2
0
)4(
33
)()(
!4
)(
)()()( xxxx
f
xHxfxR=?=
ξ
6,三次样条插值
设定义在[a,b]上的函数f(x)的型值点为(，要求[a,b]上的分段三次插值函数S(x)，使得
),(,),,(),,
1100 nn
yxyxyx L
),,1,0 nL(,)( iyxS
ii
==，且。称为在[a,b]上的三次样条插值函数。
],[)(
2
baCxS ∈ )(xS f
)
)(x
要求满足,(xS
),,1,0(,)( nkyxS
kk
L==
24
)(lim)(lim xSxS
kk
xxxx
′=′
+
→→
，)1,,2,1(?= nk L
)(lim)(lim xSxS
kk
xxxx
′′=′′
+
→→
，( )1,,2,1?= nk L
此外，实际问题对三次样条插值函数在端点处的状态也有要求，即所谓的边界条件,
第一类边界条件(夹持条件),
′=′
′=′
nn
yxS
yxS
)(
)(
00
第二类边界条件(自然边界条件),
′′=′′
′′=′′
nn
yxS
yxS
)(
)(
00
第三类边界条件(周期性条件),
n
yy =
0
，)(lim)(lim
0
xSxS
n
xxxx
′=′
+
→→
，)(lim)(lim
0
xSxS
n
xxxx
′′=′′
+
→→
。
满足上述边界条件之一的三次样条插值多项式是存在且唯一的。构造三次样条插值多项式常用的有三弯矩法和三转角法。
7,最小二乘法
(1) 最小二乘法的基本概念：设( ),
ii
yx ),,1,0( mi L=为给定的一组数据，
i
ω ),,1,0( mi L=为各点的权系数(通常要求
i
ω >0)，要求在函数类
)}(,),(),({
10
xxxspan
n
L=Φ
中，求一函数
∑
=
=
n
j
jj
xaxS
0
**
)()(?，( n )m≤
满足
∑∑
=
Φ∈
=
=?=
m
i
iii
S
m
i
iii
yxSyxS
0
2
0
2*
2
2
))((min))(( ωωδ
其中为中任意函数。
∑
=
=
n
j
jj
xaxS
0
()(? )
S
Φ
称按上述条件求函数的方法为数据拟合的最小二乘法，简称最小二乘法。并称为最小二乘解，为拟合函数。
)(
*
xS
)x
)(
*
x
(S
(2) 对于给定的函数表，),(
ii
yx ),,1,0( mi L=，在函数类
25
)}(,),(),({
10
xxxspan
n
L=Φ中存在唯一的函数，使得关系式
∑
=
=
n
j
jj
xaxS
0
**
)()(?
∑
=
m
i
iii
yxS
0
2
))((ω
),()
kj
fa?= (k
∑
Φ∈
=
=?=
S
m
i
iii
yxS
0
2*
2
2
min))((ωδ
**
1
*
0
,,,
n
aaa L
,(
0
n
j
kj
∑
=
成立；最小二乘解的系数可以通过求解以下法方程组得到,
),,1,0 nL=
=
),(
),(
),(
1
0
1
0
nn
f
f
f
a
a
a
MM
ik
x )(
i
)
=
∑
∑
∑
=
=
=
m
i
i
n
ii
m
i
iii
m
i
ii
n
yx
yx
y
a
a
a
0
0
0
1
0
ω
ω
ω
M
M
),()
),()
),()
1
11
01
nn
n
n
L
MO
L
L
∑
=
=
m
i
ijik
x
0
)()ω?
∑
=
=
m
i
kiik
xy
0
(),?ω?
),,1 nL
∑
∑
∑
=
+
=
+
=
m
i
n
ii
m
i
n
ii
m
i
n
ii
x
x
x
0
21
0
12
0
ω
ω
ω
L
MO
L
L
,(),(
,(),(
,(),(
0
101
000
nn
MM
j
,(?
f(
,0()( jxx
j
j
==?
∑∑
∑∑
∑∑
==
==
==
m
i
n
ii
m
i
n
ii
m
i
ii
m
i
ii
m
i
ii
m
i
i
xx
xx
x
00
00
00
ωω
ωω
ωω
MM
,
即
其中
当取基底为时，相应的法方程组为
注：为避免求解病态的法方程组，往往选取正交多项式为基函数。
26
3.2 典型例题
例3.1 已知函数的函数表如表所示,xy sin=
试分别用拉格朗日线性插值和抛物线插值求的近似值，并估计截断误差。 3367.0sin
x
0.32 0.34 0.36
y sin= 0.314 567 0.333 487 0.352 274 x
解 1) 线性插值,
取两个节点,34.0,32.0
10
== xx相应的,333487.0,314567.0
10
== yy由公式(3.2)
得
314567.0
32.034.0
32.0
314567.0
34.032.0
34.0
)()()(
0
01
0
0
10
1
11001
×
+×
=
+
=+=
xx
y
xx
xx
y
xx
xx
yxlyxlxL
将代入，即得 3367.0=x
330365.03367.0sin ≈,
其截断误差由公式(3.3)得
))((
2
)(
10
1
1
xxxx
M
xR≤
其中,)(max
10
1
xfM
xxx
′′=
≤≤
因,sin)(,sin)( xxfxxf?=′′=可取
,3335.0sinsinmax
11
10
≤==
≤≤
xxM
xxx
于是
5
11
1092.0
)0033.0()0167.0()3335.0(
2
1
)3367.0(3367.0sin)3367.0(
×≤
×××≤
= LR
2) 抛物线插值,取,36.0,34.0,32.0
210
=== xxx由公式(3.2)得
2
1202
10
1
2101
20
0
2010
21
2211002
))((
))((
))((
))((
))((
))((
)()()()(
y
xxxx
xxxx
y
xxxx
xxxx
y
xxxx
xxxx
yxlyxlyxlxL
+
+
=
++=
将和的数值代入，即得 3367.0=x )2,1,0(,=iyx
ii
27
,330374.03367.0sin ≈
这个结果与6位有效数字的正弦函数表完全一致，这说明查表时用二次插值精度已相当高了。其截断误差由公式(3.3)得
))()((
!3
)(
210
2
2
xxxxxx
M
xR≤
其中,828.0cos)(max
02
20
<=′′′=
≤≤
xxfM
xxx
于是
.10178.0
)0233.0()033.0()0167.0()828.0(
6
1
)3367.0(3367.0sin)3367.0(
6
22
×<
××××≤
= LR
应该指出，由插值的截断误差估计式
∏
=
+
≤
n
j
j
n
n
xx
n
M
xR
0
)(
)!1(
)(
可看出，(xR
n
的大小除了与有关外，还与因子
n
M
∏
=
n
j
j
xx
0
)(有关，为了使)(x
n
R尽量小，显然要尽可能选取靠近插值点x的插值节点。
)
在本题中，若取36.0,34.0
21
== xx为节点作线性插值，可得
330387.0)3367.0(
~
3367.0sin
1
=≈ L
其截断误差为 ))((
2
~
)(
~
21
1
2
xxxx
M
xR≤
其中,3523.0)(max
~
10
1
≤′′=
≤≤
xfM
xxx
于是
.1036.1
)0233.0)(0033.0)(3523.0(
2
1
)3367.0(
~
3367.0sin)3367.0(
~
5
11
×≤
≤
= LR
比较以上两种线性插值的节点的不同取法，显然，前者因
∏
=
n
j
j
xx
0
)(比后者小，从而
(
1
xR亦小，相应的插值效果较好。所以，选择插值节点时，应尽量选取与插值点x距离)
28
较近的一些节点。
例3.2 已知对应的函数值为,5,4,3,2,1=x,6,8,7,4,1)( =xf作4次牛顿插值多项式。
解 由给定的数据先作差商表如表3-2所示,
表 3-2 差 商 表
k
k
x )(
k
xf 一阶差商 二阶差商 三阶差商 四阶差商
0 1 1
3
1 2 4 0
3 -1/3
2 3 7 -1 1/24
1 -1/6
3 4 8 -3/2
-2
4 5 6
29
以表格中斜线上对应的数字为系数，写出所求的牛顿插值多项式为
4
43 2
1
( ) 1 ( 1) 3 ( 1)( 2) 0 ( 1)( 2)( 3) ( )
3
1
( 1)( 2)( 3)( 4)
24
1 9 83 33
1
24 12 24 12
Nx x x x x x x
xx xx
xxxx
= +?×+×+×?
+×
=?+?+
例3.3 已知的函数表如表3-3所示用牛顿向前插值公式计算cos的值，并估计余项。
xxf cos)( = 048.0
表 3-3
x
0.0 0.1 0.2 0.3 0.4 0.5 0.6
cos
1.000 00 0.995 00 0.980 07 0.955 34 0.921 06 0.877 58 0.825 34
x
解 先构造差分表(表 3-4),
表 3-4 差分表
i
x )? (
i
xf
i
f?
i
f
2
i
f
3
i
f
4
0.0 1.000 00
-0.005 00
0.1 0.995 00 -0.009 93
-0.014 93 0.000 13
0.2 0.980 07 -0.009 80 0.000 12
-0.024 73 0.000 25
0.3 0.955 34 -0.009 55 0.000 10
-0.034 28 0.000 35
0.4 0.921 06 -0.009 20 0.000 09
-0.043 48 0.000 44
0.5 0.877 58 -0.008 76
-0.052 24
0.6 0.825 34
从差分表中可见，四阶差分已接近于常数，用4次插值多项式)()(
4
xfxN ≈是恰当的。在计算cos时，因为0.048是靠近表头的值，因此用牛顿向前插值公式计算，这时可取 048.0
48.0
1.0
0.0048.0
,1.0,4.0,3.0,2.0,1.0,0.0
0
43210
=
=
=
======
h
xx
t
hxxxxx
公式中用到的各阶差分正好是表3-4中第一条斜线“——”上的对应值。
0
4
0
3
0
2
00
,,,,fffff
于是
4
23
00 0 0 0
cos 0.048 (0.048)
( 1) ( 1)( 2) ( 1)( 2)( 3)
1! 2! 3! 4!
0.48 0.48(0.48 1)
1.00000 ( 0.005 00) ( 0.009 93)
1! 2!
0.48 (0.48 1) (0.48 2)
0.000 13
3!
0.48 (0.48 1) (0.48 2) (0.48 3)
4
N
ttt tt ttt
ff f f f
≈
=+?+?+?+?
=+×? + ×?
×?×?
+×
×?×?×?
+
4
0.000 12
!
1 0.0024 0.1248 0.00993 0.06323 0.00013 ( 0.03984) 0.00012
0.99884
×
=? + × + × +? ×
≈
由余项公式得误差
)6.00.0(),(
!5
)4)(3)(2)(1(
)(
)5(5
4
<<
= ξξfh
ttttt
xR
由于 1sin)(
5
≤?= ξξf
则
65
4
1028.01.0
!5
)448.0()348.0()248.0()148.0(48.0
)(
×=×
×?×?×?×
≤xR
例3.4 已知的函数表如表3-5所示，求函数在[0，2]之间的零点近似值。 )
i
x
(xfy = )(xf
表 3-5 函数表
0 1 2
8 -7.5 -18
i
y
解 当函数的解析式未知时，而仅知其在某区间[ a ]上的函数表，那么，如何)(xf b,
30
31
b y y
x
) g
)
求其在[ ]上的零点呢？一般地，可先求在[ a ]上的插值函数，然后求的零点，把此零点就作为的近似零点。特别是，若的反函数存在，记为于是，求的零点问题就变成求函数值的问题了。若用插值法构造出，从而求得的零点的近似值，一般称其为“反插值”，但用反插值时，必须注意反插值的条件是函数有反函数，也即要求
a,
f
)(xf
y =
)(xf
)0(g
)(xf
b,
f
)(x )(x
),(yg)(xf )(x =
)(x
f
(y
(x y =单调。本题中，因为是按严格单调下降排列，所以可用反插值法求的零点的近似值。具体做法是先把原表变成反函数表(表 3-6)。
i
y
)(xf
i
y
i
x
L 0=y (xf
).(
2
yL≈
445
8
0
.0
(
(
)
≈
)5.7
)5.7
+()8
0()8
×
×
1
)185.
)18
()85.7
()80(
)18
)18
.7
.7
0 ×
+×?
×?
×
+
+
+
+
.0
=
*
≈x 445
16
2
+? xx
()
),()
11
xfx
xfx
ii
′=
)( = xx
b,
5.0
*
x
0
,x )(x
)2
f
,1,0==
H (xH
表 3-6 反函数表
8 -7.5 -18
0 1 2
由上表构造二次插值多项式，且令)(
2
y，即可得在[0，2]之间的零点近似值
)
*
x
具体计算过程为
2
1818(
0(
7(
0
0
8()5
0()5
(
2
×
+?
+
+
+
×
×
L
即 。
事实上，此题的8
2
1
3
f 在[0，2]之间的零点为。 =
应该指出，当所给函数表不满足严格单调的条件时，在利用反插值求根时会出错。
例3.5 已知是[ a ]上三个互异的节点，函数在[ a ]上具有连续的四阶导数。试求满足插值条件
21
,xx b,
)(
((
H
iH
′
的插值多项式，并估计其误差且证明是唯一的。 )
)x )
(x )
解 由给定的4个插值条件，显然可确定一个次数不超过3次的埃尔米特插值多项式
。又由应满足插值条件(H (xH ),2,1,0)(()( == ixfxH
ii
而节点上的二次插
210
,,xxx
32
)值函数也满足插值条件(
2
xN ),2,1,0)(()(
2
== ixfxN
ii
))()((
],,[))((
],[)()(
)()(()(
210
21010
1000
102
xxxxxxA
xxxfxxxx
xxfxxx
xxxxxAxN
+
故可设的形式为 )(xH
+
)(xH )( 0)( =ix
i
xH
i
)])((
,[)(],,[
10
00210
xxxx
xxfxxxxxf
(
[)
x
fx
+
[(xA?+?+=
1
x= )()(
11
xfxH ′=′
],,[)(]
210011
Axxxfxxx )(
10
xx?+?+
))((
)(][)(
2101
01101
xxxx
fxxxxfx
],
2
x
0 01 0
10110012
1012
) ( )[,] ( )(
,]( )[,,]
()()
x x fx x x x
xx x xfxxx
01
01
[,
(
fx x
xxx
xxxx
+? +?
+?
)()()( xHxfxR?=
2
)(xR
1
x
())()(()
2
10
xxxxxxkx=
)(xk
ba,)2,1,0=
)x)()(()() xtxtxktHt
))(
2
f
xxH
+
+
=
=
其中A为待定系数；显然，由上式确定的满足)2,1,(= f，常数A可由插值条件)()(
11
xfxH ′=′来确定。为此，对上式两边求导得
))(
))(],)(],(
20
2121110
xxx
xxxxxxxxH
+
+′
令，并利用插值条件x就有
)(,[)(
2101
xxxfxf?=′
于是
,[
10
xxf
A
′
=
即得所求的插值多项式为 )(xH
0 1 2
2
() ( ),]
() [
() )()
Hx fx x x x
fx f
x xx
=?
′
再来估计误差
由插值条件，显然，都是的零点 (并且是的二重零点)，故可设
10
,,xxx )(xR
)(
2
xR?
其中，为待定函数。
为了求得，暂视)(xk x为区间[ ]上任意固定的异于(ix
i
的点，并作辅助函数
()()(
2
2
10
tftg?=
33
)
b )
显然，具有下列性质,(tg
(1) 在[ ]上，存在四阶导数,a,(tg !4)()()(
)4()4(
×?= xktftg
(2) 0)()()()()(
1210
=′==== xgxgxgxgxg
所以，在(内有5个零点(其中，是二重零点，算2个)，于是类似于余项定理证明，反复应用罗尔定理，最后可得在[ ]内至少有一点
) )(tg,ba
1
x
ba,ξ，使
0)(
)4(
=ξg
而
0!4)()()(
)4()4(
=×?= xkfg ξξ
于是
)(
!4
1
)(
)4(
ξfxk =
从而，余项的表达式为
)())()((
!4
1
)(
2
2
10
)4(
xxxxxxfxR= ξ
其中，ξ介于x与之间。
i
i
x
)
x
当时，上式显然也成立。 x =
再证是唯一的（用反证法）。 (xH
设另有一个三次多项式)(
~
xH也满足插值条件
)2,1,0)(()(
~
== ixfxH
ii
和)()(
~
11
xfxH ′=′,
于是，函数)(
~
)()( xHxHx?=?也是次数不超过三次的多项式，且有
2),1,0( 0)(
~
)()( ==?=?= iffxHxHx
iiiii
0)(
~
)()(
11111
=′?′=?′=′ ffxHxHx?
即)(x?至少有4个零点，但)(x?是次数不超过三次的多项式，所以，只能0)( ≡x?，即
)(
~
)( xHxH ≡,
唯一性得证。
例3.6已知函数)(xfy =的函数表如表3-7所示，在区间[0，0.60]上求三次样条插值函数
，使满足边界条件,)(xS 0)0( =′S 64879.0)60.0(?=′S,
表3-7 函数表
i 0 1 2 3 4
0 0.15 0.30 0.45 0.60
(
i
xf 1 0.97800 0.91743 0.83160 0.73529
i
x
)
解 这是在第一种边界条件下的样条插值。若用三弯矩法求解，则确定
),4 3,2,1,0( )( ==′′ jMxS
jj
的方程组为
=
4
3
2
1
0
4
3
2
1
0
33
22
11
21
2
2
2
12
d
d
d
d
d
M
M
M
M
M
λμ
λμ
λμ
由已知，得 15.0
1
=?=
jjj
xxh
5.0
1
1
=
+
=
+jj
j
hh
h
μ，)3,2,1(
1
1
=
+
=
+
+
j
hh
h
jj
j
j
λ
为确定方程组的右端项，先构造差商表（表3-8）,)4,3,2,1,0( =jd
j
表3-8 差商表
i
i
x )(
i
xf 一阶差商 二阶差商
0
1
2
3
4
0
0,15
0.30
0.45
0.60
1
0,97800
1,0.91743
2,0.83160
3,0.73529
-0.14666
-0.40380
-0.57220
-0.64206
-0.85711
-0.56133
-0.23288
再由公式算得各],,[6
11 +?
=
jjjj
xxxfd )4,3,2,1,0( =jd
j
之值，从而方程组为
=
26893.0
39733.1
36799.3
14266.5
86666.5
21
5.025.0
5.025.0
5.025.0
12
4
3
2
1
0
M
M
M
M
M
34
解此方程组得
04452.2
0
=M，M，77762.1
1
= 13031.1
2
=M，43710.0
3
=M，M,08408.0
4
=
代入三次样条插值函数在各子区间[上的表达式,)(xS ],
1 jj
xx
),,2,1(,)
6
(
)
6
(
6
)(
6
)(
)(
1
2
2
1
1
3
1
1
3
nj
h
xxhM
y
h
xxhM
y
M
h
xx
M
h
xx
xS
j
jjj
j
j
jjj
j
j
j
j
j
j
j
j
L=
+
+
+
=
（2.28）
经整理得所求的三次样条插值函数在区间[上的表达式为 )(xS ]60.0,0
+
++?
++?
+?
=
]60.0,45.0[ 01461.107380.000035.157911.0
]45.0,30.0[ 99720.004230.025837.177023.0
]30.0,15.0[ 99857.002853.021247.171923.0
]15.0,0[ 102226.129655.0
)(
23
23
23
23
xxx
xxx
xxx
xx
xS
例3.7 已知函数)(xfy =的实测数据组如表3-9第式曲线与此数据组拟合。
二，三列所示，试用最小二乘法求多项解 （1）可先描点：把表中的数据描在坐标纸上，如图3-1所示.从图3-1中可以看出，这些
210
axaay ++=
作为已给函数)的拟合函数。
（2）建立法方程组（取
35
),(
ii
yx
点位于一条抛物线的附近。因此可取二次曲线
2
x
(xfy =
1)( =
i
xω）
∑∑
∑∑
∑∑
==
==
==
7
1
3
7
1
2
7
1
2
7
1
7
1
7
1
1
xx
xx
x
i
i
i
i
i
i
i
i
i
i
i
∑
∑
∑
=
=
=
2
1
0
7
1
4
7
1
3
7
1
2
a
a
a
x
x
x
i
i
i
i
i
i
=
∑
∑
∑
=
=
=
7
1
2
7
1
7
1
i
ii
i
ii
i
i
yx
yx
y
0 2 4 6 8
2
3
4
5
6
7
图3-1
36
i
i
i
i
x
i
x
由已知数据计算，，，，，，,
∑
=
7
1i
x
∑
=
7
1
2
i
i
x
∑
=
7
1
3
i
i
x
∑
=
7
1
4
i
i
x
∑
=
7
1i
i
y
i
i
i
yx
∑
=
7
1
∑
=
7
1
2
ii
yx
计算结果列在表3-9中。
表3-9 数值表
x
i
y
2
i
x
3
i
x
4
i
x
ii
y
i
y
2
1
2
3
4
5
6
7
1
2
3
4
5
6
7
2
3
6
7
5
3
2
1
4
9
16
36
49
64
1
8
27
64
216
343
512
1
16
81
256
1296
2401
4096
2
6
18
28
30
21
16
2
12
54
112
180
147
128
∑
31 28 179 1171 8147 121 635
将表3-9中算得的结果代入方程，即得a，，a应满足的线性代数方程组,
0 1
a
2
=+
=++
=++
63581471171179
121117117931
28179317
210
210
210
aaa
aaa
aaa
(3)求解法方程得
=
=
=
3864.0
4321.3
3185.1
2
1
0
a
a
a
故所求拟合曲线为
2
3864.04321.33185.1 xxy?+?=
它最佳地拟合了表3-9给出的数据。
例3.8 设已知一组实验数据如表3-10中的第二，三列所示.试从这组数据出发，建立变量x与之间的经验公式,y
i
i
x
表3-10 实验数据
x
i
y
i
u
2
i
x
ii
u
1
2
3
4
5
6
1
2
3
4
5
6
15.3
20.5
27.4
36.6
49.1
65.6
1.1847
1.3118
1.4378
1.5635
1.6911
1.8169
1
4
9
16
25
36
1.1847
1.6236
4.3134
6.2540
8.4555
10.9014
7
8
7
8
87.8
117.6
1.9435
2.0704
49
64
13.6045
16.5632
∑
=
8
i 1
36 419.9 13.0197 204 63.9003
解,解题步骤如前，此数据点分布图如图3-2所示，它接近一指数曲线，因而可取指数函数a，b待定常数）作为拟合函数。
可
bx
为
是，这是一个关于的非线性模型，不是我因此，首先应通过适当的变换，把化为线性模型，然后再利用最小二乘法来求解。
为此，对函数两边取常数对数得,
e y
aey =（
37
a，b
们上面讨论的线性无关函数组的线性组合，
bx
bx
aey =
aey =
bxay lglglg +=，令u lg=，aA lg=，
，则得到了线性模型,
BxAu +=
以下，可类似于例3.7的计算，先列出表3-10
3-10
求解得
ebB lg=
中第四，五，六列数据，再利用数据表
9003.63
0197.
0 2 4 6 8
0
50
100
150
图3-2
，可直接写出法方程,
=+
=+
20436
13368
BA
BA
0583.1=A，1265.0=B
b再算出 44.11=a，2913.0=
因此，经验公式为
x
ey
2912.0
44.11=
从这个例题可知，虽然原来的变量之间不是线性关系，但经过变量替换后，仍可化为线性关
（表3-11）取自于一个多项式，试确定这个多项式。
系计算。
例3.9 下列数据表3-11
x -2 -1 1 2 3 0
)x(P -43 -13 1 5 5 7
解 由差商的性质知，只要多项式的某阶差商（或差分）均为同一常数，则该多项式的次数差商
必同于该阶差商（或差分）的阶数。由此，作差商表如表3-12所示,
表3-12 差商表
x
1 2 3 4
)(xf
30 -2 -43
38
-1 -13
-8
0
1
2
3
1
5
5
7
14
4
0
2
-5
-2
1
1
1
1
0
0
由此表可看出，函数的三阶差商为常数，应是一个三次多项式，故可从已知数据中任取4组数据用牛顿插值或拉格朗日插值构造出该三次多项式为
例3.10 考虑下列插值问题：求一个2次多项式使得
)(xf )(xf
185)(
23
++?= xxxxP,
)(xp
22100
)(,)(,)( yxpmxpyxp ==′=
其中
1020
,,,ymyxx ≠
2
为已知数据．试给出使这一问题的解存在惟一的条件,
解 则设,)(
2
cxbxaxp ++= cxbxp 2)( +=′．要存在惟一的2次多项式p
插值条件(1)，当且仅当下列方程组,
存在惟一解,
阵的行列式为
)(x满足
=++
=+
=++
2
2
22
1
0
2
00
2
ycxbxa
mcxb
ycxbxa
方程组的系数矩
)2)((
0
210
1
1
210
1
10202
2
0
2
202
1
2
00
13
2
22
1
2
00
xxxxx
xxxx
x
xx
rr
xx
x
xx
D
+?=
=
由于,
20
≠≠ Dxx 0的充分必要条件为)(
2
1
201
xxx +≠．所以插值问题存在惟一解的充分必要条件为)(
20
xx +
2
1
1
x ≠,
例3.11 ]，求证 设,[)(
2
baCxf ∈
)(max)(
8
1
)](
)()(
)([)(max
2
xfabax
ab
afbf
afxf
bxabxa
′′?≤?
+?
≤≤≤≤
证 b两点为插值节点作的一次插值多项式，利用插值多项式的余项定理，以a,)(xf
得到
)),max(),,(min(),)()((
2
1
)](
)()(
)([)(
bxaxbxaxf
ax
ab
afbf
afxf
∈′′=
+?
ξξ
当时,],[ bax∈
2
)(
4
1
))((max,)(max)( abbxaxxff
bxabxa
≤′′≤′′
≤≤≤≤
ξ
得到
],[,)(max)(
8
1
))((max)(max
2
1
)](
)()(
)([)(
2
baxxfab
bxaxxf
ax
ab
afbf
afxf
bxa
bxabxa
∈′′?≤
′′≤
+?
≤≤
≤≤≤≤
从而
)(max)(
8
1
)](
)()(
)([)(max
2
xfab
ax
ab
afbf
afxf
bxa
bxa
′′?≤
+?
≤≤
≤≤
例3.12 已知插值节点
3210
xxxx <<<。证明当3,2,1,
1
==?
ihxx
ii
时
1) 二次插值多项式的误差界为
)(max
27
3
)()(max
2020
3
2
xfhxLxf
xxxxxx
′′′≤?
≤≤≤≤
2) 三次插值多项式的误差界为
)(max
24
1
)()(max
)4(4
3
3030
xfhxLxf
xxxxxx ≤≤≤≤
≤?
证 令 thxx +=
0
1)
202102
),)()()((
6
1
)()( xxxxxxxxfxLxf <<′′′=? ξξ
39
)(max)2)(1(max
6
1
))()((max)(max
6
1
)()(max
20
2020
20
20
3
210
2
xfttth
xxxxxxxf
xLxf
xxxt
xxxxxx
xxx
′′′=
′′′≤
≤≤≤≤
≤≤≤≤
≤≤
考虑函数
20),2)(1()( ≤≤= ttttt?
由0)( =′ t?即3解得 026
2
=+? tt
3
1
1,
3
1
1
21
+=?= tt
由0)2(,
9
32
)(,
9
32
)(,0)0(
21
=?=== tt知
9
32
)(max
20
=
≤≤
t
t
得
)(max
27
3
)()(max
2020
3
2
xfhxLxf
xxxxxx
′′′≤?
≤≤≤≤
2) ))()()((
!4
)(
)()(
3210
)4(
3
xxxxxxxx
f
xLxf=?
ξ
)(max)3)(2)(1(max
24
))()()((max)(max
24
1
)()(max
)4(
30
4
3210
)4(
3
30
3030
30
xftttt
h
xxxxxxxxxf
xLxf
xxxt
xxxxxx
xxx
≤≤≤≤
≤≤≤≤
≤≤
=
≤
作变换,
2
3
= ts则
2
3
3,
2
1
2,
2
1
1,
2
3
==?+=?+= stststst,
因而
40
1)
4
9
)(
4
1
(max)
4
9
)(
4
1
(max
)
4
9
)(
4
1
(max
)
2
3
)(
2
1
)(
2
1
)(
2
3
(max
)3)(2)(1(max
4
9
0
22
2
3
0
22
2
3
2
3
2
3
2
3
30
===
=
++=
≤≤≤≤
≤≤?
≤≤?
≤≤
zzss
ss
ssss
tttt
ss
s
s
t
)(max
24
)()(max
)4(
4
3
3030
xf
h
xLxf
xxxxxx ≤≤≤≤
≤?
例 3.13 设为个互异的节点，l为拉格朗日基本插值多项式，即
n
xxx,,,
10
L 1+n )(x
i
∏
≠
=
=
n
ij
j
i
xl )(
ji
j
xx
xx
0
，试证,
1) nkxxlx
k
n
i
i
k
i
,,1,0,)(
0
L==
∑
=
2) nkxlxx
n
i
j
k
i
,,2,1,0)()(
0
L==?
∑
=
3) 设为任一个最高次项系数为1的)(xp 1+n次多项式，则
∏∑
==
=?
n
i
i
n
i
ii
xxxlxpxp
00
)()()()(
证 1) 记．则nkxx
k
k
,,1,0,)( L==? )(x
k
以为插值节点的拉格朗日插值多项式为且有
n
xxx,,,
10
L
∑
=
n
i
iik
xlx
0
)()(?
0)(
)!1(
)(
)()()(
0
)1(
0
=?
+
=?
∏∑
=
+
=
n
i
i
nn
i
iikk
xx
n
xlxx
ξ?
即
0)(
0
=?
∑
=
n
i
i
k
i
k
xlxx
因而
41
nkxxlx
k
n
i
i
k
i
,,1,0,)(
0
L==
∑
=
2) 记则,,,2,1,)()( nktxx
k
k
L=?=ψ )(x
k
ψ以为插值节点
n
xxx,,,
10
L
的拉格朗日插值多项式为且有
∑
=
n
i
iik
xlx
0
)()(ψ
0)(
)!1(
)(
)()()(
0
)1(
0
=?
+
=?
∏∑
=
+
=
n
i
i
nn
i
iikk
xx
n
xlxx
ξψ
ψψ
即
0)()()(
0
=
∑
=
n
i
i
k
i
k
xltxtx
因而
nktxxltx
k
n
i
i
k
i
,,2,1,)()()(
0
L=?=?
∑
=
在上式中令xt =，即得
nkxlxx
n
i
i
k
i
,,2,1,0)()(
0
L==?
∑
=
3) 由插值余项公式有
∏∏
∏∑
==
=
+
=
=?
+
+
=
+
=?
n
i
i
n
i
i
n
i
i
nn
i
ii
xxxx
n
n
xx
n
p
xlxpxp
00
0
)1(
0
)()(
)!1(
)!1(
)(
)!1(
)(
)()()(
ξ
例3.14 若有个互异的零点，则有
nn
nn
axaxaxaxf ++++=
1
1
10
)( L n
n
xxx,,,
21
L
=
≤≤
=
′
∑
=
1,
1
20,0
)(
0
1
nk
a
nk
xf
x
n
j j
k
j
证 由于是的n个互异的零点，所以
n
xxx,,,
21
L )(xf
∏
=
==
n
i
in
xxaxxxxxxaxf
1
0210
)()())(()( L
42
∏
≠
=
=′
n
ji
i
ijj
xxaxf
1
0
)()(
∑
∏
∑
∏
∑
=
≠
=
=
≠
=
=
=
=
′
n
j
n
ji
i
ij
k
j
n
j
n
ji
i
ij
k
j
n
j j
k
j
xx
x
a
xxa
x
xf
x
1
1
01
1
0
1
)(
1
)(
)(
记，则
k
k
xxg =)(
=?
≤≤
=
1,)!1(
20,0
)(
)1(
nkn
nk
xg
n
于是得
=
≤≤
=
==
=
′
=
≠
=
=
∑
∏
∑
1,
1
20,0
)!1(
)(1
],,,[
1
)(
)(
1
)(
0
)1(
0
21
0
1
1
0
1
nk
a
nk
n
g
a
xxxg
a
xx
xg
axf
x
n
k
nk
n
j
n
ji
i
ij
jk
n
j
j
k
j
ξ
L
例3.15 给出的等距节点函数表，如用线性插值计算的近似值，使其截断误差为
xxf cos)( = )(xf
5
10
×
2
1
，则函数表的步长应取多大？
解 由于是一个周期为xxf cos)( = π2的函数，只要给出一个周期内的数据即可,
设给定如下等距节点,
M
hMiihx
i
π2
,0,=≤≤=
则任给]2,0[ π∈x，一定存在i使得
],[
1+
∈
ii
xxx
以为节点作的一次插值多项式
1+
,
ii
xx )(xf
ii
i
i
ii
i
i
xx
xx
xf
xx
xx
xfxL
+
=
+
+
+
+
1
1
1
1
1
)()()(
43
44
) )并以作为的近似值．由插值余项定理得 (
1
xL (xf
),(),)()((
2
1
)()(
111 ++
∈′′=?
iiiiii
xxxxxxfxLxf ξξ
注意到
2
1
4
1
))((max,1)(max
11
hxxxxxf
ii
xxxxxx
iiii
≤≤′′
+
≤≤≤≤
++
可得
2
11
8
1
))((max)(max
2
1
)()(
11
hxxxxxfxLxf
ii
xxxxxx
iiii
≤′′≤?
+
≤≤≤≤
++
故只要
52
10
2
1
8
1
×≤h
即有
5
1
10
2
1
)()(
×≤? xLxf
解得
3
10102
×≤h
即当步长
3
10102
×≤h时截断误差不超过
5
10
2
1
×,
例3.16 求超定方程组
=+
=+
=?
=+
72
62
353
1142
21
21
21
21
xx
xx
xx
xx
的最小二乘解，并求误差平方和,
解 方程组写成矩阵形成为
=
7
6
3
11
12
21
53
42
2
1
x
x
正规方程组为
=
7
6
3
11
1254
2132
12
21
53
42
1254
2132
2
1
x
x
即
=
48
51
463
318
2
1
x
x
解得
2418.1,0403.3
21
== xx
误差平方和
34066.0
)3224.0(4761.00881.0)0478.0(
)27()26()533()4211(
2222
2
21
2
21
2
21
2
21
=
+++?=
+++= xxxxxxxxE
45
第四章 数值积分与微分
4.1 内容提要
1,插值型求积公式
(1) 考虑积分
∫
=
b
a
dxxffI )()( (4.1)
给定一组节点bxxxa
n
≤<<<≤ L
10
.作的n次插值多项式 )(xf
∑
=
=
n
i
iin
xlxfxL
0
)()()( (4.2)
用代替(4.1)右端的被积函数得到 )(xL
n
)(xf
∑∑
∫∫
==
===
n
i
iii
n
i
b
a
i
b
a
nn
xfAxfdxxldxxLfI
00
)()(])([)()( (4.3)
其中
nidxxlA
b
a
ii
,,1,0,)( L==
∫
称(4.3)为插值型求积公式。
(2) 若取,1,,
10
=== nbxax则得到求积公式
)]()([
2
)( bfaf
ab
fT +
=
称T为梯形公式。 )( f
若取,2,,
2
,
210
==
+
== nbx
ba
xax则得到求积公式
)]()
2
(4)([
6
)( bf
ba
faf
ab
fS +
+
+
=
称为辛卜生公式。 )( fS
关于它们的截断误差有公式
),(),()
2
(
180
)()(
),(),(
12
)(
)()(
)4(4
3
baf
abab
fSfI
baf
ab
fTfI
∈
=?
∈′′
=?
ηη
ηη
47
(3) 如果一个求积公式
∑
=
=
n
i
iin
xfAfI
0
)()( (A为常数，i
i
nL,1,0= ) (4.4)
对于次数≤的多项式均能准确成立，但至少对一个m 1+m次多项式不准确成立，则称该求积公式具有次代数精度。 m
m
定理 6.1 求积公式(6.4)具有次代精度的充分必要条件为该分式对精确成立，而对不精确成立。
m
= xxxxf,,,1)( L=
1
)(
+m
xf
2,复化求积公式
(1) 将积分区间[作等分，记],ba n,,,1,0,,niihax
n
ab
h
i
L=+=
=则
∑
∫∫
=
+
==
1
0
1
)()()(
n
i
x
x
b
a
i
i
dxxfdxxffI
对每一个积分用梯形公式得到
∫
+1
)(
i
i
x
x
dxxf
∑
=
+
+=
1
0
1
)]()([
2
)(
n
i
iin
xfxf
h
fT
称T为复化梯形公式；若对每一个积分用辛卜生公式，得到 )( f
n
∫
+1
)(
i
i
x
x
dxxf
∑
=
+
+
++=
1
0
1
2
1
)]()(4)([
6
)(
n
i
i
i
in
xfxfxf
h
fS
其中).(
2
1
1
2
1 +
+
+=
ii
i
xxx称为复化辛卜生公式。从计算考虑可将上面两式改写为 )( fS
n
6/])(4)(2)()([)(
])())()((
2
1
[)(
1
0
2
1
1
1
0
1
1
0
∑∑
∑
=
+
=
=
+++=
++=
n
i
i
n
i
inn
n
i
inn
xfxfxfxfhfS
xfxfxfhfT
(2) 复代梯形公式和复化辛卜生公式的截断误差分别为
),(),()
2
(
180
)()(
),(),(
12
)()(
)4(4
2
baf
hab
fSfI
bafh
ab
fTfI
n
n
∈
=?
∈′′
=?
ηη
ηη
(3) 如果一个复化求积公式有估计式 )( fI
n
48
)()()(
p
n
hOfIfI =?
49
)
)
则称是p阶的。 ( fI
n
(4) T和T之间有如下关系,( f
n
)(
2
f
n
n
ab
hxfhfTfT
n
i
i
nn
=+=
∑
=
+
],)()([
2
1
)(
1
0
2
12
(5) 和T之间有如下关系,) )( fS
n
( f
n
)]()(4[
3
1
)(
2
fTfTfS
nnn
=
3,外推算法
由
)]()([
12
1
)(
12
1)()(
1
0
2
bfaffh
h
fTfI
n
i
i
n
′?′→′′?=
∑
=
η
知当适当小时有 h
2
)()( chfTfI
n
≈?
于是
)]()([
4
1
)
2
()()(
2
2
fTfI
h
cfTfI
nn
≈≈?
或者
)]()([
3
1
)()(
22
fTfTfTfI
nnn
≈? (4.5)
由此可见，只要当二分前后两个积分值T和T相当接近，就可以保证T的误差很小，且大致等于
)( f
n
)(
2
f
n
)(
2
f
n
)].()([
3
1
2
fTfT
nn
(4.5)是一个误差的事后估计式。于是如果用这个误差作为T的一种修正,)(
2
f
n
可望得到的结果
)(
3
1
)(
3
4
)]()([
3
1
)()(
~
222
fTfTfTfTfTfT
nnnnnn
=?+=
会更加精确。事实上),()(
~
fSf
nn
=T即
)(
3
1
)(
3
4
)(
2
fTfTfS
nnn
= (4.6)
这种由低阶精度公式T和T作线性组合，得到一个高阶精度公式的方法称)(
2
f
n
)( f
n
)( f
n
S
为外推方法。
类似地，由四阶精度公式和作线性组合，可得六阶精度公式 )(
2
fS
n
)( fS
n
)(
15
1
)(
15
16
)(
2
fSfSfC
nnn
= (4.7)
由六阶精度公式和C作线性组合，可得八阶精度公式 )(
2
fC
n
)( f
n
)(
63
1
)(
63
64
)(
2
fCfCfR
nnn
= (4.8)
称和分别为柯特斯公式和龙贝格公式。 ) )
S
)
)
( fC
n
( fR
n
应用(4.6)，(4.7)和(4.8)将粗糙的梯形值逐步加工成精度较高的辛卜生值，
柯特斯值C和龙贝格值的方法，称为龙贝格求积法。
)( fT
n
)( f
n
( f
n
)( fR
n
4,数值微分
根据函数在一些离散点上的函数值推算它在某点处导数值的近似值的方法称为数值微分。
给定函数,并设为的一个n次插值多项式，称 (xf )(xp
n
)(xf
)()( xpxf
n
′≈′
为插值型求导公式。特别，称
)]()([
2
1
),(
000
hxfhxf
h
hxD+=
为求的中点公式，并有误差估计式 )(
0
xf ′
),(),(
6
),()(
00
2
00
hxhxf
h
hxDxf +?∈′′′?=?′ ξξ
4.2 典型例题
例4.1用梯形公式和辛卜生公式计算积分
∫
1
0
dxe
x
并估计误差,
解:记则．,)(,1,0
x
exfba
===
xxxx
exfexfexfexf
=?=′′′=′′?=′ )(,)(,)(,)(
)4(
∫
=
1
0
)( dxefI
x
50
6839.0)(
2
01
)]()([
2
)(
10
=+
=+
=
eebfaf
ab
fT
由
)1,0(,
12
1
)(
12
)(
)()(
3
∈?=′′
=?
ηη
η
ef
ab
fTfI
得
08333.0
12
1
)()( =≤? fTfI
6323.0)4(
6
01
)]()
2
(4)([
6
)(
15.00
=++
=
+
+
+
=
eee
bf
ba
faf
ab
fS
由
)1,0(,)
2
1
(
180
1
)()
2
(
180
)()(
4)4(4
∈?=
=?
ηη
η
ef
abab
fSfI
得
0003472.0
2880
1
)
2
1
(
180
1
)()(
4
==×≤? fSfI
例 4.2 求近似公式
)]
4
3
(2)
2
1
()
4
1
(2[
3
1
)(
1
0
fffdxxf +?≈
∫
的代数精度,
解 当时,1)( =xf
左边＝1,右边＝1)212(
3
1
=+?，左边＝右边；
当时,xxf =)(
左边＝
2
1
，右边＝,
2
1
)
4
3
2
2
1
4
1
2(
3
1
=×+?×左边＝右边；
当时,
2
)( xxf =
左边＝
3
1
，右边＝,
3
1
])
4
3
(2)
2
1
()
4
1
(2[
3
1
222
=×+?×左边＝右边；
当时,
3
)( xxf =
左边＝
4
1
，右边＝,
4
1
])
4
3
(2)
2
1
()
4
1
(2[
3
1
333
=×+?×左边＝右边；
当时,
4
)( xxf =
51
左边＝
5
1
，右边＝,
192
37
])
4
3
(2)
2
1
()
4
1
(2[
3
1
444
=×+?×左边≠右边；
因而所给求积公式具有3次代数精度,
例 4.3 给定求积公式
(1) )()0()()(
2
2
hCfBfhAfdxxf
h
h
++?≈
∫
试决定使它的代数精度尽可能地高,C
1
BA,,
解 当时，左边，右边＝)( =xf hdx
h
h
41
2
2
==
∫
CBA ++；
当时，左边,xxf =)( 0
2
2
==
∫
h
h
xdx
右边＝hCAhCBhA )(0)( +?=×+×+?×；
当时，左边
2
)( xxf =
3
2
2
2
3
16
hdxx
h
h
==
∫
,
右边＝,
2222
)(0)( hCAhCBhA +=×+×+?×
要使求积公式至少具有2次代数精度，其充分必要条件为满足如下方程组,CBA,,
=+
=+?
=++
32
3
16
)(
0)(
4
hhCA
hCA
hCBA
解得
hChBhA
3
8
,
3
4
,
3
8
=?==
代入(1)得
)](2)0()(2[
3
4
)(
2
2
hffhf
h
dxxf
h
+≈
∫
(2)
当时，(2)的左边＝，(2)的右边＝
3
)( xxf = 0
2
2
3
=
∫
h
h
dxx
0]2(2[
3
4
3
=×× h
h
0)
33
+h，左边＝右边；
当时，(2)的左边＝
4
)( xxf =
5
2
2
4
5
64
hdxx
h
h
=
∫
，(2)的右边＝
54
3
16
](2[
3
4
hh
h
=××
44
20)h +，左边≠右边,
所以当求积公式(1)中求积系数取为hChBhA
3
8
,
3
4
,
3
8
=?==时得到求积公式(2)，其
52
代数精度取到最高，此时代数精度为3,
例 4.4 求三个不同的节点和常数c使求积公式
321
,,xxx
)]()()([)(
210
1
1
xfxfxfcdxxf ++≈
∫
(1)
具有尽可能高的代数精度,
解 不妨设
210
xxx <<,
当时，左边＝2,右边＝c1)( =xf c3)111( =++；
当时，左边＝0,右边＝cxxf =)( )(
210
xxx ++；
当时，左边＝
2
)( xxf =
3
2
,右边＝； )(
2
2
2
1
2
0
xxxc ++
当时，左边＝0,右边＝c,
3
)( xxf = )(
3
2
3
1
3
0
xxx ++
要使求积公式至少具有3次代数精度，当且仅当满足如下线性代数方程组
210
,,,xxxc
=++
=++
=++
=
)5(0)(
)4(
3
2
)(
)3(0)(
)2(23
3
2
3
1
3
0
2
2
2
1
2
0
210
xxxc
xxxc
xxxc
c
由(2)得
3
2
=c
于是(3)-(5)变为
=++
=++
=++
)8(0
)7(1
)6(0
3
2
3
1
3
0
2
2
2
1
2
0
210
xxx
xxx
xxx
由(6)得
)(
201
xxx +?= (9)
将其代入(7)和(8)得
=++?
=+++
0)(
1)(
3
2
3
20
3
0
2
2
2
20
2
0
xxxx
xxxx
即
53
2
1
2
220
2
0
=++ xxxx (10)
0)(
2020
=+ xxxx (11)
由(11)得，或或0
0
=x 0
2
=x 0
20
=+ xx,
当时，由(10)得0
0
=x
2
1
2
=x，再由(9)得
2
1
1
=x，与
10
xx <矛盾；
当时，由(10)得0
2
=x
2
1
0
=x，再由(9)得
2
1
1
=x，与
21
xx <矛盾；
当时，由(10)得0
20
=+ xx
2
1
2
=x，因而
2
1
0
=x，再由(9)得0
1
=x,
综上我们得到要使求积公式(1)至少具有3次代数精度当且仅当取
2
1
,0,
2
1
,
3
2
210
==?== xxxc
此时求积公式为
)]
2
1
()0()
2
1
([
3
2
)(
1
1
fffdxxf ++?≈
∫
(12)
当时，(12)的左边
4
)( xxf =
5
2
1
1
4
==
∫
dxx，(12)的右边,
3
1
])
2
1
(0)
2
1
[(
3
2
444
=++?=
左边右边，所以求积公式(12)的代数精度为3,≠
例4.6设，且],[)( baCxf ∈ )(af ′与)(bf ′存在．是满足的三次样条插值函数．记
)(xS
)() bfb ′=(),()( SafaS ′′=′
n
b
hn
a
iihax
i
=,,,2,1 L=+=,0,．试证
)]()([
12
])(2)()([
2
)(
21
1
bfaf
h
xfbfaf
h
dxxS
n
i
i
b
a
′?′+++=
∑
∫
=
证 的三次样条插值函数的分段表达式为 )(xf
54
1,,2,1,0),,(
)
6
1
()
6
1
(
)(
6
1
)(
6
1
)(
1
2
11
12
3
1
3
1
=∈
+
+
+?=
+
++
+
++
nixxx
h
xx
hMy
h
xx
hMy
xxM
h
xxM
h
xS
ii
i
ii
i
ii
iiii
L
(1)
其中为下列方程组的解,
n
MMM,,,
10
L
=
n
n
n
n
n
n
d
d
d
d
d
d
M
M
M
M
M
M
1
2
2
1
0
1
2
2
1
0
21
2
1
2
2
1
2
1
2
2
1
2
1
2
2
1
2
1
2
2
1
12
MMOOO (2)
1,,2,1,
2
],[],[
6
11
=
=
+
ni
h
xxfxxf
d
iiii
i
L
h
xxfbf
d
h
afxxf
d
nn
n
],[)(
6,
)(],[
6
110
0
′
=
′?
=
nixfy
ii
,,2,1,0),( L==．将(2)中的第一个方程和最后一个方程同乘以
2
1
，并将结果和
(2)中其它各方程相加，得
)]()([
3
2
3
3333
2
3
12210
afbf
h
MMMMMM
nnn
′?′=
++++++
L
即
)]()([
2
2
1
1
0
afbf
h
MMM
n
n
i
i
′?′=++
∑
=
(3)
对(1)的两边在区间[ ]积分得
1
,
+ii
xx
55
1,,2,1,0,)(
24
1
)(
2
1
)
6
(
2
1
)
6
(
2
1
24
1
24
1
)
6
()
6
(
)(
6
)(
6
)(
3
11
21
1
23
1
3
21
1
12
313
1
11
111
=+?+=
+?++=
+
+
+?=
++
+
++
+
+
+
+
+
∫∫
∫∫∫
++
+++
nihMMhyy
hh
M
yhh
M
yhMhM
dx
h
xx
h
M
ydx
h
xx
h
M
y
dxxx
h
M
dxxx
h
M
dxxS
iiii
i
i
i
iii
x
x
ii
i
x
x
ii
i
x
x
i
i
x
x
i
i
x
x
i
i
i
i
i
i
i
i
i
i
L
将上式对i从0到1?n求和，得
)2(
24
1
)2(
2
1
)(
24
1
)(
2
1
)()(
1
1
0
3
1
1
0
1
0
1
3
1
0
1
1
0
1
n
n
i
in
n
i
i
n
i
ii
n
i
ii
n
i
x
x
b
a
MMMhyyyh
MMhyyhdxxSdxxS
i
i
++?++=
+?+==
∑∑
∑∑∑
∫∫
=
=
=
+
=
+
=
+
应用(3)得
)]()([
12
1
])(2)()([
2
)]()([
12
1
)]()(2)([
2
1
)(
2
1
2
1
bfafhxfbfaf
h
afbfhbfxfafhdxxS
n
i
i
n
i
i
b
a
′?′+++=
′?′?++=
∑
∑
∫
=
=
例 4.7 用复化梯形公式计算积分
∫
+
=
1
0
)1( xx
dx
I
精确至3位有效数,
解 所给积分I的被积函数
xx
xf
)1(
1
)(
+
=
在的附近具有奇性．首先应用变量代换对进行变形,0=x )( fI
∫∫∫
+
=
+
=
+
=
1
0
1
0
2
1
0
1
2
1
2
)1(
dt
tx
xd
xx
dx
I
现在被积函数变为
2
1
2
)(
t
tg
+
=
应用复化梯形公式计算得
56
55.1)
5.01
2
5.1(
2
1
)]5.0(1[
2
1
5.1)
11
2
01
2
(
2
1
)]1()0([
2
01
2
12
22
1
=
+
+=×+=
=
+
+
+
=+
=
gTT
ggT
5695.1
)]
875.01
2
625.01
2
375.01
2
125.01
2
(
4
1
5656.1[
2
1
))]875.0()625.0()375.0()125.0((
4
1
[
2
1
5656.1)]
75.01
2
25.01
2
(
2
1
55.1[
2
1
))]75.0()25.0((
2
1
[
2
1
2222
48
22
24
=
+
+
+
+
+
+
+
+=
+++×+=
=
+
+
+
×+=
+×+=
ggggTT
ggTT
因为
2
48
10
2
1
×≤?TT
所以
57.1≈I
例 4.8 给定积分．当要求误差小于10时，用复化梯形公式及复化辛卜生公式计算时所需节点数及步距分别为多少？
∫
3
1
sin xdxe
x 6?
解 设．设将区间
［1,3］作n等分，
∫
=====
b
a
x
dxxffIbaxexf )()(,10,3,1,sin)(
6
ε
n
h
2
=．则
)(
1
3
2
)()
2
(
12
13
)(
12
)()(
2
22
η
ηη
f
n
f
n
fh
ab
fTfI
n
′′××?=
′′
=′′
=?
(1)
)(
1
90
1
)()
1
(
180
13
)()
2
(
180
)()(
)4(
4
)4(4)4(4
η
ηη
f
n
f
n
f
hab
fSfI
n
×?=
=
=?
(2)
对求导得 )(xf
)
4
sin(2)cos(sin)(
π
+=+=′ xexxexf
xx
57
xexexxexf
xxx
cos2)
2
sin(2)]
4
cos()
4
[sin(2)( =+=+++=′′
πππ
)
4
3
sin()2(
)]
2
cos()
2
[sin(2)(
3
π
ππ
+=
+++=′′′
xe
xxexf
x
x
xexe
xxexf
xx
x
sin4)sin(4
)]
4
3
cos()
4
3
[sin()2()(
3)4(
=+=
+++=
π
ππ
因而
3
3131
2cos2max)(max exexf
x
xx
≤=′′
≤≤≤≤
(3)
3
31
)4(
31
4sin4max)(max exexf
x
xx
≤?=
≤≤≤≤
(4)
由(1)和 (3)，得
2
3
3
2
3
4
2
1
3
2
)()(
n
e
e
n
fTfI
n
=××≤?
当
ε<
2
3
3
4
n
e
(5)
时
ε<? )()( fTfI
n
由(5)解得
04.492610
3
2
3
4
3
3
=×=>
e
e
e
n
ε
因而取，即取4928个节点，步长为4927=n
4927
2
时，用复化梯形公式计算误差小于10,
6?
由(2)和(4)，得
4
3
3
4
45
2
4
1
90
1
)()(
n
e
e
n
fSfI
n
=××≤?
当
ε<
4
3
45
2
n
e
(6)
时
58
ε<? )()( fSfI
n
由(6)解得
74.30
45
2
4
3
=>
ε
e
n
因而取，即取32个节点，步长为31=n
31
2
时，用复化辛卜生公式计算误差小于10,
6?
例 4.11 用龙贝格求积法求积分
∫
+
=
1
0
2
1
4
dx
x
I
的近似值，要求误差不超过
5
10
2
1
×,
解 记1,0,
1
4
)(
2
==
+
= ba
x
xf,
2
11
4
)1(,4
01
4
)0(
22
=
+
==
+
= ff
3)]1()0([
2
01
)]()([
2
1
=+
=+
= ffbfaf
ab
T
2.3
5.01
4
)5.0(
2
=
+
=f
1.3)2.33(
2
1
)]5.0([
2
1
12
=+=+= fTT
3133333.3)4(
3
1
121
=?= TTS
5600000.2
75.01
4
)75.0(,7647059.3
25.01
4
)25.0(
22
=
+
==
+
= ff
1311765.3))]75.0()25.0((
2
1
[
2
1
24
=+×+= ffTT
3
12242
10549.0
15
1
,1415687.3)4(
3
1
×=?=?= SSTTS
1421179.3)16(
15
1
121
=?= SSC
9384615.3
125.01
4
)125.0(
2
=
+
=f
5068493.3
375.01
4
)375.0(
2
=
+
=f
8764045.2
625.01
4
)625.0(
2
=
+
=f
2654867.2
875.01
4
)875.0(
2
=
+
=f
59
1389885.3))]875.0()625.0()375.0()125.0((
4
1
[
2
1
48
=+++×+= ffffTT
1415925.3)4(
3
1
484
=?= TTS
55
24
10
2
1
101587.0
15
1
×<×=? SS
所以
14159.3≈I
其误差不超过
5
10
2
1
×,
例 4.14 用中点差商公式
h
hxfhxf
xf
2
)()(
)(
00
0
+
≈′ (1)
计算xxf =)(
.0,05.0,5
在处的一阶导数的近似值．取
，观测计算结果，解释所发生的现象．设计高精度的算法，并重新计算,
2=x
0005.0 00005.0,,005.0=h
解 对于xxf =)(，中点差商公式为
h
hh
f
2
22
)2(
+
≈′ (2)
在计算过程中保留6位有效数字．计算结果列于下表．准确值(6位有效数字) 353553.0)2( =′f
h h+2 h?2
h
hh
2
22+
h
hh
f
2
22
)2(
+
′
0.5 1.58114 1.22474 0.356400 -0.002847
0.05 1.43178 1.39642 0.353600 -0.000047
0.005 1.41598 1.41244 0.354000 -0.000447
0.0005 1.41439 1.41404 0.350000 0.003553
0.00005 1.41423 1.41420 0.300000 0.053553
由表可知当05.0=h时计算效果最佳．当05.0<h时，(2)的分子是相近的两数相减，
由于舍入误差的影响，精度降低．越小，计算效果越差,h
为了减少误差，把中点差商公式(2)改写为
60
hh
f
++
=′
22
1
)2( (3)
下表列出了用公式(3)计算的结果．可见当h越小时，精度越高,
h
h+2 h?2
hh?++ 22
1
hh
f
++
′
22
1
)2(
0.5 1.58114 1.22474 0.356394 -0.002864
0.05 1.43178 1.39642 0.353582 -0.000029
0.005 1.41598 1.41244 0.353554 -0.000001
0.0005 1.41439 1.41404 0.353553 0.000000
0.00005 1.41423 1.41420 0.353555 0.000000
61
第五章 常微分方程数值解法
5.1 内容提要
本章主要介绍一阶常微分方程初值问题
0
(,) (5.1)
( ) (5.2)
dy
fxy
dx
ya y
=
=
在区间
[ ]
,ab上的数值解法。
本章讨论的数值方法，它不是求方程的解的 解析表达式或近似表达式，而是求在一系列离散节点( 1,2,......,)
i
x = n
0
上的函数值的近似值。数值解法的优点在于用它可以直接求出具有任意精确度的数值解。
()
i
yx
i
y
1,欧拉（Euler）方法及改进的欧拉方法
称
1
(,),( 0,1,2,...) (5.3)
nn nn
yyhfxyn
+
=+ =
为欧拉公式，它是一个显式的差分方程，由于初值已知，从式（5.3）
0
y
可得
10
(,)
n
yyhfxy= +
相继地可以求得利用欧拉公式求初值问题式（5.1）~式(5.2)的方法就称为欧拉方法。
23
,,...,,....
n
yy y
称预测—校正格式
)
n
x
1n
1
111
(,) (5.4)
(,) (,),
2
nnnn
nn nn nn
yyhfxy
h
yy fxyfxy
+
+++
=+
=+ +
为改进的欧拉公式，用它求初值方程式(5.1)~式(5.2)的数值的方法称为改进的欧拉方法。
2,截断误差
设则称误差(
n
yy=
11
()
nn
yx yε
++
=?
+
n
为方法的局部截断误差。
称在节点x处的初值问题的精确解与方法的数值解的误差()
nnn
yx y=?l的 方法为整体的截断误差。
63
定义 5.1 若某种微分方程数值公式的局部截断误差为
1
1
()
p
n
ohε
+
+
=
则称这种方法具有p阶精度或称为p阶方法这里，p为非负整数。通常，p越大，h越小，
数值方法越精确。
3,龙格－库塔(Runge－Kutta)方法
龙格－库塔方法是一类应用度比较广的高精度单步法，简称为R－K方法,
显式R—K方法的一般形式为
1
1
1
1
1
(,) (5.5)
(,)
r
nn i
i
nn
i
inin ij
j
yyhk
kfxy
kfx hyh k
ω
αβ
+
=
=
=+
=
=+ +
∑
∑
( i=2,3,…,r )
64
ij
其中，参数,,
ii
ω αβ是既与方程式 (5.1)右端(,)f xy无关，又与步数n无关的常数，选择它们的原则是要求式(5.5)的右端在
( )
,nn
x y处作泰勒展开后，按h的升幂排列重新整理得式
23
12 3
11
....,(5.6)
2! 3!
yy h h hγγ γ=+ + + +
再与微分方程式(5.1)的精确解y(x)在点
n
x x=处的泰勒展开式
'2'
1
1
( ) ( ) ( ) ( ) ( ),..,( ) ( )
2! !
p
pp
nn nn n n
h
yx yx h yx hy x hy x y x h
p
ο
+
+
=+= + + ++ +
1
相比较,使其有尽可能多的项重合，例如，要求
'' '' ( 1)
12 3
,,,...,
p
nnnpn
rfrfrf r f
=== =
就得到p个方程，从而定出参数,,
iiij
ω αβ，再代入的表达式，就可得到计算微分方程问题式(5.1)、(5.2)的数值计算公式,
12
,,...,
r
kk k
1
1
(5.5)
k
nn i
i
yyhkω
+
=
=+
∑
式(5.5)就称为r段的 R－K方法的计算公式,
若式(5.6)与的 泰勒展开式的前p+1项完全一致,即使局部截断误差达到
1
()
n
yx
+
1
(
p
hο
+
)
ij
,则称公式(5.5)为p阶r段的R—K方法,
通常用得比较多的是r取到2,3,4在这些情况下,为了计算方便,可适当选取
,,
ii
ω αβ,使p=r(即阶数=段数),
例如,称二阶二段的R—K计算公式
)
12
1
221
(,) (5.7)
11
(,
22
nn
nn
n
yyhk
kfxy
kfx hy hk
+
=+
=
=+ +
为中点公式,
称三阶三层的R—K计算公式
()
)
1123
1
221
4
6
(,) (5.8)
11
(,
22
nn
nn
n
h
yy kkk
kfxy
kfx hy hk
+
=+ + +
=
=+ +
为库塔公式,
在实际应用中常用的四阶R—K公式，例如
(1)标准四阶R—K公式
11234
1
21
32
43
(22 )
6
(,)
1
(5.9)(,)
22
1
(,)
22
(,)
nn
nn
nn
nn
nn
h
yy kkkk
kfxy
h
kfx y hk
h
kfx y hk
kfxhyhk
+
=+ + + +
=
=++
=++
=++
(2)吉尔(Gill)公式
11 2 34
1
21
312
423
(2 2) (2 2)
6
(,)
1
(,)
(5.10)22
21 2 2
(,)
22 2
222
(,)
nn
nn
nn
nn
nn
h
yy k k kk
kfxy
h
kfx y hk
h
k f x y hk hk
kfxhy hk hk
+
=+ +? ++ +
=
=++?
=++ +
+
=+? +
其局部截断误差均为O( ),
5
h
4,收敛性与稳定性
收敛性与稳定性从不同角度描述了数值方法的可靠性，只有既收敛又稳定的方法才能提供比
65
较可靠的计算结果。下面只给出单步法的收敛性与稳定性,
定义 5.2 若数值方法对任意固定的节点
0
( 0,1,2,.....)
n
xxnhn= +=，当且时有，则称方法是收敛的,
0(→h n
)
n
x
)→∞
(
n
yy→
由此可见，数值方法的收敛性并不涉及计算过程的舍入误差，而只是与方法的截断误差有关。下面给出关于估计方法的截断误差及判别单步法收敛的充分定理。
定理 5.1 若(1)单步法
1
(,,)
nn nn
yyhxyhφ
+
=+中的增量函数(,,)x yhφ在区域上连续,并且关于y满足李普希茨条件,
0
,0hh≤≤,axb y≤≤?∞<<+∞
1212
(,,) (,,),( 0);xy h xy h Ly y Lφφ?≤?>
(2)方法的局部截断误差为
1
( ),( 1,2,....)
p
n
oh nε
+
==
(3)初始值是精确的，即
0
()ya y=
则：(1)方法的整体截断误差为，(2)当p时，方法是收敛的。 1≥
数值稳定性有各式各样的定义，下面只简单的介绍"绝对稳定性"的概念,
定义 5.3 设用某一种数值方法计算时，所得到的实际计算结果为且由误差引起以后各节点处的误差为
n
y
n>
%
y
n
y )
,
n
null
nn
yδ =? (
m
ym
m
δ，如果
mn
δ δ≤，则称该数值方法绝对稳定的。
数值稳定性的分析相当复杂，它不仅与方法本身有关，而且总跟方程的右端(,)f xy以及步长有关，所以，对一般的微分方程，研究这个问题很困难但事实上，研究数值方是否稳定，
不可能也不必对每个不同的右端函数进行讨论，为简单起见，通常只对试验方程(也称为模型方程)
),( yxf
y
dy
dx
λ= （5.11）
(其中λ为常数，当λ复数时，0)Re( <λ )进行讨论，即研究将数值方法用于解方程(5.11)
得到的差分方程是否数值稳定,
在把某一数值方法应用于实验方程，为了保证数值稳定，步长和h λ都要受到一定的限制，故有
定义5.4 一个数值方法用于解试验方程(5.11)，若在hλμ =平面中的某个区域Ｒ中方法都是绝对稳定的，而在域R外，方法是不稳定的，则称区域R是该数值方法的绝对稳定域,
66
67
h
显然，绝对稳定性区域越大，该方法的绝对稳定性越好，若某个方法的绝对稳定性区域包含λ复平面的左半平面，则称此数值方法是Ａ-稳定的,
对于＂标准的四阶R-K方法＂可得其绝对稳定性区域为1
!4!3!2
1
432
<++++
μμμ
μ，
若设0<λ，则可得到当
λ
78.2
<< h
y
0时，四阶R-K方法绝对稳定．且当方程(5.1)的右端函数在其定义域上关于满足李普希兹条件时，四阶的R-K方法也是收敛的,),y(xf
5,线性多步法
其一般形式为
)(
110111101 rnrnnnrnrnnn
ffffhyayayay
++
++++++++= ββββ LL (5.12)
可简写为 (5.13)
∑∑
=
=
+
+=
r
k
knk
r
k
knkn
fhyay
10
1
β
这里，f
kkjjj
nrnrnjyxf βα,),1,,1,(),,( ++== L均为常数．若0
1
≠
β，则式(5.12)
是隐式;如0
1
=
β，则式(5.12)是显式,
构造线性多步法公式常用泰勒展开方法和数值积分的方法,
例如，常用的阿达姆斯四步外推显式公式
)9375955(
24
3211+
+?+=
nnnnnn
ffff
h
yy (5.14)
其局部截断误差为
),(
720
251
13
55
+?
∈=
nn
xxyhR ηη）（
）（
所以它是一个四阶公式，也称为四阶阿达姆斯外推公式,
又如四阶的阿达姆斯隐式公式
)5199(
24
2111++
+?++=
nnnnnn
ffff
h
yy (5.15)
其局部截断误差为
),(),(
720
19
12
)5(5
+?
∈?=
nn
xxyhR ηη
由于积分区间在插值区间[ ]内，故阿达姆斯隐式公式又称为阿达姆斯内插公式．它是三步隐式公式,
12
,
+? nn
xx
由阿达姆斯计算公式(5.14)与公式(5.15)可知，阿达姆斯方法每一步仅调用一次函数
,因而计算量比R-K方法少．实际计算时，常常把阿达姆斯外推公式和内插公式联合使用，构成阿达姆斯预测－校正系统，经常使用的是
),( yxf
3=r的情形，即
+?++=
+?+=
+++
+
]519),(9[
24
)9375955(
24
21111
3211
nnnnnnn
nnnnnn
fffyxf
h
yy
ffff
h
yy
校正预测
这是一个四步四阶公式，至于初始值，除以外，其他的通常用四阶R-K公式计算,
0
y
321
,,yyy
5.2 典型例题
例5.１对常微分方程初值问题
≤≤=
=
)10(1)0( xy
y
dx
dy
取步长，分别用欧拉方法、改进的欧拉方法及标准的四阶R-K方法作数值计算，并将计算结果与方程的精确解进行比较,
1.0=h
x
exy
=)(
解 由欧拉公式(5.3)改进的欧拉公式(5.4)及标准的四阶R-K公式(5.9)，可得三种方法的具体计算公式分别如下,
(1) 欧拉方法,
nnnn
yyyy 9.0)(1.0
1
=?+=
+
(2) 改进的欧拉方法,
=?+?+=
=?+=
++
+
nnnnn
nnnn
yyyyy
yyyy
905.0)]()[(
2
1.0
9.0)(1.0
11
1
(3) 标准的四阶R-K方法,
=+?=
=+?=
=+?=
=
=++++=
+
nn
nn
nn
n
nnn
ykyk
ykyk
ykyk
yk
ykkkkyy
90475.0)1.0(
9525.0)05.0(
95.0)05.0(
9048375.0)22(
6
1.0
34
23
12
1
43211
计算结果与精确解的比较如表5-1所示,
68
表 5-1
欧拉方法
改进的
欧拉方法
标准四阶
R-K方法
精确解法
x
n
y
n
y
n
y )(
n
x
0.1 0.900000 0.905000 0.904837 0.904837
0.2 0.810000 0.819025 0.818730 0.818730
0.3 0.729000 0.741218 0.740818 0.740818
0.4 0.656100 0.670802 0.670320 0.670320
0.5 0.590490 0.607076 0.606530 0.606530
0.6 0.531441 0.549404 0.548812 0.548811
0.7 0.478297 0.497210 0.496585 0.496585
0.8 0.430467 0.449975 0.449329 0.449328
0.9 0.387421 0.407228 0.406569 0.406569
1.0 0.348679 0.368541 0.367879 0.367879
n
y
可见，标准四阶R-K方法计算的精确度最高,
例5.2 用欧拉方法解初值问题
=
+=
0)0(y
bax
dx
dy
证明：其整体截断误差为
,
2
1
)(
2
anhyxy
nnn
=?=e
这里，是欧拉方法的数值解．,ynhx
n
=
n
证 显然，该初值问题的精确解为
bxaxxy +=
2
2
1
)(，则有
nnn
bxaxxy +=
2
2
1
)(,
由欧拉公式(6.3)，有
),2,1(),(
11111
L=++=
nhfyyxhfyy
nnnnnn
记
今将已知,代入上述递推式，有 baxyxf +=),(
∑∑
=
=
++=+=
+++=
++=+=
1
0
1
0
0
1233
12211
0
)(
)(
n
i
i
n
i
i
nnnn
nnnnnn
baxhfhy
hfhfhfy
hfhfyhfyy
L
69
nnn
n
n
i
i
ahxaxbxn
n
ahnbh
nhhhahnbh
xxxahnbhxhanbh
2
1
2
1
)1(
2
)1(1
)20(
)(
2
2
110
1
0
+=
+
+=
+++++=
++++=+=
=
∑
L
L
而
nnn
bxaxxy +=
2
2
1
)(
所以
2
22
2
1
)(
2
1
2
1
)
2
1
2
1
(
2
1
)(
anhnhahahx
ahxaxbxbxax
yxye
n
nnnnn
nnn
===
+?+=
=
例5.3 用泰勒展开方法求解初值问题
=
+=
00
22
)( yxy
yx
dx
dy
解 泰勒展开的方法也是微分方程离散化的一个重要方法．这个方法的具体做法是，假设初值问题式(5.1)-式(5.2)的精确解在a)(xy bx ≤≤上存在1+p阶连续导数，那么,
)
n
(
1+n
xy在处作泰勒展开，有 xx =
1
23
()
()( )
() () () ()
2! 3!
() ()(5.17)
!
nn
nn n
p
p
np
yx yx h
hh
yx hy x y x y x
h
yx Rx
p
+
=+
′ ′′ ′′′=+ + +
++ +L
n
其中截断误差为
)()()(
)!1(
)(
1
1)1(
1
+
++
+
<<=
+
=
nn
pp
p
p
xxhOy
p
h
xR ξξ
在式(6.17)中略去截断误差，并以近似值代替，则得
)(k
n
y )(
)(
n
k
xy
)(
2
1
!!2
p
n
p
nnnn
y
p
h
y
h
yhyy ++′′+′+=
+
L (5.18)
用公式(5.18)解初值问题式(5.1)-式(5.2)的数值方法称为泰勒展开法，由于它的截断误差是
，故称为)(
1+p
hO p阶方法,
70
71
2
) y
对于本题，用泰勒展开法的关键是先求出的各阶导数，再代入泰勒展开式(5.17)，由已知，直接求其各阶导数，得
y
2
yxy +=′
M
)32)((8)53(44
62422
)3)((242
)(222
)(2222
222222
)4(
2222
2
22
yxyxyyxxy
yyyyyyyyyyy
yxyxxy
yyyy
yxyxyyxy
+++++=
′′′+′′′=′′′+′′′+′′′=
++++=
′+′′+=′′′
++=′+=′′
当然，可以继续求导下去，但随着求导次数增大，工作量越来越大，如果要求得到公式的局部截断误差为O，则把上列求得的各阶导数代入泰勒展开式(5.17)中，并设,
我们就构造了一个求解微分方程
(
4
h
nn
yx =)(
22
yx
dx
dy
+=的数值公式,
″
+
′
++=
+ nnnnn
fhfhhfyy
22
1
6
1
2
1
其中，表示对
)k(
n
f ),( yxf x的阶编导数在k
n
xx =点的值。由于这个公式的局部截断误差是，故称为三阶方法。 )(
4
hO
从理论上讲，只要方程的解充分光滑，用泰勒展开方法可以构造任意有限阶的数值计算公式，但由于求复合函数的高阶导数往往是很繁琐的，因此，泰勒展开法一般只用于求
“表头”，即开头几点的的数值解，如等。
)(xy
4321
,,,yyyy
例 5.4 取步长h，写出用标准四阶R-K方法求解初值问题 4.0=
≤≤=
+=
)91(0)1(
)sin(
xy
yxx
dx
dy
的计算公式。
解,0,4.014.01
00
=+=×+=+= ynnnhxx
n
取,202,1,0 L=n其标准四阶R-K
方法的计算公式为
=
++++=
++++=
++++=
+++=
++++=
+
)19,2,1,0(
)4.04.04.1sin()4.04.1(
)2.04.02.1sin()4.02.1(
)2.04.02.1sin()4.02.1(
)4.01sin()4.01(
)22(
6
4.0
34
23
12
1
43211
Ln
kynnk
kynnk
kynnk
ynnk
kkkkyy
n
n
n
n
nn
例 5.5 用泰勒展开的构造方法构造一个以为节点的线性三步公式,
112
,,,
+ nnnn
xxxx
)(
11011221101?++
+++++=
nnnnnnn
fffhyayayay βββ (5.19)
解 利用泰勒展开的方法构造线性多步法的计算公式的基本思想类同于R-K方法的推导，它是先将线性多步法公式(5.11)的右端均在
n
xx =处作泰勒展开，然后按的幂次作升幂排列重新整理得
h
L+?++?
+
+
+?+=
+
==
+
+
==
==
+
∑∑
∑∑∑∑
1
11
1
1
1
)(
1
1
10
1
])()1()([
)!1(
]])()([
!
)(
p
n
r
k
k
p
k
r
k
p
p
p
j
j
n
r
k
k
j
r
k
k
j
j
n
r
k
kn
ykpk
p
h
ykjk
j
h
yy
βα
βαα
(5.20)
于是，欲使公式(5.12)成为p阶的，即局部截断误差为O只要对展开式(5.20)与在处的泰勒展开式
),(
1+p
h
)
1 n
(
+n
xy x
L+
+
+=
+
+
=
+ ∑
)1(
1
0
)(
1
)!1(!
)(
p
n
pp
j
j
n
j
n
y
p
h
y
j
h
xy （5.21）
相比较，能符合到项，为此令
p
,
h
（5.22）
==?+?
=
∑∑
∑
=
=
=
),,2,1(1)()(
1
1
1
1
0
pjkjk
r
k
k
j
r
k
k
j
r
k
k
Lβα
α
解出系数,
kk
βα代入式(5.13)，即得p阶的线性多步法公式，式（5.21）与式（5.20）相减，
可得方法的局部截断误差为
L+?+
+
=?
+
==
+
+
++ ∑∑
)1(
11
1
1
11
])()1()(1[
)!1(
)(
p
n
r
k
k
p
r
k
k
p
p
nn
ykpk
p
h
yxy βα （5.23）
72
对本题，显然，关键只要确定公式中的系数
101210
,,,,,βββααα
，设
)1,,1()(),(
),1,2()(
+?=′=
==
nnnixyyxf
nnnixyy
iii
ii
然后将函数)(),(),(),(
1112 +
′′
nnnn
xyxyxyxy在点
n
xx =处作泰勒展开后，代入式（5.19），
再按h的升幂排列整理得
L+
+++
+++
′′′+++′′?+++
′++++++=
+
5)5(
1121
4)4(
1121
3
1121
2
1121
101212101
)()
24
1
24
1
120
32
120
1
(
)()
6
1
6
1
3
2
24
1
(
)()
2
1
2
1
3
4
6
1
()()2
2
1
(
)()2()()(
hxy
hxy
hxyhxy
hxyxyy
n
n
nn
nnn
ββαα
ββαα
ββααββαα
βββααααα
(5.24)
另一方面，将在处作泰勒展开，有 )()(
1
hxyxy
nn
+=
+ n
x
L++
+′′′+′′+′+=
+
5)5(
4)4(32
1
)(
!5
1
)(
!4
1
)(
!3
1
)(
!2
1
)()()(
hxy
hxyhxyhxyhxyxyxy
n
nnnnnn
（5.25）
显然，要使公式成为五阶，也即局部截断误差为，则要使公式（5.24）与式（5.25）中的前六项关于的同次幂系数完全相同，并由的任意性可得六个方程,
)()(
6
11
hOyxy
nn
=?
++
)(xyh
=++
=?++
=++
=?++
=+++
=++
L
120
1
24
1
24
1
120
32
120
1
:
24
1
6
1
6
1
3
2
24
1
:
6
1
2
1
2
1
3
4
6
1
:
2
1
2
2
1
:
12:
1:
1121
5
1121
4
1121
3
1121
2
10121
1
210
0
ββαα
ββαα
ββαα
ββαα
βββαα
ααα
h
h
h
h
h
h
(5.26)
73
从式(5.26)的前六个方程中解出6个未知元,,,,,,
101210
βββααα
再代入式(5.19)中，即得到一个五阶的数值计算公式。当然，我们也可以只要求前面几个方程成立，譬如说，我们只要求前面四个方程成立，则从前面四个方程中解6个未知元，这样，就有两个自由参数，
如令,0,0
0
==
2
αα这时，从前面四个方程中可解得,
3
1
,
3
4
1
=β,
3
1
,1
011
===
ββα代入式
（5.19），得到计算公式,
)4(
3
1111?+?+
+++=
nnnnn
fff
h
yy （5.27）
称为辛甫生公式，此公式也可用数值积分得到。
从上面的推导过程知此公式的局部截断误差是O，但这时，式（5.26）中第五个方程恰巧也满足，所以阶数可再提高一阶，是一个4阶公式，但此公式的稳定性较差，哈明用
)(
4
h
的不同数值进行试验发现，若令0
1
=α，并要求式（5.26）中的前5个方程成立，则可得
8
3
=,
8
6
,
8
1
,
8
9
1120
=?==
ββαα,
8
3
0
= β
1
α
由此得稳定性较好的哈明(Hamming)公式,
)]2(39[
8
1
1121?+?+
++?=
nnnnnn
fffhyyy （5.28）
此式是隐式三步法，这个公式不能用数值积分方法推得，可证得其局部截断误差为
)(
40
1
)(
)5(5
11
ξyhyxyR
nn
=?=
++
例5.6 对初值问题
=
=
1)0(y
y
dx
dy
证明：用梯形法求得的数值解为
n
n
h
h
y
+
=
2
2
，并证明当步长时，y收敛于初值问题的精确解。
0→h
n
x
exy
=)(
证 梯形法的计算公式为
)],(),([
2
111 +++
++=
nnnnnn
yxfyxf
h
yy
本题中，,),( yyxf?=代入梯形公式中，得计算式,
),(
2
11 ++
+=
nnnn
yy
h
yy
则有 )(
2
11 nnnn
yy
h
yy+=
74
经整理得
1
2
2
+
=
nn
y
h
h
y
从而有
02
2
1
2
2
)
2
2
(
2
2
y
h
h
y
h
h
y
h
h
y
n
nnn
+
==
+
=
+
=
L
由已知，故证明了用梯形法求得的数值解为 1
0
=y
n
n
h
h
y
+
=
2
2
以下证明当时，有；此处，0→h
x
n
ey
→ nhx =，计算极限
x
x
x
h
h
h
h
x
h
h
x
h
n
h
e
e
h
h
h
h
h
h
h
h
h
h
+
→
→→→
==
+
+
=
+
=
+
=
+
1
]
2
2
1
2
2
1[
1
lim
2
2
1
1
lim
2
2
1lim
2
2
lim
2
1
2
2
0
000
这说明，对于本题的初值问题，梯形法是收敛的。
例5.7 讨论梯形法的绝对稳定性区域。
解 由定义5.4知，把梯形法的计算公式用于试验方程y
dx
dy
λ=（其中λ为常数，当λ
是复数时，0)( <λ
e
R）得其计算公式为
()
][
2
],),([
2
1
111
+
+++
+=
++=
nnn
nnnnnn
yy
h
y
yxfyxf
h
yy
λλ
若的实际计算值为
n
y
n
y
~
，则由误差
nnn
yy
~
=δ引起的误差为
1+n
y
nn
h
h
δ
λ
λ
δ
+
=
+
2
2
1
。
由此可见，要使的误差
n
y
n
δ逐步削弱，hλμ =应满足,1
2
2
<
+
μ
μ
显然，对于任意的步长都满足，也即对步长不需作任何限制，梯形法对0>h h +∞<< h0都是绝对稳定的，
并称它为A-稳定。
例5.8 证明：对初值问题
75
≤≤=
=
)()(
),(
0
bxayay
yxf
dx
dy
当方程的右端函数在其定义区域上关于满足李普希兹条件时，标准R-K方法是收敛的。
),( yxf y
证 由标准R-K方法的计算公式（5.8）和定理5.1知,
当假设在区域),( yxf ∞+<<?∞≤≤ ybxa,上关于满足李普希兹条件，即 y
2121
),(),( yyLyxfyxf?≤?
很容易验证标准四阶R-K方法的增量函数),,( hyxφ关于也满足李普希兹条件，从而标准四阶R-K方法收敛。
y
显然，R-K方法的增量函数
)],,(),,(2),,(2),,([
6
1
),,(
4321
hyxkhyxkhyxkhyxkhyx +++=φ
其中
)),,(,(),,(
)
2
),,(
,
2
(),,(
)
2
),,(
,
2
(),,(
),,(),,(
34
2
3
1
2
1
hyxhkyhxfhyxk
hyxhk
y
h
xfhyxk
hyxhk
y
h
xfhyxk
yxfhyxk
++=
++=
++=
=
由于关于满足李普希兹条件，从而 ),( yxf y
*
**
*
1*1
*
1*1
*
22
**
11
2
1
1
)
2
1
(
)
2
),,(
(
2
),,(
2
),,(
,
22
),,(
,
2
),,(),,(
),,(),,(
yyhLL
yyhLyyL
hyxhk
y
hyxhk
yL
hyxhk
y
h
xf
hyxhk
y
h
xf
hyxkhyxk
yyLhyxkhyxk
+=
+?≤
+?+≤
++
++=
≤?
同理，可得以下不等式,
76
*32
*
33
*2
*
33
])(
4
1
)(
2
1
1[
),,(),,(
,])(
4
1
2
1
1[
),,(),,(
yyhLhLhLL
yyxkhyxk
yyhLhLL
hyxkhyxk
+++≤
+
+≤
于是，作为线性组合的)4,3,2,1( =ik
i
),,( hyxφ在
0
0,,hhybxa ≤≤+∞<<?∞≤≤上有
**3
0
2
00
*
])(
24
1
)(
6
1
2
1
1[
),,(),,(
yyLyyLhLhLhL
hyxhyx
=+++≤
φφ
即),,( hyxφ关于满足李普希兹条件，其李普希兹常数为 y
])(
24
1
)(
6
1
2
1
1[
3
0
2
00
LhLhLhLL +++=
所以，标准四阶R-K方法收敛。
例5.9 分别用四阶阿达姆斯显式和隐式公式求初值问题
≤≤=
++?=
)10(1)0(
1
xy
xy
dx
dy
的数值解，取步长1.0=h。
解 本题的精确解为，根据题意，有 xexy
x
+=
)(
1,1.0 ++?===
nnnn
xyfnnhx
由四阶阿达姆斯显式公式（5.13），有
)9,,4,3(
]52.224.09.07.39.55.18[
24
1
)9375955(
24
321
3211
L=
+++?+=
+?+=
+
n
nyyyy
ffff
h
yy
nnnn
nnnnnn
由四阶阿达姆斯隐式公式（5.14），有
]52.224.01.05.01.229.0[
24
1
)5199(
24
211
2111
++?++?=
+?++=
+
++
nyyyy
ffff
h
yy
nnnn
nnnnnn
77
由上式可解出
)9,,3,2(
]52.224.01.05.01.22[
9.24
1
211
L=
++?+=
+
n
nyyyy
nnnn
利用精确解求出起步值（一般可用四阶R-K公式计算起步值）后，按上面的公式计算，其结果如表5-2所示。
,)( xexy
x
+=
表 5-2
阿达姆斯显式法 阿达姆斯隐式法
n
x
n
y
nn
yxy?)(
n
y
0.3 初值 1.04081800
7
10146.
×2
0.4 1.07032292
6
102.874
×
1.07031966
7
10846.
×3
0.5 1.10653547
6
10816.4
×
1.10653013
7
10213.
×5
0.6 1.14881840
6
106.772
×
1.14881100
7
10285.
×6
0.7 1.19659339
6
10090.8
×
1.19658459
7
10106.
×7
0.8 1.24933815
6
10192.9
×
1.24932819
7
10714.
×7
0.9 1.30657961
6
10954.9
×
1.30656884
7
10141.
×8
1.0 1.36788995
6
10052.1
×
1.36787859
7
10418.
×8
nn
yxy?)(
从表5-2中可以看出，四阶阿达姆斯隐式法比显式法的精度高，一般地，同阶的隐式法比显式法精确，而且数值稳定性也好。但在隐式公式中，通常很难解出，需要用迭代法求解，
这样，又增加了计算量。因此，实际计算时，很少单独使用显式公式或隐式公式，而是将它们联合使用：先用显式公式求出的预测值，记作
1+n
y
)(
1+n
xy
1+n
y，再用隐式公式对预测值进行校正，求出的近似值，例如)(
1+n
xy 3=r时的公式（5.15）。
78
第六章 逐次逼近法
6.1 内容提要
1．根的概念
给定方程() 0=xf。如果有使得
*
x ( ) 0
*
=xf，则称为
*
x ( ) 0=xf的根，或的零点。设有正整数m使得
() 0=xf
0
() )()(
*
xgxxxf
m
=
且，则当时，称为)(
*
≠xg 2≥m
*
x ( ) 0=xf的m重根；当1=m时，称为的单重根。
*
x f
) ]
() 0=x
2．求根步骤
（1）确定所给方程存在多少根，找出每个根所在区间。所找区间要越小越好。
（2）求出含根区间中的根。
3．二分法
设在上连续，(xf [ba,( ) ( ) 0<bfaf，则( ) 0=xf在( )ba,内至少有一根。记
。一般地，设bba
0
,a ==
0
( ) ( ) 0<
nn
bfaf，则( )
nn
ax
*
∈ b,。取[的中点]
nn
ba,
(
n
x += )
n
b
n
a
2
1
，则有
2
nn
n
ab
x
*
x
≤?。如果ε≤
2
n
a
n
b
（ε为给定精度），则取
，停止计算。否则计算
n
xx ≈
*
()0=
n
xf x
*
，则，停止计算。否则和两式中有且仅有一式成立。若
n
x= ()()
nn
xf 0<af
()()
nn
bfxf 0< ( ) ( ) 0<
n
x
n
faf
nnnn
bbx
，令a；
若，令a
n
x=
nn
ba =
+1
,
n+1
()
n
fxf ()0<
n
b ==
++ 11
,。我们有( )( )
1
<
+n
bf 0
1+n
af且
(
n
b?
2
1
)
0 ]
n
a
nn
ab =?
++ 11
。重复以上过程。
4．迭代法
（1）迭代格式的构造及其敛散性条件
设求在区间[内的根。将其改写成等价形式 ()=xf ba,
( )xx?= （6.1）
构造迭代格式
79
(),...2,1,0,
1
==
+
kxx
kk
（6.2）
给定，可得序列
0
x { }
∞
=0kk
x。称(x)?为迭代函数，称{ }
∞
=0kk
x为迭代序列。
定理6.1 设( )x?在上具有一阶连续导数，且满足如下2个条件：①，[ba,] x
]
[]ba,∈
() [bax,∈?；②存在正常数L，使得对任意[ ]bax,∈有( ) 1<≤′ Lx?。则
1）方程（6.1）在上有唯一根； [ba,]
]
*
x
2）对任意[bax,
0
∈，迭代格式（6.2）收敛，且；
*
lim xx
k
k
=
∞→
3）,...;2,1,
1
1
*
=?
≤?
kxx
L
L
xx
kkk
4）,...;2,1,
1
01
*
=?
≤? kxx
L
L
xx
k
k
5）().
*
*
*
1
x
xx
xx
k
k
k
′=
+
∞→
lim
定理6.2 设方程（6.1）在区间[ ]ba,内有根，且当
*
x [ ]bax,∈时() 1≥′ x?，则对任意初值且，迭代格式（6.2）发散。 []bax,
0
∈
*
0
xx ≠
（2）迭代法的局部收敛性
对于方程（6.1），若在的某个邻域
*
x { }δ≤?=
*
xxxS内，对任意初值，迭代格式（6.2）收敛，则称迭代格式（6.2）是局部收敛的。
S∈
0
x
定理6.3 设方程（6.1）有根，且在的某个邻域
*
x
*
x { }δ≤?=
*
xxxS内()x?存在一阶连续的导数，则
1）当()1
*
<′ x?时，迭代格式（6.2）局部收敛；
2）当()1
0
*
>′ x?时，迭代格式（6.2）发散；
（3）迭代法的收敛速度
设序列{}收敛于，记e，如果存在非零常数c和正数
p，使得
∞
=kk
x
*
x,...2,1,0,
*
=?= kxx
kk
c
e
e
p
k
k
k
=
+
∞→
1
lim
80
81
0
则称序列{}是p阶收敛的。
∞
=kk
x
当p=1且10 << c时，称为线性收敛，当时，称为超线性收敛。特别，当p=2
时，称为平方收敛。p越大，收敛越快。
1>p
定理6.4 若(x)?在附近的某个邻域内有阶连续导数，且
*
x )1( ≥pp
( ) ( )
()
( )
( )
( ) 0,0,,0,
**1***
≠==′=
xxxxx
pp
L则对任意靠近于的初值，迭代格式（6.2）是p阶收敛的，且有
*
x x
0
()
( )
( )
!
lim
*
*
*
1
p
x
xx
xx
p
p
k
k
k
=
+
∞→
如果p=1，要求( ) 1
*
<′ x?
（4）埃特金加速法
设已有迭代格式（6.2），且是收敛的，则下列算法称为埃特金加速法,
() L,2,1,0,
1
==
+
kxx
kk
φ （2.3）
其中
()
( )( ) ( )
()()()xxx
xxx
x
+?
=
φ
2
2
定理6.5 设在附近，
*
x (x)?有p+1阶导数，则对一个充分靠近的初值，有
*
x x
k
0
1）如果（6.2）是线性收敛，则（6.3）是2阶收敛；
2）如果（6.2）是阶收敛，则（6.3）是2p-1阶收敛。 )2( ≥pp
5．牛顿法
（1）牛顿迭代公式局部收敛性
设为x () 0=xf的一个近似值，在点( )( )
kk
xfx,作( )xf的切线
() ( )( )
kkk
xxxfxfy?′+=
该切线与x轴交点的横坐标作为新的根近似值，即
1+k
x
( )
()
k
k
kk
xf
xf
xx
′
=
+1
（6.4）
（6.4）称为牛顿公式，它是局部收敛的。
当为单根时有
*
x
( )
()
*
*
2*
*
1
2)(
lim
xf
xf
xx
xx
k
k
k
′
′′
=
+
∞→
当是重根时有
*
)x (2≥mm
mxx
xx
k
k
k
1
1lim
*
*
1
=
+
∞→
（2）求重根的修正牛顿公式
①重数m已知：
( )
()
k
k
kk
xf
xf
mxx
′
=
+1
,
②重数m未知：
( )
()
( )
()xf
xf
xu
xu
xu
xx
k
k
kk
′
=
′
=
+
)(,
1
。
6．割线法
设是方程的两个近似值。以,x
1?k
x
k
() 0=xf ( )( )
11
,
kk
xfx和()( )
kk
xfx,两点作的割线 (f )x
()
() ( )
()
k
kk
kk
k
xx
xx
xfxf
xfy?
+=
1
1
该割线与x轴交点的横坐标作为新的近似根，即
()
(
1
1
1
)()(
+
+=
kk
kk
k
kk
xx
xfxf
xf
xx ) （6.5）
（6.5）称为割线公式。它也是局部收敛的。
7,向量范数和矩阵范数
(1) 常用的三种向量范数
向量的1-范数：
∑
=
=
n
i
i
xx
1
1
向量的2-范数：
∑
=
=
n
i
i
xx
1
2
2
向量的∞-范数：
i
ni
xx
≤≤
∞
=
1
max
(2) 常用的三种矩阵范数
矩阵的1-范数：
∑
=
≤≤
=
n
i
ij
nj
aA
1
1
1
max
82
矩阵的2-范数：)(
2
AAA
T
ρ=
矩阵的∞-范数：
∑
=
≤≤
∞
=
n
j
ij
ni
aA
1
1
max
其中)(Bρ为矩阵B的谱半径，即}0|max{)( =?= BEB λλρ
矩阵范数和矩阵的谱半径之间满足AA ≤)(ρ,
8,迭代法
(1) 迭代格式及其收敛的充要条件
将线性代数方程组
bAx = (6.1)
改写为
fBxx +=
构造迭代格式
(6.2) L,2,1,0,
)()1(
=+=
+
kfBxx
kk
给定，由上式可得向量序列{,如果该向量序列收敛，则称迭代格式(6.2)收敛，
称矩阵B为迭代矩阵。
)0(
x
∞
=0
)(
}
k
k
x
定理6.6 迭代格式(6.2)对任意的初值均收敛的充分必要条件为
)0(
x,1)( <Bρ
(2) 雅可比迭代格式
+=
+=
+=
+
+
+
)(
1
)(
1
)(
1
)(
11,
)(
22
)(
11
)1(
2
)(
2
)(
323
)(
121
22
)1(
2
1
)(
1
)(
313
)(
212
11
)1(
1
n
k
nnn
k
n
k
n
nn
k
n
k
nn
kkk
k
nn
kkk
bxaxaxa
a
x
bxaxaxa
a
x
bxaxaxa
a
x
L
LLLL
L
L
迭代矩阵J的特征方程为
0=++ UDL λ
其中
83
0
0
0
0
0
1,321
3,12,11,1
3231
21
=
nnnnn
nnn
aaaa
aaa
aa
a
L
L
L
OMMM
=
=
0
0
0
0
0
,
,1
31,3
21,223
11,11312
33
22
11
nn
nn
nn
nn
nn
a
aa
aaa
aaaa
U
a
a
a
a
D
MMO
L
L
L
O
如果A是严格对角占优矩阵，则雅可比迭代法收敛。
(3) 高斯－赛德尔迭代格式
+=
+=
+=
+
+++
++
+
)(
1
)(
1
)(
1
)1(
11,
)1(
22
)1(
11
)1(
2
)(
2
)(
323
)1(
121
22
)1(
2
1
)(
1
)(
313
)(
212
11
)1(
1
n
k
nnn
k
n
k
n
nn
k
n
k
nn
kkk
k
nn
kkk
bxaxaxa
a
x
bxaxaxa
a
x
bxaxaxa
a
x
L
LLLL
L
L
迭代矩阵G的特征方程为
0=++ UDL λλ
如果A是严格对角占优矩阵或者A是对称正定矩阵，则高斯－赛德尔迭代法收敛。
6.2 典型例题
例6.1判断下列方程有几个实根，并求出其隔根区间,
1） 035
3
= xx
2）
x
ex
= 2
解：1）设，则() 35
3
= xxxf ()
=?=′
3
5
353
22
xxxf。
84
当
3
5
<x时，( ) xfxf,0< ()′为减函数；当
3
5
>x时，( )(xfxf,0> )′为增函数。
)2(,03 =?<? f,03
3
5
>? f
() 93 >=f。根据以上信息，可作( )xf的图象。可见( ) 0=xf有三个根
3
∈
∈,
3
5
,0,
3
5
),
3
5
*
3
*
2
xx∈,2(
*
1
x。
,01)0(,03
3
5
3
10
3
5
3
10
3
5
<?=<=
=
ff
0
2）将原方程改写为
x
ex
=?2
记，做函数() ()
x
exfxxf
=?=
21
,2 ( )xf
1
和( )xf
2
的图象。由图可知和有两个交点，其横坐标
()x
1
f f ()x
2
( )2,1),1
*
2
∈x,2(∈
*
1
x。因而所给方程有两个根
( )2,1),1,2(
*
2
∈ x
*
1
∈x。
例6.2 设具有m阶连续导数，证明是(xf )
*
x ( )xf的m重零点的充分必要条件为
( ) ( )
()
( )
( )
( ) 0,0,0
*1**
≠===
fxfxxf
mm
L,0 ′f
*
x
证 必要性。设是的重零点，则
1
m
x )(xf m
)()()( xgxxxf
m?
=，且 0)( ≠
xg
∑
∑
=
=
+=
=
k
i
ikimi
k
k
i
ikimi
k
k
xgxximmmc
xgxxcxf
0
)(
0
)()()(
)())(1()1(
)(])[()(
L
当时 0?≤≤ mk
0)())(1()1()(
0
)()(
=?+=
∑
=
k
i
ikimi
k
k
xgxximmmcxf L
当时 k =
0)(!)())(1()1()(
0
)()(
≠=?+=
=
∑
xgmxgxximmmcxf
m
i
imimi
m
k
L
充分性，设使得 x
85
86
0)( =
xf，，L，,0)( =′
xf 0)(
)1(
=
xf
m
0)(
)(
≠
xf
m
由泰勒公式得
m
m
m
m
m
m
xx
m
xxxf
xx
m
xxxf
xx
m
xf
xxxfxfxf
)(
!
))((
)(
!
))((
)(
)!1(
)(
))(()()(
)(
)(
1
)1(
+
=
+
+
++?′+=
θ
θ
L
其中0 1<<θ,令
))((
!
1
)(
)(
+= xxxf
m
xg
m
θ
则有
)()()( xgxxxf
m?
=且0)(
!
1
)(
)(
≠=
xf
m
xg
m
根据定义为的重零点,
)
0
x (xf m
例6.3设在[上有根.根的第次近似值为，且)( =xf ],ba
x k
k
x ],[ bax
k
∈，证明
m
xf
xx
k
)(
≤?
其中)(min
],[
xfm
bax
′=
∈
,
证 由泰勒展开式得
0))(()()( =?′+=
kkk
xxfxfxf ξ
其中介于和之间.由上式解得
k
x x
k
ξ
)(
)(
k
k
k
f
xf
xx
ξ′
=?
两边取绝对值得到
m
xf
f
xf
xx
k
k
k
)(
)(
)(
≤
′
=?
ξ
例6.4 证明方程
0210 =?+ xe
x
存在惟一实根.用二分法求出此根，要求误差不超过)1,0(∈
x
2
10
2
1
×,
解 记，则对任意210)(?+= xexf
x
Rx∈，,010)( >+=′
x
exf
因而是严格单调的，有惟一实根,)(xf 0)( =xf )1,0(∈
x
用二分法求解，要使
2
10
2
1
×≤? xx
k
，只要
2
1
10
2
1
2
01
+
×≤
k
解的64.6
2lg
2
=≥k,取.所以只要二等分7次，即可求得满足精度要求的根，计算过程列表如下表,
7=k
k
b f
k
a（的符号）)(
k
af
k
x（的符号）)(
k
xf
k
（的符号）)(
k
b
0
1
2
3
4
5
6
7
)(0?
)(0?
)(0?
)(0?
)(0625.0?
)(0625.0?
)(078125.0?
)(0859375.0?
)(5.0 +
)(25.0 +
)(125.0 +
)(0625.0?
)(09375.0 +
)(078125.0?
)(0859375.0?
)+(1
0
0
0
0
0
0
0
)(5,+
)(25,+
)(125,+
)(125,+
)(09375,+
)(09375,+
)(09375,+ (1 +)
所以，09.0)09375.00859375.0(
2
1
≈+≈
x,
例6.4 用迭代法求方程
04 =? xe
x
的根，精确至三位有效数
解 记，，xexf
x
4)(?=
x
exf =)(
1
xxf 4)(
2
=
作，的图象，见图。由图可知两曲线有两个交点，其横坐标分别为
，,
) )
0
1
(
1
xf
)1,0(∈
(
2
xf
3,2(
2
∈
1
x )x
因而有两个根和,)( =xf
1
x
2
x
1） 求.将在区间[0，1]内改写成等价形式,x 0)( =xf
x
ex
4
1
=
87
记
x
ex
4
1
)( =?，则
x
ex
4
1
)( =′?，当)1,0(∈x时]1,0[]
4
1
,
4
1
[)]1(),0([)( ∈=∈ ex，
1
4
)1()( <=′≤′
e
x，因而迭代法格式
L,2,1,0,
4
1
1
==
+
kex
k
x
k
对任意均收敛。取，计算得]1,0[
0
∈x 5.0
0
=x 4122.0
1
=x，3775.0
2
=x，x，
，
3675.0
3
=
3600.0=
4
x 3583.0
5
=x，3577.0
6
=x，3575.0
7
=x,所以,358.0
1
=
x
2） 求.将在区间[2，3]内改写成等价形式,
2
x 0)( =xf
)4ln( xx =
记)4ln()( xx =?，则
x
x
1
)( =′?,当]3,2[∈x时]3,2[]12ln,8[ln)]3(),2([)( cx =∈,
1因而迭代法格式
对任意]均收敛，取，计算得
2
1
)( <≤′ x?
L,2,1,0),4ln(
1
==
+
kxx
kk
3,2[∈x 5.2
0
=x 303.2
1
=x，221.2
2
=x，184.2
3
=，x
。所以例6.5 证明对任何初始值，由迭代公式
167.2=，
4
x 156.2
5
=x，.2
7
=x 155 16.2
2
=
x,
Rx ∈
0
L,2,1,0,cos
1
==
+
kxx
kk
所产生的序列
0
都收敛于方程
∞
=
}{
kk
x xx cos=的根。
证 x记x cos)( =?，则xx sin)(?=′?。
1） 先考虑区间]。当时，1,1[? ]1,1[?∈x ]1,1[cos)(?∈= xx?，11sin)1()( <=′≤′ x，
=1
都收敛于方程故对任意初值由迭代公式]1,1?，[∈
0
x（1）产生的序列}{
k
x
∞
k
xx cos=,
2） 对任意初值Rx ∈
0
，有]1,1[cos
01
∈= xx
∞
=1
}{
kk
x都收敛于方程
，将此看成新的迭代初值，则由1）可知，
由迭代公式（）产生的序列
1
x
1 xx cos=的根。
例6.6 已知方程
08.0
23
= xx
88
89
5附近有一个根。将此方程写成如下2个等价形式,在.1
0
=x
32
8.0 xx +=，8.0
3
= xx
构造如下两个迭代格式,
1）L,2,1,0,8.0
3 2
=+= kxx
2） （2）
（1）
1+k
L,2,1,0,8.0
3
1
=?=
+
kxx
k
值判断这两个迭代格式是否收敛，选一种收敛较快的迭代格式，求出具有4位有效数字的近似
,
1解
2x
）记
3 2
8.0)( xx +=?则
3/23
)8.0(3
)(
x
x
+
=′?
14755.0
05.3
1
)5.18.0(
1
)5.18.0(3
5.12
)5.1(
3/23/223/22
<==
+
=
+
×
=′?
所以迭代格式（1）是局部收敛的
2）8.0)(
3
= xx?，则)( =′ x?
8.02
3
3
2
x
x
,
1103.2
8.05.12
5.13
)5.1(
3
2
>=
×
=′?
所以迭代格式（2）是发散的,
3）5，利用迭代格式（1）计算，得 取.1=x 4502.1
1
=x，4265.1
2
=x，4153.1
3
=，x
41.1= 00
4
x，4075.1
5
=x，4063.1
6
=x，4057.1
7
=x，4054.1
8
=x,
,
0
所以405.1≈
x
例6.7 试用简单迭代法和埃特金加速法求方程
5附近的根，精确至四位有效数,
解
x
.简单迭代法公式
x
x
= e
在.0=x
记ex
=)(?
1
=
+
x
k
L,2,1,0),( =kx
k
计算得
k 0 1 2 3 4 5 6
90
k
x 0.5 0.60653 0.54524 0.57970 0.56006 0.57117 0.56486
k 7 8 9 10 11 12 13
k
x 0.55844 0.56641 0.56756 0.56691 0.56728 0.56707 0.56719
k 14 15 16 17 18
k
x 0.56712 0.56716 0.56714 0.56715 0.56714
所以5671.0≈
x,
埃特金算法公式
kkk
kkk
kkk
kkk
k
xxx
xxx
xxx
xxx
x
+?
=
+?
=
+
2)(2))((
)(((
22
1
,
))
其中( )
~
kk
xx?=，)
~
(
kk
xx?=。计算得
k
k
x
k
x
~
k
1 0.56762 0.56687 0.56730
2 0.56714 0.56714
x
0 0.5 0.60653 0.54524
所以5671.0≈
x,
例6.8 )设(x?在,[ ba ]上连续可微，且1)(0 <′< x?，)(xx?=在],ba[上有根，但x
则
≠ x
0
，
由
L,2,1,0),(
1
==
+
kxx
kk
产生的迭代序列}单调收敛于
证 则由
1
{
k
x
x
设bxx ≤<
0
))(()()(
0001
′=?=? xxxxxx ξ
及)(0 <′< x?知
于是
<?< xxxx
01
0
01
xxx <<
91
x
k
由归纳法可得，若，则.因而{单调下降并以为下界。
因而存在。在（1）的两边取极限得
bxx
k
≤<
kk
xxx <<
+
1
∞
=0
}
kk
x
k
x
∞→
lim
)lim(lim
k
k
k
k
xx
∞→∞→
=?
因而为方程
k
k
x
∞→
lim
)(xx?= （2）
的根。由于1)(0 <′< x?知方程（2）在[内的根是惟一的，因而 ],ba
∞→
= xx
k
k
lim
类似地，如果，可证{是单调上升收敛于,
<≤ xxa
0
∞
=0
}
kk
x
x
例6.9 设)(x?满足
1）],[,)(
000
baxxx ∈>?；
2）],[,0)(
0
bxxx ∈≥′?；
3）)(xx?=在[上有根,]
0
,
0
bx
则由出发，由 x
L,2,1,0),(
1
==
+
kxx
kk
（1）
产生迭代序列单调上升收敛于)(xx?=在[上的最小根,],
0
bx
证 记)(xx?=在[上的最小正根为.由 ],
0
bx
x
0))(()()(
0001
≤?′=?=?
xxxxxx η
及条件1）得
≤≤ xxx
10
一般地设，则由
≤≤ xxx
kk 1
0))(()()(
111
≥?′=?=?
+ kkkkkkk
xxxxxx η
有
1+
≤
kk
xx （2）
由
0))(()()(
1
≤?′=?=?
+
xxxxxx
kkkk
有
+
≤ xx
k 1
（3）
综合（2）和（3）知
+
≤≤ xxx
kk 1
92
} x由归纳原理{为单调上升序列且以为上界.故存在，记
k
x
x
k
k
x
∞→
lim
k
k
x
∞→
= lim
~
，则
x≤< xx
~
0
，在（1）的两边取极限得)
~
(
~
xx?=，即x
~
为)(xx?=的根，又因为是
x
)(xx?=在上的最小正根，所以],[
0
x b
x=x
~
。综上{单调上升收敛于}
k
x )(xx?=在上的最小正根,],[
0
bx
例6.10 以知在[内仅有一个根，而当)(xFx = ],ba ],[ bax∈时，1)( >≥′ LxF，为常数。试问如何将化为适合迭代的形式。将
L
)x(Fx = tgxx =化为适合迭代的形式，并求方程最小正根，即求（弧度）附近的根，计算结果精确至6位有效数,5.4=x
解 注意到
)(
1
])([
1
xF
xF
′
=′
则当时,],[ bax∈
1
1
])([
1
<≤′
L
xF
所以将方程写成等价形式构造迭代格式 )(xFx = )(
1
xFx
=
L,2,1,0),(
1
1
==
+
kxFx
kk
可望收敛。上面几式中表示的反函数 )(
1
xF
)(xF
对于方程
tgxx =
在4.5（弧度）附近写成等价形式,
arctgxx +=π
构造如下迭代格式,
L,2,1,0,
1
=+=
+
karctgxx
kk
π （1）
此时迭代函数为arctgxx +=π? )(。求导得
2
1
1
(
k+
=′?。
1
5.41
1
)5.4(
2
<
+
=′?。因而迭代格式（1）是局部收敛的。取5.4=x，计算结果如下,
k 0 1 2 3 4
4.5 4.49372 4.49342 4.49641 4.49341
k
x
所以tgxx =的最小正根为4.49341,
例6.11 证明迭代格式
L,2,1,0,
3
)3(
2
2
1
=
+
+
=
+
k
ax
axx
x
k
kk
k
是计算ax =
的三阶方法，并求
3
1
)(
lim
k
k
x
xa
xa
+
∞→
假设初值充分靠近,
0
x
x
解 记
ax
axx
x
+
+
=
2
2
3
)3(
)(?
则
a
aa
aaa
a =
+
+
=
3
)3(
)(? （1）
axxxaa 3)()3(
32
+=+? （2）
对（2）两边求一阶导数，得
axxxxax 33)(6)()3(
22
+=+′+
令ax =，并利用（1）得
aaaaaaa 336)()3( +=+′+?
因而
0)( =′ a? （3）
93
对（2）两边求二阶导数，得
xxxaxxax 6)(6)()3(2)()3(
22
=+′′++′′+
令ax =，并利用（1）和（3）得
0)( =′′ a? （4）
对（2）两边求三阶导数，得
6)()3(3)()3(3)()3(
222
=′′′++′′′++′′′+ xaxxaxxax
令ax =，并利用（3）和（4）得
0
2
3
)( ≠=′′′
a
a? （5）
由（1），（3）~（5）知迭代格式)(
1 kk
xx?=
+
三阶收敛，且有
a
a
ax
ax
xa
xa
k
k
k
k
k
k
4
1
)(
6
1
)(
lim
)(
lim
3
1
3
1
=′′′=
=
+
∞→
+
∞→
例6.12 Leonardo于1225年研究了方程
020102
23
=?++ xxx
并得出了是该方程的一个根。无人知道他是用什么方法得出的，在当时这是一个非常著名的结果。试用牛顿方法求此结果,
368808107.1=x
解 记则20102)(
23
++= xxxxf
3
26
)
3
2
(31043)(
22
++=++=′ xxxxf。当时，又Rx∈ 0)( >′ xf 07)1( <?=f，016)2( >=f 0)( =xf，所以有惟一实根
。用牛顿迭代格式 )2,1(∈
x
=
=
′
=
+
5.1
,2,1,0,
)(
)(
0
1
x
k
xf
xf
xx
k
k
kk
L
k
f
k
x )(
k
xf )(
k
x′
0 1.5 2.875 75.22
1 31.37362637 0.101788669 15505372.21
2 91.36881481
4
104158.1
×
09622130.21
94
3 71.36880810
8
106.1 ×?
09613922.21
4 71.36880810
并改写，10)43()( ++=′ xxxf，求得所以
,368808107.1≈
x
20)10)2(()(?++= xxxxf
例6.13 证明：对于 。，牛顿法仅为线性收敛的多重根
*
0)( xxf =
证 设重根（m>2），则 mxfx的为0)(
*
=
)()()(
*
xgxxxf
m
=
且牛顿迭代法的迭代函数为,0)(
*
≠xg
)()(()(
)()(
)()()()(
)()(
)(
)(
)(
'*
*
'*1*
*
'
xgxxxmg
xgxx
x
xgxxxgxxm
xgxx
x
xf
xf
xx
mm
m
+
=
+?
=?=
*
)(lim
*
xx
xx
=
→
令
**
)( xx ≡?
)(
1
1]
)()()(
)(
1[lim
)()()(
)()(
lim
)()(
lim
*'
'*
*
'
'*
*
*
*
*
**
x
mxgxxxmg
xg
xx
x
xgxxxmg
xgxx
x
xx
xx
xx
xxxx
=?=
+
=
+
=
→
→→
由于所以当是,1)(0
*'
<< x?
*
x 0)( =xf的多重根时，牛顿法是线性收敛的。
例6.14 是的几重根？取0=
x 0221)(
22
== xxexf
x
5.0
0
=x，分别用牛顿公式与求重根的修正牛顿公式计算此根的近似值，精确至
4
10
)( ≤
k
xf,
解,
22
221)( xxexf
x
= 0)0( =f
,xexf
x
422)(
2
=′ 0)0( =′f
,44)(
2
=′′
x
exf 0)0( =′′f
95
96
x
,exf
2
8)( =′′′ 08)0( ≠=′′′f
因而是0
*
=x 0)( =xf的三重根。
牛顿迭代公式为
L,2,1,0,
)(
)(
1
=
′
=
+
k
xf
xf
xx
k
k
kk
(1)
计算得
k 0 1 2 3 4
0.5 0.34805 0.23905 0.16264 0.10993
(
k
xf 0.21828 0.067537 0.020617
3
102347.6
×
3
10837.
×
k
x
)
1
k 5 6 7
0.074659 0.050085 0.033530
)(
k
xf
4
107621.5
×
4
107180.1
×
5
101116.5
×
k
x
求重根牛顿迭代公式为
L,2,1,0,
)(
)(
3
1
=
′
=
+
k
xf
xf
xx
k
k
kk
(2)
计算得
k 0 1 2
0.5 0.044161 0.00032687
(
k
xf 0.21828
4
101741.1
×
10
1012006.3
×
k
x
)
容易看出公式(2)比公式(1)收敛快得多。
例6.15 设在[上有高阶导数，是) ](xf,ba ),(
*
bax ∈ ( ) 0=xf的m重根，且牛顿法收敛，证明牛顿迭代序列[有下列极限关系,
)2(≥
∞
=0
]
kk
x
m
xxx
xx
kkk
kk
k
=
+?
+?
∞→
11
1
2
lim
这个结论有什么用处？
证 牛顿迭代格式为
L,2,1,0),(
1
==
+
kxx
kk
其中
)(
)(
)(
xf
xf
xx
′
=?
且有
**
)(lim)(
*
xxx
xx
==
→
mh
xhx
x
h
1
1
)()(
lim)(
**
0
*
=
+
=′
→
(1)
由
))(()()(
111+
′=?=?
kkkkkkk
xxxxxx ξ
其中介于与之间，得到
1?k
x
k
x
k
ξ
)(1
1
))((
2
11
1
11
1
11
1
kkkkkk
kk
kkkk
kk
kkk
kk
xxxx
xx
xxxx
xx
xxx
xx
ξ?ξ? ′?
=
′+?
=
+?
=
+?
+?
+?
对上式两边取极限，并利用及(1)，得到
*
lim x
k
k
=
∞→
ξ
m
m
xxx
xx
kkk
kk
k
=
=
+?
+?
∞→
)
1
1(1
1
2
lim
11
1
由于牛顿法对重根是一阶局部收敛的。当适当大时，k
11
1
2
+?
+?
kkk
kk
xxx
xx
接近于重数。
求出之后，再应用求重根的修正牛顿迭代法进行求解，而后者具有二阶收敛速度。
m
m
例6.16 用割线法求方程在区间[内的一个实根，精确至5位有效数。 093
23
=+ xxx ]2,1
解 记
9)1)3((93)(
23
+=+= xxxxxxxf
割线法公式为
97
L,3,2,1),(
)()(
)(
1
1
1
=?
=
+
kxx
xfxf
xf
xx
kk
kk
k
kk
取，计算得 6.1,4.1
10
== xx
k 0 1 2 3 4 5 6
1.4 1.6 1.52967 1.51069 1.52417 1.52511 1.52510
(
k
xf -2.168 1.176 0.0692609 -0.216464 -0.0140970
4
1017173.
×
5
1041118.3
×
k
x
)
1?
所以。 5251.1
*
≈x
例6.17 设为),(],,[)(
*2
baxbaCxf ∈∈ 0)( =xf的单根，取),(
0
bax ∈。在牛顿法中，若用
0
0
)()
x
xf
k
(
k
xf
(
x
xf
k
近似代替导数)′，得到单点割线法
L,3,2,1),(
)()(
0
0
1
=
=
+
kxf
xfxf
xx
xx
k
k
k
kk
证明单点割线法是局部收敛的，且收敛阶一般为1,
证 迭代函数为
)()(
)()(
)(
0
0
xfxf
xfxx
xx
=?
2
0
000
)]()([
)()()()]()()][()()([
1)(
xfxf
xfxfxxxfxfxfxxxf
x
′′?+
=′?
)(
)()()(
)]([
)]()[()(
1)(
0
*
0
*
0
2
0
0
*
0
*
*
xf
xfxxxf
xf
xfxfxx
x
′?+
=
′?
=′?
由泰勒展开
2*
0
*
0
**
0
))((
2
1
))(()()( xxfxxxfxfxf?′′+?′+= ξ
其中ξ介于与之间，得
0
*
x x
98
))(()(2
))((
))((
2
1
))((
))((
2
1
)(
*
0
*
*
0
2*
0
*
0
*
2*
0
*
xxfxf
xxf
xxfxxxf
xxf
x
′′+′
′′
=
′′+?′
′′
=′
ξ
ξ
ξ
ξ
由于在〔a,b〕在有界，且，当适当靠近时，有 )(xf ′′ 0)('
*
≠xf
0
x
*
x
1)(
*
<′ x?
由局部收敛性定理知单点割线法是局部收敛的，且收敛阶一般为1,
例6.18 设，求
T
x )8513(?=,,,
21
xxx
∞
解
1139985)1(3
8}8,5,1,3max{
178513
2222
2
1
==++?+=
=?=
=++?+=
∞
x
x
x
例6.19 设，证明
n
Rx∈
1）
∞∞
≤≤ xnxx
1
2)
∞∞
≤≤ xnxx
2
证 1) 因为
∞
≤≤
=
≤≤
=
∞
≤≤
=
==≤=
=≥=
∑∑
∑
xnxnxxx
xxxx
j
nj
n
i
j
nj
n
i
i
i
ni
n
i
i
1
1
1
1
1
1
1
1
maxmax
max
所以
∞∞
≤≤ xnxx
1
3） 因为
∞
≤≤
=
≤≤
=
∞
=
==≤=
=≥=
∑∑
∑
xnxnxxx
xxxx
j
nj
n
i
j
nj
n
i
i
i
n
i
i
2
1
1
2
1
1
2
2
2
1
2
2
)max()max(
}max{
所以
99
∞∞
≤≤ xnxx
2
例6.20 已知,求
=
61
34
A
21
,,AAA
∞
解
{ }
{}761,34max
963,14max
1
=++=
=++=
∞
A
A
=
=
4518
1817
61
34
63
14
AA
T
018)45)(17(
4518
1817
2
==
=? λλ
λ
λ
λ AAE
T
即 044162
2
=+? λλ
解得 ( ),130231,130231,130231
21
+=?=+= AA
T
ρλλ
() 130231
2
+== AAA
T
ρ
例6.21 若是对称矩阵，求证A
2
)( AA =ρ
证 因为A是对称的，所以A有个实特征值n
n
λλλ,,,
21
L和一组标准正交特征向量系使得
n
xxx L,,
21
nixxAx
iiii
,,2,1,1,
2
L===λ
不妨设
n
λλλ ≥≥≥ L
21
则
1
)( λρ =A
任一向量x可表示为
∑
=
=
n
i
ii
xx
1
α
则
100
∑∑
∑∑
==
==
==
==
n
i
ii
n
i
i
n
i
iii
n
i
ii
Axx
xAxAx
1
22
2
1
2
2
11
,λαα
λαα
设,1
2
=x则
1
1
2
1
1
2
1
2
1
22
2
λαλλαλα ==≤=
∑∑∑
===
n
i
i
n
i
i
n
i
ii
Ax
又因为
1
2
11
2
11
2
1
λλλ =?== xxAx
所以
)(max
1
2
1
2
2
AAxA
x
ρλ ===
=
例6.22 用雅可比迭代法和高斯－赛德尔迭代法解线性方程组
=+
=+
12
23
21
21
xx
xx
精确至3位有效数。
解 1) 雅可比迭代格式为
=?=
=
+
+
L,2,1,0,2/)1(
,3/)2(
)(
1
)1(
2
)(
2
)1(
1
kxx
xx
kk
kk
取计算结果如下：,
k
0,0
)0(
2
)0(
1
== xx
0 1 2 3 4
(
1
x 0 0.66667 0.50000 0.61111 0.58333
(
2
x 0 0.50000 0.16667 0.25000 0.19445
)k
)k
101
k 5 6 7 8 9
(
1
x 0.60185 0.59722 0.60031 0.59954 0.60005
(
2
x 0.20833 0.19908 0.20139 0.19985 0.20023
)k
)k
因而,200.0,600.0
*
2
*
1
≈≈ xx
2) 高斯－赛德尔迭代格式为
L,2,1,0
,2/)1(
,3/)2(
)1(
1
)1(
2
)(
2
)1(
1
=
=
=
++
+
k
xx
xx
kk
kk
102
,取计算结果如下,0,0
)0(
2
)0(
1
== xx
k 0 1 2 3 4 5
(
1
x 0 0.66667 0.61111 0.60185 0.60031 0.60005
(
2
x 0 0.16667 0.19445 0.19908 0.19985 0.19975
)k
)k
因而 。 200.0,600.0
*
2
*
1
≈≈ xx
例6.23 给定线性方程组
=++?
=+?+
=?+
=?+
17722
23823
1138
7510
4321
4321
321
431
xxxx
xxxx
xxx
xxx
1) 写出雅可比迭代格式和高斯－赛德尔迭代格式；
2) 考查雅可比迭代格式和高斯－赛德尔迭代格式的敛散性。
解 1) 雅可比迭代格式为
+?=
=
+?=
+=
+
+
+
+
7/)2217(
)8/()2323(
8/)311(
10/)57(
)(
3
)(
2
)(
1
)1(
4
)(
4
)(
2
)(
1
)1(
3
)(
3
)(
1
)1(
2
)(
4
)(
3
)1(
1
kkkk
kkkk
kkk
kkk
xxxx
xxxx
xxx
xxx
高斯－赛德尔迭代格为
+?=
=
+?=
+=
++++
+++
++
+
7/)2217(
)8/()2323(
8/)311(
10/)57(
)1(
3
)1(
2
)1(
1
)1(
4
)(
4
)1(
2
)1(
1
)1(
3
)(
3
)1(
1
)1(
2
)(
4
)(
3
)1(
1
kkkk
kkkk
kkk
kkk
xxxx
xxxx
xxx
xxx
3) 由于所给线性方程组的系数矩阵
=
7221
1823
0381
51010
A
是严格对角占优的，所以雅可比迭代格式和高斯－赛德尔迭代格式均是收敛的
例6.24 给定线性方程组,其中，证明雅可比迭代法收敛，而高斯－赛德尔迭代法发散。
bAx =
=
321
011
101
A
证 1) 雅可比迭代矩阵J的特征方程为
0
321
01
10
=
λ
λ
λ
展开得
023
3
=++λλ (1)
记，则。又因为23)(
3
++= λλλf 019)(
2
>+=′ λλf
01( >?f
8
9
)
2
1
(,02) =?<?= f。所以（1）有惟一实根),
2
1
,1(
*
1
∈λ另二个根为共轭复根，即
*
3
*
2
,λλ
*
3
*
2
λλ =。由韦达定理有
3
2
*
3
*
2
*
1
=λλλ (2)
用牛顿迭代法 L,2,1,0,
)(
)(
1
=
′
=
+
k
f
f
k
k
kk
λ
λ
λλ
求,
*
1
k
λ
0 1 2 3 4
-1 -0.8 -0.75030 -0.74742 -0.747415
(
k
f λ -2 -0.336 -0.017444
-2.8629×
5
10
1.5094×
6
10
(
k
f λ′ 10 6.76 6.0666 6.0277
k
λ
)
)
所以。 747415.0
*
1
=λ
将代入(2)并利用
*
1
λ
*
2
*
3
λλ =得
891963.0
3
2
*
1
2
*
2
=
=
λ
λ
103
因而 94444.0891963.0
*
3
*
2
=== λλ
综上我们得到{ },194444.0,,max)(
*
3
*
2
*
1
<== λλλρ J所以雅可比迭代法收敛。
2) 高斯－赛德尔迭代矩阵G的特征方程为
0
32
0
10
=
λλλ
λλ
λ
展开得 0)1(3
2
=+λλ
解得 因而.1,0,0
*
3
*
2
*
1
=== λλλ 1)( =Gρ。所以高斯－赛德尔迭代法发散。
104
