浙江大学研究生学位课程
《实用数值计算方法,1
第六章 常微分方程及方程组的解法
6.1 常微分方程及其求解概述
6.2 初值问题解法
6.3 边值问题解法浙江大学研究生学位课程
《实用数值计算方法,2
6.1 常微分方程及其求解概述
6.2 初值问题解法
1.Euler方法
2.线性多步法
3.Runge--Kutta法
4.方程组及刚性问题的 Gear方法
6.3 边值问题解法
1.Shooting(试射法 )
2.差分法
Eq u a t io n salD if f e r e n t iO r d in a r yO D E s
浙江大学研究生学位课程
《实用数值计算方法,3
6.1 常微分方程及求解概述
（ Ordinary Differential Equations,ODE)
6.1.1 基本概念描述自由落体的 ODE:
的运动轨迹。
唯一确定了那么时时或设：
时若设：
tssO D E
st
st
sst
mg
dt
sd
16
36
01
00
261,00
16
1
0
'
0
2
2
浙江大学研究生学位课程
《实用数值计算方法,4
只有一个自变量的微分方程为 ODE,
否则称为偏微分方程 PDE。
方程中未知函数导数的最高阶数称为方程的阶。 （ 6-4）是二阶的
方程中关于未知函数及其各阶导数均是一次的，则称为线性微分方程。
4622 xrdxdyxqdx yd如：
和（ 6-1）都是线性二阶 ODE。
(6-2),(6-3)是 (6-1)的 初始条件 。亦称定解条件。 (6-1)(6-2)叫做 初值问题 。
6.1.1
浙江大学研究生学位课程
《实用数值计算方法,5
6.1.1
（ 6-1），（ 6-3）叫做 边值问题 。
在没有给定解条件时。方程一般有一族解曲线 y(x,c) 。如：
为任意常数。有解 ccecxy
y
dx
dy
x?
,
对任意的 n阶 ODE，如果能写成：
56,,,,1'nn yyyxfy?
则称该方程为 显式 的。方程（ 6-4）是显式的。而下面方程是 隐式 的。
xyyy '' e x ps in
浙江大学研究生学位课程
《实用数值计算方法,6
对于高阶显式方程。通过定义 n-1个新变量，可以写成 n维一阶方程组。
即令：
向量形式：一般我们把方程组写成组：可以改写为如下的方程那么
76
,,,,
56
66
1,,1,
110
'
1
1
'
2
2
'
1
1
'
0
0
1
nn
nn
i
i
yyyxfy
yy
yy
yy
xyxy
ni
dx
xdy
xy
6.1.1
浙江大学研究生学位课程
《实用数值计算方法,7
yxfy,'?
Tn
T
n
ffff
yyyy
110
110
,,,
,,,
这里：
在讨论初值问题时，我们从一阶方程开始：
00
',
yxy
yxfy
然后毫不费力地套用来解方程组。
当 f(x,y)与 y无关时,f(x,y)=g(x)
广泛的意义。故求解微分方程具有更变为积分问题则初值问题
00
00
'
,yxycdxxgy
yxy
xgy
6.1.1
浙江大学研究生学位课程
《实用数值计算方法,8
6.1.2 数值解及其重要性
难以求积。而的解是例如但仍然得不到精确解。
量的显示形式，有的解可以表示为自变故得不到精确解。不能表示为初等函数，
例如方程不定积分的组合表示。
及其的解很少能用初等函数
x
t
x
tx
dte
dteexy
xyy
yxy
O D E
0
0
'
22'
2
22
21
浙江大学研究生学位课程
《实用数值计算方法,9
6.1.3 ODE数值解的基本思想和方法特点基本思想有两点
1,离散化用 Taylor级数，数值积分和差商逼近导数等手段，把 ODE转化为离散的代数方程 (称差分方程 )。
2,递推化在具有唯一解的条件下，通过步进法逐步计算出解在一系列离散点上的值。从而得到原 ODE的数值近似解。
浙江大学研究生学位课程
《实用数值计算方法,10
6.2 初值问题解法我们讨论一阶 ODE，而高阶可能化为一阶 ODEs。一阶初值问题可以一般地写成：
86,,
00
0?
yxy
Xxxyxf
dx
dy
6.2.1 欧拉（ Euler） 方法
Euler方法是求解 (6-8)最 简单 方法，
但 精度差,故 不实用 。然而对理论分析很有用。
浙江大学研究生学位课程
《实用数值计算方法,11
6.2.1.1 方法原理及推导设初值问题 (6-8)满足：
2121
21
0
,,
,,,,,0
)
,:)
yyLyxfyxf
RyxyxL
L i ps c hi t zyfb
yXxxRfa
对常数条件：满足关于内连续。
在区域
内有界。在条件来代替：
验的可用更强的，但便于检就存在唯一的解。那么
R
y
f
ba
)),
86
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,12
为变步长。则称上述划分是任意的：
称为定步长。若其中等分，得到节点分成把
1
11
0
0
0
,,2,1
,,1,0
:,
i
iii
i
i
h
nihxx
nxXh
niihxx
xnXx
11110 xxxxx?
0y
0y
0y
6.2.1
图 6.1 常微分方程初值问题的数值解浙江大学研究生学位课程
《实用数值计算方法,13
公式。称为的近似。作为于是，我们取故
。其中级数构造：用
E ul e rxy
niyxy
yxhfyy
ni
xyxhfxy
xhyxyxy
xx
y
h
xhyxyhxy
T ay lori
i
iiii
iii
iii
ii
i
iii
96
96
1,,1,0
,
1,,1,0
,
2
)
1
00
1
'
1
1
''
2
'
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,14
欧拉方法的几何意义：
nx
Xxxx 210
0y
Xy
ny
1xy
2xy
1y
2y
h步长
6.2.1
图 6.2 Euler方法的几何意义浙江大学研究生学位课程
《实用数值计算方法,15
余项，得：
分，舍去若用梯形公式来计算积
。公式得即并使得舍去余项分，得并用左矩形公式计算积令改写成积分形式把用数值积分法构造
96
,,
,
,
,
86
)
1
E ul e r
yhxyR
Rxyxhfxyhxy
xx
dxxyxfxyhxy
ii
iii
iiiii
i
hx
x
106
1,,1,0
,,
2
00
111?
niyxy
yxfyxf
h
yy iiiiii
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,16
,改为以下预估一校正式一般地，把隐式格式隐式。半隐法：可化为显式的公式隐式：如改进的公式显式：如全部算出可计算单步法：由公式。称为改进的这个公式
106
)
106)
96)
)
)96(
1
e
E u le rd
E u le rc
yyya
E u le r
iii
1161,,1,0
,,
2
,
00
0
111
0
1
niyxy
yxfyxf
h
yy
yxhfyy
iiiiii
iiii
校正予估
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,17
6.2.1
0x 2?ix 1?ix?2?ixix1?ix nx
多步法单步法自动起步显式隐式半隐式图 6.3 ODE求解方法的类型浙江大学研究生学位课程
《实用数值计算方法,18
hvss
hvv
s
v
E ul e r
gm
v
s
v
dt
ds
mg
dt
dv
ts
ts
mg
dt
sd
mmm
mm
11
1
0
0
0
'
0
2
2
10
0
10
.10,1
100
00
10
0
方法：则：由加速度令质量初值问题：
20
0
5
10
ttvts
tvtv
真解：
m t v
m
s
m
m t v
m
s
m
s (t )
0
1
2
3
4
0.
1.
2.
3.
4.
10.
0.
-1 0,
-2 0,
-3 0,
0.
10.
10.
0.
-2 0,
0
1
2
3
4
5
6
7
8
0.
0,5
0.
-5,0
-1 0,
-1 5,
-2 0,
-2 5,
-3 0,
10.
5.
0.
-5,
-1 0,
-1 5,
-2 0,
-2 5,
-3 0,
0.
5.
7,5
7,5
5.
0.
-7,5
-1 7,5
-3 0,
0.
3,7 5
5.
3,7 5
0.
-6,2 5
-1 5,
2 6,2 5
-4 0,
0.1?h 5.0?h
表 6.1 自由落体运动方程的 Euler公式求解
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,19
t4321 t4321
40
30
20
10
0
10
40
30
20
10
0
10
v s
图 6.4 运动轨迹
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,20
图 6.5 Euler公式的误差
6.2.1.2 Euler方法的误差估计一般其它方法的误差估计也类似。
这里误差是指 截断误差 (算法理论误差 )
而不是 舍入误差 。后者由计算机字长等决定，属于稳定性问题。
i) 几何分析
ii yxy?
iy
ix 1?ix
1?ixy
1?iy
xy
iiii yxhfyy,1
浙江大学研究生学位课程
《实用数值计算方法,21
6.2.1
证明作为习题）
这里：
或有方法而对改进的或则有若满足生。称为该误差由当前一步所产有是精确的，则对设定量分析
(
,,
,
2
106
,
2
,,
!2
96
)
0
'''
3
3
2
2
0
''
''
2
11
XxxMxy
hORM
h
R
E u le r
hORM
h
R
XxxMxy
Ry
h
yxy
xyy
ii
ii
ii
iiii
ii
局部截断误差，
浙江大学研究生学位课程
《实用数值计算方法,22
6.2.1
整体截断误差,设 ym是 Euler公式（ 6-9）
精确解，而 y(x) 是初值问题（ 6-8）的解。
则 整体截断误差定义为它是局部截断误差的积累。
定理,若 f(x,y)关于 y满足 Lipschitz条件。
则有估计式：
)126()1(01 LXLXm ehLRe
的步长的区间上限局部截断误差上界常数
xh
xX
R
L ips c hi t zL
mmm xyy
浙江大学研究生学位课程
《实用数值计算方法,23
6.2.1
hL
hL
RhL
hLRhL
hLRhL
RhL
Rxyxfyxfh
Rxyxhfxyhxy
yxhfyy
m
m
m
j
jm
m
mm
mmmmmmm
mmmmm
mmmm
11
1
11
111
1
,,
,
,
1
0
1
0
0
1
1
2
1
1
1
递推相减得：
证明：
浙江大学研究生学位课程
《实用数值计算方法,24
6.2.1
。方法，对改进同阶。整体截断误差与方法的说明可记为上式得：
代入，将若于是得：
注意到：
2
1
2
000
01
1
1
1
,
1
2
2
0
,2,1,0
1
11
1
hOE ul e r
h
E ul e rhO
e
L
hM
M
h
Rxyy
m
e
hL
R
e
ehLhL
Xhmx
m
m
LX
m
LXLX
m
LX
X
h
m
m
ex
x
x
11lim
浙江大学研究生学位课程
《实用数值计算方法,25
6.2.1
立是差分方程的解，若成初值问题的解是收敛性：设收敛性、稳定性讨论使用上述结果。是不知道的。因此难以的估计式中的实际上，
误差容易估计。
决定，但局部方法的精度由整体误差误差低一阶。一般来说，整体比局部说明：
m
y
O D Exy
MMxyiii
ii
i
,
.3.1.2.6
''
0m a xlim 10 mmnmh xyy
浙江大学研究生学位课程
《实用数值计算方法,26
6.2.1
况。严格定义为：的差分方程解的变化情初值的改变导致公式不能无限缩小。所以实际上对对解的影响。误差的积累可能扩大了时离散点增多。舍入当方法是收敛的，且阶为
。则称差分方程是收敛的
96
0
h
h
hOE ul e r
注意,
稳定性,
,2,1,
,0
,
,,,,
00
0
000
mvucvu
Xmhhh
vu
vuhc
mm
mm
对当满足：
它们的解使对定义：
浙江大学研究生学位课程
《实用数值计算方法,27
6.2.1
证毕。
证明：记方法是稳定的。则条件，的事实上，若方法是稳定的。则称
0
0
1
1
1
1
1
1
,,
,
,
,
LX
m
m
mmmmmm
mmmm
mmmm
mmm
e
hL
hL
vxfuxfh
vxhfvv
uxhfuu
vu
E ul e r
L ips c hi zyxf
E ul e r
浙江大学研究生学位课程
《实用数值计算方法,28
6.2.1
的区域范围。的存在稳定实际上，对给定的问题为常数。这里可推知和初值，由对给定的步长
。而不是总是适当大小，限状况。实际计算中，
一种极对以后解的影响。它是始误差时，初稳定性实际上描述了当
0
0
,
0
0
yh
c
mncxyy
xyy
h
h
h
nn
mm
绝对稳定：
浙江大学研究生学位课程
《实用数值计算方法,29
6.2.1
时可提高一阶精度。
故取
：的一次项为组合后，可使时，有而步长为法的近似解：
外推法外推技巧
xyxyxy
h
h
hchcxyxy
E ul e r
R ic har ds on
h
h
h
2
2
21
2
0
2
.4.1.2.6
2212 4121 hchcxyxy h
222 hOxyxyxy hh
浙江大学研究生学位课程
《实用数值计算方法,30
6.2.1
阶精度。它具有故可取则阶方法计算若用
1
122
1462
136
:,
1
2
1
2
1
2
p
hOxyxyxy
hOhcxyxy
hOhcxyxy
xyxyp
pph
h
p
Pp
p
p
h
pp
p
h
h
h
156212 1 2 xyxyxy hhpp
浙江大学研究生学位课程
《实用数值计算方法,31
6.2.1
12
21
146136
pp
p
p
h
h
hOhCxyxy
：，由估计。外推法用于实用的误差
。和误差估计的有效方法外推法是提高精度因此故
R ic har ds on
xyxyhC
h
h
p
p
p
p
166
12
2
2
浙江大学研究生学位课程
《实用数值计算方法,32
6.2.1
。果见表预估校正方法。比较结用的数值解求用其真解为例：
2.6
1.0,21
10
1,0,
2
E ul e r
hxy
y
x
y
x
y
dx
dy
x
i
y
i
y( x
i
)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.0
1.0959
1.1841
1.2662
1.3434
1.4164
1.4860
1.5525
1.6165
1.6782
1.7379
1.0
1.095 4
1.18 32
1.26 49
1.34 16
1.41 42
1.48 32
1.5 492
1.61 25
1.67 33
1.73 21
表 6.2 予估校正求解结果对比浙江大学研究生学位课程
《实用数值计算方法,33
6.2.1
时的解及误差估计：技巧解此方程，
方法及外推。用其真解为例：
1
10
1,0
'
x
E ul e rey
y
xyy
x
1
1 / 2
1 /4
1 /8
1 /1 6
1 /3 2
,0 0 0 0 0
,2 5 0 0 0
,3 1 6 4 1
,3 4 3 6 1
,3 5 6 0 7
,3 6 2 0 6
,5 0 0 0 0
,3 8 2 8 1
,3 7 0 8 1
,3 6 8 5 4
,3 6 8 0 4
,3 6 7 9 2
-,3 6 7 8 8
-,1 1 7 8 8
-,0 5 1 4 1
-,0 2 4 2 7
-,0 1 1 8 1
-,0 0 5 8 2
-,5 0 0 0
-,1 3 2 8 1
-,0 5 4 4 1
-,0 2 4 9 3
-,0 1 1 9 6
-,0 0 5 8 6
,1 3 2 1 2
,0 1 4 9 3
,0 0 2 9 3
,0 0 0 6 6
,0 0 0 1 6
,0 0 0 0 4
h hy hh yy?22 hy? 估计? 外推?
表 6.3 Euler法与外推结果的比较浙江大学研究生学位课程
《实用数值计算方法,34
6.2.2 线性多步法
1
1
11
,1
1
1
1
00
10
:,
,
,,
1.2.2.6
,,1,
,,,
m
m
m
m
m
m
m
m
x
x
kmmm
mkmkm
x
x
mm
x
x
x
x
m
ii
m
dxxLyy
yxf
xxxL agr ange
dxyxfxyxy
dxyxfdyyxf
dx
dy
y
mixyyxyy
xxx
的插值多项式近似上，，，在用插值法欲求处已求得：假定在
外插法：Ad a m s
外插浙江大学研究生学位课程
《实用数值计算方法,35
6.2.2
k
i
imkimm
imimim
im
kmmm
imk
k
i imimk
imk
imkm
fbhyy
A dam s
kiyxff
ki
xx
xxxxxx
xD
xD
xD
fxL
0
1
1
,
0,
,
,
,,1,0,
,,1,0
外插值公式：由此得这里：
(?)
方法。时，线性多步法就是特别，当
。其值见表其中
E u le rk
dx
xD
xD
h
b
m
m
x
x
imimk
imk
ki
0
4.6
1 1
,
,
(?)
浙江大学研究生学位课程
《实用数值计算方法,36
6.2.2
i 0 1 2 3 4 5
b
0 i
2 b
1 I
12 b
2 I
24 b
3 I
720 b
4 I
1440 b
5 i
1
3
23
55
1901
4277
-1
-1 6
-5 9
-2 7 7 4
-7 9 2 3
5
37
2616
9982
-9
-1 2 7 4
-7 2 9 8
251
2877 -4 7 5
211 51623
12
2
mmmmm fff
hyy
k 时，线性多步法为例如：
Adams 外插法 (k=2) 3阶 3步 显式
2?m
1?m
m
表 6.4 外插系数 bki值图 6.6 3阶 3步外插法
浙江大学研究生学位课程
《实用数值计算方法,37
6.2.2
，舍入误差也小。内插法较外插法更精确最常用，其插值公式为法。时，即得到改进当
。它是隐式或半隐式格式
：多项式的取插值节点为内插法：
1
1,1
,,,,
,
p
E ul e rkp
xLL agr ang e
xxx
A dam s
p
km
pmmpkm
1
,1
m
m
x
x
p
kmmm dxxLyy
。见表其中 5.6
0
11
i
i
k
k
i
imkmm
b
fbhyy
浙江大学研究生学位课程
《实用数值计算方法,38
6.2.2
i 0 1 2 3 4 5
b 0i
*
2 b 1i
*
1 2 b 2i
*
2 4 b 3i
*
7 2 0 b 4i
*
1 4 4 0 b 5i
*
1
1
5
9
251
475
1
8
19
646
1427
-1
-5
-2 6 4
-7 9 8
1
106
482
-1 9
-1 7 3 27
值内插系数表?kib5.6
111 85
12
12
mmmmm fff
hyy
pk 时，线性多步法为，例如：
Adams 外插法 (k=2) 4阶 3步
1?m1?m m
图 6.7 4步 3阶 Adams内插公式
浙江大学研究生学位课程
《实用数值计算方法,39
6.2.2
外插法：不稳定。
稳定内插法：舍入误差小，
稳定性：
单。外插法：显式，计算简
，复杂内插法：隐式或半隐式迭代格式：
阶方法外插法：
阶方法内插法：
局部截断误差：
)
)
1
2
)
2
3
iii
ii
khO
khO
i
k
k
1?p
浙江大学研究生学位课程
《实用数值计算方法,40
6.2.2
k
i
ii
k
ii
ihxyhihxyhxy
k
k
0
'
00
0,
,0
,
2.2.2.6
；
定义算符：
步法不同时为步法线性多步法是常数一般的插值型公式待定系数法
k
i
k
i
imiimi fhy
0 0
1?k?取浙江大学研究生学位课程
《实用数值计算方法,41
6.2.2
0x
1?mx mx 1?mx
kmx?
my
nk
k
nk
k
k
f
f
h
y
y
y
0
0
0
1
0
0
0
0
mmm yxff
hkN
,
,1
个方程。等步长图 6.8 一般化插值形式浙江大学研究生学位课程
《实用数值计算方法,42
6.2.2
k
rr
k
rr
r
k
i
ii
k
rr
r
k
r
k
r
C
iC
C
xyhC
xhyCxyChxy
xy
ih
xihyxyihxy
xy
ih
xihyxyihxy
T ay lor
1
2
1
1
21
0
1
100
'
10
'''
2
''''
''
2
'
2
!1
1
2
!
1
!2
!2
其中：；
代入得：
展开右端：
浙江大学研究生学位课程
《实用数值计算方法,43
6.2.2
步方法。
阶阶及外插及内插法为隐式显式由此得于是：
充分大，使次连续可微，则可取有因此，若
k
kkA dam s
xyhCfhy
CCCk
rxy
k
k
k
i
rr
rihxii
r
1
0
0
.0
1
0
11
1
10
k
i
k
i
imiimi fhy
kr
0 0
步方法：阶线性浙江大学研究生学位课程
《实用数值计算方法,44
6.2.2
ak 02,1,2
2
记步法：例，考虑
04
!2
1
8
!3
1
024
!2
1
02
01
2113
2112
21011
10
C
C
C
aC
解之得,
a
a
a
a
5
12
1
1
3
2
,51
12
1
,1
2
1
0
1
浙江大学研究生学位课程
《实用数值计算方法,45
6.2.2
方法型：它是以下阶方法时，
阶的方法是时，
故同时，可以算出：
故一般二步法为：
eMiSi m ps on
CCa
Ca
ln
40,01
3,01
54
4
mmm
mmm
fafafah
ayyay
51185
12
1
12
12
aC
aC
1317
!53
1
1
!4
1
5
4
mmmmm fffhyy 122 43
浙江大学研究生学位课程
《实用数值计算方法,46
6.2.2
1,1
01
0,2
11
1
3.2.2.6
11
1
'
即特征根的特征方程：
步长方法：例对典型方程：
绝对稳定性
h
h
hh
h
yh
yhy
hfyy
hE u le r
yy
m
mm
mmm
来讨论绝对稳定性）
程一般只限于这个典型方(
浙江大学研究生学位课程
《实用数值计算方法,47
6.2.2
稳定可以任意取，步长时，恒有当其特征根为隐式方法例，后退
A
h
h
y
h
y
yhy
hyyy
E ul e r
mm
mm
mmm
10Re
1
1
1
1
1
1
1
1
'
11
图 6.9
浙江大学研究生学位课程
《实用数值计算方法,48
6.2.2
线性多步法的绝对稳定性：
次多项式的其特征方程为多步法为典型方程则常数设
k
hh
yhy
yy
y
f
k
i
imi
k
i
imi
00
'
k
i
im
i
k
i
im
i h
00
0
浙江大学研究生学位课程
《实用数值计算方法,49
6.2.2
定义,
，称为，取值范围绝对稳定的则称线性多步法关于此
，按模小于，若特征方程的根给定
h
h
h i 1
绝对稳定 。
绝对稳定区域。
.0,,0,
01
01
1
0
1
0
1
EI
k
i
i
k
k
k
i
i
k
k
i
i
bh
bh
A d am s
其稳定区域为程：内插及外插法的特征方
浙江大学研究生学位课程
《实用数值计算方法,50
6.2.2
的取值表 EI,5.6
k 1 2 3 4
E
2? 1
I
6? 3
11
6?
49
90?
10
3?
5
3
1
2
5
1
2
22
0
3
1
1
3
4
3
1
14
3
1
hOe
hOe
hh
h
h
h
h
它的根：
方法的特征方程
Milne
浙江大学研究生学位课程
《实用数值计算方法,51
6.2.2
：见表及误差，得到的解对不同的方法计算例，用对稳定。
方法不绝，因此，对任何
6.6
10
10
ln
ln
'
yh
y
yy
eMi
eMih
h u 误 差
0,1
0,0 1
0,0 0 1
0,0 0 0 1
0,0 0 0 0 1
,4 8 6 E - 4
,4 8 2 E - 4
,3 8 6 E - 4
,3 2 9 E - 4
,5 7 8 E - 4
-,3 2 3 E - 5
-,2 8 1 E - 5
,6 7 6 E - 5
-,3 7 5 E - 4
-,1 2 4 E - 4
表 6.6 计算结果浙江大学研究生学位课程
《实用数值计算方法,52
6.2.2
6.2.2.4 预估 --校正方法
（ Predictor--Corrector Method）
迭代法。这是个非线性方程已知。，，，，，其中阶隐式方法：但计算困难。现考虑其精度高，稳定性好。尤隐式方法：
110
0
kify
k
y
f
imim
176
,
1
0
1
0
k
i
imi
k
i
kmkmkimikm
fh
yxfhyy
，，
，
10
186
1
0
1
0
1
n
fh
yxfhyy
k
i
imi
k
i
n
kmkimi
n
km km
浙江大学研究生学位课程
《实用数值计算方法,53
6.2.2
。此时应减小步长
。太多的校正不合适。可满足二，三步校正。，一般只要对给定算法。通称为算式。称为校正式，即算式。称为预估式，即构成的算法，时，
选为显式算法而
h
yy
PC
CC or r e c t or
Pe di c t or
y
n
km
n
km
km
1
0
0
196
Pr186
196186
1961
0
1
0
0
k
i
k
i
imiimikm fhyy
浙江大学研究生学位课程
《实用数值计算方法,54
6.2.2
注意,一步校正的计算量? 预估计算量。
所以要适当选取 h才能发挥 PC的优点。
设绝对稳定区域：
达到精度的校正次数为 N，则 h的选取，
应满足：
否则可用 N步显式算法稳定达到目的 。
0,
0,
c
p
SC
SP
算式：
算式：
p
cp
SNh
ShS
1.2
.1
hmx 1?mx
pS
cS
2?N
浙江大学研究生学位课程
《实用数值计算方法,55
6.2.2
Adams四步四阶 预估校正算法
123434
12334
937599
24
:54
9375955
24
:44
mmmmmm
mmmmmm
ffff
h
yyC
ffff
h
yyP
阶步阶步
0
,,
270
19
0
4
1
4
0
44
1
44
0
4
1
4
55
5
，
稳定性判别可用很大的麻烦但变步长会带来可用来控制步长方法误差：
C
mm
mmmm
mm
S
y
f
hh
yy
yxfyxf
y
f
h
yyyhC
浙江大学研究生学位课程
《实用数值计算方法,56
6.2.3 Runge--Kutta 方法
，考虑为构造高精度的单步法
。后者比前者精度高一阶型预估一校正法：
方法：
12
1
211
1
11
,
,
2
1
2
1
,
KyhxhfK
yxhfK
KKyy
E ul e r
yxhfK
Kyy
E ul e r
mm
mm
mm
mm
mm
12
21
,mm
mm
yhxhfk
kyy
E u le r 法后退浙江大学研究生学位课程
《实用数值计算方法,57
6.2.3
方法的级正整数，称为 KRN
Niba
KbhyhaxfK
yxfK
KChyx
hyxhyy
i
l
ili
i
l
lilmimi
mm
N
i
iimm
mmmm
,,3,2
,
,
,,
,,
1
1
1
1
1
1
1
上的线性组合：在某些点是而 1,?mm xxf?
其中：
阶次项一致），，，即项完全一致。级数。取的展开，并比较进行把
PP
PT ay lor
xyT ay lorK mi
10(
1
1
浙江大学研究生学位课程
《实用数值计算方法,58
6.2.3
时，有当公式的阶数。称为这里的就可以得到系数
3?
N
KRP1,,1;,,1, ilNibC ili
232132333
1222
1
3
1
1
,
,
,
KhbKbahyhaxfK
KhayhaxfK
fyxfK
KChyy
mm
mm
mmm
i
iimm
浙江大学研究生学位课程
《实用数值计算方法,59
6.2.3
的系数得：比较而其中：
展开：在将
3,2,1,0
6
1
2
2
2
2
,,
43
2
2
''
32
33
2
223223
2
3
1
3322321
32
ih
hOGFfh
F
h
hfxyhxy
fffffG
xyfffF
hOGaCaCFfbaC
h
FaCaChfCCCKC
T ay loryxKK
i
y
mm
yyxyxx
yx
y
i
ii
mm
浙江大学研究生学位课程
《实用数值计算方法,60
6.2.3
6
1
3
1
2
1
1
3223
2
33
2
22
3322
321
baC
aCaC
aCaC
CCC
六个未知数，二个自由,故可取
21
2
1
6
1
3
2
6
1
3232
321
baa
CCC
，，
，，
213
12
1
3211
2
2
1
2
1
,
4
6
1
hKhKyhxfK
hKyhxfK
yxfK
KKKhyy
mm
mm
mm
mm
，
，
故：
浙江大学研究生学位课程
《实用数值计算方法,61
6.2.3
N=4,四级四阶 R--K方法
34
23
12
1
43211
,
2
1
,
2
1
2
1
,
2
1
,
22
6
hKyhxfK
hKyhxfK
hKyhxfK
yxfK
KKKK
h
yy
mm
mm
mm
mm
mm
----最常用的古典 Runge-Kutta方法。
增加计算函数值的次数（级）与提高精度（阶）的关系见下表,级 N 2 3 4 5 6 7 n? 8
阶 P 2 3 4 4 5 6 n? 2
表 6.7 Runge-Kutta方法中 级 与 阶 的关系浙江大学研究生学位课程
《实用数值计算方法,62
6.2.3
各种解法精度比较：
例高精度方法的主要优缺点
10
22
).1(
2'
y
yxxy
KR
X
m
E u l er 改进 E u l er R— K ( 4 阶 ) 精确解
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.0 000
1.2 000
1.4 420
1,7384
2,1041
2,5569
3,1 183
3.8 139
4.6 747
5.7 376
7.0 472
1.0 000
1.2 2 10
1.4 9 23
1.8 284
2.2 466
2,7680
3,4174
4,2257
5,2288
6,4704
8,0032
1.0 000
1.2 221
1.4 977
1.8 432
2.2 783
2.8 274
3.5 20 1
4.3 92 7
5.4 89 4
6.8 64 3
8.5 83 4
1.0 000
1.2 221
1.4 977
1.8 432
2.2 783
2.8 274
3.5 202
4.3 928
5.4 895
6.8 645
8.5 836
表 6.8 各种解法在例题中的结果比较浙江大学研究生学位课程
《实用数值计算方法,63
6.2.3
（ 2）,单步法，自动起步
（ 3）,易改为变步长
（ 4）,绝对稳定区域较同阶线性多步法大
（ 5）,计算工作量较大，有时大于隐式方法
（ 6）,估计误差不易绝对稳定性讨论
mm
mm
y
hh
hy
f
h
f
h
hfyy
yyN
3
3
2
2
''
3
'
2
1
'
!3!2
!3!2
3
此时：
为例。取典型方程以浙江大学研究生学位课程
《实用数值计算方法,64
6.2.3
D
KR
h
hhh
hhh
y
y
hh
hy
mm
im
im
mm
区域方法的绝对稳定一般地，各级时，，当其根为：
故特征方程：
1051.2
6
1
2
1
1
0
6
1
2
1
1
!3!2
1
1
32
1
32
1
32
1
级? 1 D
1
2
3
4
(-2,0)
(-2,0)
(-2.51,0)
(-2.78,0)
432
32
2
24
1
6
1
2
1
1
6
1
2
1
1
2
1
1
1
hhhh
hhh
hh
h
表 6.9 各级 R-K方法的绝对稳定区域浙江大学研究生学位课程
《实用数值计算方法,65
6.2.3
0,78.24,2.0
0,78.222.0
10
20
2
2
1
1
'
hh
hh
y
yy
KR
h
hKR
，
取方法计算用四级经典例绝对稳定区域。
使方法时，要取在使用
x h
1
时的误差 h
2
时的误差
0.0
0.2
0.4
0.6
0.8
1.0
0
-0.092795
-0.012010
-0.001366
-0.000152
-0.000017
0
4.98
25.0
125.0
625.0
3125.0
表 6.10 不同步长对精度的影响浙江大学研究生学位课程
《实用数值计算方法,66
6.2.3
。可能导致舍入误差积累的计算量，步长太小：增加不必要非绝对稳定。步长太大：精度不够，
实际步长选取：
外推技巧。需用选取误差估计和实际步长的属于绝对稳定区域。
近似确定步长。使可以用常数的问题对于
R ic har ds on
xhxh
x
x
y
f
mm
浙江大学研究生学位课程
《实用数值计算方法,67
6.2.3
是适当的。，则认为，若对给定误差估计方法：对四阶方法自选步长的
2
12
2
2
1
1
2
2
11
4
4
4
54
4
2
11
5
4
1
2
1
54
11
h
yych
hOchyy
hO
h
cxyy
hOchxyy
KR
KR
h
m
h
m
h
m
h
m
m
h
m
m
h
m
浙江大学研究生学位课程
《实用数值计算方法,68
P阶图 6.10 变步长 Runge--Kutta方法框图
00,0,,,xyym PhX给定
211 211m a x,hmhm hmhm yy yy 计算
hmm
mm
yy
hxx
11
1
hh?2
Xxm1
121 P
,2,1
,
m
yx mm输出
hh?2
1mm
Y
Y
N
N
Y
N
hh
25.0
9.0?
hh?4
hh 2.09.0
5
49.0
6.2.3
浙江大学研究生学位课程
《实用数值计算方法,69
展开：处在 T a y lo ryxKK mm,,32
xym
yy
m
xx
ymxm
m
m
m
m
m
m
m
m
m
fKbfbaah
f
Kbfbah
f
ha
fKbfbahfhafK
yx
f
fha
y
ffha
x
fha
y
f
fha
x
f
hafK
2323233
2
2
232323
22
3
23232333
2
2
22
22
2
2
22
2
222
!2!2
!2!2
1
1
1
1
1
1
1
1
1
!
,,3,2
,
,
P
P
m
P
mm
i
j
iij
i
j
jijmimi
mmm
N
i
iimm
hO
P
h
xyxyxy
ab
Ni
KbhyhaxfK
fyxfK
KChyy
6.2.3
浙江大学研究生学位课程
《实用数值计算方法,70
1.0
00
0
!
0
!2
0
!1
0
)
)
,,,
,,,
2'
1
1
''
2
'
110
11
h
y
yxy
T a y l or
R
Ry
s
x
y
x
y
x
yxy
T a y l orxyb
KR
a
yyy
kyyyy
s
s
s
s
k
kmkmmm
前三个点的值级数求初值问题例：用不影响所需精度为止。直至展开式：的用精确解计算所需点上的值。
方法，或用小步长高精度用单步法求出方法取得。
的值可以通过两种出发值步方法
6.2.4 出发值的计算浙江大学研究生学位课程
《实用数值计算方法,71
00499 950 00.0000500 050 00.0
160
1
20
1
2
1
0000000
0
721352
615102
5102
432
32
2
21
,,
101
852
764''''
76'5''4'''8
65'4''2'''7
54''''''6
4''''2''5
''''''4
''2''''
'''
2'
2
yyy
xxxxy
yyyyyy
x
yyyyyyyyxy
yyyyyyyxy
yyyyyyxy
yyyyyxy
yyyyxy
yyyxy
yyxy
yxxy
yxyxf
，，
于是，
时，故所以解：由于
25206010 85'' yyy,，
得，，，取前三点 1.00 101 hhxxhx
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,72
质量控制 Runge--Kutta 步进程序
SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)
PARAMETER (NMAX=10,PGROW=-0.20,PSHRINK=-0.25,FCOR=1./15.
ONE=1.0,SAFETY=0.9,ERRCON=6.E-4
EXTERNAL DERIVS
DIMENSION Y(N),DYSX(N),YSCAL(N),YTEMP(NMAX),YSAV(NMAX),
DYSAV(NMAX)
XSAV=X 保留初值
DO 11 I=1,N
YSAV(I)=Y(I)
DYSAV(I)=DYDX(I)
H=HTRY
HH=0.5?H
CALL RK4 (YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)
X=XSAV+HH
CALL DERIVS (X,YTEMP,DYDX)
CALL RK4 (YTMP,DYDX,N,X,HH,Y,DERIVS)
X=XSAV+H
IF (X,EQ,XSAV) PAUS E ‘步长无意义’
CALL RK4 (YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)
ERRMAX=0.
DO 12 I=1,N
YTEMP(I)=Y(I)-YTEMP(I)
ERRMAX = MAX(ERRMAX,ABS(YTEMP(I) / YSCAL(I)))
ERRMAX = ERRMAX / EPS
误差计算一个单步计算两个半步长计算置初始值
11
1
12
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,73
PARAMETER (NVAR=4)
DIMENSION VS(NVAR)
COMMON /PATH/ KMAX,KOUNT,DXSAV,X(200),Y(10,200)
EXTERNAL DERIVS RKQC
X1=1.0
X2=10.0
VS(1)=BESSO(X1)
VS(2)=BESSI(X1)
VS(3)=BESSJ(2,X1)
VS(4)=BESSJ(3,X1)
EPS=1.0E-4
HI=1.0
HMIN=0.0
KMAX=100
DXSAV=(X1-X2)/20.0
CALL ODEINT (VS,NVAR,X1,X2,EPS,HI,HMIN,NOK
NBAD,NOK,NBAD,DERIVS,RKQC)
WRITE (?,?)
END
SUBROUTINE DERIVS (X,Y,DYDX)
DIMENSION Y(1),DYDX(1)
DYDX(1) = -Y(2)
DYDX(2) = Y(1)-(1.0/X)?Y(2)
DYDX(3) = Y(2)-(2.0/X)?Y(3)
DYDX(4) = Y(3)-(3.0/X)?Y(4)
RETURN
END
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,74
IF (ERRMAX.GT.ONE) THEN 不满足要求
H = SAFETY?H?(ERRMAXPSHRINK) 估计下次步长（缩小）
GOTO 1
ELSE 满足要求
HDID = H
IF (ERRMAX,GT,ERRCON) THEN
HNEXT = SAFETY?H? (ERRMAXPGROW) 估计可放大的步长
ELSE
HNEXT =4,?H 最多可放大四倍
ENDIF
ENDIF
DO 13 I=1,N 完成五阶截断误差
Y(I)=Y(I) + YTEMP(I)?FCOR
CONTINUE
RETURN
END
四阶 Runge--Kutta 步进程序（续）
13
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,75
计算结果
NOK=30
NBAD=0
KOUNT=15
X I n t eg ral B E SS J ( 3,x )
1,0 0 0 0
1,6 3 2 9
2,1 0 3 6
2,6 8 1 3
3,3 7 6 9
4,1 0 3 9
4,7 9 1 0
5,5 1 4 2
6,1 4 2 5
6,8 4 1 3
7,4 9 5 9
8,2 1 1 2
8,9 1 7 4
9,5 8 5 5
1 0,0 0 0 0
0,0 1 9 5 6 3
0,0 7 6 5 6 4
0,1 4 5 8 8 2
0,2 5 0 5 4 4
0,3 7 0 1 1 5
0,4 3 3 3 9 8
0,3 9 6 3 8 4
0,2 5 2 4 6 0
0,0 7 1 7 1 0
-0,1 2 9 1 0 4
-0,2 5 7 5 3 6
-0,2 8 6 3 9 3
-0,1 9 7 3 8 9
-0,0 4 3 9 6 8
0,0 5 8 3 8 0
0,1 4 5 8 8 1
0,4 3 3 3 9 7
0,3 9 6 3 8 3
-0,1 2 9 1 0 3
-0,2 5 7 5 3 5
-0,2 8 6 3 9 1
-0,1 9 7 3 8 8
0,0 5 8 3 7 9
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,76
四阶 Runge--Kutta 方法：
SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
PARAMETER (NMAX=10)
DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX)
DYT(NMAX),DYM(NMAX)
HH=H?0.5
H6=H/6.0
XH=X+HH
DO 11 I=1,N
YT(I)=Y(I)+HH?DYDX(I) 第一步
CALL DERIVS(XH,YT,DYT)
DO 12 I=1,N 第二步
YT(I)=Y(I)+HH?DYT(I)
CALL DERIVS(XH,YT,DYM)
DO 13 I=1,N 第三步
YT(I)=Y(I)+H? DYM(I)
DYM(I)=DYT(I) + DYM(I)
CALL DERIVS(X+H,YT,DYT)
DO 14 I=1,N
YOUT(I)=DYDX(I)+DYT(I)+2,?DYM(I)
YOUT(I)=Y(I) + H6?YOUT(I)
CONTINUE
RETURN
END
14
13
12
11
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,77
用四阶 RUNGE--KUTTA公式，等步长计算 NSTEP步
SUBROUTINE RKDUMB(VSTART,NVAR,X1,X2,NSTEP,DERIVS)
PARAMETER (NMAX=10) 函数最大个数
COMMON /PATH/ XX(200),Y(10,200)
DIMENSION VSTART(NVAR),V(NMAX),DV(NMAX)
置初值
DO 13 K=1,NSTEP
CALL DERIVS(X,V,DV)
CALL RK4(V,DV，,NVAR,X,H,V,DERIVS)
X=X+H
XX(K+1)=X 存贮中间步骤
DO 12 I=1,NVAR
Y(I,K+1)=V(I)
CONTINUE
CONTINUE
结束返回
12
13
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,78
带自适应步长控制的四阶 Runge--Kutta 驱动程序
SUBROUTINE ODEINT(YSTART,NVAR,X1,X2,EPS,H1
HMIN,NOK,NBAD,DERIVS,,RKQC)
PARAMETER (MAXSTEP=10000,NMAX=10,TWO=2.0,
ZERO=0.0,TINY=1.0E-30)
COMMON /PATH/ KMAX,KOUNT,DXSAV,XP(200),YP(10,200)
DIMENSION YSTART(NVAR),YSCAL(NMAX),Y(NMAX),DYDX(NMAX)
X=X1
H=SIGH(H1,X2-X1)
NOK=0
NBAD=0
KOUNT=0
DO 11 I=1,NVAR
Y(I)=YSTART(I)
XSAV=X-DXSAV?TWO ----保证第一步存贮
DO 16 NSTEP=1,MAXSTEP ----最多取 MAXSTEP步数
CALL DERIVS(X,Y,DYDX)
DO 12 I=1,NVAR 改变比例因子控制准确性
YSCAL(I)=ABS(Y(I))+ABS(H?DYDX(I)) + TINY
IF(KMAX.GT.O.AND.ABS(X-XSAV).GT.ABS(DXSAV)
KOUNT,LT,KMAX-1) THEN
KOUNT=KOUNT+1
XP(KOUNT)=X
DO 13 I=1,NVAR
YP(I,KOUNT)=Y(I)
XSAV=X
ENDIF
存贮中间结果
13
12
11
置初值
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,79
IF ((X+H-X2)?(X+H-X1),GT,ZERO) H=X2-X1 越界处理
CALL RKQC(Y,DYDX,NVAR,X,H,EPS,YSCAL,HDID,HNEXT,DERIVS)
IF (HDID,EQ,H) THEN
NOK=NOK+1 好步长计算
ELSE
NBAD=NBAD+1 坏步长计算
ENDIF
IF ((X-X2)?(X2-X1),GE,ZERO) THEN 进行完毕
DO 14 I=1,NVAR
YSTART(I)=Y(I)
IF (KMAX,NE,O) THEN 保留最终步结果
KOUNT=KOUNT+1
XP(KOUNT)=X
DO 15 I=1,NVAR
YP(I,KOUNT)=Y(I)
ENDIF
RERURN
ENDIF
IF (ABS(HNEXT),LT,HMIN) PAUSE '步长太小 '
H=HNEXT
CONTINUE
PAUSE '步数太多，越界 '
RETURN
END
16
15
14
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,80
i
i
h
mm
h
m
h
m
m
h
m
m
h
m
y s c a l e
sh
sh
hhh
h
h
hOyxy
hOchyy
hO
h
cxyy
hOchxyy
K ut t aR unge
k ut t aR unge
m a x
12
2
2
1
2
0
10
25.0
1
0
1
10
2.0
1
0
1
0
2.0
1
0
10
00
11
5
4
2
1
54
5
2
11
5
4
1
2
1
54
11
对方程组调整步长的依据。
当当实用上应是误差的期望则有对提高精度（五阶方法）
例：
记方法对四阶方法的自选步长进行误差控制
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,81
6.2.5 方程组与刚性问题
（ Stiff Sets of Equation）
1,考虑 2阶常微分方程组的初值问题,
2
0
2
1
0
1
211
2
211
1
0
0
,,
,,
yy
yy
yyxf
dx
dy
yyxf
dx
dy
解此方程组的 Euler方法
2221222
1
1121111
1
,,
,,
mmmmmm
mmmmmm
hfyyyxhfyy
hfyyyxhfyy
浙江大学研究生学位课程
《实用数值计算方法,82
6.2.5
也一样，但也完全类似。理论结果方法线性多步法及是向量。其中，
：用向量形式记法，则有
K u t t aR u nge
fyfy
h
yx
dx
d
mm
mmm
,,,
0
,
1
0
fyy
yy
f
y
函数的绝对值 向量范数偏微商
y
f
2
2
1
2
2
1
1
1
y
f
y
f
y
f
y
f
浙江大学研究生学位课程
《实用数值计算方法,83
6.2.5
满足当其特征根为步长。，其中为线性多步法的特征方程则有类似方程的结果：
的特征值为又设方程组为矩阵时，则常系数线性为常数矩阵设
ji
j
j
j
hhh
njJ
xgJy
dx
dy
J ac obiJ
y
f
,,2,1
)(
可能是复数
0
00
k
i
im
ij
K
i
im
i h
kinjji,,2,1,,2,11 ；?
浙江大学研究生学位课程
《实用数值计算方法,84
6.2.5
是复平面上的区域。
它构成了绝对稳定区域这样的使线性多步法时，则
,Rh
h
j
j
绝对收敛 。
Adams内插法的绝对稳定区域,
.3
.2
.1
.2?
.3?
0,1
.2?.4?.6?.7?
2?k
4?k
3?k
5?k
图 6.11 Adams内插法的绝对稳定区域浙江大学研究生学位课程
《实用数值计算方法,85
6.2.5
.1
.1?
0
.2?
2?k
4?k
3?k1?k 5?k
6?k
图 6.12 Adams外推法的绝对稳定区域
.3?,1?
.1
.2
.3
图 6.13 R-K法的绝对稳定区域
4?k
3?k
2?k
1?k
.2?
浙江大学研究生学位课程
《实用数值计算方法,86
6.2.5
2,刚性方程组 （ Stiff Equations）
T
y
JJyy
h
J
y
f
J ac obi
hh
h
2,1,20
.120.70.0
0.50.0
07.491.0
,
'
，
例如的选取就有困难。时，
但当的特征值。矩阵取遍里属于绝对稳定区域。这使得选步长为使算法绝对稳定，应
jj m inm a x
浙江大学研究生学位课程
《实用数值计算方法,87
6.2.5
然而，由稳定性要求部分很快可以不计。的，即对应适当大时，故当其解为：
的特征值：
xy
eex
J
i
xx
21
12050
3
2
1
.0,
1.0
.50
.120
xx
x
xx
eexy
exy
eexy
1 2 050
3
50
2
501.0
1
78.24
2
120
120
CKR
CE ul e r
C
h
Ch
Chj
：阶方法：
浙江大学研究生学位课程
《实用数值计算方法,88
6.2.5
导致解的面目全非。
差积累将非绝对稳定，故舍入误则因为来取步长如果按增。只能取很小，工作量激故
,3 hCh
h
6
10
10
Rem inRem a x
Rem inRem a x
,,2,10Re
O
O
nj
gJy
dx
dy
j
j
j
j
j
j
j
j
j
且如果称为方程组
定义,
刚性方程,
一般坏条件问题严重坏条件问题刚性比,
浙江大学研究生学位课程
《实用数值计算方法,89
6.2.5
0Re
.3
,
0
h
Xxxx
y
f
j
半平面：
区域是左绝对稳定该算法的刚性方程组的解法。
也称为刚性方程组。
上满足上述条件。在的特征值对非线性方程组，
稳定：
没有限制。其数值方法应当对
A
h
图 6.14 A-稳定区域浙江大学研究生学位课程
《实用数值计算方法,90
6.2.5
一稳定：?A
证明了：
，
是绝对稳定区域该数值方法的
W idl und
h
hW
o
900
a r g
隐式方法。
定是一稳定的线性多步法必?A
图 6.15一稳定区域?A
浙江大学研究生学位课程
《实用数值计算方法,91
6.2.5
k
i
kmkimi
mmmm
mmm
fhy
G e arb
yy
h
yy
hfyy
E ul e r
St if fAa
0
''
11
11
)
2
)
方法：
重要因素但精度和计算工作量是梯形公式：
公式：例如：后退方程组的求解。
一稳定的算法均可用于浙江大学研究生学位课程
《实用数值计算方法,92
6.2.5
稳定）方法系数表（表AG e a r11.6
K?
k
6
5
4
3
2
1
0
1
2
3
4
5
6
1
1
1
1
1
1
1 -1 90
0
90
0
88
0
73
0
51
0
18
0
147360? 147
450
137
300?
147
400
137
300
25
48
1118?
25
36
137200?
147225
43?
119
2516?
13775
14772?
31
112?
253
13712?
14710
32
116
25
12
13760
14760
浙江大学研究生学位课程
《实用数值计算方法,93
6.2.5
Ni
KbhyhaxfK
KChyy
AK ut t aR ungec
G e ar
N e w t on
h
y
f
h
y
N
j
jijmimi
N
i
iimm
k
km
,,2,1
,
)
1
1
1
1
稳定方法隐式程组。方法也适用于非刚性方法。一般采用修正受到限制。步长时，要求用迭代法求解
浙江大学研究生学位课程
《实用数值计算方法,94
6.2.4
6
3
4
1
6
3
2
1
4
1
4
1
2
1
3
2
1
2
1
2
22
1
1
212
211
211
212
1
211
1
ba
hKbhKyahxfK
bhKhKyahxfK
KKhyy
KKhyhxfK
yxfK
KKhyy
hK
yhxfK
hKyy
mm
mm
mm
mm
mm
mm
mm
mm
，
，
，
二级四阶
，
，
二级二阶
，
中点法（一级二阶）
浙江大学研究生学位课程
《实用数值计算方法,95
mx 2?mx1?mx
Euler法一阶
mx 2?mx1?mx
mx 2?mx1?mx
中点法二阶
(Runge--Kutta)
mx 1?mx
四级四阶
RungeKutta
hxm 21? hxm 211
Euler 预估 --校正 法二阶
hxm 21
my
1?my
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,96
6.3 边值问题解法
000
0000
10
'
10
'
''
'''
，，
第三边值条件：
，
第二边值条件：
，
第一边值条件：
，，，
问题：二阶微分方程
byby
ayay
byay
byay
baxyyxfy
浙江大学研究生学位课程
《实用数值计算方法,97
2212
1211
21
,,2,10,,,,
,,2,10,,,,
,,2,1,
,,,,
nkyyybB
njyyyaB
Nibax
yyyxf
dx
xdy
Nk
Nj
Ni
i
一般地，
6.3.1 试射法 （ Shooting）
y
01?jB x
02?kB
图 6.16 试射法求解示意浙江大学研究生学位课程
《实用数值计算方法,98
6.3.1
的函数是显然直到满足条件。
问题。值，继续求解新的初值修改则满足要求，否则若设求得
，
，，
，则可解初值问题设
，
，，
考虑第一边值问题
11
1
1
1
1
'
'''
1
'
'''
~
uuFby
u
v
vby
uayay
yyxfy
uay
byay
yyxfy
浙江大学研究生学位课程
《实用数值计算方法,99
6.3.1
uGu
uGubyby
uGbyuGby
uayuay
uay
uFu
uu
uu
FuuF
使要找的函数是于是设解为
，
：则问题变为解初值问题对第三边值问题，设使得直到找到
：不知道的。用线性近似是但的要求
,
~
,
~
,
.
,
0
'
2
'
1
01
'
1
12
12
1
11
浙江大学研究生学位课程
《实用数值计算方法,100
6.3.2 差分方法 （ Difference Methods）
(Relaxation Methods)
0x 1x 2x
Mx
01?jB
02?kB
iy
True solution
2nd iteration
1st iteration
x
212
1101
021
,,10,,,
,,10,,,
,,,,,,
nkyyxB
njyyxB
xxxyyyxf
dx
xdy
NMk
Nj
MNi
i
图 6.17 差分方法示意浙江大学研究生学位课程
《实用数值计算方法,101
6.3.2
：个方程由边值条件提供另外个。上述方程共有未知数。共有；：要求的值在即其中，
还有别的方程可用用差分方程代替：
维向量是，
N
MN
MN
MkNjy
Mk
xyyyyy
Nfyyxf
dx
dy
jk
k
T
Nkkkk
1
,,1,0,,1
,,2,1
,,,,
,,
21
02,2 1111 kkkkkkkkk yyxxfxxyyE
NnnnyxBE
nyxBE
MMM
21
22
'
10010
0,
0,
个个浙江大学研究生学位课程
《实用数值计算方法,102
6.3.2
2
'
1
'
10
1
0
0
0
11
1
1
1 1
1
1
111
1
,,1
,,1
,,2,1;,,2,1
,0
,,
:,
njEy
y
E
njEy
y
E
MkNj
Ey
y
E
y
y
E
yyE
y
y
E
y
y
E
yyEyyyyE
yy
jM
N
i
iM
iM
jM
j
N
i
i
i
j
jk
N
i
ik
ik
jk
N
i
ik
ik
jk
N
i
N
i
ik
ik
k
ik
ik
k
kkkkkkkk
kk
另外，由边值条件得则有设每次迭代有修正值浙江大学研究生学位课程
《实用数值计算方法,103
6.3.2
很大。
很大时，矩阵的阶个特殊稀疏问题。当它是一可解下面方程组求得。的修正量时每次迭代当
NM
M
y
nnMN
jk
1
2,3,3,5
21
23
'
13
'
53
43
33
23
12
52
42
32
22
12
51
41
31
21
11
30
20
10
53
43
33
23
13
52
42
32
22
12
51
41
31
21
11
50
40
30
20
10
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
图 6.18 修正值矩阵浙江大学研究生学位课程
《实用数值计算方法,104
块对角线性代数方程组的快速求解利用矩阵的特殊结构使总运算次数达到极小，
并使矩阵系数的存贮达到极小
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
图 6.19 Gauss 消去法的目标结构
6.3.2
浙江大学研究生学位课程
《实用数值计算方法,105
D
D
D
1
0
0
1
0
0
Z
Z
Z
Z
Z
1
0
0
0
0
0
0
0
D
D
D
0
1
0
0
1
0
Z
Z
Z
Z
Z
0
1
0
0
0
0
0
0
D
D
D
0
0
1
0
0
1
Z
Z
Z
Z
Z
0
0
1
0
0
0
0
0
A
A
A
S
S
S
S
S
S
D
D
D
D
D
S
S
S
1
0
0
0
0
A
A
A
S
S
S
S
S
S
D
D
D
D
D
S
S
S
0
1
0
0
0
D
D
D
D
D
0
0
1
0
0
V
V
V
V
V
V
D
D
D
D
D
0
0
0
1
0
D
D
D
D
D
0
0
0
0
1
A
A
A
A
A
S
S
S
S
S
A
A
A
S
S
S
A
A
A
A
A
S
S
S
S
S
V
V
V
V
V
V
V
V
V
V
V
V
V
V
V
V
S
S
S
A
A
A
A
A
S
S
S
S
S
S
S
S
D--有待于对角化
A--将要改动
S--需要存贮
Z--将要变为 0
a)
b)
a)
b)
6.3.2
图 6.20 特殊矩阵结构 1
浙江大学研究生学位课程
《实用数值计算方法,106
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
1
0
0
0
0
0
1
0
0
0
0
1
0
0
0
0
0
1
0
0
Z
Z
0
0
1
0
0
0
0
0
0
0
1
0
Z
Z
0
0
0
1
0
0
0
0
0
0
0
1
Z
Z
0
0
0
0
1
0
0
S
S
S
S
S
D
D
S
S
S
S
S
1
0
S
S
S
S
S
D
D
S
S
S
S
S
0
1
V
V
V
V
V
V
V
V
V
V
V
V
V
V
S
S
S
S
S
A
A
S
S
S
S
S
S
S
a)
b)
包括四步运算
1,使用前一步的结果，简化部分元素为零
2,消去剩下子块正方结构为单位阵
3,存储以后要使用的非零系数
4,回代
6.3.2
图 6.22 特殊矩阵结构 2
浙江大学研究生学位课程
《实用数值计算方法,107
6.3.3 样条函数在两点边值问题中的应用考察 二阶线性微分方程 的边值问题
226
216
206
'
0
'
00
'''
nnn rbyby
rayay
xryxqxyxpxy
可以证明：
，个约束条件：
时，当上的：是定义在节点而
，
来近似：用三次样条函数可以的解假定微分方程边值问题的区间和节点为：设
2,1,0,,3,2,1)1(3
,,2,1
226206
1
1
23
10
jnin
xSxS
xxxni
dxcxbxaxSxS
xxS
baxxSxy
xS
xy
bxxxa
x
k
j
ik
j
i
ii
iiiii
k
n
1
0
32
00 6
1
2
1 n
k
kk xxhxxgxxfexS
bax,?
浙江大学研究生学位课程
《实用数值计算方法,108
6.3.3
此时为，而边值条件其中则有
，，并令代入微分方程把易知：
个待定常数。是而
，
其中，
226216
,,
206,,
,
2
1
3,,,,,,
,0
'''
1
0
''
1
0
2
0
'
110
iiiiii
i
n
k
kk
n
k
kk
n
k
kk
k
xrrxqqxpp
xxxSxSxS
xxhgxS
xxhxxgfxS
nhhhgfe
xx
xxxx
xx
236,,2,1,0,
6
1
2
1
1
0
32
00
1
0
1
0
2
0
nir
xxhxxgxxfeq
xxhxxgfpxxhg
i
n
k
kikiii
n
k
n
k
kikiikik
246000 rfe
浙江大学研究生学位课程
《实用数值计算方法,109
6.3.3
256
2
1
6
1
2
1
1
0
2
0
1
0
32
00
n
n
k
knknn
n
k
knknnn
rxxhxxgf
xxhxxgxxfe
3
2
3
1
3
0
2
333
000
''
110
3
2
3
1
0246
001
00,1
110
1
3
2
3
1
010
01
00
1,001
,,,,,,
33256246236
xhxhxhgxfxxS
e
r
r
rqp
y
y
xxyxy
xShhhgfe
nn
n
现设得由
，，
，
，，
有：
，，，，三等分。有节点：，把区间问题例：给定微分方程边值求得。的值，故未知数个个方程可解出共，，于是联立
浙江大学研究生学位课程
《实用数值计算方法,110
6.3.3
1,0
3
2
061224.0
3
1
061224.0
061224.05.0540814.0
~
~
061224.0061224.0061224.0
0.05.0540814.0
12
2795755
27186611655
272758118911655
02727278
3
2
6
3
1
662
3
2
3
3
1
332
33
32
210
0
01
012
012
210
''
2
2
2
1
2
0
'
x
xx
xxxxy
xyxy
hhh
egf
g
fgh
fghh
fghhh
fghhh
xhxhxhgxS
xhxhxhgxfxS
可写为的近似解故这组方程的解为条件有代入原微分方程和边界有浙江大学研究生学位课程
《实用数值计算方法,111
的实数均为大于时的捕食系统考虑系数为两种群间的交叉亲疏为本种群的亲疏系数系数为两个种群的自然增长模型例，生态系统的
0
,,,
0
,
,
,
2122
211211
21
12
21
00
2221102
1221101
xbxkx
bakkxaxkx
ba
ba
ba
ba
xxbxbbx
xxaxaax
V ol t e r r a
竞争系统，
捕食系统，
共栖系统，
时，当
00)
00)
00)
0
12
12
12
12
baiii
baii
bai
ba
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,112
为实数，
函数：阶为第一类其中，
，
，
精确解为
，积分区间
，
，
例，求解初值问题
v
kvk
x
xxJ
B e s s e lvvxJ
xJyxJy
xJyxJy
xx
xJxyxJxy
xJxyxJxy
y
x
y
dx
dy
y
x
y
dx
dy
y
x
y
dx
dy
y
dx
dy
k
k
v
v
v
0
2
34123
12101
21
13141213
11121011
43
4
32
3
21
2
2
1
1!
4
1
2
1
3,2,1,0
0.100.1,
3
2
1
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,113
相平面上的轨迹为
0201,xx
2x
1x
ak1
62k
因。而食用鱼反而减少的原期间捕鱼量减少，该模型解释了一次大战周期解的平均值平衡点
a
k
b
k 12
,
年度 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923
比例 1 1,9 2 1,4 2 2,1 2 1,2 3 6,4 2 7,3 1 6,0 1 5,9 1 4,8 1 0,7
1914年 --1923年捕获的鱼中，食用鱼所占的比例 (%)
图 6.23 相平面上的轨迹图
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,114
为扩散系数
，
，
统考虑密度分布的捕食系少。食用鱼增加，食肉鱼减此时，平衡点为
，则模型改为系数因此，我们可控制捕捞
21
2122
2
2
2
2
2
2111
2
1
2
1
1
2
1
12
2122
1211
,
,
DD
SbSSk
x
S
D
t
S
SaSSk
x
S
D
t
S
txS
txS
a
k
b
k
xbxkx
xaxkx
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,115
隐式。半隐式：可化为显式的新节点值组才能解出隐式格式：应通过方程节点值式可显式解出新显式格式：通过递推公每一个初值即可计算自动起步：由即可计算出单步法：仅由方法类型与特点：
.5
,,,,,
.4
,,,,,,
.3
,2,1
.2
.1
00111
0011
00
1
hyxyxFy
hxyxyxFy
iy
yxy
yy
iii
iiii
i
ii
公式简单，但精度不好编程简单，使用方便计算量少，但稳定性差计算量大，但稳定性好
6.3.3
总 结浙江大学研究生学位课程
《实用数值计算方法,116
.1,
2
1
0
2
1
1
111
1,11
.4
1.0
10
.3
06
,,2,4
6
1
.2
0 02.00 01.03.005.0
1001 000
.1
2''
5
'
111
21
21
2122
211211
，，，取节点为值问题用样条函数方法求解边之值。项。并计算取至初值问题：级数展开求解微分方程用
，的绝对稳定区域为算法为三阶的，并证明上式所定义的隐式单步公式试证明由
，，，
，
轨迹图，并讨论的相，作出捕食系统方法求解初值问题用四阶
yy
xyxxy
yx
y
xyyx
dx
dy
Tay lor
yxhfyxfyxfhyy
bakk
xx
xbxkx
txtxxaxkx
K ut t aR un ge
nnnnnnnn
第六章 习题浙江大学研究生学位课程
《实用数值计算方法,117
)(
1
,,,,,
)2(;//1
)(
)(
.5
211221
2
1122
2211222
1221111
O D E I N T
K ut t aR u ng e
aabakk
x
akak
xbxxakx
xxaaxkx
果。
）而导致的结条件（以上结论，并解释由于方法验证用带精度控制的四阶和初值。试选取适当的系数
，最后趋于灭绝。某一最大值后开始减少开始时增大，直至达到食肉种群量则
。初始时刻食饵十分丰富
）（
满足捕食系统
《实用数值计算方法,1
第六章 常微分方程及方程组的解法
6.1 常微分方程及其求解概述
6.2 初值问题解法
6.3 边值问题解法浙江大学研究生学位课程
《实用数值计算方法,2
6.1 常微分方程及其求解概述
6.2 初值问题解法
1.Euler方法
2.线性多步法
3.Runge--Kutta法
4.方程组及刚性问题的 Gear方法
6.3 边值问题解法
1.Shooting(试射法 )
2.差分法
Eq u a t io n salD if f e r e n t iO r d in a r yO D E s
浙江大学研究生学位课程
《实用数值计算方法,3
6.1 常微分方程及求解概述
（ Ordinary Differential Equations,ODE)
6.1.1 基本概念描述自由落体的 ODE:
的运动轨迹。
唯一确定了那么时时或设：
时若设：
tssO D E
st
st
sst
mg
dt
sd
16
36
01
00
261,00
16
1
0
'
0
2
2
浙江大学研究生学位课程
《实用数值计算方法,4
只有一个自变量的微分方程为 ODE,
否则称为偏微分方程 PDE。
方程中未知函数导数的最高阶数称为方程的阶。 （ 6-4）是二阶的
方程中关于未知函数及其各阶导数均是一次的，则称为线性微分方程。
4622 xrdxdyxqdx yd如：
和（ 6-1）都是线性二阶 ODE。
(6-2),(6-3)是 (6-1)的 初始条件 。亦称定解条件。 (6-1)(6-2)叫做 初值问题 。
6.1.1
浙江大学研究生学位课程
《实用数值计算方法,5
6.1.1
（ 6-1），（ 6-3）叫做 边值问题 。
在没有给定解条件时。方程一般有一族解曲线 y(x,c) 。如：
为任意常数。有解 ccecxy
y
dx
dy
x?
,
对任意的 n阶 ODE，如果能写成：
56,,,,1'nn yyyxfy?
则称该方程为 显式 的。方程（ 6-4）是显式的。而下面方程是 隐式 的。
xyyy '' e x ps in
浙江大学研究生学位课程
《实用数值计算方法,6
对于高阶显式方程。通过定义 n-1个新变量，可以写成 n维一阶方程组。
即令：
向量形式：一般我们把方程组写成组：可以改写为如下的方程那么
76
,,,,
56
66
1,,1,
110
'
1
1
'
2
2
'
1
1
'
0
0
1
nn
nn
i
i
yyyxfy
yy
yy
yy
xyxy
ni
dx
xdy
xy
6.1.1
浙江大学研究生学位课程
《实用数值计算方法,7
yxfy,'?
Tn
T
n
ffff
yyyy
110
110
,,,
,,,
这里：
在讨论初值问题时，我们从一阶方程开始：
00
',
yxy
yxfy
然后毫不费力地套用来解方程组。
当 f(x,y)与 y无关时,f(x,y)=g(x)
广泛的意义。故求解微分方程具有更变为积分问题则初值问题
00
00
'
,yxycdxxgy
yxy
xgy
6.1.1
浙江大学研究生学位课程
《实用数值计算方法,8
6.1.2 数值解及其重要性
难以求积。而的解是例如但仍然得不到精确解。
量的显示形式，有的解可以表示为自变故得不到精确解。不能表示为初等函数，
例如方程不定积分的组合表示。
及其的解很少能用初等函数
x
t
x
tx
dte
dteexy
xyy
yxy
O D E
0
0
'
22'
2
22
21
浙江大学研究生学位课程
《实用数值计算方法,9
6.1.3 ODE数值解的基本思想和方法特点基本思想有两点
1,离散化用 Taylor级数，数值积分和差商逼近导数等手段，把 ODE转化为离散的代数方程 (称差分方程 )。
2,递推化在具有唯一解的条件下，通过步进法逐步计算出解在一系列离散点上的值。从而得到原 ODE的数值近似解。
浙江大学研究生学位课程
《实用数值计算方法,10
6.2 初值问题解法我们讨论一阶 ODE，而高阶可能化为一阶 ODEs。一阶初值问题可以一般地写成：
86,,
00
0?
yxy
Xxxyxf
dx
dy
6.2.1 欧拉（ Euler） 方法
Euler方法是求解 (6-8)最 简单 方法，
但 精度差,故 不实用 。然而对理论分析很有用。
浙江大学研究生学位课程
《实用数值计算方法,11
6.2.1.1 方法原理及推导设初值问题 (6-8)满足：
2121
21
0
,,
,,,,,0
)
,:)
yyLyxfyxf
RyxyxL
L i ps c hi t zyfb
yXxxRfa
对常数条件：满足关于内连续。
在区域
内有界。在条件来代替：
验的可用更强的，但便于检就存在唯一的解。那么
R
y
f
ba
)),
86
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,12
为变步长。则称上述划分是任意的：
称为定步长。若其中等分，得到节点分成把
1
11
0
0
0
,,2,1
,,1,0
:,
i
iii
i
i
h
nihxx
nxXh
niihxx
xnXx
11110 xxxxx?
0y
0y
0y
6.2.1
图 6.1 常微分方程初值问题的数值解浙江大学研究生学位课程
《实用数值计算方法,13
公式。称为的近似。作为于是，我们取故
。其中级数构造：用
E ul e rxy
niyxy
yxhfyy
ni
xyxhfxy
xhyxyxy
xx
y
h
xhyxyhxy
T ay lori
i
iiii
iii
iii
ii
i
iii
96
96
1,,1,0
,
1,,1,0
,
2
)
1
00
1
'
1
1
''
2
'
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,14
欧拉方法的几何意义：
nx
Xxxx 210
0y
Xy
ny
1xy
2xy
1y
2y
h步长
6.2.1
图 6.2 Euler方法的几何意义浙江大学研究生学位课程
《实用数值计算方法,15
余项，得：
分，舍去若用梯形公式来计算积
。公式得即并使得舍去余项分，得并用左矩形公式计算积令改写成积分形式把用数值积分法构造
96
,,
,
,
,
86
)
1
E ul e r
yhxyR
Rxyxhfxyhxy
xx
dxxyxfxyhxy
ii
iii
iiiii
i
hx
x
106
1,,1,0
,,
2
00
111?
niyxy
yxfyxf
h
yy iiiiii
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,16
,改为以下预估一校正式一般地，把隐式格式隐式。半隐法：可化为显式的公式隐式：如改进的公式显式：如全部算出可计算单步法：由公式。称为改进的这个公式
106
)
106)
96)
)
)96(
1
e
E u le rd
E u le rc
yyya
E u le r
iii
1161,,1,0
,,
2
,
00
0
111
0
1
niyxy
yxfyxf
h
yy
yxhfyy
iiiiii
iiii
校正予估
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,17
6.2.1
0x 2?ix 1?ix?2?ixix1?ix nx
多步法单步法自动起步显式隐式半隐式图 6.3 ODE求解方法的类型浙江大学研究生学位课程
《实用数值计算方法,18
hvss
hvv
s
v
E ul e r
gm
v
s
v
dt
ds
mg
dt
dv
ts
ts
mg
dt
sd
mmm
mm
11
1
0
0
0
'
0
2
2
10
0
10
.10,1
100
00
10
0
方法：则：由加速度令质量初值问题：
20
0
5
10
ttvts
tvtv
真解：
m t v
m
s
m
m t v
m
s
m
s (t )
0
1
2
3
4
0.
1.
2.
3.
4.
10.
0.
-1 0,
-2 0,
-3 0,
0.
10.
10.
0.
-2 0,
0
1
2
3
4
5
6
7
8
0.
0,5
0.
-5,0
-1 0,
-1 5,
-2 0,
-2 5,
-3 0,
10.
5.
0.
-5,
-1 0,
-1 5,
-2 0,
-2 5,
-3 0,
0.
5.
7,5
7,5
5.
0.
-7,5
-1 7,5
-3 0,
0.
3,7 5
5.
3,7 5
0.
-6,2 5
-1 5,
2 6,2 5
-4 0,
0.1?h 5.0?h
表 6.1 自由落体运动方程的 Euler公式求解
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,19
t4321 t4321
40
30
20
10
0
10
40
30
20
10
0
10
v s
图 6.4 运动轨迹
6.2.1
浙江大学研究生学位课程
《实用数值计算方法,20
图 6.5 Euler公式的误差
6.2.1.2 Euler方法的误差估计一般其它方法的误差估计也类似。
这里误差是指 截断误差 (算法理论误差 )
而不是 舍入误差 。后者由计算机字长等决定，属于稳定性问题。
i) 几何分析
ii yxy?
iy
ix 1?ix
1?ixy
1?iy
xy
iiii yxhfyy,1
浙江大学研究生学位课程
《实用数值计算方法,21
6.2.1
证明作为习题）
这里：
或有方法而对改进的或则有若满足生。称为该误差由当前一步所产有是精确的，则对设定量分析
(
,,
,
2
106
,
2
,,
!2
96
)
0
'''
3
3
2
2
0
''
''
2
11
XxxMxy
hORM
h
R
E u le r
hORM
h
R
XxxMxy
Ry
h
yxy
xyy
ii
ii
ii
iiii
ii
局部截断误差，
浙江大学研究生学位课程
《实用数值计算方法,22
6.2.1
整体截断误差,设 ym是 Euler公式（ 6-9）
精确解，而 y(x) 是初值问题（ 6-8）的解。
则 整体截断误差定义为它是局部截断误差的积累。
定理,若 f(x,y)关于 y满足 Lipschitz条件。
则有估计式：
)126()1(01 LXLXm ehLRe
的步长的区间上限局部截断误差上界常数
xh
xX
R
L ips c hi t zL
mmm xyy
浙江大学研究生学位课程
《实用数值计算方法,23
6.2.1
hL
hL
RhL
hLRhL
hLRhL
RhL
Rxyxfyxfh
Rxyxhfxyhxy
yxhfyy
m
m
m
j
jm
m
mm
mmmmmmm
mmmmm
mmmm
11
1
11
111
1
,,
,
,
1
0
1
0
0
1
1
2
1
1
1
递推相减得：
证明：
浙江大学研究生学位课程
《实用数值计算方法,24
6.2.1
。方法，对改进同阶。整体截断误差与方法的说明可记为上式得：
代入，将若于是得：
注意到：
2
1
2
000
01
1
1
1
,
1
2
2
0
,2,1,0
1
11
1
hOE ul e r
h
E ul e rhO
e
L
hM
M
h
Rxyy
m
e
hL
R
e
ehLhL
Xhmx
m
m
LX
m
LXLX
m
LX
X
h
m
m
ex
x
x
11lim
浙江大学研究生学位课程
《实用数值计算方法,25
6.2.1
立是差分方程的解，若成初值问题的解是收敛性：设收敛性、稳定性讨论使用上述结果。是不知道的。因此难以的估计式中的实际上，
误差容易估计。
决定，但局部方法的精度由整体误差误差低一阶。一般来说，整体比局部说明：
m
y
O D Exy
MMxyiii
ii
i
,
.3.1.2.6
''
0m a xlim 10 mmnmh xyy
浙江大学研究生学位课程
《实用数值计算方法,26
6.2.1
况。严格定义为：的差分方程解的变化情初值的改变导致公式不能无限缩小。所以实际上对对解的影响。误差的积累可能扩大了时离散点增多。舍入当方法是收敛的，且阶为
。则称差分方程是收敛的
96
0
h
h
hOE ul e r
注意,
稳定性,
,2,1,
,0
,
,,,,
00
0
000
mvucvu
Xmhhh
vu
vuhc
mm
mm
对当满足：
它们的解使对定义：
浙江大学研究生学位课程
《实用数值计算方法,27
6.2.1
证毕。
证明：记方法是稳定的。则条件，的事实上，若方法是稳定的。则称
0
0
1
1
1
1
1
1
,,
,
,
,
LX
m
m
mmmmmm
mmmm
mmmm
mmm
e
hL
hL
vxfuxfh
vxhfvv
uxhfuu
vu
E ul e r
L ips c hi zyxf
E ul e r
浙江大学研究生学位课程
《实用数值计算方法,28
6.2.1
的区域范围。的存在稳定实际上，对给定的问题为常数。这里可推知和初值，由对给定的步长
。而不是总是适当大小，限状况。实际计算中，
一种极对以后解的影响。它是始误差时，初稳定性实际上描述了当
0
0
,
0
0
yh
c
mncxyy
xyy
h
h
h
nn
mm
绝对稳定：
浙江大学研究生学位课程
《实用数值计算方法,29
6.2.1
时可提高一阶精度。
故取
：的一次项为组合后，可使时，有而步长为法的近似解：
外推法外推技巧
xyxyxy
h
h
hchcxyxy
E ul e r
R ic har ds on
h
h
h
2
2
21
2
0
2
.4.1.2.6
2212 4121 hchcxyxy h
222 hOxyxyxy hh
浙江大学研究生学位课程
《实用数值计算方法,30
6.2.1
阶精度。它具有故可取则阶方法计算若用
1
122
1462
136
:,
1
2
1
2
1
2
p
hOxyxyxy
hOhcxyxy
hOhcxyxy
xyxyp
pph
h
p
Pp
p
p
h
pp
p
h
h
h
156212 1 2 xyxyxy hhpp
浙江大学研究生学位课程
《实用数值计算方法,31
6.2.1
12
21
146136
pp
p
p
h
h
hOhCxyxy
：，由估计。外推法用于实用的误差
。和误差估计的有效方法外推法是提高精度因此故
R ic har ds on
xyxyhC
h
h
p
p
p
p
166
12
2
2
浙江大学研究生学位课程
《实用数值计算方法,32
6.2.1
。果见表预估校正方法。比较结用的数值解求用其真解为例：
2.6
1.0,21
10
1,0,
2
E ul e r
hxy
y
x
y
x
y
dx
dy
x
i
y
i
y( x
i
)
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.0
1.0959
1.1841
1.2662
1.3434
1.4164
1.4860
1.5525
1.6165
1.6782
1.7379
1.0
1.095 4
1.18 32
1.26 49
1.34 16
1.41 42
1.48 32
1.5 492
1.61 25
1.67 33
1.73 21
表 6.2 予估校正求解结果对比浙江大学研究生学位课程
《实用数值计算方法,33
6.2.1
时的解及误差估计：技巧解此方程，
方法及外推。用其真解为例：
1
10
1,0
'
x
E ul e rey
y
xyy
x
1
1 / 2
1 /4
1 /8
1 /1 6
1 /3 2
,0 0 0 0 0
,2 5 0 0 0
,3 1 6 4 1
,3 4 3 6 1
,3 5 6 0 7
,3 6 2 0 6
,5 0 0 0 0
,3 8 2 8 1
,3 7 0 8 1
,3 6 8 5 4
,3 6 8 0 4
,3 6 7 9 2
-,3 6 7 8 8
-,1 1 7 8 8
-,0 5 1 4 1
-,0 2 4 2 7
-,0 1 1 8 1
-,0 0 5 8 2
-,5 0 0 0
-,1 3 2 8 1
-,0 5 4 4 1
-,0 2 4 9 3
-,0 1 1 9 6
-,0 0 5 8 6
,1 3 2 1 2
,0 1 4 9 3
,0 0 2 9 3
,0 0 0 6 6
,0 0 0 1 6
,0 0 0 0 4
h hy hh yy?22 hy? 估计? 外推?
表 6.3 Euler法与外推结果的比较浙江大学研究生学位课程
《实用数值计算方法,34
6.2.2 线性多步法
1
1
11
,1
1
1
1
00
10
:,
,
,,
1.2.2.6
,,1,
,,,
m
m
m
m
m
m
m
m
x
x
kmmm
mkmkm
x
x
mm
x
x
x
x
m
ii
m
dxxLyy
yxf
xxxL agr ange
dxyxfxyxy
dxyxfdyyxf
dx
dy
y
mixyyxyy
xxx
的插值多项式近似上，，，在用插值法欲求处已求得：假定在
外插法：Ad a m s
外插浙江大学研究生学位课程
《实用数值计算方法,35
6.2.2
k
i
imkimm
imimim
im
kmmm
imk
k
i imimk
imk
imkm
fbhyy
A dam s
kiyxff
ki
xx
xxxxxx
xD
xD
xD
fxL
0
1
1
,
0,
,
,
,,1,0,
,,1,0
外插值公式：由此得这里：
(?)
方法。时，线性多步法就是特别，当
。其值见表其中
E u le rk
dx
xD
xD
h
b
m
m
x
x
imimk
imk
ki
0
4.6
1 1
,
,
(?)
浙江大学研究生学位课程
《实用数值计算方法,36
6.2.2
i 0 1 2 3 4 5
b
0 i
2 b
1 I
12 b
2 I
24 b
3 I
720 b
4 I
1440 b
5 i
1
3
23
55
1901
4277
-1
-1 6
-5 9
-2 7 7 4
-7 9 2 3
5
37
2616
9982
-9
-1 2 7 4
-7 2 9 8
251
2877 -4 7 5
211 51623
12
2
mmmmm fff
hyy
k 时，线性多步法为例如：
Adams 外插法 (k=2) 3阶 3步 显式
2?m
1?m
m
表 6.4 外插系数 bki值图 6.6 3阶 3步外插法
浙江大学研究生学位课程
《实用数值计算方法,37
6.2.2
，舍入误差也小。内插法较外插法更精确最常用，其插值公式为法。时，即得到改进当
。它是隐式或半隐式格式
：多项式的取插值节点为内插法：
1
1,1
,,,,
,
p
E ul e rkp
xLL agr ang e
xxx
A dam s
p
km
pmmpkm
1
,1
m
m
x
x
p
kmmm dxxLyy
。见表其中 5.6
0
11
i
i
k
k
i
imkmm
b
fbhyy
浙江大学研究生学位课程
《实用数值计算方法,38
6.2.2
i 0 1 2 3 4 5
b 0i
*
2 b 1i
*
1 2 b 2i
*
2 4 b 3i
*
7 2 0 b 4i
*
1 4 4 0 b 5i
*
1
1
5
9
251
475
1
8
19
646
1427
-1
-5
-2 6 4
-7 9 8
1
106
482
-1 9
-1 7 3 27
值内插系数表?kib5.6
111 85
12
12
mmmmm fff
hyy
pk 时，线性多步法为，例如：
Adams 外插法 (k=2) 4阶 3步
1?m1?m m
图 6.7 4步 3阶 Adams内插公式
浙江大学研究生学位课程
《实用数值计算方法,39
6.2.2
外插法：不稳定。
稳定内插法：舍入误差小，
稳定性：
单。外插法：显式，计算简
，复杂内插法：隐式或半隐式迭代格式：
阶方法外插法：
阶方法内插法：
局部截断误差：
)
)
1
2
)
2
3
iii
ii
khO
khO
i
k
k
1?p
浙江大学研究生学位课程
《实用数值计算方法,40
6.2.2
k
i
ii
k
ii
ihxyhihxyhxy
k
k
0
'
00
0,
,0
,
2.2.2.6
；
定义算符：
步法不同时为步法线性多步法是常数一般的插值型公式待定系数法
k
i
k
i
imiimi fhy
0 0
1?k?取浙江大学研究生学位课程
《实用数值计算方法,41
6.2.2
0x
1?mx mx 1?mx
kmx?
my
nk
k
nk
k
k
f
f
h
y
y
y
0
0
0
1
0
0
0
0
mmm yxff
hkN
,
,1
个方程。等步长图 6.8 一般化插值形式浙江大学研究生学位课程
《实用数值计算方法,42
6.2.2
k
rr
k
rr
r
k
i
ii
k
rr
r
k
r
k
r
C
iC
C
xyhC
xhyCxyChxy
xy
ih
xihyxyihxy
xy
ih
xihyxyihxy
T ay lor
1
2
1
1
21
0
1
100
'
10
'''
2
''''
''
2
'
2
!1
1
2
!
1
!2
!2
其中：；
代入得：
展开右端：
浙江大学研究生学位课程
《实用数值计算方法,43
6.2.2
步方法。
阶阶及外插及内插法为隐式显式由此得于是：
充分大，使次连续可微，则可取有因此，若
k
kkA dam s
xyhCfhy
CCCk
rxy
k
k
k
i
rr
rihxii
r
1
0
0
.0
1
0
11
1
10
k
i
k
i
imiimi fhy
kr
0 0
步方法：阶线性浙江大学研究生学位课程
《实用数值计算方法,44
6.2.2
ak 02,1,2
2
记步法：例，考虑
04
!2
1
8
!3
1
024
!2
1
02
01
2113
2112
21011
10
C
C
C
aC
解之得,
a
a
a
a
5
12
1
1
3
2
,51
12
1
,1
2
1
0
1
浙江大学研究生学位课程
《实用数值计算方法,45
6.2.2
方法型：它是以下阶方法时，
阶的方法是时，
故同时，可以算出：
故一般二步法为：
eMiSi m ps on
CCa
Ca
ln
40,01
3,01
54
4
mmm
mmm
fafafah
ayyay
51185
12
1
12
12
aC
aC
1317
!53
1
1
!4
1
5
4
mmmmm fffhyy 122 43
浙江大学研究生学位课程
《实用数值计算方法,46
6.2.2
1,1
01
0,2
11
1
3.2.2.6
11
1
'
即特征根的特征方程：
步长方法：例对典型方程：
绝对稳定性
h
h
hh
h
yh
yhy
hfyy
hE u le r
yy
m
mm
mmm
来讨论绝对稳定性）
程一般只限于这个典型方(
浙江大学研究生学位课程
《实用数值计算方法,47
6.2.2
稳定可以任意取，步长时，恒有当其特征根为隐式方法例，后退
A
h
h
y
h
y
yhy
hyyy
E ul e r
mm
mm
mmm
10Re
1
1
1
1
1
1
1
1
'
11
图 6.9
浙江大学研究生学位课程
《实用数值计算方法,48
6.2.2
线性多步法的绝对稳定性：
次多项式的其特征方程为多步法为典型方程则常数设
k
hh
yhy
yy
y
f
k
i
imi
k
i
imi
00
'
k
i
im
i
k
i
im
i h
00
0
浙江大学研究生学位课程
《实用数值计算方法,49
6.2.2
定义,
，称为，取值范围绝对稳定的则称线性多步法关于此
，按模小于，若特征方程的根给定
h
h
h i 1
绝对稳定 。
绝对稳定区域。
.0,,0,
01
01
1
0
1
0
1
EI
k
i
i
k
k
k
i
i
k
k
i
i
bh
bh
A d am s
其稳定区域为程：内插及外插法的特征方
浙江大学研究生学位课程
《实用数值计算方法,50
6.2.2
的取值表 EI,5.6
k 1 2 3 4
E
2? 1
I
6? 3
11
6?
49
90?
10
3?
5
3
1
2
5
1
2
22
0
3
1
1
3
4
3
1
14
3
1
hOe
hOe
hh
h
h
h
h
它的根：
方法的特征方程
Milne
浙江大学研究生学位课程
《实用数值计算方法,51
6.2.2
：见表及误差，得到的解对不同的方法计算例，用对稳定。
方法不绝，因此，对任何
6.6
10
10
ln
ln
'
yh
y
yy
eMi
eMih
h u 误 差
0,1
0,0 1
0,0 0 1
0,0 0 0 1
0,0 0 0 0 1
,4 8 6 E - 4
,4 8 2 E - 4
,3 8 6 E - 4
,3 2 9 E - 4
,5 7 8 E - 4
-,3 2 3 E - 5
-,2 8 1 E - 5
,6 7 6 E - 5
-,3 7 5 E - 4
-,1 2 4 E - 4
表 6.6 计算结果浙江大学研究生学位课程
《实用数值计算方法,52
6.2.2
6.2.2.4 预估 --校正方法
（ Predictor--Corrector Method）
迭代法。这是个非线性方程已知。，，，，，其中阶隐式方法：但计算困难。现考虑其精度高，稳定性好。尤隐式方法：
110
0
kify
k
y
f
imim
176
,
1
0
1
0
k
i
imi
k
i
kmkmkimikm
fh
yxfhyy
，，
，
10
186
1
0
1
0
1
n
fh
yxfhyy
k
i
imi
k
i
n
kmkimi
n
km km
浙江大学研究生学位课程
《实用数值计算方法,53
6.2.2
。此时应减小步长
。太多的校正不合适。可满足二，三步校正。，一般只要对给定算法。通称为算式。称为校正式，即算式。称为预估式，即构成的算法，时，
选为显式算法而
h
yy
PC
CC or r e c t or
Pe di c t or
y
n
km
n
km
km
1
0
0
196
Pr186
196186
1961
0
1
0
0
k
i
k
i
imiimikm fhyy
浙江大学研究生学位课程
《实用数值计算方法,54
6.2.2
注意,一步校正的计算量? 预估计算量。
所以要适当选取 h才能发挥 PC的优点。
设绝对稳定区域：
达到精度的校正次数为 N，则 h的选取，
应满足：
否则可用 N步显式算法稳定达到目的 。
0,
0,
c
p
SC
SP
算式：
算式：
p
cp
SNh
ShS
1.2
.1
hmx 1?mx
pS
cS
2?N
浙江大学研究生学位课程
《实用数值计算方法,55
6.2.2
Adams四步四阶 预估校正算法
123434
12334
937599
24
:54
9375955
24
:44
mmmmmm
mmmmmm
ffff
h
yyC
ffff
h
yyP
阶步阶步
0
,,
270
19
0
4
1
4
0
44
1
44
0
4
1
4
55
5
，
稳定性判别可用很大的麻烦但变步长会带来可用来控制步长方法误差：
C
mm
mmmm
mm
S
y
f
hh
yy
yxfyxf
y
f
h
yyyhC
浙江大学研究生学位课程
《实用数值计算方法,56
6.2.3 Runge--Kutta 方法
，考虑为构造高精度的单步法
。后者比前者精度高一阶型预估一校正法：
方法：
12
1
211
1
11
,
,
2
1
2
1
,
KyhxhfK
yxhfK
KKyy
E ul e r
yxhfK
Kyy
E ul e r
mm
mm
mm
mm
mm
12
21
,mm
mm
yhxhfk
kyy
E u le r 法后退浙江大学研究生学位课程
《实用数值计算方法,57
6.2.3
方法的级正整数，称为 KRN
Niba
KbhyhaxfK
yxfK
KChyx
hyxhyy
i
l
ili
i
l
lilmimi
mm
N
i
iimm
mmmm
,,3,2
,
,
,,
,,
1
1
1
1
1
1
1
上的线性组合：在某些点是而 1,?mm xxf?
其中：
阶次项一致），，，即项完全一致。级数。取的展开，并比较进行把
PP
PT ay lor
xyT ay lorK mi
10(
1
1
浙江大学研究生学位课程
《实用数值计算方法,58
6.2.3
时，有当公式的阶数。称为这里的就可以得到系数
3?
N
KRP1,,1;,,1, ilNibC ili
232132333
1222
1
3
1
1
,
,
,
KhbKbahyhaxfK
KhayhaxfK
fyxfK
KChyy
mm
mm
mmm
i
iimm
浙江大学研究生学位课程
《实用数值计算方法,59
6.2.3
的系数得：比较而其中：
展开：在将
3,2,1,0
6
1
2
2
2
2
,,
43
2
2
''
32
33
2
223223
2
3
1
3322321
32
ih
hOGFfh
F
h
hfxyhxy
fffffG
xyfffF
hOGaCaCFfbaC
h
FaCaChfCCCKC
T ay loryxKK
i
y
mm
yyxyxx
yx
y
i
ii
mm
浙江大学研究生学位课程
《实用数值计算方法,60
6.2.3
6
1
3
1
2
1
1
3223
2
33
2
22
3322
321
baC
aCaC
aCaC
CCC
六个未知数，二个自由,故可取
21
2
1
6
1
3
2
6
1
3232
321
baa
CCC
，，
，，
213
12
1
3211
2
2
1
2
1
,
4
6
1
hKhKyhxfK
hKyhxfK
yxfK
KKKhyy
mm
mm
mm
mm
，
，
故：
浙江大学研究生学位课程
《实用数值计算方法,61
6.2.3
N=4,四级四阶 R--K方法
34
23
12
1
43211
,
2
1
,
2
1
2
1
,
2
1
,
22
6
hKyhxfK
hKyhxfK
hKyhxfK
yxfK
KKKK
h
yy
mm
mm
mm
mm
mm
----最常用的古典 Runge-Kutta方法。
增加计算函数值的次数（级）与提高精度（阶）的关系见下表,级 N 2 3 4 5 6 7 n? 8
阶 P 2 3 4 4 5 6 n? 2
表 6.7 Runge-Kutta方法中 级 与 阶 的关系浙江大学研究生学位课程
《实用数值计算方法,62
6.2.3
各种解法精度比较：
例高精度方法的主要优缺点
10
22
).1(
2'
y
yxxy
KR
X
m
E u l er 改进 E u l er R— K ( 4 阶 ) 精确解
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.0 000
1.2 000
1.4 420
1,7384
2,1041
2,5569
3,1 183
3.8 139
4.6 747
5.7 376
7.0 472
1.0 000
1.2 2 10
1.4 9 23
1.8 284
2.2 466
2,7680
3,4174
4,2257
5,2288
6,4704
8,0032
1.0 000
1.2 221
1.4 977
1.8 432
2.2 783
2.8 274
3.5 20 1
4.3 92 7
5.4 89 4
6.8 64 3
8.5 83 4
1.0 000
1.2 221
1.4 977
1.8 432
2.2 783
2.8 274
3.5 202
4.3 928
5.4 895
6.8 645
8.5 836
表 6.8 各种解法在例题中的结果比较浙江大学研究生学位课程
《实用数值计算方法,63
6.2.3
（ 2）,单步法，自动起步
（ 3）,易改为变步长
（ 4）,绝对稳定区域较同阶线性多步法大
（ 5）,计算工作量较大，有时大于隐式方法
（ 6）,估计误差不易绝对稳定性讨论
mm
mm
y
hh
hy
f
h
f
h
hfyy
yyN
3
3
2
2
''
3
'
2
1
'
!3!2
!3!2
3
此时：
为例。取典型方程以浙江大学研究生学位课程
《实用数值计算方法,64
6.2.3
D
KR
h
hhh
hhh
y
y
hh
hy
mm
im
im
mm
区域方法的绝对稳定一般地，各级时，，当其根为：
故特征方程：
1051.2
6
1
2
1
1
0
6
1
2
1
1
!3!2
1
1
32
1
32
1
32
1
级? 1 D
1
2
3
4
(-2,0)
(-2,0)
(-2.51,0)
(-2.78,0)
432
32
2
24
1
6
1
2
1
1
6
1
2
1
1
2
1
1
1
hhhh
hhh
hh
h
表 6.9 各级 R-K方法的绝对稳定区域浙江大学研究生学位课程
《实用数值计算方法,65
6.2.3
0,78.24,2.0
0,78.222.0
10
20
2
2
1
1
'
hh
hh
y
yy
KR
h
hKR
，
取方法计算用四级经典例绝对稳定区域。
使方法时，要取在使用
x h
1
时的误差 h
2
时的误差
0.0
0.2
0.4
0.6
0.8
1.0
0
-0.092795
-0.012010
-0.001366
-0.000152
-0.000017
0
4.98
25.0
125.0
625.0
3125.0
表 6.10 不同步长对精度的影响浙江大学研究生学位课程
《实用数值计算方法,66
6.2.3
。可能导致舍入误差积累的计算量，步长太小：增加不必要非绝对稳定。步长太大：精度不够，
实际步长选取：
外推技巧。需用选取误差估计和实际步长的属于绝对稳定区域。
近似确定步长。使可以用常数的问题对于
R ic har ds on
xhxh
x
x
y
f
mm
浙江大学研究生学位课程
《实用数值计算方法,67
6.2.3
是适当的。，则认为，若对给定误差估计方法：对四阶方法自选步长的
2
12
2
2
1
1
2
2
11
4
4
4
54
4
2
11
5
4
1
2
1
54
11
h
yych
hOchyy
hO
h
cxyy
hOchxyy
KR
KR
h
m
h
m
h
m
h
m
m
h
m
m
h
m
浙江大学研究生学位课程
《实用数值计算方法,68
P阶图 6.10 变步长 Runge--Kutta方法框图
00,0,,,xyym PhX给定
211 211m a x,hmhm hmhm yy yy 计算
hmm
mm
yy
hxx
11
1
hh?2
Xxm1
121 P
,2,1
,
m
yx mm输出
hh?2
1mm
Y
Y
N
N
Y
N
hh
25.0
9.0?
hh?4
hh 2.09.0
5
49.0
6.2.3
浙江大学研究生学位课程
《实用数值计算方法,69
展开：处在 T a y lo ryxKK mm,,32
xym
yy
m
xx
ymxm
m
m
m
m
m
m
m
m
m
fKbfbaah
f
Kbfbah
f
ha
fKbfbahfhafK
yx
f
fha
y
ffha
x
fha
y
f
fha
x
f
hafK
2323233
2
2
232323
22
3
23232333
2
2
22
22
2
2
22
2
222
!2!2
!2!2
1
1
1
1
1
1
1
1
1
!
,,3,2
,
,
P
P
m
P
mm
i
j
iij
i
j
jijmimi
mmm
N
i
iimm
hO
P
h
xyxyxy
ab
Ni
KbhyhaxfK
fyxfK
KChyy
6.2.3
浙江大学研究生学位课程
《实用数值计算方法,70
1.0
00
0
!
0
!2
0
!1
0
)
)
,,,
,,,
2'
1
1
''
2
'
110
11
h
y
yxy
T a y l or
R
Ry
s
x
y
x
y
x
yxy
T a y l orxyb
KR
a
yyy
kyyyy
s
s
s
s
k
kmkmmm
前三个点的值级数求初值问题例：用不影响所需精度为止。直至展开式：的用精确解计算所需点上的值。
方法，或用小步长高精度用单步法求出方法取得。
的值可以通过两种出发值步方法
6.2.4 出发值的计算浙江大学研究生学位课程
《实用数值计算方法,71
00499 950 00.0000500 050 00.0
160
1
20
1
2
1
0000000
0
721352
615102
5102
432
32
2
21
,,
101
852
764''''
76'5''4'''8
65'4''2'''7
54''''''6
4''''2''5
''''''4
''2''''
'''
2'
2
yyy
xxxxy
yyyyyy
x
yyyyyyyyxy
yyyyyyyxy
yyyyyyxy
yyyyyxy
yyyyxy
yyyxy
yyxy
yxxy
yxyxf
，，
于是，
时，故所以解：由于
25206010 85'' yyy,，
得，，，取前三点 1.00 101 hhxxhx
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,72
质量控制 Runge--Kutta 步进程序
SUBROUTINE RKQC(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)
PARAMETER (NMAX=10,PGROW=-0.20,PSHRINK=-0.25,FCOR=1./15.
ONE=1.0,SAFETY=0.9,ERRCON=6.E-4
EXTERNAL DERIVS
DIMENSION Y(N),DYSX(N),YSCAL(N),YTEMP(NMAX),YSAV(NMAX),
DYSAV(NMAX)
XSAV=X 保留初值
DO 11 I=1,N
YSAV(I)=Y(I)
DYSAV(I)=DYDX(I)
H=HTRY
HH=0.5?H
CALL RK4 (YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)
X=XSAV+HH
CALL DERIVS (X,YTEMP,DYDX)
CALL RK4 (YTMP,DYDX,N,X,HH,Y,DERIVS)
X=XSAV+H
IF (X,EQ,XSAV) PAUS E ‘步长无意义’
CALL RK4 (YSAV,DYSAV,N,XSAV,HH,YTEMP,DERIVS)
ERRMAX=0.
DO 12 I=1,N
YTEMP(I)=Y(I)-YTEMP(I)
ERRMAX = MAX(ERRMAX,ABS(YTEMP(I) / YSCAL(I)))
ERRMAX = ERRMAX / EPS
误差计算一个单步计算两个半步长计算置初始值
11
1
12
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,73
PARAMETER (NVAR=4)
DIMENSION VS(NVAR)
COMMON /PATH/ KMAX,KOUNT,DXSAV,X(200),Y(10,200)
EXTERNAL DERIVS RKQC
X1=1.0
X2=10.0
VS(1)=BESSO(X1)
VS(2)=BESSI(X1)
VS(3)=BESSJ(2,X1)
VS(4)=BESSJ(3,X1)
EPS=1.0E-4
HI=1.0
HMIN=0.0
KMAX=100
DXSAV=(X1-X2)/20.0
CALL ODEINT (VS,NVAR,X1,X2,EPS,HI,HMIN,NOK
NBAD,NOK,NBAD,DERIVS,RKQC)
WRITE (?,?)
END
SUBROUTINE DERIVS (X,Y,DYDX)
DIMENSION Y(1),DYDX(1)
DYDX(1) = -Y(2)
DYDX(2) = Y(1)-(1.0/X)?Y(2)
DYDX(3) = Y(2)-(2.0/X)?Y(3)
DYDX(4) = Y(3)-(3.0/X)?Y(4)
RETURN
END
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,74
IF (ERRMAX.GT.ONE) THEN 不满足要求
H = SAFETY?H?(ERRMAXPSHRINK) 估计下次步长（缩小）
GOTO 1
ELSE 满足要求
HDID = H
IF (ERRMAX,GT,ERRCON) THEN
HNEXT = SAFETY?H? (ERRMAXPGROW) 估计可放大的步长
ELSE
HNEXT =4,?H 最多可放大四倍
ENDIF
ENDIF
DO 13 I=1,N 完成五阶截断误差
Y(I)=Y(I) + YTEMP(I)?FCOR
CONTINUE
RETURN
END
四阶 Runge--Kutta 步进程序（续）
13
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,75
计算结果
NOK=30
NBAD=0
KOUNT=15
X I n t eg ral B E SS J ( 3,x )
1,0 0 0 0
1,6 3 2 9
2,1 0 3 6
2,6 8 1 3
3,3 7 6 9
4,1 0 3 9
4,7 9 1 0
5,5 1 4 2
6,1 4 2 5
6,8 4 1 3
7,4 9 5 9
8,2 1 1 2
8,9 1 7 4
9,5 8 5 5
1 0,0 0 0 0
0,0 1 9 5 6 3
0,0 7 6 5 6 4
0,1 4 5 8 8 2
0,2 5 0 5 4 4
0,3 7 0 1 1 5
0,4 3 3 3 9 8
0,3 9 6 3 8 4
0,2 5 2 4 6 0
0,0 7 1 7 1 0
-0,1 2 9 1 0 4
-0,2 5 7 5 3 6
-0,2 8 6 3 9 3
-0,1 9 7 3 8 9
-0,0 4 3 9 6 8
0,0 5 8 3 8 0
0,1 4 5 8 8 1
0,4 3 3 3 9 7
0,3 9 6 3 8 3
-0,1 2 9 1 0 3
-0,2 5 7 5 3 5
-0,2 8 6 3 9 1
-0,1 9 7 3 8 8
0,0 5 8 3 7 9
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,76
四阶 Runge--Kutta 方法：
SUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS)
PARAMETER (NMAX=10)
DIMENSION Y(N),DYDX(N),YOUT(N),YT(NMAX)
DYT(NMAX),DYM(NMAX)
HH=H?0.5
H6=H/6.0
XH=X+HH
DO 11 I=1,N
YT(I)=Y(I)+HH?DYDX(I) 第一步
CALL DERIVS(XH,YT,DYT)
DO 12 I=1,N 第二步
YT(I)=Y(I)+HH?DYT(I)
CALL DERIVS(XH,YT,DYM)
DO 13 I=1,N 第三步
YT(I)=Y(I)+H? DYM(I)
DYM(I)=DYT(I) + DYM(I)
CALL DERIVS(X+H,YT,DYT)
DO 14 I=1,N
YOUT(I)=DYDX(I)+DYT(I)+2,?DYM(I)
YOUT(I)=Y(I) + H6?YOUT(I)
CONTINUE
RETURN
END
14
13
12
11
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,77
用四阶 RUNGE--KUTTA公式，等步长计算 NSTEP步
SUBROUTINE RKDUMB(VSTART,NVAR,X1,X2,NSTEP,DERIVS)
PARAMETER (NMAX=10) 函数最大个数
COMMON /PATH/ XX(200),Y(10,200)
DIMENSION VSTART(NVAR),V(NMAX),DV(NMAX)
置初值
DO 13 K=1,NSTEP
CALL DERIVS(X,V,DV)
CALL RK4(V,DV，,NVAR,X,H,V,DERIVS)
X=X+H
XX(K+1)=X 存贮中间步骤
DO 12 I=1,NVAR
Y(I,K+1)=V(I)
CONTINUE
CONTINUE
结束返回
12
13
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,78
带自适应步长控制的四阶 Runge--Kutta 驱动程序
SUBROUTINE ODEINT(YSTART,NVAR,X1,X2,EPS,H1
HMIN,NOK,NBAD,DERIVS,,RKQC)
PARAMETER (MAXSTEP=10000,NMAX=10,TWO=2.0,
ZERO=0.0,TINY=1.0E-30)
COMMON /PATH/ KMAX,KOUNT,DXSAV,XP(200),YP(10,200)
DIMENSION YSTART(NVAR),YSCAL(NMAX),Y(NMAX),DYDX(NMAX)
X=X1
H=SIGH(H1,X2-X1)
NOK=0
NBAD=0
KOUNT=0
DO 11 I=1,NVAR
Y(I)=YSTART(I)
XSAV=X-DXSAV?TWO ----保证第一步存贮
DO 16 NSTEP=1,MAXSTEP ----最多取 MAXSTEP步数
CALL DERIVS(X,Y,DYDX)
DO 12 I=1,NVAR 改变比例因子控制准确性
YSCAL(I)=ABS(Y(I))+ABS(H?DYDX(I)) + TINY
IF(KMAX.GT.O.AND.ABS(X-XSAV).GT.ABS(DXSAV)
KOUNT,LT,KMAX-1) THEN
KOUNT=KOUNT+1
XP(KOUNT)=X
DO 13 I=1,NVAR
YP(I,KOUNT)=Y(I)
XSAV=X
ENDIF
存贮中间结果
13
12
11
置初值
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,79
IF ((X+H-X2)?(X+H-X1),GT,ZERO) H=X2-X1 越界处理
CALL RKQC(Y,DYDX,NVAR,X,H,EPS,YSCAL,HDID,HNEXT,DERIVS)
IF (HDID,EQ,H) THEN
NOK=NOK+1 好步长计算
ELSE
NBAD=NBAD+1 坏步长计算
ENDIF
IF ((X-X2)?(X2-X1),GE,ZERO) THEN 进行完毕
DO 14 I=1,NVAR
YSTART(I)=Y(I)
IF (KMAX,NE,O) THEN 保留最终步结果
KOUNT=KOUNT+1
XP(KOUNT)=X
DO 15 I=1,NVAR
YP(I,KOUNT)=Y(I)
ENDIF
RERURN
ENDIF
IF (ABS(HNEXT),LT,HMIN) PAUSE '步长太小 '
H=HNEXT
CONTINUE
PAUSE '步数太多，越界 '
RETURN
END
16
15
14
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,80
i
i
h
mm
h
m
h
m
m
h
m
m
h
m
y s c a l e
sh
sh
hhh
h
h
hOyxy
hOchyy
hO
h
cxyy
hOchxyy
K ut t aR unge
k ut t aR unge
m a x
12
2
2
1
2
0
10
25.0
1
0
1
10
2.0
1
0
1
0
2.0
1
0
10
00
11
5
4
2
1
54
5
2
11
5
4
1
2
1
54
11
对方程组调整步长的依据。
当当实用上应是误差的期望则有对提高精度（五阶方法）
例：
记方法对四阶方法的自选步长进行误差控制
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,81
6.2.5 方程组与刚性问题
（ Stiff Sets of Equation）
1,考虑 2阶常微分方程组的初值问题,
2
0
2
1
0
1
211
2
211
1
0
0
,,
,,
yy
yy
yyxf
dx
dy
yyxf
dx
dy
解此方程组的 Euler方法
2221222
1
1121111
1
,,
,,
mmmmmm
mmmmmm
hfyyyxhfyy
hfyyyxhfyy
浙江大学研究生学位课程
《实用数值计算方法,82
6.2.5
也一样，但也完全类似。理论结果方法线性多步法及是向量。其中，
：用向量形式记法，则有
K u t t aR u nge
fyfy
h
yx
dx
d
mm
mmm
,,,
0
,
1
0
fyy
yy
f
y
函数的绝对值 向量范数偏微商
y
f
2
2
1
2
2
1
1
1
y
f
y
f
y
f
y
f
浙江大学研究生学位课程
《实用数值计算方法,83
6.2.5
满足当其特征根为步长。，其中为线性多步法的特征方程则有类似方程的结果：
的特征值为又设方程组为矩阵时，则常系数线性为常数矩阵设
ji
j
j
j
hhh
njJ
xgJy
dx
dy
J ac obiJ
y
f
,,2,1
)(
可能是复数
0
00
k
i
im
ij
K
i
im
i h
kinjji,,2,1,,2,11 ；?
浙江大学研究生学位课程
《实用数值计算方法,84
6.2.5
是复平面上的区域。
它构成了绝对稳定区域这样的使线性多步法时，则
,Rh
h
j
j
绝对收敛 。
Adams内插法的绝对稳定区域,
.3
.2
.1
.2?
.3?
0,1
.2?.4?.6?.7?
2?k
4?k
3?k
5?k
图 6.11 Adams内插法的绝对稳定区域浙江大学研究生学位课程
《实用数值计算方法,85
6.2.5
.1
.1?
0
.2?
2?k
4?k
3?k1?k 5?k
6?k
图 6.12 Adams外推法的绝对稳定区域
.3?,1?
.1
.2
.3
图 6.13 R-K法的绝对稳定区域
4?k
3?k
2?k
1?k
.2?
浙江大学研究生学位课程
《实用数值计算方法,86
6.2.5
2,刚性方程组 （ Stiff Equations）
T
y
JJyy
h
J
y
f
J ac obi
hh
h
2,1,20
.120.70.0
0.50.0
07.491.0
,
'
，
例如的选取就有困难。时，
但当的特征值。矩阵取遍里属于绝对稳定区域。这使得选步长为使算法绝对稳定，应
jj m inm a x
浙江大学研究生学位课程
《实用数值计算方法,87
6.2.5
然而，由稳定性要求部分很快可以不计。的，即对应适当大时，故当其解为：
的特征值：
xy
eex
J
i
xx
21
12050
3
2
1
.0,
1.0
.50
.120
xx
x
xx
eexy
exy
eexy
1 2 050
3
50
2
501.0
1
78.24
2
120
120
CKR
CE ul e r
C
h
Ch
Chj
：阶方法：
浙江大学研究生学位课程
《实用数值计算方法,88
6.2.5
导致解的面目全非。
差积累将非绝对稳定，故舍入误则因为来取步长如果按增。只能取很小，工作量激故
,3 hCh
h
6
10
10
Rem inRem a x
Rem inRem a x
,,2,10Re
O
O
nj
gJy
dx
dy
j
j
j
j
j
j
j
j
j
且如果称为方程组
定义,
刚性方程,
一般坏条件问题严重坏条件问题刚性比,
浙江大学研究生学位课程
《实用数值计算方法,89
6.2.5
0Re
.3
,
0
h
Xxxx
y
f
j
半平面：
区域是左绝对稳定该算法的刚性方程组的解法。
也称为刚性方程组。
上满足上述条件。在的特征值对非线性方程组，
稳定：
没有限制。其数值方法应当对
A
h
图 6.14 A-稳定区域浙江大学研究生学位课程
《实用数值计算方法,90
6.2.5
一稳定：?A
证明了：
，
是绝对稳定区域该数值方法的
W idl und
h
hW
o
900
a r g
隐式方法。
定是一稳定的线性多步法必?A
图 6.15一稳定区域?A
浙江大学研究生学位课程
《实用数值计算方法,91
6.2.5
k
i
kmkimi
mmmm
mmm
fhy
G e arb
yy
h
yy
hfyy
E ul e r
St if fAa
0
''
11
11
)
2
)
方法：
重要因素但精度和计算工作量是梯形公式：
公式：例如：后退方程组的求解。
一稳定的算法均可用于浙江大学研究生学位课程
《实用数值计算方法,92
6.2.5
稳定）方法系数表（表AG e a r11.6
K?
k
6
5
4
3
2
1
0
1
2
3
4
5
6
1
1
1
1
1
1
1 -1 90
0
90
0
88
0
73
0
51
0
18
0
147360? 147
450
137
300?
147
400
137
300
25
48
1118?
25
36
137200?
147225
43?
119
2516?
13775
14772?
31
112?
253
13712?
14710
32
116
25
12
13760
14760
浙江大学研究生学位课程
《实用数值计算方法,93
6.2.5
Ni
KbhyhaxfK
KChyy
AK ut t aR ungec
G e ar
N e w t on
h
y
f
h
y
N
j
jijmimi
N
i
iimm
k
km
,,2,1
,
)
1
1
1
1
稳定方法隐式程组。方法也适用于非刚性方法。一般采用修正受到限制。步长时，要求用迭代法求解
浙江大学研究生学位课程
《实用数值计算方法,94
6.2.4
6
3
4
1
6
3
2
1
4
1
4
1
2
1
3
2
1
2
1
2
22
1
1
212
211
211
212
1
211
1
ba
hKbhKyahxfK
bhKhKyahxfK
KKhyy
KKhyhxfK
yxfK
KKhyy
hK
yhxfK
hKyy
mm
mm
mm
mm
mm
mm
mm
mm
，
，
，
二级四阶
，
，
二级二阶
，
中点法（一级二阶）
浙江大学研究生学位课程
《实用数值计算方法,95
mx 2?mx1?mx
Euler法一阶
mx 2?mx1?mx
mx 2?mx1?mx
中点法二阶
(Runge--Kutta)
mx 1?mx
四级四阶
RungeKutta
hxm 21? hxm 211
Euler 预估 --校正 法二阶
hxm 21
my
1?my
6.2.4
浙江大学研究生学位课程
《实用数值计算方法,96
6.3 边值问题解法
000
0000
10
'
10
'
''
'''
，，
第三边值条件：
，
第二边值条件：
，
第一边值条件：
，，，
问题：二阶微分方程
byby
ayay
byay
byay
baxyyxfy
浙江大学研究生学位课程
《实用数值计算方法,97
2212
1211
21
,,2,10,,,,
,,2,10,,,,
,,2,1,
,,,,
nkyyybB
njyyyaB
Nibax
yyyxf
dx
xdy
Nk
Nj
Ni
i
一般地，
6.3.1 试射法 （ Shooting）
y
01?jB x
02?kB
图 6.16 试射法求解示意浙江大学研究生学位课程
《实用数值计算方法,98
6.3.1
的函数是显然直到满足条件。
问题。值，继续求解新的初值修改则满足要求，否则若设求得
，
，，
，则可解初值问题设
，
，，
考虑第一边值问题
11
1
1
1
1
'
'''
1
'
'''
~
uuFby
u
v
vby
uayay
yyxfy
uay
byay
yyxfy
浙江大学研究生学位课程
《实用数值计算方法,99
6.3.1
uGu
uGubyby
uGbyuGby
uayuay
uay
uFu
uu
uu
FuuF
使要找的函数是于是设解为
，
：则问题变为解初值问题对第三边值问题，设使得直到找到
：不知道的。用线性近似是但的要求
,
~
,
~
,
.
,
0
'
2
'
1
01
'
1
12
12
1
11
浙江大学研究生学位课程
《实用数值计算方法,100
6.3.2 差分方法 （ Difference Methods）
(Relaxation Methods)
0x 1x 2x
Mx
01?jB
02?kB
iy
True solution
2nd iteration
1st iteration
x
212
1101
021
,,10,,,
,,10,,,
,,,,,,
nkyyxB
njyyxB
xxxyyyxf
dx
xdy
NMk
Nj
MNi
i
图 6.17 差分方法示意浙江大学研究生学位课程
《实用数值计算方法,101
6.3.2
：个方程由边值条件提供另外个。上述方程共有未知数。共有；：要求的值在即其中，
还有别的方程可用用差分方程代替：
维向量是，
N
MN
MN
MkNjy
Mk
xyyyyy
Nfyyxf
dx
dy
jk
k
T
Nkkkk
1
,,1,0,,1
,,2,1
,,,,
,,
21
02,2 1111 kkkkkkkkk yyxxfxxyyE
NnnnyxBE
nyxBE
MMM
21
22
'
10010
0,
0,
个个浙江大学研究生学位课程
《实用数值计算方法,102
6.3.2
2
'
1
'
10
1
0
0
0
11
1
1
1 1
1
1
111
1
,,1
,,1
,,2,1;,,2,1
,0
,,
:,
njEy
y
E
njEy
y
E
MkNj
Ey
y
E
y
y
E
yyE
y
y
E
y
y
E
yyEyyyyE
yy
jM
N
i
iM
iM
jM
j
N
i
i
i
j
jk
N
i
ik
ik
jk
N
i
ik
ik
jk
N
i
N
i
ik
ik
k
ik
ik
k
kkkkkkkk
kk
另外，由边值条件得则有设每次迭代有修正值浙江大学研究生学位课程
《实用数值计算方法,103
6.3.2
很大。
很大时，矩阵的阶个特殊稀疏问题。当它是一可解下面方程组求得。的修正量时每次迭代当
NM
M
y
nnMN
jk
1
2,3,3,5
21
23
'
13
'
53
43
33
23
12
52
42
32
22
12
51
41
31
21
11
30
20
10
53
43
33
23
13
52
42
32
22
12
51
41
31
21
11
50
40
30
20
10
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
E
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
y
图 6.18 修正值矩阵浙江大学研究生学位课程
《实用数值计算方法,104
块对角线性代数方程组的快速求解利用矩阵的特殊结构使总运算次数达到极小，
并使矩阵系数的存贮达到极小
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
BV
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
图 6.19 Gauss 消去法的目标结构
6.3.2
浙江大学研究生学位课程
《实用数值计算方法,105
D
D
D
1
0
0
1
0
0
Z
Z
Z
Z
Z
1
0
0
0
0
0
0
0
D
D
D
0
1
0
0
1
0
Z
Z
Z
Z
Z
0
1
0
0
0
0
0
0
D
D
D
0
0
1
0
0
1
Z
Z
Z
Z
Z
0
0
1
0
0
0
0
0
A
A
A
S
S
S
S
S
S
D
D
D
D
D
S
S
S
1
0
0
0
0
A
A
A
S
S
S
S
S
S
D
D
D
D
D
S
S
S
0
1
0
0
0
D
D
D
D
D
0
0
1
0
0
V
V
V
V
V
V
D
D
D
D
D
0
0
0
1
0
D
D
D
D
D
0
0
0
0
1
A
A
A
A
A
S
S
S
S
S
A
A
A
S
S
S
A
A
A
A
A
S
S
S
S
S
V
V
V
V
V
V
V
V
V
V
V
V
V
V
V
V
S
S
S
A
A
A
A
A
S
S
S
S
S
S
S
S
D--有待于对角化
A--将要改动
S--需要存贮
Z--将要变为 0
a)
b)
a)
b)
6.3.2
图 6.20 特殊矩阵结构 1
浙江大学研究生学位课程
《实用数值计算方法,106
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
1
0
0
0
0
0
1
0
0
0
0
1
0
0
0
0
0
1
0
0
Z
Z
0
0
1
0
0
0
0
0
0
0
1
0
Z
Z
0
0
0
1
0
0
0
0
0
0
0
1
Z
Z
0
0
0
0
1
0
0
S
S
S
S
S
D
D
S
S
S
S
S
1
0
S
S
S
S
S
D
D
S
S
S
S
S
0
1
V
V
V
V
V
V
V
V
V
V
V
V
V
V
S
S
S
S
S
A
A
S
S
S
S
S
S
S
a)
b)
包括四步运算
1,使用前一步的结果，简化部分元素为零
2,消去剩下子块正方结构为单位阵
3,存储以后要使用的非零系数
4,回代
6.3.2
图 6.22 特殊矩阵结构 2
浙江大学研究生学位课程
《实用数值计算方法,107
6.3.3 样条函数在两点边值问题中的应用考察 二阶线性微分方程 的边值问题
226
216
206
'
0
'
00
'''
nnn rbyby
rayay
xryxqxyxpxy
可以证明：
，个约束条件：
时，当上的：是定义在节点而
，
来近似：用三次样条函数可以的解假定微分方程边值问题的区间和节点为：设
2,1,0,,3,2,1)1(3
,,2,1
226206
1
1
23
10
jnin
xSxS
xxxni
dxcxbxaxSxS
xxS
baxxSxy
xS
xy
bxxxa
x
k
j
ik
j
i
ii
iiiii
k
n
1
0
32
00 6
1
2
1 n
k
kk xxhxxgxxfexS
bax,?
浙江大学研究生学位课程
《实用数值计算方法,108
6.3.3
此时为，而边值条件其中则有
，，并令代入微分方程把易知：
个待定常数。是而
，
其中，
226216
,,
206,,
,
2
1
3,,,,,,
,0
'''
1
0
''
1
0
2
0
'
110
iiiiii
i
n
k
kk
n
k
kk
n
k
kk
k
xrrxqqxpp
xxxSxSxS
xxhgxS
xxhxxgfxS
nhhhgfe
xx
xxxx
xx
236,,2,1,0,
6
1
2
1
1
0
32
00
1
0
1
0
2
0
nir
xxhxxgxxfeq
xxhxxgfpxxhg
i
n
k
kikiii
n
k
n
k
kikiikik
246000 rfe
浙江大学研究生学位课程
《实用数值计算方法,109
6.3.3
256
2
1
6
1
2
1
1
0
2
0
1
0
32
00
n
n
k
knknn
n
k
knknnn
rxxhxxgf
xxhxxgxxfe
3
2
3
1
3
0
2
333
000
''
110
3
2
3
1
0246
001
00,1
110
1
3
2
3
1
010
01
00
1,001
,,,,,,
33256246236
xhxhxhgxfxxS
e
r
r
rqp
y
y
xxyxy
xShhhgfe
nn
n
现设得由
，，
，
，，
有：
，，，，三等分。有节点：，把区间问题例：给定微分方程边值求得。的值，故未知数个个方程可解出共，，于是联立
浙江大学研究生学位课程
《实用数值计算方法,110
6.3.3
1,0
3
2
061224.0
3
1
061224.0
061224.05.0540814.0
~
~
061224.0061224.0061224.0
0.05.0540814.0
12
2795755
27186611655
272758118911655
02727278
3
2
6
3
1
662
3
2
3
3
1
332
33
32
210
0
01
012
012
210
''
2
2
2
1
2
0
'
x
xx
xxxxy
xyxy
hhh
egf
g
fgh
fghh
fghhh
fghhh
xhxhxhgxS
xhxhxhgxfxS
可写为的近似解故这组方程的解为条件有代入原微分方程和边界有浙江大学研究生学位课程
《实用数值计算方法,111
的实数均为大于时的捕食系统考虑系数为两种群间的交叉亲疏为本种群的亲疏系数系数为两个种群的自然增长模型例，生态系统的
0
,,,
0
,
,
,
2122
211211
21
12
21
00
2221102
1221101
xbxkx
bakkxaxkx
ba
ba
ba
ba
xxbxbbx
xxaxaax
V ol t e r r a
竞争系统，
捕食系统，
共栖系统，
时，当
00)
00)
00)
0
12
12
12
12
baiii
baii
bai
ba
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,112
为实数，
函数：阶为第一类其中，
，
，
精确解为
，积分区间
，
，
例，求解初值问题
v
kvk
x
xxJ
B e s s e lvvxJ
xJyxJy
xJyxJy
xx
xJxyxJxy
xJxyxJxy
y
x
y
dx
dy
y
x
y
dx
dy
y
x
y
dx
dy
y
dx
dy
k
k
v
v
v
0
2
34123
12101
21
13141213
11121011
43
4
32
3
21
2
2
1
1!
4
1
2
1
3,2,1,0
0.100.1,
3
2
1
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,113
相平面上的轨迹为
0201,xx
2x
1x
ak1
62k
因。而食用鱼反而减少的原期间捕鱼量减少，该模型解释了一次大战周期解的平均值平衡点
a
k
b
k 12
,
年度 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923
比例 1 1,9 2 1,4 2 2,1 2 1,2 3 6,4 2 7,3 1 6,0 1 5,9 1 4,8 1 0,7
1914年 --1923年捕获的鱼中，食用鱼所占的比例 (%)
图 6.23 相平面上的轨迹图
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,114
为扩散系数
，
，
统考虑密度分布的捕食系少。食用鱼增加，食肉鱼减此时，平衡点为
，则模型改为系数因此，我们可控制捕捞
21
2122
2
2
2
2
2
2111
2
1
2
1
1
2
1
12
2122
1211
,
,
DD
SbSSk
x
S
D
t
S
SaSSk
x
S
D
t
S
txS
txS
a
k
b
k
xbxkx
xaxkx
6.3.3
浙江大学研究生学位课程
《实用数值计算方法,115
隐式。半隐式：可化为显式的新节点值组才能解出隐式格式：应通过方程节点值式可显式解出新显式格式：通过递推公每一个初值即可计算自动起步：由即可计算出单步法：仅由方法类型与特点：
.5
,,,,,
.4
,,,,,,
.3
,2,1
.2
.1
00111
0011
00
1
hyxyxFy
hxyxyxFy
iy
yxy
yy
iii
iiii
i
ii
公式简单，但精度不好编程简单，使用方便计算量少，但稳定性差计算量大，但稳定性好
6.3.3
总 结浙江大学研究生学位课程
《实用数值计算方法,116
.1,
2
1
0
2
1
1
111
1,11
.4
1.0
10
.3
06
,,2,4
6
1
.2
0 02.00 01.03.005.0
1001 000
.1
2''
5
'
111
21
21
2122
211211
，，，取节点为值问题用样条函数方法求解边之值。项。并计算取至初值问题：级数展开求解微分方程用
，的绝对稳定区域为算法为三阶的，并证明上式所定义的隐式单步公式试证明由
，，，
，
轨迹图，并讨论的相，作出捕食系统方法求解初值问题用四阶
yy
xyxxy
yx
y
xyyx
dx
dy
Tay lor
yxhfyxfyxfhyy
bakk
xx
xbxkx
txtxxaxkx
K ut t aR un ge
nnnnnnnn
第六章 习题浙江大学研究生学位课程
《实用数值计算方法,117
)(
1
,,,,,
)2(;//1
)(
)(
.5
211221
2
1122
2211222
1221111
O D E I N T
K ut t aR u ng e
aabakk
x
akak
xbxxakx
xxaaxkx
果。
）而导致的结条件（以上结论，并解释由于方法验证用带精度控制的四阶和初值。试选取适当的系数
，最后趋于灭绝。某一最大值后开始减少开始时增大，直至达到食肉种群量则
。初始时刻食饵十分丰富
）（
满足捕食系统
