2004-12-13 1
第九章常微分方程初值问题的数值解法
§9.1 引言
§9.2 Euler方法
§9.3 Runge –Kutta公式
§9.4单步法的进一步讨论
§9.5线性多步法数值算例
2004-12-13 2
§9.1 引言一、问题本章着重讨论一阶常微分方程初值问题
=
==
′
00
yxy
yxf
dx
dy
y
）（
），（
在区间[a,b]上的数值解法。
这些问题多数情况下求不出解析解，只能用近似的方法求解。常用的近似方法有两类。一类称为近似解析法，如级数解法，逐次逼近法等。另一类称为数值解法，它可以给出解在一些离散点上的近似值。
2004-12-13 3
二、解的存在唯一性若f（x,y）在区域D=,上连续，且关于y满足李普希兹（Lipschitz）条件,即存在常数L,使
{ }Rybxa ∈≤≤，
，
2121
),(),( yyLyxfyxf?≤?
G
21
,yy
对内任两个均成立，其中L是与x,y无关的常数,则上面的初值问题存在唯一解，且解是连续可微的。
Remark：在f(x,y)对y可微的情况下，若偏导数有界，
则可取，此时有之间。与介于
21212121
y,)(
),(
),(),( yyyLyy
y
xf
yxfyxf ξ
ξ
≤?
=?
y
yxf
L
Dyx
=
∈
),(
max
),(
此时Lipschitz连续条件显然成立。这是验证该条件的最简便的方法。
2004-12-13 4
三、解的适定性解的适定性是指解的存在唯一性以及数值稳定性。
此处主要是指解对于右端项以及初值的扰动的适应性。关于适定性有如下的结论：
定理：若f（x,y）在区域D= 上满足Lipschitz连续条件，则初值问题是适定的。
{ }Rybxa ∈≤≤，
四、等价的积分方程若y(x)是初值问题的解，对方程两边同时积分，
利用初始条件可得：
∫
+=
x
x
dttytfxyxy
0
))(,()()(
0
该方程为与初值问题同解的积分方程，我们可以从积分方程出发去构造初值问题的求解公式。
2004-12-13 5
五、数值解法常微分方程初值问题的数值解是求上述初值问题的解y(x)在存在区间[a，b]中的点列
），，（nihxx
iii
L10
1
=+=
上的近似值。称为步长，i
y
i
h
一般情况下我们取等步长，记为h。
)(
n
xy
初值问题的解析解（理论解）用表示,数值解法的精确解用表示，并记f
n
=f(x
n
,y
n
),
而。
n
y
))(,()(
nnn
xyxfxy =
′
求初值问题的数值解一般是逐步进行的，即计算出y
n
之后计算y
n+1
。
2004-12-13 6
数值解法的分类数值解法一般分为：
（1）单步法：在计算时，只用到,
和，即前一步的值。
n
x
1+n
y
n
y
1+n
x
（2）多步法：计算时，除用到,和以外，还要用和,即前
k步的值。
1+n
y
pn
y
1+n
x
n
x
n
y
pn
x
)0;,2,1( >= kkp L
单步法和多步法都有显式和隐式方法之分。显式和隐式的单步法可以分别写成：
），（hyxhyy
nnnn
,
1
φ+=
+
），（hyyxhyy
nnnnn
,,
11 ++
+= φ
对多步法来说，显式和隐式方法具有相同的意义。
Back
2004-12-13 7
§9.2 Euler方法一、显式欧拉（Euler）方法设节点为，欧拉方法的计算公式
），，（nnnhxx
n
,,210
0
L=+=
），，，）（，（L210
1
=+=
+
nyxhfyy
nnnn
这是一种最简单的显式单步方法，该方法可以通过不同的途径获得。
),2,1,0)(,(
1
L==
+
nyxf
h
yy
nn
nn
1、差商法用两点差商公式代替导数,再用表示的近似值,则得到
nn
nn
xx
xyxy
+
+
1
1
)()(
n
y
)(
n
xy
)(
n
xy
′
2004-12-13 8
显式欧拉（Euler）方法（续）
假设在附近把y（x）展成Tayler级数
n
x
L+++=+ )(''
2
)(')()(
2
nnnn
xy
h
xhyxyhxy
取h的线性部分,并用表示的近似值,得
n
y )(
n
xy
),2,1,0)(,(
1
L=+=
+
nyxhfyy
nnnn
2、Taylor展开法
3、数值积分法
∫
+
+=
+
1
))(,()()(
1
n
n
x
x
nn
dttytfxyxy
对微分方程两端从x
n
到x
n+1
积分，得等价的积分方程对右端的积分部分采用左矩形公式近似即得Euler公式。
Remark：Taylor展开法与数值积分法是构造微分方程数值解的两类主要的方法。
2004-12-13 9
欧拉（Euler）方法的几何意义
Euler方法有明显的几何意义。
如右图所示，一阶常微分方程初值问题的解曲线y(x)过点P
0
(x
0
,y
0
)。从P
0
出发以
f(x
0
,y
0
)为斜率作一直线段，
与x=x
1
相交于点P
1
(x
1
,y
1
)，
显然有y
1
=y
0
+hf(x
0
,y
0
)。
同理，再由P
1
出发以f(x
1
,y
1
)为斜率作一直线段，与x=x
2
相交于点P
2
(x
2
,y
2
)，显然有y
2
=y
1
+hf(x
1
,y
1
)。这样一直做下去，得到一条折线P
0
P
1
P
2
……，作为y(x)的近似曲线。
因此，Euler方法又称为Euler折线法。
2004-12-13 10
二、隐式Euler方法和Euler方法的改进
dttytfxyxy
n
n
x
x
nn
），（）（）
∫
+
+=
+
1
)((
1
中将积分用右矩形公式
)),((,)(
2
)('
2
baab
f
bfabdxxf
b
a
∈≈
∫
η
η
）（）（）（
代入，有
K,2,1,0)),(,()()(
111
=+≈
+++
nxyxhfxyxy
nnnn
若在等价积分方程
K,2,1,0),,(
111
=+≈
+++
nyxhfyy
nnnn
从而得到上式是一个隐式的单步方法，称为隐式欧拉法或后退的欧拉法。利用此方法，每一步都要把上式作为y
n+1
的一个方程来求解。从数值积分的误差分析，很难期望隐式欧拉法比显式欧拉法更精确。
2004-12-13 11
隐式Euler方法和Euler方法的改进（续）
为了得到更精确的方法，在等价的积分方程中用梯形公式
),(),(''
12
][
2
3
baf
ab
bfaf
ab
dxxf
b
a
∈
+
=
∫
ηη
）（
）（）（）（
L,2,1,0)],,(),([
2
111
=++=
+++
nyxfyxf
h
yy
nnnnnn
近似积分项，再分别用y
n
，y
n+1
代替y(x
n
)和y(x
n＋1
)，即可得到梯形方法：
该方法也是一种隐式的单步方法。对该方法，
)],(),([
2
1
),,,(
111 +++
+=
nnnnnnn
yxfyxfhyyxφ
从n＝0开始计算，每步都要求解y
n+1
的一个方程。一般来说，这是一非线性方程，可迭代计算如下：
2004-12-13 12
隐式Euler方法和Euler方法的改进（续）
L,2,1,0,
)],(),([
2
),(
][
11
]1[
1
]0[
1
=
++=
+=
++
+
+
+
s
yxfyxf
h
yy
yxhfyy
s
nnnnn
s
n
nnnn
使用上式时，先用第一式计算出y
n+1
的近似值，再用第二式反复进行迭代，得到数列，用来控制是否继续进行迭代，其中ε为允许误差。把满足要求的作为y(x
n+1
)的近似值y
n+1
，类似地可得出y
n+2
，
y
n+3
，…。
]0[
1+n
y
{}
∞
=
+
0
][
1
s
s
n
y
ε≤?
+
+
+
][
1
]1[
1
s
n
s
n
yy
]1[
1
+
+
s
n
y
当f(x,y)关于y满足Lipschitz条件时，且步长h满足时，上述迭代过程是收敛的。这是因为：
1
2
1
<hL
2004-12-13 13
]0[
11
][
11
][
11
][
11
][
1111
]1[
11
2
)(
)(
2
)
2
++++++++
++++
+
++
<<?<?≤?
=
=?
nn
s
nn
s
nn
s
nn
s
nnnn
s
nn
yyyyyy
hL
yy
y
fh
yxfyxf
h
yy
L
ξ
，（），（
)(0)
2
(
]0[
11
1]1[
11
∞→→?≤≤?
++
++
++
syy
hL
yy
nn
ss
nn
L
实用中，h取得较小时,梯形方法第二式只迭代一次就结束，得到Euler预估—校正方法(改进的Euler方法)：
++=
+=
+++
+
)],(),([
2
),(
]0[
111
]0[
1
nnnnnn
nnnn
yxfyxf
h
yy
yxhfyy
其中第一式称为预估算式，第二式称为校正算式。
隐式Euler方法和Euler方法的改进（续）
其中ξ介于与之间。
1+n
y
][
1
s
n
y
+
1
][
1
,
++
→∞
n
s
n
yy
即为s趋于。
2004-12-13 14
隐式Euler方法和Euler方法的改进（续）
若将Euler预-校方法中的第一式带入第二式，得
))],(,(),([
2
11 nnnnnnnn
yxhfyxfyxf
h
yy +++=
++
Remark：这是一种显示的单步方法。有时为了计算方便，常将上式改写成：
++=
=
++=
+
),(
),(
)(
2
12
1
211
hKyhxfK
yxfK
KK
h
yy
nn
nn
nn
2004-12-13 15
三、单步法的局部截断误差和阶设一般的单步法为：
显式公式：
），，（hyxhyy
nnnn
φ+=
+1
隐式公式：
），，，（hyyxhyy
nnnnn 11 ++
+= φ
设为数值方法的精确值，y（）为微分方程的精确解。
n
x
n
y
nnn
yxye?=）（
n
x
定义1：为某一数值方法在处的整体截断误差。
Remark：整体截断误差不仅与这步的计算有关，而且与前面所有点的计算的误差累计有关。为了简化误差分析，我们着重分析计算中的某一步。对一般的显式单步法，有如下定义：
n
x
2004-12-13 16
单步法的局部截断误差和阶（续）
定义2：对单步法，在的假设下，称为在处的局部截断误差。
）（
nn
xyy =
)),(,()()(
1,
hxyxhxyxyR
nnnnhn
φ=
+
1+n
x
）（
nn
xyy =
若设，即第n步及以前各步都没有误差，
则由显示单步法计算一步所得之与之差为：
1+n
y
)(
1+n
xy
)),(,()()(
)],,([)()(
1
111
hxyxhxyxy
hyxhyxyyxy
nnnn
nnnnnn
φ
φ
=
+?=?
+
+++
即在的假设下，。
11,
)(
++
=
nnhn
yxyR
）（
nn
xyy =
这就是上面定义中称R
n,h
为“局部”的含义，我们一般用该式作为定义。这里应该注意，R
n,h
和整体截断误差
e
n+1
是不同的。
2004-12-13 17
单步法的局部截断误差和阶（续）
定义3：若一个单步方法的局部截断误差为O(h
p+1
)，即则称该方法为p阶方法（其中p为正整数）。
)()),(,()()(
1+
=+
p
hOhxyxhxyhxy φ
Remark：由前面的定义可知，若某个单步方法是一种
p阶方法，则有R
n,h
=O(h
p+1
)，即p阶方法的局部截断误差为h的p＋1阶。我们往往比较关心R
n,h
按h展开式的第一项。
定义4：若一个单步方法是一种p阶方法，其局部截断误差可以写成：
)())(,(
21
,
++
+=
pp
nnhn
hOhxyxR ψ
则Ψ(x
n
,y(x
n
))h
p+1
称为方法的主局部截断误差，或局部截断误差的主项。
2004-12-13 18
单步法的局部截断误差和阶（续）
例1：求Euler方法的局部截断误差。
由Taylor展式，有
)()('''
!3
)(''
!2
)()()(
))(,()()()),(,()()(
2
32
hOxy
h
xy
h
xyhxyhxy
xyxhfxyhxyhxyxhxyhxy
=++=
′
+=
+=+
L
φ
故Euler方法是一阶方法,局部截断误差为：
)()(''
2
3
2
,
hOxy
h
R
nhn
+=
主局部截断误差为：。
)(''
2
2
n
xy
h
2004-12-13 19
单步法的局部截断误差和阶（续）
例2：求Euler预估-校正方法
++=
+=
+++
+
)],(),([
2
),(
]0[
111
]0[
1
nnnnnn
nnnn
yxfyxf
h
yy
yxhfyy
的局部截断误差。
改写上述公式为：
++=
=
++=
+
),(
),(
)(
2
12
1
211
hkyhxfk
yxfk
kk
h
yy
nn
nn
nn
)('))(,(),(
1 nnnnn
xyxyxfyxfk ===
2004-12-13 20
单步法的局部截断误差和阶（续）
),(
12
hkyhxfk
nn
++=
L+
+
+
+
+
+=
)],(),(2),([
2
1
)],(),([),(
2
2
2
1
2
2
1
2
2
2
2
1
nnnnnn
nnnnnn
yx
y
f
khyx
yx
f
khyx
x
f
h
yx
y
f
hkyx
x
f
hyxf
)(),(),(),(),(
2
hOyxfyxhfyxhfyxf
nnynnnnxnn
+++=
)(][),(
2
))(,(
hOfffhyxf
nn
xyx
yxnn
+++=
)()('')('
2
hOxhyxy
nn
++=
2004-12-13 21
单步法的局部截断误差和阶（续）
代入，得
)()]('')(')('[
3
2
1
hOxhyxyxyyy
nnn
h
nn
++++=
+
而
L++
′
+=
+
)(''
2
)()(
2
1 nnnn
xy
h
xyhxyxy（）
因此有
)()(
3
11,
hOyxyR
nnhn
=?=
++
故Euler预估-校正方法为二阶方法。
2004-12-13 22
单步法的局部截断误差和阶（续）
例3：求梯形公式的局部截断误差。
))(,(''
12
))](())(([
2
)()(
3
111
εε yf
h
xyxfxyxf
h
xyxy
nnnnnn
++=
+++
，，
由微分中值定理
))((),())(,(
11),(1111
1
++++++
=?
+
nnxnnnn
yxy
y
f
yxfxyxf
n
η
代入上式得
)(
12
))()(
2
1(
3
11),(
1
ε
η
y
h
yxy
y
fh
nnx
n
′′′
=?
++
+
)]()([
2
111 +++
++=
nnnnnn
yxfyxf
h
yy，，
上面两式相减，有
)(
12
)]())(([
2
)(
3
11111
εy
h
yxfxyxf
h
yxy
nnnnnn
′′′
=?
+++++
，，
2004-12-13 23
单步法的局部截断误差和阶（续）
故
)()]('''
12
][
2
1[)(
3
3
),(11
1
hOy
h
y
fh
yxy
n
xnn
=?+
+=?
+
++
ε
η
L
1
2
),(
1
<
+
η
n
x
y
fh
当h充分小时，若，则可以选择h，s.t,
,从而有
1<≤
L
y
f
L+
+=
+
+
),(
1
1
2
1
2
1
1
η
η
n
n
x
x
y
fh
y
fh
），（
即梯形公式为二阶公式。
Back
2004-12-13 24
§9.3 Runge –Kutta公式一、Taylor 方法
Taylor展开法与数值积分法是推导高阶方法的常用手段。
本节以Taylor展开法为基础，介绍如何推导高阶单步法。
假定初值问题的解y(x)以及函数f(x,y)是足够光滑的，有
)()(
!
)(''
!2
)(')(
1)(
2
1
+
+
+++++=
p
n
p
p
nnnn
hOxy
p
h
xy
h
xhyxyxy L）（
)())(,(
!
))(,('
!2
))(,()(
1)1(
2
+?
+++++=
p
nn
p
p
nnnnn
hOxyxf
p
h
xyxf
h
xyxhfxy L
)(
2
1
!!2
p
n
p
nnnn
y
p
h
y
h
yhyy ++
′′
+
′
+=
+
L
当h充分小时，略去余项，将y
(k)
(x
n
)用来代替，
则有p阶Taylor方法：
)(
1+p
hO
)(k
n
y
)(k
n
y
其中(k=1,2,…,p)根据求导法则，其计算公式为：
2004-12-13 25
Taylor 方法（续）
),('
nnn
yxfy =
),(),(),(''
nnynnnnxn
yxfyxfyxfy +=
M
),](2['''
22
nnyyxyyxyxxn
yxfffffffffy ++++=
Remark1：显然，p阶Taylor方法的局部截断误差为
。当p=1时，
Taylor方法就是Euler方法。当p≥2时，需要计算公式中的高阶导数。
)()(
)!1(
1
2)1(1
,
+++
+
+
=
p
n
pp
hn
hOxyh
p
R
Remark2：显然，Taylor方法可以得到任意阶精度的方法。但在实际计算中，Taylor方法往往相当困难，因为公式中的高阶导数会很复杂。故Taylor方法很少单独使用，但常用它来启发思路。
2004-12-13 26
二、Runge—Kutta方法
Euler公式：
）（
），（
2
,
1
11
hOR
yxfk
hkyy
hn
nn
nn
=
=
+=
+
改进的Euler预估校正格式：
)(
)
2
3
,
12
1
211
hOR
hkyhxfk
yxfk
kk
h
yy
hn
nn
nn
nn
=
++=
=
++=
+
），（
），（
（
基本思想：用不同点的函数值作线性组合，构造高阶单步的近似公式，把近似公式和解的Taylor展开式比较，使前面尽可能多的项完全相同。这种方法间接应用Taylor展开的思想，避免了高阶导数计算的困难。
2004-12-13 27
一般的Runge-Kutta方法的形式为
=
++=
=
=
+=
∑
∑
=
=
+
riKhyhxfk
yxfk
Kchyx
hyxhyy
i
j
jijnini
nn
r
i
iinn
nnnn
,,3,2,
),,(
),,(
1
1
1
1
1
Lβα
φ
φ
，
），（
其中，为常数，选取这些常数的原则是，
要求第一式的右端在处泰勒展开后，按h的幂次重新整理，得到
iiji
C，，βα
），
nn
yx(
Runge—Kutta方法（续）
L++++=
+
3
3
2
211
3
1
2
1
hVhVhVyy
nn
！！
2004-12-13 28
Runge—Kutta方法（续）
L+
′′′
+
′′
+
′
+=
+
)(
3
)(
2
)()()(
32
1 nnnnn
xy
h
xy
h
xyhxyxy
！！
有若干项重合，即要求与微分方程的解y(x
n+1
)在x
n
处的展开式
L，，，)()()(
321 nnn
xyVxyVxyV
′′′
=
′′
=
′
=
上述公式叫做r级Runge-Kutta方法(c
r
≠0)。
当r＝1时，上式中φ(x,y,h)=c
1
f(x,y)。和泰勒展开式中
p=1时的φ(x,y,h)=f(x,y)比较知，取c
1
＝1即可。这就是
Euler方法。
下面用Taylor展开方法对r(≥2)级的Runge-Kutta方法确定参数，使其成为p (≥2)阶的方法。
iiji
C，，βα
2004-12-13 29
Runge—Kutta方法（续）
下面以二级Runge—Kutta公式为例进行具体推导。
++=
=
+=
+=
+
），（
），（
12122
1
2211
1
),,(
),,(
KhyhxfK
yxfK
KcKchyx
hyxhyy
nn
nn
nn
nnnn
βα
φ
φ
对要求适当选取系数，使当时，
上式的局部截断误差为，即成为二阶方法。
21221
,,,βαcc）（
nn
xyy =
）
3
(hO
做法一：按照定义，根据局部截断误差的定义，将上式代入，并按h展开，得：
))](,()([)()(
),,()()(
21221
xyhyhxfcxychxyhxy
hyxhxyhxy
′
+++
′
+=
+
βα
φ
2004-12-13 30
Runge—Kutta方法（续）
L+
′′′
+
′′
+
′
=+ )(
3
1
)(
2
1
)()()(
32
xyhxyhxyhxyhxy
！！
＋
)(]))(()(2)[(
!2
1
])([
!1
1
),())(,(
32
21212
2
2
212212
hOfxyhfxyhhfh
fxyhhfyxfxyhyhxf
yyxyxx
yx
+
′
+
′
++
′
++=
′
++
ββαα
βαβα
将上两式代入，利用Taylor展式，并按h的幂次整理，得：
)()](
6
1
)
2
1
6
1
(
)
3
1
()
2
1
6
1
[(
])
2
1
()
2
1
[()()1(
),,()()(
422
212
2122
2
22
3
21222
2
21
hOffffffc
ffcfch
ffcfchxyhcc
hyxhxyhxy
yxyyy
xyxx
yx
+++?
+?+?+
+?+
′
=
+
β
βαα
βα
φ
2004-12-13 31
Runge—Kutta方法（续）
==
=+
2
1
1
21222
21
βα cc
cc
要使，则需
)(
3
,
hOR
hn
=
上述方程有无穷多组解，所确定的解都能使二级
Runge-Kutta方法成为一个二阶方法。
由于二级方法，方程也可写为
0
2
≠c
==
=
2
212
21
2
1
1
c
cc
βα
在上述方程中，选取不同的c
2
，即可获得不同的二级二阶Runge-Kutta方法。
2004-12-13 32
Runge—Kutta方法（续）
)()],()),((
),(),(2),()[(
2
1
)],(),(),([
!1
1
),(
)),(,(
32
21
212
2
2
212
2122
hOyxfyxfh
yxfyxfhhyxfh
yxfyxfhyxhfyxf
yxfhyhxfK
nnyynn
nnxynnnnxx
nnynnnnxnn
nnnn
+
+?+
++=
++=
β
βαα
βα
βα
！
＋
做法二：按照定义，将K
2
在(x
n
,y
n
)处展开，有将K
1
,K
2
代入y
n+1
,有，
)()],(),(),(),(2),([
2
)],(),(),([)()()(
),(
422
21212
2
2
3
2
212
2
221
2211
hOyxfyxfyxfyxfyxf
hc
yxfyxfyxfhcxyhccxy
Khcyxfhcyy
nnyynnnnxynnnnxx
nnynnnnxnn
nnnn
++++
++
′
++=
++=
+
ββαα
βα
2004-12-13 33
Runge—Kutta方法（续）
),()(' yxfxy =
),(),(),()('' yxfyxfyxfxy
yx
+=
由
M
),](2[)('''
22
yxfffffffffxy
yyxyyxyxx
++++=
)(),)]((
6
1
)
2
1
6
1
(
)
3
1
()
2
1
6
1
[(),]()
2
1
(
)
2
1
[()()1()(
422
212
2122
2
22
3
212
22
2
2111,
hOyxffffffc
ffcfchyxffc
fchxyhccyxyR
nnyxyyy
xyxxnny
xnnnhn
+++?+
+
+?+
′
=?=
++
β
βααβ
α
＋
故同第一种做法一样，可以获得无穷多组二级二阶
Runge-Kutta方法。常用的二阶Runge-Kutta方法有：
2004-12-13 34
Runge—Kutta方法（续）
1,
2
1
21
221
====
β
αc
++=
=
++=
+
),(
),(
)(
2
12
1
211
hKyhxfK
yxfK
KK
h
yy
nn
nn
nn
++=
=
+=
+
)
2
,
2
(
),(
12
1
21
K
h
y
h
xfK
yxfK
hKyy
nn
nn
nn
2
1
,1,0
21
221
====
β
αcc
②
c
①
Euler预估-校正格式中（间）点公式
2004-12-13 35
Runge—Kutta方法（续）
3
2
,
4
3
,
4
1
21221
==== βαcc
③
++=
=
++=
+
)
3
2
,
3
2
(
),(
)3(
4
12
1
211
hKyhxfK
yxfK
KK
h
yy
nn
nn
nn
此时可以论证：
,0
2
1
6
1
2
2
2
=?
cα
0
3
1
2
21
2
=?
c
β
α
二阶Heun（休恩）方法此时，即二阶Heun
方法是局部截断误差项数最少的方法且二级Runge-
Kutta方法不可能达到三阶。
)()(
6
4
),(
3
,
hOffff
h
R
nn
yxyxyhn
++=
2004-12-13 36
Runge—Kutta方法（续）
用类似的方法可以研究其它级的Runge-Kutta方法。
常用的其它级的Runge-Kutta方法有：
++=
++=
=
++=
+
)
3
2
,
3
2
(
)
3
1
,
3
(
),(
)3(
4
23
12
1
311
hKy
h
xfK
hKy
h
xfK
yxfK
KK
h
yy
nn
nn
nn
nn
+?+=
++=
=
+++=
+
)2,(
)
2
1
,
2
(
),(
)4(
6
213
12
1
3211
hKhKyhxfK
hKy
h
xfK
yxfK
KKK
h
yy
nn
nn
nn
nn
三阶Heun（休恩）公式三阶Kutta公式
2004-12-13 37
标准（经典）四阶Runge-Kutta公式
Runge—Kutta方法（续）
++=
++=
++=
=
++++=
+
),(
)
2
1
,
2
(
)
2
1
,
2
(
),(
)22(
6
34
23
12
1
43211
hKyhxfK
hKy
h
xfK
hKy
h
xfK
yxfK
KKKK
h
yy
nn
nn
nn
nn
nn
++?+=
+
++=
++=
=
+++?++=
+
))
2
2
1(
2
2
,(
)
2
22
2
12
,
2
(
)
2
1
,
2
(
),(
))22()22([
6
324
213
12
1
43211
hKhKyhxfK
hKhKy
h
xfK
hKy
h
xfK
yxfK
KKKK
h
yy
nn
nn
nn
nn
nn
基尔（Gill）公式
Back
2004-12-13 38
§9.4单步法的进一步讨论一、收敛性定义：
定义：对于初值问题，如果一个单步显式方法产生近似解对于任一固定的,均有则称该单步法是收敛的。
nhxxbxx +=∈
00
],,[ )(lim
0
xyy
n
h
=
→
Remark：该定义可类似地应用于单步隐式方法以及后面的线性多步法。从定义可知，若格式收敛，整体截断误差e
n
=y(x
n
)-y
n
必然趋于零。
定理（整体截断误差与局部截断误差的关系）：若初值问题的一个单步方法之局部截断误差为且单步法中函数关于y满足lipschitz条件，则
)1(),(
1
,
≥=
+
phOR
p
hn
),,( hyxΦ
)()(
111
p
nnn
hOyxye =?=
+++
2004-12-13 39
整体截断误差与局部截断误差的关系（续）
)],,([)(
~
)(
111,
hyxhyxyyxyR
nnnnnnhn
φ+?=?=
+++
证明：根据局部截断误差的定义，有
)()),(,()()(
1
1
+
+
==
p
nnnn
hOhxyxhxyxy φ
其中为在假定下计算的处的数值解。
1
~
+n
y
)(
nn
xyy =
1+n
x
即
)),(,()(
~
1
hxyxhxyy
nnnn
φ+=
+
1
11
~
)(
+
++
<?
p
nn
chyxy
则有常数c>0，使得由于且关于满足Lipschitz条件，得
),,,(
1
hyxhyy
nnnn
φ+=
+
),,( hyxφ
y
11
~
++
nn
yy
),,()),(,()( hyxhxyxhyxy
nnnnnn
φφ?+?≤
nn
yxyhL?+≤ )()1(
2004-12-13 40
整体截断误差与局部截断误差的关系（续）
1111
~~
)(
++++
+?≤
nnnn
yyyxy
111
)(
+++
=
nnn
yxye故
n
p
ehLCh )1(
1
++≤
+
从而有递推关系：
])1()[1(
1
11
1?
++
+
++++≤
n
pp
n
ehLchhLche
1
21
)1()]1(1[
+
++++=
n
p
ehLhLch
0
1
)1( ehL
n+
++
])1()1()1(1[
21 np
hLhLhLch +++++++≤≤
+
LK
0
1
1
1
)1(
1)1(
1)1(
ehL
hL
hL
ch
n
n
p +
+
+
++
+
+
=
2004-12-13 41
若则且
,)(
00
yxy =
0
0
=e
整体截断误差与局部截断误差的关系（续）
hL
ehLhLhL =+++≤+< L
2
)(
2
1
110
得：
hLnn
ehL
)1(1
)1(0
++
≤+<
[]1
1
)1(
)1(
1
1
=
≤
+
+
+
+
hLnP
hLn
p
n
eh
L
c
hL
e
che
故当固定时，1+
=
n
xx
abxxhn
n
≤?=+
+ 01
)1(
[ ]
pLabp
n
hceh
L
c
e
1
)(
1
1 =?≤
+
所以证毕
2004-12-13 42
整体截断误差与局部截断误差的关系（续）
Remark1：该定理说明，数值方法的整体截断误差比局部截断误差低一阶。收敛的方法至少是一阶方法。
Remark2：在该定义的条件下，欧拉方法是一阶方法，欧拉预估-校正方法是二阶方法。当f(x,y)
关于y满足利普希茨条件时，r级龙格-库塔方法中的φ关于y也满足利普希茨条件，所以定理中的条件得到满足，解的收敛性得到了保证。
2004-12-13 43
二.相容性设单步法为，
),,(
1
hyxhyy
nnnn
φ+=
+
)],,([)(
1
hyxhyxy
nnnn
φ+?=
+
)),(,()()(
1
hxyxhxyxy
nnnn
φ=
+
11,
)(
++
=
nnhn
yxyR
由于,且为任意点，故该式相当于用近似方程
）（00
,
→→ hR
hn
n
x
)),(,(
)()(
hxyx
h
xyhxy
φ≈
+
))(,( xyxf
dx
dy
=来替代通过处求解近似方程而获得原方程的近似解。因此，必须要求当时，近似方程应逼近于原微分方程。
n
xx =
0→h
2004-12-13 44
相容性（续）
因此要使时，近似方程的极限状态成为原微分方程，需且只需极限
0→h
)(
)()(
lim
0
xy
h
xyhxy
h
′
=
+
→
由于
))(,()),(,(lim
0
xyxfhxyx
h
=
→
φ
成立。
)),(,( hxyxφ
由于假设为连续函数，因而上式可以表示为
))(,()0),(,( xyxfxyx =φ
定义：如果当时，近似方程能逼近微分方程，则称数值公式与原微分方程相容。
0→h
),,(
1
hyxhyy
nnnn
φ+=
+
相容的充要条件为
))(,()0),(,( xyxfxyx =φ
2004-12-13 45
相容性（续）
Remark：可以证明，若单步法的阶大于或等于1，则单步法与微分方程相容；反之，如果单步法与微分方程相容，且关于h满足Lipschitz条件，则单步法至少为一阶方法。事实上：
),,( hyxφ
⑴若单步法的阶大于或等于1，由
0)),(,()()(
1,
→=
+
hxyxhxyxyR
nnnnhn
φ
（当）知
0→h
)),(,(lim
)()(
lim
00
hxyx
h
xyhxy
hh
φ
→→
=
+
)0),(,())(,()(' xyxxyxfxy φ==
故有即单步法与微分方程相容。
2004-12-13 46
相容性（续）
⑵若单步法与微分方程相容，且关于
h满足Lipschitz条件，则
),,( hyxφ
),,()()(
111,
hyxhyxyyxyR
nnnnnnhn
φ=?=
+++
)),(,()()(
1
hxyxhxyxy
nnnn
φ=
+
)()()()(
2
nnn
xyhOxyhxy?+
′
+=
)]()0),(,([ hOxyxh
nn
+? φ
)())(,()(
2
hOxyxhfxyh
nnn
+?
′
=
)(
2
hO=
即与微分方程相容的单步法至少为一阶方法。
2004-12-13 47
相容性（续）
通过上面的讨论可以得到下面的定理：
定理：设增量函数在区域中连续，并对变量y满足利普希茨条件，
则单步法收敛的充要条件为相容性条件成立。
),,( hyxφ
}0,,),{(
0
hhRybxayxD ≤≤∈≤≤=
Remark：在满足定理的条件下，Euler方法，
Euler预估-校正格式，Runge-Kutta方法等都与原微分方程相容。
2004-12-13 48
三、稳定性定义1 *用一个数值方法求解微分方程初值问题时，对给定的步长h>0，若在计算时引入误差（也称扰动），
由此引起计算后面的时的误差按绝对值均不增加，则称这个数值方法是绝对稳定的。
n
y
n
δ
),2,1( L=
+
ky
kn
关于单步法收敛性的概念和收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的。但实际计算时通常会有舍入误差及其累积。数值求解微分方程的过程是一个递推公式，必须考虑误差积累是否得到控制。
设f(x,y)关于y满足Lipschitz条件，我们总是针对模型方程进行讨论，其中为复常数。
λ
yy λ=
′
0Re <）（λ
为保证微分方程本身的稳定性，这里假定。
2004-12-13 49
稳定性（续）
定义1设步长为h>0的单步法用于求解时，
yy λ=
′
中由引起的误差满足，则称
n
δ
m
δ
)( nm
nm
>< δδ
单步法对于所用步长h和复数是绝对稳定的。
λ
若在计算时有误差,但在计算后面的
)( nmy
m
>
n
δ
n
y
Remark1：在上面的定义中，可以取小于或等于关系符。取小于号是为了和线性多步法相一致。
Remark2：从上面的定义可知，单步法是否稳定，
与模型方程中的复数λ以及所用步长h有关。若对复平面上的某个区域G，当时，单步法绝对稳定，则称G为单步法的绝对稳定区域，G与实轴的交集为绝对稳定区间。
Gh∈λ
2004-12-13 50
稳定性（续）
①Euler显式公式将Euler显式公式用于模型方程，则有
yy λ=
′
nnnn
yhyhyy )1(
1
λλ +=+=
+
下面考察几个常用公式的稳定性。
设有误差，
n
y
n
δ nnn
yy δ+=
~
参与运算的量为，
误差为，则实际得到近似的量为
1+n
y
1+n
δ
由此引起的要求误差不增加，即，nn
δδ <
+1
1+n
y
111
~
+++
+=
nnn
yy δ
，即。与前面公式
))(1(
11 nnnn
yhy δλδ ++=+
++
nn
h δλδ )1(
1
+=
+
相减，有。
必须。
11
1
<+=
+
λ
δ
δ
h
n
n
是保证绝对稳定性对11 <+ λh
步长h加的限制。当λ为实数时，可以得到用hλ表示的绝对稳定的区间(-2,0)。
2004-12-13 51
稳定性（续）
②隐式Euler公式后退的Euler公式用到上，有，
yy λ=
′
11 ++
+=
nnn
yhyy λ
故
11 ++
+=
nnn
hδλδδ
即
nn
h
δ
λ
δ
=
+
1
1
1
令
1
1
1
1
<
=
+
λδ
δ
h
n
n
11 >? λh
得绝对稳定区域为。
可见，若取λ为实数，则对于任意的h都是绝对稳定的。
2004-12-13 52
稳定性（续）
③梯形公式梯形公式用到上，有yy λ=
′
)(
2
11 ++
++=
nnnn
yy
h
yy
λ
故)(
2
11 ++
++=
nnnn
h
δδ
λ
δδ
即
nn
h
h
δ
λ
λ
δ
+
=
+
2
1
2
1
1
设
iyxh +=λ
2004-12-13 53
稳定性（续）
iyx
iyx
++
=
2
2
2
1
2
1
λ
λ
h
h
+
=
n
n
δ
δ
1+
由于
2
1
22
22
}
)2(
)2(
{
yx
yx
+?
++
=
2
1
22
22
}
44
44
{
yxx
yxx
++?
+++
=
当x=Re( )<0,上式右端总小于1，故梯形法的绝对稳
hλ
hλ
定区域为Re( )<0,即左半平面。
2004-12-13 54
稳定性（续）
④四阶经典R-K方法四阶经典R-K方法用到上，有
yy λ=
′
n
yk λ=
1
)
2
1
()
2
1
(
2
12
λλλ hyhkyk
nn
+=+=
)]
2
1
(
2
1
[)
2
1
(
2
23
hhyhkyk
nn
λλλλλ ++=+=
)]
4
1
2
1
([)(
322
34
λλλλλλ hhhyhkyk
nn
+++=+=
)}
4
1
2
1
()
4
1
2
1
(2
43322322
λλλλλλλ hhhhh +++++++
)22(
6
43211
kkkk
h
yy
nn
++++=
+
)
2
1
(2{
6
2
λλλ hy
h
y
nn
+++=
})(
24
1
)(
6
1
)(
2
1
1{
432
λλλλ hhhhy
n
++++=
因此
2004-12-13 55
稳定性（续）
扰动满足
})(
24
1
)(
6
1
)(
2
1
1{
432
1
λλλλδδ hhhh
nn
++++=
+
令,得四阶经典R－K方法的绝对稳定区域为
1
1
<
+
n
n
δ
δ
1)(
24
1
)(
6
1
)(
2
1
1
432
<++++ λλλλ hhhh
Back
2004-12-13 56
§9.5线性多步法一、线性多步法的一般理论步的信息来预测，则可以期望获得较高的精度。
kn
y
+
线性多步法的基本思想：如果充分利用前面多
)(
n
xy
记的近似值为,并记，
则k步线性多步方法一般形式为
nhxxy
nn
+=
0
,),(
nnn
yxff =
)(
110
0
knknn
k
j
jnj
fffhfh
++
=
+
+++==
∑
ββββ L
knknnnn
k
j
jnj
yyyyyy
++++
=
+
+++++=
∑
αααααα L
3322110
0
jj
βα,
00
,,0 βαα ≠
k
其中为常数，不全为零。
2004-12-13 57
线性多步法的一般理论（续）
因为上式两边可同乘一常数，故可规定。
1=
k
α
则上式可写成
)()(
11011110 knknnknknnkn
fffhyyyy
++?+?++
+++++++?= βββααα LL
若则为隐式方法，若则为显式方法
0≠
k
β 0=
k
β
=
+=
∑
=
++
+++
1
0
)(
),(
k
j
jnjjnj
knknkkn
yfhg
gyxfhy
αβ
β
或写成上式称为线性k步法。
Remark：前面的R－K方法是增加一些非节点处的函数值的计算来提高单步法的精度的，这样使计算量增加了许多。而线性多步法每步只需要计算一个函数值。
2004-12-13 58
线性多步法的一般理论（续）
对于隐式情形（）的公式，由于f(x,y)一般是非线性函数，故难以求解到y
n+k
的显示表达式，
故常用迭代法求解：
0≠
k
β
gyxfhy
s
knknk
s
kn
+=
++
+
+
),(
)()1(
β
其中任意给出，s＝0,1,2,…,迭代到满足给定的精度要求为止。
)0(
kn
y
+
可以证明，当f(x,y)满足Lipschitz条件（或）
时，只要，迭代关系式就是收敛的。
L
y
f
≤
k
L
h
β
1
<
2004-12-13 59
线性多步法的一般理论（续）
定义：对于线性多步法定义处的局部截断误差为kn
x
+
)]()()([
110 knknn
xyxyxyh
++
′
++
′
+
′
βββ L
)()()(
110 knknn
xyxyxy
++
+++= ααα L
∑
=
+++
′
=
k
j
jnjjnjkn
xyhxyR
0
)]()([ βα
)(
1103322110 knknnknknnnn
fffhyyyyy
++++++
+++=+++++ βββααααα LL
kn
R
+按h展开的首项称主局部截断误差。
)1,,2,1,0(),(?==
++
kjxyy
jnjn
L
knknkn
yxyR
+++
= )(
kn
x
+
或在的假定下，定义处的局部截断误差为。
2004-12-13 60
线性多步法的一般理论（续）
若,则称线性多步法为p阶方法。
)(
1+
+
=
p
kn
hOR
若线性多步法为p阶方法，则
)()(
2)1(1
1
+++
++
+=
p
n
pp
pkn
hOxyhCR
设则有,1,,1,0),(?==
++
kjxyy
jnjn
L
knknknn
Rxyxyxyh
+++
+
′
++
′
+
′
= ))()()((
110
βββ L
∑
=
+++
+=
k
j
knjnjnj
Rxyxfh
0
))(,(β
)()()()(
110
0
knnn
k
j
jnj
xyxyxyxy
++
=
+
+++=
∑
Lααα
1+p
C其中称为主局部截断误差系数。
)(
)1(1
1 n
pp
p
xyhC
++
+
即主局部截断误差为
2004-12-13 61
线性多步法的一般理论（续）
将上式与有,1,,1,0),(?==
++
kjxyy
jnjn
L
相减，并注意
∑∑
=
++
=
+
=
k
j
jnjnj
k
j
jnj
yxfhy
00
),(βα
knknknknknk
Ryxfxyxfh
+++++
+?= )],())(,([β
knkn
yxy
++
)(
利用微分中值定理
))((
),(
knkn
knkn
yxy
y
xf
++
++
=
η
),())(,(
knknknkn
yxfxyxf
++++
其中介于与之间。
kn+
η kn
y
+ )(
kn
xy
+
))(](
),(
1[
knkn
knkn
kkn
yxy
y
xf
hR
++
++
+
=
η
β
故
2004-12-13 62
线性多步法的一般理论（续）
故在的假定下，若
1,,1,0),(?==
++
kjxyy
jnjn
L
0=
k
β
（显示公式），
knknkn
yxyR
+++
= )(
则：
若,且
（隐式公式）0≠
k
β
)()(
2)1(1
1
+++
++
+=
p
n
pP
pkn
hOxyhCR
即为p阶方法，则当y(x)充分可微时，
)]()(][
),(
1[
2)1(1
1
+++
+
++
++
+=
p
n
pp
p
knkn
k
hOxyhC
y
xf
h L
η
β
)]()([
),(
1
1
2)1(1
1
+++
+
++
+
=
p
n
pp
p
knkn
k
hOxyhC
y
xf
h
η
β
knkn
yxy
++
)(
)()(
2)1(1
1
+++
+
+=
p
n
pp
p
hOxyhC
2004-12-13 63
线性多步法的一般理论（续）
即的首项与的首项相同，因此两种局部截断误差的定义是相同的。
knkn
yxy
++
)(
kn
R
+
Remark1：可以利用此处的截断误差定义分析前面的单步隐式方法。对于Euler方法，其主局部截断误差为，而对于梯形方法，
其主局部截断误差为。
)(
2
1
2
n
xyh
′′
)(
12
1
3
n
xyh
′′′
Remark2：可以证明，显式线性多步法的整体截断误差比局部截断误差低一阶。
2004-12-13 64
二、基于数值积分的构造方法将方程两端从积分得))(,( xyxfy =
′
knjn
xx
+?
到
=
+=
∫
+
+
00
)(
))(,()()(
yxy
dttytfxyxy
kn
jn
x
x
jnkn
（*）
对，取等距插值节点，对应))(,( tytf
pnnn
xxx
,,,
1
L
))(,(,)),(,(
pnpnnn
xyxfxyxf
L
的函数值为，
构造p次Lagrange插值多项式：
))(,()()(
0
inin
p
i
ip
xyxftltL
=
∑
=
2004-12-13 65
基于数值积分的构造方法（续）
代入（*）式有
∑
=
+=
p
i
ininpijn
xyxfhxy
0
))(,()( β
∫
∑
+
=
+=
kn
jn
x
x
i
p
i
ininjn
dttlxyxfxy )())(,()(
0
∫
+
+≈
+
kn
jn
x
x
pjnkn
dttLxyxy )()()(
其中
，
ii
xxh?=
+1
∫
+
=
kn
jn
x
x
ipi
dttl
h
)(
1
β
),2,1,0(
)(
)( pi
xx
xt
tl
p
il
il
lnin
ln
i
L=
=
∏
≠
=
2004-12-13 66
基于数值积分的构造方法（续）
若作变量代换
。则],[],[,kjxxshxt
knjnn
+=
+?
hlslhxshxxshxxt
nnlnnln
)()( +=+=?+=?
hlilhxihxxx
nnlnin
)()()( +?==?
代入
∫
+
=
kn
jn
x
x
ipi
dttl
h
)(
1
β
则
∫
∏
≠
=
=
+?
+
=
k
j
p
il
l
pi
pids
li
ls
0
),,2,1,0( Lβ
，表示代替用),(),(
ininininin
yxffxyy
∑
=
+
+=
p
i
inpijnkn
fhyy
0
β
对k，j和p选择不同的值，即可以得到不同类型的具体公式。
则可得线性多步法显式公式
2004-12-13 67
基于数值积分的构造方法（续）
取可得到Adams显式公式
0,1 == jk
)(
1101 pnppnpnpnn
fffhyy
+
++++= βββ L
其中。最常用的是p＝3，
∫
∏
≠
=
+?
+
=
1
0
0
p
il
l
pi
ds
li
ls
β
]9375955[
24
3211+
+?+=
nnnnnn
ffff
h
yy
取k=0,j=1可得到Adams隐式公式
)(
*
1
*
1
*
01 pnppnpnpnn
fffhyy
++++= βββ L
再用n＋1代替n，得到
)(
1
**
11
*
01 +?++
++++=
pnppnpnpnn
fffhyy βββ L
其中
∫
∏
≠
=
+?
+
=
0
1
0
*
p
il
l
pi
ds
li
ls
β
2004-12-13 68
基于数值积分的构造方法（续）
最常用的是p＝3的情况：
]5199[
24
2111++
+?++=
nnnnnn
ffff
h
yy
取k＝1，j=1，得到Nystr?m显式公式：
=
+?
+
=
++++=
∫
∏
≠
=
+
1
1
0
11011
,,2,1,0,
)(
p
il
l
pi
pnppnpnpnn
pids
li
ls
fffhyy
L
L
β
βββ
取k=0，j=2，得到Milne隐式公式（用n＋1代替n）：
=
+?
+
=
++++=
∫
∏
≠
=
+?+?+
0
2
0
111011
,,2,1,0,
)(
p
il
l
pi
pnppnpnpnn
pids
li
ls
fffhyy
L
L
β
βββ
2004-12-13 69
基于数值积分的构造方法（续）
∏
=
+
+
=
p
i
in
p
p
xtyf
0
)1(
)!1(
1
)())(,( ξξ
)())(,( tLtytf
p
由
),max(),min(
npn
xtxt <<
ξ
于是线性多步法
==
+=
∫
∑
+
=
+
),2,1,0()(
1
0
pidttl
h
fhyy
kn
jn
x
x
ipi
p
i
inpijnkn
Lβ
β
的局部截断误差为：
2004-12-13 70
基于数值积分的构造方法（续）
knknkn
yxyR
+++
= )(
∫∫
+
+
+=
jn
jn
kn
n
x
x
p
x
x
jnjn
dttLydttytfxy )())(,()(
j－
∫
+
=
kn
jn
x
x
p
dttLtytf ))())(,((
dtxtyf
p
R
in
x
x
p
i
p
kn
kn
jn
)())(,(
)!1(
1
0
)1(
=
+
+
Π
+
=
∫
+
ξξ
即
∫
=
+
+
+Π
+
=
k
j
p
i
p
p
dsissysf
p
h
)()))((),((
)!1(
0
)1(
2
ξξ
2004-12-13 71
基于数值积分的构造方法（续）
对于Adams显式公式与隐式公式，由于)(
0
is
p
i
+Π
=
积分中值定理在[-1，0] (j=1，k=0，隐式公式)上恒负，故根据在[0，1]（j=0，k=1，显式公式）上恒正，
∫
+Π
+
=
=
+
+
+
1
0 0
)1(
2
,
)())(,(
)!1(
dsisyf
p
h
R
p
i
p
p
Ekn
ηη
∫
+Π
+
=
=
+
+
1
0 0
)2(
2
)()(
)!1(
dsisy
p
h
p
i
p
p
η
∫
=
+
+
+
+Π
+
=
0
1 0
)2(
2
,
)()(
)!1(
dsisy
p
h
R
p
i
p
p
Ikn
η
为某中间点，E表示Explicit，I表示Implicit。
η
2004-12-13 72
基于数值积分的构造方法（续）
当p=3时,
)()(
720
251
)(
720
251
6)5(5)5(5
,1
hOxyhyhR
nEn
+==
+
η
)()(
720
19
)(
720
19
6)5(5)5(5
,
hOxyhyhR
nIkn
+?=?=
+
η
由局部截断误差的分析可知，在f(x,y(x))具有p+1
阶连续偏导数的条件下，p+1步Adams显式方法与p步Adams隐式方法的局部截断误差与h
p+2
同阶，
即它们都是p＋1阶方法。特别地，当p＝3时，
Adams显、隐方法都是四阶的。
2004-12-13 73
三、基于Taylor 展开的构造方法
Taylor展开法更具一般性。下面构造四步方其为四阶方法，并求其局部截断误差的主项。
法，即k＝4。用Taylor展开法构造下述公式，使
33221104 ++++
+++=
nnnnn
yyyyy αααα
)(
443322110 ++++
+++++
nnnnn
fffffh βββββ
)3,2,1,(),( +++== nnnnixyy
ii
在上式右边近似认为
。)4,3,2,1,(),())(,( ++++=
′
== nnnnnixyxyxff
iiii
以及
)3,2,01()( =
+
ixy
in
)4,3,2,01()( =
′
+
ixy
in
n
x
然后将函数及在点按 的幂展开：
h
2004-12-13 74
基于Taylor 展开的构造方法（续）
])(
!5
)(
!4
)(
!3
)(
2
)()([)(
)5(
5
)4(
432
104
L+++
′′′
+
′′
+
′
++=
+ nnnnnnnn
xy
h
xy
h
xy
h
xy
h
xyhxyxyy
！
αα
])(
!4
)(
!3
)(
2
)()([)(
)5(
4
)4(
32
10
L+++
′′′
+
′′
+
′
+
′
+
nnnnnn
xy
h
xy
h
xy
h
xyhxyhxyh
！
ββ
])(
!5
)3(
)(
!4
)3(
)(
!3
)3(
)(
2
)3(
)(3)([
)5(
5
)4(
432
3
L
nnnnnn
xy
h
xy
h
xy
h
xy
h
xyhxy +++
′′′
+
′′
+
′
++
！
α
])(
!5
)2(
)(
!4
)2(
)(
!3
)2(
)(
2
)2(
)(2)([
)5(
5
)4(
432
2
L+++
′′′
+
′′
+
′
++
nnnnnn
xy
h
xy
h
xy
h
xy
h
xyhxy
！
α
])(
!4
)4(
)(
!3
)4(
)(
2
)4(
)(4)([
)5(
4
)4(
32
4
L+++
′′′
+
′′
+
′
+
nnnnn
xy
h
xy
h
xy
h
xyhxyh
！
β
])(
!4
)3(
)(
!3
)3(
)(
2
)3(
)(3)([
)5(
4
)4(
32
3
L+++
′′′
+
′′
+
′
+
nnnnn
xy
h
xy
h
xy
h
xyhxyh
！
β
])(
!4
)2(
)(
!3
)2(
)(
2
)2(
)(2)([
)5(
4
)4(
32
2
L+++
′′′
+
′′
+
′
+
nnnnn
xy
h
xy
h
xy
h
xyhxyh
！
β
2004-12-13 75
基于Taylor 展开的构造方法（续）
)()(
3210 n
xyαααα +++=
hxy
n
)()32(
43210321
′
++++++++ βββββααα
2
4321321
)()432
2
9
2
2
1
( hxy
n
′′
+++++++ ββββααα
3)3(
432
1
321
)()8
2
9
2
22
9
3
4
6
1
( hxy
n
βββ
β
ααα +++++++
(1)
4)4(
4321321
)()
3
32
2
9
3
4
6
1
8
27
3
2
24
1
( hxy
n
ββββααα +++++++
L++++++++
5)5(
4321321
)()
3
32
8
27
3
2
24
1
40
81
15
4
120
1
( hxy
n
ββββααα
2004-12-13 76
基于Taylor 展开的构造方法（续）
另一方面有
)()4(
!3
1
)()4(
!2
1
)('4)()(
32
4 nnnnn
xyhxyhxhyxyxy
′′′
+
′′
++=
+
（2）
L+++ )()4(
!5
1
)()4(
!4
1
)5(5)4(4
nn
xyhxyh
比较（1），（2）两式中h的同次幂，由于
)(xy
的任意性，可得：
2004-12-13 77
基于Taylor 展开的构造方法（续）
0
h
1
3210
=+++ αααα
432
43210321
=+++++++ βββββααα
1
h
8432
2
9
2
2
1
4321321
=++++++ ββββααα
2
h
3
32
8
2
9
2
22
9
3
4
6
1
432
1
321
=++++++ βββ
β
ααα
3
h
3
32
3
32
2
9
3
4
6
1
8
27
3
2
24
1
4321321
=++++++ ββββααα
4
h
15
128
3
32
8
27
3
2
24
1
40
81
15
4
120
1
4321321
=++++++ ββββααα
5
h
2004-12-13 78
基于Taylor 展开的构造方法（续）
如果要求前六个方程成立，得到的公式的局部截断误差为。也可以只要求前面的几个方程成立。
)(
6
hO
1、Simpson公式要求前四个方程成立。令,，
0
310
=== ααα
0
10
==ββ
.
3
1
,
3
4
,
3
1
,1
4322
==== βββα
]4[
3
23424 +++++
+++=
nnnnn
fff
h
yy
]4[
3
122 nnnnn
fff
h
yy +++=
+++
可得，
公式为即
2004-12-13 79
基于Taylor 展开的构造方法（续）
第五个方程恰好也成立。
将得
)1()2(?
44,4
)(
+++
=
nnhn
yxyR
)()()]
3
32
8
27
3
2
24
1
40
81
15
4
120
1
(
!5
4
[
65)5(
4321321
5
hOhxy
n
+++++++?= ββββααα
)()(
90
1
65)5(
hOhxy
n
+?=
2004-12-13 80
基于Taylor 展开的构造方法（续）
2.隐式的Adams公式若要求前五个方程成立，令0,0
0210
==== βααα
.
24
9
,
24
19
,
24
5
,
24
1
,1
43213
==?=== ββββα
得隐式的Adams公式为
]5199[
24
12334 nnnnnn
ffff
h
yy +?++=
+++++
将得：
)1()2(?
44,4
)(
+++
=
nnhn
yxyR
)()(
720
19
6)5(5
hOxyh
n
+
=
2004-12-13 81
基于Taylor 展开的构造方法（续）
3.显式的Adams公式要求前五个方程成立，令
0,0
4210
==== βααα
24
55
,
24
59
,
24
37
,
24
9
,1
32103
=?==?== ββββα
得显式Adams公式：
)9375955(
24
12334 nnnnnn
ffff
h
yy?+?+=
+++++
解得：
44,4
)(
+++
=
nnhn
yxyR
21321
5
3
2
24
1
40
81
15
4
120
1
(
!5
4
[ ββααα ++++?=
将得：
)1()2(?
)()(
720
251
65)5(
hOhxy
n
+=
)()()]
3
32
8
27
65)5(
43
hOhxy
n
+++ ββ
2004-12-13 82
基于Taylor 展开的构造方法（续）
4.Milne公式（四步四阶）
要求前四个方程成立，令
0,0
40210
===== ββααα
此时第五方程恰好成立。
得Milne公式：
)22(
3
4
1234 ++++
+?+=
nnnnn
fff
h
yy
44,4
)(
+++
=
nnhn
yxyR
)()()]
3
32
8
27
3
2
24
1
40
81
15
4
120
1
(
!5
4
[
65)5(
43
21321
5
hOhxy
n
+++
++++?=
ββ
ββααα
)()(
45
14
6)5(5
hOxyh
n
+=
解得
3
4
,
8
3
,1
2310
==== βββα
将得
)1()2(?
2004-12-13 83
基于Taylor 展开的构造方法（续）
5,Hamming公式要求前五个方程成立,令
0,0
1020
==== ββαα
可得Hamming公式：
解得
8
3
,
8
6
,
8
3
,
8
9
,
8
1
43232
==?=== βββαα
)2(
8
3
)9(
8
1
234134 ++++++
++?=
nnnnnn
fff
h
yyy
)2(
8
3
)9(
8
1
12323 +++++
++?=
nnnnnn
fff
h
yyy
将得：
)1()2(?
44,4
)(
+++
=
nnhn
yxyR
321
5
40
81
15
4
120
1
(
!5
4
[ ααα ++?=
)()(
40
1
65)5(
hOhxy
n
+?=)()()]
3
32
8
27
3
2
24
1
65)5(
4321
hOhxy
n
+++++ ββββ
2004-12-13 84
基于Taylor 展开的构造方法（续）
Remark1：用Taylor展开法比用数值积分法推导线性多步法更加灵活，推导出的方法也比积分法的更多。用数值积分法可以得到的方法，用Taylor
展开法都可以得到，但用Taylor展开法得到的方法，用数值积分法却不一定能够得到。
Back
Remark2：应用线性多步法求解初值问题时，开头几点处的数值要用别的方法先计算出来。一般选用与多步法同阶或更高阶的单步法，如Runge-Kutta法、
Taylor展开法等。
Remark3：对线性隐式多步方法，除开头几点的函数值需要单独计算外，还需要迭代求解或采用预估-校正方法求解。
2004-12-13 85
数值算例求解常微分方程初值问题
=
=
1)0(
2
y
y
x
y
dx
dy
在[0,1]上的解，取步长h＝0.1。
计算结果如下图所示：
2004-12-13 86
2004-12-13 87
第九章常微分方程初值问题的数值解法
§9.1 引言
§9.2 Euler方法
§9.3 Runge –Kutta公式
§9.4单步法的进一步讨论
§9.5线性多步法数值算例
2004-12-13 2
§9.1 引言一、问题本章着重讨论一阶常微分方程初值问题
=
==
′
00
yxy
yxf
dx
dy
y
）（
），（
在区间[a,b]上的数值解法。
这些问题多数情况下求不出解析解，只能用近似的方法求解。常用的近似方法有两类。一类称为近似解析法，如级数解法，逐次逼近法等。另一类称为数值解法，它可以给出解在一些离散点上的近似值。
2004-12-13 3
二、解的存在唯一性若f（x,y）在区域D=,上连续，且关于y满足李普希兹（Lipschitz）条件,即存在常数L,使
{ }Rybxa ∈≤≤，
，
2121
),(),( yyLyxfyxf?≤?
G
21
,yy
对内任两个均成立，其中L是与x,y无关的常数,则上面的初值问题存在唯一解，且解是连续可微的。
Remark：在f(x,y)对y可微的情况下，若偏导数有界，
则可取，此时有之间。与介于
21212121
y,)(
),(
),(),( yyyLyy
y
xf
yxfyxf ξ
ξ
≤?
=?
y
yxf
L
Dyx
=
∈
),(
max
),(
此时Lipschitz连续条件显然成立。这是验证该条件的最简便的方法。
2004-12-13 4
三、解的适定性解的适定性是指解的存在唯一性以及数值稳定性。
此处主要是指解对于右端项以及初值的扰动的适应性。关于适定性有如下的结论：
定理：若f（x,y）在区域D= 上满足Lipschitz连续条件，则初值问题是适定的。
{ }Rybxa ∈≤≤，
四、等价的积分方程若y(x)是初值问题的解，对方程两边同时积分，
利用初始条件可得：
∫
+=
x
x
dttytfxyxy
0
))(,()()(
0
该方程为与初值问题同解的积分方程，我们可以从积分方程出发去构造初值问题的求解公式。
2004-12-13 5
五、数值解法常微分方程初值问题的数值解是求上述初值问题的解y(x)在存在区间[a，b]中的点列
），，（nihxx
iii
L10
1
=+=
上的近似值。称为步长，i
y
i
h
一般情况下我们取等步长，记为h。
)(
n
xy
初值问题的解析解（理论解）用表示,数值解法的精确解用表示，并记f
n
=f(x
n
,y
n
),
而。
n
y
))(,()(
nnn
xyxfxy =
′
求初值问题的数值解一般是逐步进行的，即计算出y
n
之后计算y
n+1
。
2004-12-13 6
数值解法的分类数值解法一般分为：
（1）单步法：在计算时，只用到,
和，即前一步的值。
n
x
1+n
y
n
y
1+n
x
（2）多步法：计算时，除用到,和以外，还要用和,即前
k步的值。
1+n
y
pn
y
1+n
x
n
x
n
y
pn
x
)0;,2,1( >= kkp L
单步法和多步法都有显式和隐式方法之分。显式和隐式的单步法可以分别写成：
），（hyxhyy
nnnn
,
1
φ+=
+
），（hyyxhyy
nnnnn
,,
11 ++
+= φ
对多步法来说，显式和隐式方法具有相同的意义。
Back
2004-12-13 7
§9.2 Euler方法一、显式欧拉（Euler）方法设节点为，欧拉方法的计算公式
），，（nnnhxx
n
,,210
0
L=+=
），，，）（，（L210
1
=+=
+
nyxhfyy
nnnn
这是一种最简单的显式单步方法，该方法可以通过不同的途径获得。
),2,1,0)(,(
1
L==
+
nyxf
h
yy
nn
nn
1、差商法用两点差商公式代替导数,再用表示的近似值,则得到
nn
nn
xx
xyxy
+
+
1
1
)()(
n
y
)(
n
xy
)(
n
xy
′
2004-12-13 8
显式欧拉（Euler）方法（续）
假设在附近把y（x）展成Tayler级数
n
x
L+++=+ )(''
2
)(')()(
2
nnnn
xy
h
xhyxyhxy
取h的线性部分,并用表示的近似值,得
n
y )(
n
xy
),2,1,0)(,(
1
L=+=
+
nyxhfyy
nnnn
2、Taylor展开法
3、数值积分法
∫
+
+=
+
1
))(,()()(
1
n
n
x
x
nn
dttytfxyxy
对微分方程两端从x
n
到x
n+1
积分，得等价的积分方程对右端的积分部分采用左矩形公式近似即得Euler公式。
Remark：Taylor展开法与数值积分法是构造微分方程数值解的两类主要的方法。
2004-12-13 9
欧拉（Euler）方法的几何意义
Euler方法有明显的几何意义。
如右图所示，一阶常微分方程初值问题的解曲线y(x)过点P
0
(x
0
,y
0
)。从P
0
出发以
f(x
0
,y
0
)为斜率作一直线段，
与x=x
1
相交于点P
1
(x
1
,y
1
)，
显然有y
1
=y
0
+hf(x
0
,y
0
)。
同理，再由P
1
出发以f(x
1
,y
1
)为斜率作一直线段，与x=x
2
相交于点P
2
(x
2
,y
2
)，显然有y
2
=y
1
+hf(x
1
,y
1
)。这样一直做下去，得到一条折线P
0
P
1
P
2
……，作为y(x)的近似曲线。
因此，Euler方法又称为Euler折线法。
2004-12-13 10
二、隐式Euler方法和Euler方法的改进
dttytfxyxy
n
n
x
x
nn
），（）（）
∫
+
+=
+
1
)((
1
中将积分用右矩形公式
)),((,)(
2
)('
2
baab
f
bfabdxxf
b
a
∈≈
∫
η
η
）（）（）（
代入，有
K,2,1,0)),(,()()(
111
=+≈
+++
nxyxhfxyxy
nnnn
若在等价积分方程
K,2,1,0),,(
111
=+≈
+++
nyxhfyy
nnnn
从而得到上式是一个隐式的单步方法，称为隐式欧拉法或后退的欧拉法。利用此方法，每一步都要把上式作为y
n+1
的一个方程来求解。从数值积分的误差分析，很难期望隐式欧拉法比显式欧拉法更精确。
2004-12-13 11
隐式Euler方法和Euler方法的改进（续）
为了得到更精确的方法，在等价的积分方程中用梯形公式
),(),(''
12
][
2
3
baf
ab
bfaf
ab
dxxf
b
a
∈
+
=
∫
ηη
）（
）（）（）（
L,2,1,0)],,(),([
2
111
=++=
+++
nyxfyxf
h
yy
nnnnnn
近似积分项，再分别用y
n
，y
n+1
代替y(x
n
)和y(x
n＋1
)，即可得到梯形方法：
该方法也是一种隐式的单步方法。对该方法，
)],(),([
2
1
),,,(
111 +++
+=
nnnnnnn
yxfyxfhyyxφ
从n＝0开始计算，每步都要求解y
n+1
的一个方程。一般来说，这是一非线性方程，可迭代计算如下：
2004-12-13 12
隐式Euler方法和Euler方法的改进（续）
L,2,1,0,
)],(),([
2
),(
][
11
]1[
1
]0[
1
=
++=
+=
++
+
+
+
s
yxfyxf
h
yy
yxhfyy
s
nnnnn
s
n
nnnn
使用上式时，先用第一式计算出y
n+1
的近似值，再用第二式反复进行迭代，得到数列，用来控制是否继续进行迭代，其中ε为允许误差。把满足要求的作为y(x
n+1
)的近似值y
n+1
，类似地可得出y
n+2
，
y
n+3
，…。
]0[
1+n
y
{}
∞
=
+
0
][
1
s
s
n
y
ε≤?
+
+
+
][
1
]1[
1
s
n
s
n
yy
]1[
1
+
+
s
n
y
当f(x,y)关于y满足Lipschitz条件时，且步长h满足时，上述迭代过程是收敛的。这是因为：
1
2
1
<hL
2004-12-13 13
]0[
11
][
11
][
11
][
11
][
1111
]1[
11
2
)(
)(
2
)
2
++++++++
++++
+
++
<<?<?≤?
=
=?
nn
s
nn
s
nn
s
nn
s
nnnn
s
nn
yyyyyy
hL
yy
y
fh
yxfyxf
h
yy
L
ξ
，（），（
)(0)
2
(
]0[
11
1]1[
11
∞→→?≤≤?
++
++
++
syy
hL
yy
nn
ss
nn
L
实用中，h取得较小时,梯形方法第二式只迭代一次就结束，得到Euler预估—校正方法(改进的Euler方法)：
++=
+=
+++
+
)],(),([
2
),(
]0[
111
]0[
1
nnnnnn
nnnn
yxfyxf
h
yy
yxhfyy
其中第一式称为预估算式，第二式称为校正算式。
隐式Euler方法和Euler方法的改进（续）
其中ξ介于与之间。
1+n
y
][
1
s
n
y
+
1
][
1
,
++
→∞
n
s
n
yy
即为s趋于。
2004-12-13 14
隐式Euler方法和Euler方法的改进（续）
若将Euler预-校方法中的第一式带入第二式，得
))],(,(),([
2
11 nnnnnnnn
yxhfyxfyxf
h
yy +++=
++
Remark：这是一种显示的单步方法。有时为了计算方便，常将上式改写成：
++=
=
++=
+
),(
),(
)(
2
12
1
211
hKyhxfK
yxfK
KK
h
yy
nn
nn
nn
2004-12-13 15
三、单步法的局部截断误差和阶设一般的单步法为：
显式公式：
），，（hyxhyy
nnnn
φ+=
+1
隐式公式：
），，，（hyyxhyy
nnnnn 11 ++
+= φ
设为数值方法的精确值，y（）为微分方程的精确解。
n
x
n
y
nnn
yxye?=）（
n
x
定义1：为某一数值方法在处的整体截断误差。
Remark：整体截断误差不仅与这步的计算有关，而且与前面所有点的计算的误差累计有关。为了简化误差分析，我们着重分析计算中的某一步。对一般的显式单步法，有如下定义：
n
x
2004-12-13 16
单步法的局部截断误差和阶（续）
定义2：对单步法，在的假设下，称为在处的局部截断误差。
）（
nn
xyy =
)),(,()()(
1,
hxyxhxyxyR
nnnnhn
φ=
+
1+n
x
）（
nn
xyy =
若设，即第n步及以前各步都没有误差，
则由显示单步法计算一步所得之与之差为：
1+n
y
)(
1+n
xy
)),(,()()(
)],,([)()(
1
111
hxyxhxyxy
hyxhyxyyxy
nnnn
nnnnnn
φ
φ
=
+?=?
+
+++
即在的假设下，。
11,
)(
++
=
nnhn
yxyR
）（
nn
xyy =
这就是上面定义中称R
n,h
为“局部”的含义，我们一般用该式作为定义。这里应该注意，R
n,h
和整体截断误差
e
n+1
是不同的。
2004-12-13 17
单步法的局部截断误差和阶（续）
定义3：若一个单步方法的局部截断误差为O(h
p+1
)，即则称该方法为p阶方法（其中p为正整数）。
)()),(,()()(
1+
=+
p
hOhxyxhxyhxy φ
Remark：由前面的定义可知，若某个单步方法是一种
p阶方法，则有R
n,h
=O(h
p+1
)，即p阶方法的局部截断误差为h的p＋1阶。我们往往比较关心R
n,h
按h展开式的第一项。
定义4：若一个单步方法是一种p阶方法，其局部截断误差可以写成：
)())(,(
21
,
++
+=
pp
nnhn
hOhxyxR ψ
则Ψ(x
n
,y(x
n
))h
p+1
称为方法的主局部截断误差，或局部截断误差的主项。
2004-12-13 18
单步法的局部截断误差和阶（续）
例1：求Euler方法的局部截断误差。
由Taylor展式，有
)()('''
!3
)(''
!2
)()()(
))(,()()()),(,()()(
2
32
hOxy
h
xy
h
xyhxyhxy
xyxhfxyhxyhxyxhxyhxy
=++=
′
+=
+=+
L
φ
故Euler方法是一阶方法,局部截断误差为：
)()(''
2
3
2
,
hOxy
h
R
nhn
+=
主局部截断误差为：。
)(''
2
2
n
xy
h
2004-12-13 19
单步法的局部截断误差和阶（续）
例2：求Euler预估-校正方法
++=
+=
+++
+
)],(),([
2
),(
]0[
111
]0[
1
nnnnnn
nnnn
yxfyxf
h
yy
yxhfyy
的局部截断误差。
改写上述公式为：
++=
=
++=
+
),(
),(
)(
2
12
1
211
hkyhxfk
yxfk
kk
h
yy
nn
nn
nn
)('))(,(),(
1 nnnnn
xyxyxfyxfk ===
2004-12-13 20
单步法的局部截断误差和阶（续）
),(
12
hkyhxfk
nn
++=
L+
+
+
+
+
+=
)],(),(2),([
2
1
)],(),([),(
2
2
2
1
2
2
1
2
2
2
2
1
nnnnnn
nnnnnn
yx
y
f
khyx
yx
f
khyx
x
f
h
yx
y
f
hkyx
x
f
hyxf
)(),(),(),(),(
2
hOyxfyxhfyxhfyxf
nnynnnnxnn
+++=
)(][),(
2
))(,(
hOfffhyxf
nn
xyx
yxnn
+++=
)()('')('
2
hOxhyxy
nn
++=
2004-12-13 21
单步法的局部截断误差和阶（续）
代入，得
)()]('')(')('[
3
2
1
hOxhyxyxyyy
nnn
h
nn
++++=
+
而
L++
′
+=
+
)(''
2
)()(
2
1 nnnn
xy
h
xyhxyxy（）
因此有
)()(
3
11,
hOyxyR
nnhn
=?=
++
故Euler预估-校正方法为二阶方法。
2004-12-13 22
单步法的局部截断误差和阶（续）
例3：求梯形公式的局部截断误差。
))(,(''
12
))](())(([
2
)()(
3
111
εε yf
h
xyxfxyxf
h
xyxy
nnnnnn
++=
+++
，，
由微分中值定理
))((),())(,(
11),(1111
1
++++++
=?
+
nnxnnnn
yxy
y
f
yxfxyxf
n
η
代入上式得
)(
12
))()(
2
1(
3
11),(
1
ε
η
y
h
yxy
y
fh
nnx
n
′′′
=?
++
+
)]()([
2
111 +++
++=
nnnnnn
yxfyxf
h
yy，，
上面两式相减，有
)(
12
)]())(([
2
)(
3
11111
εy
h
yxfxyxf
h
yxy
nnnnnn
′′′
=?
+++++
，，
2004-12-13 23
单步法的局部截断误差和阶（续）
故
)()]('''
12
][
2
1[)(
3
3
),(11
1
hOy
h
y
fh
yxy
n
xnn
=?+
+=?
+
++
ε
η
L
1
2
),(
1
<
+
η
n
x
y
fh
当h充分小时，若，则可以选择h，s.t,
,从而有
1<≤
L
y
f
L+
+=
+
+
),(
1
1
2
1
2
1
1
η
η
n
n
x
x
y
fh
y
fh
），（
即梯形公式为二阶公式。
Back
2004-12-13 24
§9.3 Runge –Kutta公式一、Taylor 方法
Taylor展开法与数值积分法是推导高阶方法的常用手段。
本节以Taylor展开法为基础，介绍如何推导高阶单步法。
假定初值问题的解y(x)以及函数f(x,y)是足够光滑的，有
)()(
!
)(''
!2
)(')(
1)(
2
1
+
+
+++++=
p
n
p
p
nnnn
hOxy
p
h
xy
h
xhyxyxy L）（
)())(,(
!
))(,('
!2
))(,()(
1)1(
2
+?
+++++=
p
nn
p
p
nnnnn
hOxyxf
p
h
xyxf
h
xyxhfxy L
)(
2
1
!!2
p
n
p
nnnn
y
p
h
y
h
yhyy ++
′′
+
′
+=
+
L
当h充分小时，略去余项，将y
(k)
(x
n
)用来代替，
则有p阶Taylor方法：
)(
1+p
hO
)(k
n
y
)(k
n
y
其中(k=1,2,…,p)根据求导法则，其计算公式为：
2004-12-13 25
Taylor 方法（续）
),('
nnn
yxfy =
),(),(),(''
nnynnnnxn
yxfyxfyxfy +=
M
),](2['''
22
nnyyxyyxyxxn
yxfffffffffy ++++=
Remark1：显然，p阶Taylor方法的局部截断误差为
。当p=1时，
Taylor方法就是Euler方法。当p≥2时，需要计算公式中的高阶导数。
)()(
)!1(
1
2)1(1
,
+++
+
+
=
p
n
pp
hn
hOxyh
p
R
Remark2：显然，Taylor方法可以得到任意阶精度的方法。但在实际计算中，Taylor方法往往相当困难，因为公式中的高阶导数会很复杂。故Taylor方法很少单独使用，但常用它来启发思路。
2004-12-13 26
二、Runge—Kutta方法
Euler公式：
）（
），（
2
,
1
11
hOR
yxfk
hkyy
hn
nn
nn
=
=
+=
+
改进的Euler预估校正格式：
)(
)
2
3
,
12
1
211
hOR
hkyhxfk
yxfk
kk
h
yy
hn
nn
nn
nn
=
++=
=
++=
+
），（
），（
（
基本思想：用不同点的函数值作线性组合，构造高阶单步的近似公式，把近似公式和解的Taylor展开式比较，使前面尽可能多的项完全相同。这种方法间接应用Taylor展开的思想，避免了高阶导数计算的困难。
2004-12-13 27
一般的Runge-Kutta方法的形式为
=
++=
=
=
+=
∑
∑
=
=
+
riKhyhxfk
yxfk
Kchyx
hyxhyy
i
j
jijnini
nn
r
i
iinn
nnnn
,,3,2,
),,(
),,(
1
1
1
1
1
Lβα
φ
φ
，
），（
其中，为常数，选取这些常数的原则是，
要求第一式的右端在处泰勒展开后，按h的幂次重新整理，得到
iiji
C，，βα
），
nn
yx(
Runge—Kutta方法（续）
L++++=
+
3
3
2
211
3
1
2
1
hVhVhVyy
nn
！！
2004-12-13 28
Runge—Kutta方法（续）
L+
′′′
+
′′
+
′
+=
+
)(
3
)(
2
)()()(
32
1 nnnnn
xy
h
xy
h
xyhxyxy
！！
有若干项重合，即要求与微分方程的解y(x
n+1
)在x
n
处的展开式
L，，，)()()(
321 nnn
xyVxyVxyV
′′′
=
′′
=
′
=
上述公式叫做r级Runge-Kutta方法(c
r
≠0)。
当r＝1时，上式中φ(x,y,h)=c
1
f(x,y)。和泰勒展开式中
p=1时的φ(x,y,h)=f(x,y)比较知，取c
1
＝1即可。这就是
Euler方法。
下面用Taylor展开方法对r(≥2)级的Runge-Kutta方法确定参数，使其成为p (≥2)阶的方法。
iiji
C，，βα
2004-12-13 29
Runge—Kutta方法（续）
下面以二级Runge—Kutta公式为例进行具体推导。
++=
=
+=
+=
+
），（
），（
12122
1
2211
1
),,(
),,(
KhyhxfK
yxfK
KcKchyx
hyxhyy
nn
nn
nn
nnnn
βα
φ
φ
对要求适当选取系数，使当时，
上式的局部截断误差为，即成为二阶方法。
21221
,,,βαcc）（
nn
xyy =
）
3
(hO
做法一：按照定义，根据局部截断误差的定义，将上式代入，并按h展开，得：
))](,()([)()(
),,()()(
21221
xyhyhxfcxychxyhxy
hyxhxyhxy
′
+++
′
+=
+
βα
φ
2004-12-13 30
Runge—Kutta方法（续）
L+
′′′
+
′′
+
′
=+ )(
3
1
)(
2
1
)()()(
32
xyhxyhxyhxyhxy
！！
＋
)(]))(()(2)[(
!2
1
])([
!1
1
),())(,(
32
21212
2
2
212212
hOfxyhfxyhhfh
fxyhhfyxfxyhyhxf
yyxyxx
yx
+
′
+
′
++
′
++=
′
++
ββαα
βαβα
将上两式代入，利用Taylor展式，并按h的幂次整理，得：
)()](
6
1
)
2
1
6
1
(
)
3
1
()
2
1
6
1
[(
])
2
1
()
2
1
[()()1(
),,()()(
422
212
2122
2
22
3
21222
2
21
hOffffffc
ffcfch
ffcfchxyhcc
hyxhxyhxy
yxyyy
xyxx
yx
+++?
+?+?+
+?+
′
=
+
β
βαα
βα
φ
2004-12-13 31
Runge—Kutta方法（续）
==
=+
2
1
1
21222
21
βα cc
cc
要使，则需
)(
3
,
hOR
hn
=
上述方程有无穷多组解，所确定的解都能使二级
Runge-Kutta方法成为一个二阶方法。
由于二级方法，方程也可写为
0
2
≠c
==
=
2
212
21
2
1
1
c
cc
βα
在上述方程中，选取不同的c
2
，即可获得不同的二级二阶Runge-Kutta方法。
2004-12-13 32
Runge—Kutta方法（续）
)()],()),((
),(),(2),()[(
2
1
)],(),(),([
!1
1
),(
)),(,(
32
21
212
2
2
212
2122
hOyxfyxfh
yxfyxfhhyxfh
yxfyxfhyxhfyxf
yxfhyhxfK
nnyynn
nnxynnnnxx
nnynnnnxnn
nnnn
+
+?+
++=
++=
β
βαα
βα
βα
！
＋
做法二：按照定义，将K
2
在(x
n
,y
n
)处展开，有将K
1
,K
2
代入y
n+1
,有，
)()],(),(),(),(2),([
2
)],(),(),([)()()(
),(
422
21212
2
2
3
2
212
2
221
2211
hOyxfyxfyxfyxfyxf
hc
yxfyxfyxfhcxyhccxy
Khcyxfhcyy
nnyynnnnxynnnnxx
nnynnnnxnn
nnnn
++++
++
′
++=
++=
+
ββαα
βα
2004-12-13 33
Runge—Kutta方法（续）
),()(' yxfxy =
),(),(),()('' yxfyxfyxfxy
yx
+=
由
M
),](2[)('''
22
yxfffffffffxy
yyxyyxyxx
++++=
)(),)]((
6
1
)
2
1
6
1
(
)
3
1
()
2
1
6
1
[(),]()
2
1
(
)
2
1
[()()1()(
422
212
2122
2
22
3
212
22
2
2111,
hOyxffffffc
ffcfchyxffc
fchxyhccyxyR
nnyxyyy
xyxxnny
xnnnhn
+++?+
+
+?+
′
=?=
++
β
βααβ
α
＋
故同第一种做法一样，可以获得无穷多组二级二阶
Runge-Kutta方法。常用的二阶Runge-Kutta方法有：
2004-12-13 34
Runge—Kutta方法（续）
1,
2
1
21
221
====
β
αc
++=
=
++=
+
),(
),(
)(
2
12
1
211
hKyhxfK
yxfK
KK
h
yy
nn
nn
nn
++=
=
+=
+
)
2
,
2
(
),(
12
1
21
K
h
y
h
xfK
yxfK
hKyy
nn
nn
nn
2
1
,1,0
21
221
====
β
αcc
②
c
①
Euler预估-校正格式中（间）点公式
2004-12-13 35
Runge—Kutta方法（续）
3
2
,
4
3
,
4
1
21221
==== βαcc
③
++=
=
++=
+
)
3
2
,
3
2
(
),(
)3(
4
12
1
211
hKyhxfK
yxfK
KK
h
yy
nn
nn
nn
此时可以论证：
,0
2
1
6
1
2
2
2
=?
cα
0
3
1
2
21
2
=?
c
β
α
二阶Heun（休恩）方法此时，即二阶Heun
方法是局部截断误差项数最少的方法且二级Runge-
Kutta方法不可能达到三阶。
)()(
6
4
),(
3
,
hOffff
h
R
nn
yxyxyhn
++=
2004-12-13 36
Runge—Kutta方法（续）
用类似的方法可以研究其它级的Runge-Kutta方法。
常用的其它级的Runge-Kutta方法有：
++=
++=
=
++=
+
)
3
2
,
3
2
(
)
3
1
,
3
(
),(
)3(
4
23
12
1
311
hKy
h
xfK
hKy
h
xfK
yxfK
KK
h
yy
nn
nn
nn
nn
+?+=
++=
=
+++=
+
)2,(
)
2
1
,
2
(
),(
)4(
6
213
12
1
3211
hKhKyhxfK
hKy
h
xfK
yxfK
KKK
h
yy
nn
nn
nn
nn
三阶Heun（休恩）公式三阶Kutta公式
2004-12-13 37
标准（经典）四阶Runge-Kutta公式
Runge—Kutta方法（续）
++=
++=
++=
=
++++=
+
),(
)
2
1
,
2
(
)
2
1
,
2
(
),(
)22(
6
34
23
12
1
43211
hKyhxfK
hKy
h
xfK
hKy
h
xfK
yxfK
KKKK
h
yy
nn
nn
nn
nn
nn
++?+=
+
++=
++=
=
+++?++=
+
))
2
2
1(
2
2
,(
)
2
22
2
12
,
2
(
)
2
1
,
2
(
),(
))22()22([
6
324
213
12
1
43211
hKhKyhxfK
hKhKy
h
xfK
hKy
h
xfK
yxfK
KKKK
h
yy
nn
nn
nn
nn
nn
基尔（Gill）公式
Back
2004-12-13 38
§9.4单步法的进一步讨论一、收敛性定义：
定义：对于初值问题，如果一个单步显式方法产生近似解对于任一固定的,均有则称该单步法是收敛的。
nhxxbxx +=∈
00
],,[ )(lim
0
xyy
n
h
=
→
Remark：该定义可类似地应用于单步隐式方法以及后面的线性多步法。从定义可知，若格式收敛，整体截断误差e
n
=y(x
n
)-y
n
必然趋于零。
定理（整体截断误差与局部截断误差的关系）：若初值问题的一个单步方法之局部截断误差为且单步法中函数关于y满足lipschitz条件，则
)1(),(
1
,
≥=
+
phOR
p
hn
),,( hyxΦ
)()(
111
p
nnn
hOyxye =?=
+++
2004-12-13 39
整体截断误差与局部截断误差的关系（续）
)],,([)(
~
)(
111,
hyxhyxyyxyR
nnnnnnhn
φ+?=?=
+++
证明：根据局部截断误差的定义，有
)()),(,()()(
1
1
+
+
==
p
nnnn
hOhxyxhxyxy φ
其中为在假定下计算的处的数值解。
1
~
+n
y
)(
nn
xyy =
1+n
x
即
)),(,()(
~
1
hxyxhxyy
nnnn
φ+=
+
1
11
~
)(
+
++
<?
p
nn
chyxy
则有常数c>0，使得由于且关于满足Lipschitz条件，得
),,,(
1
hyxhyy
nnnn
φ+=
+
),,( hyxφ
y
11
~
++
nn
yy
),,()),(,()( hyxhxyxhyxy
nnnnnn
φφ?+?≤
nn
yxyhL?+≤ )()1(
2004-12-13 40
整体截断误差与局部截断误差的关系（续）
1111
~~
)(
++++
+?≤
nnnn
yyyxy
111
)(
+++
=
nnn
yxye故
n
p
ehLCh )1(
1
++≤
+
从而有递推关系：
])1()[1(
1
11
1?
++
+
++++≤
n
pp
n
ehLchhLche
1
21
)1()]1(1[
+
++++=
n
p
ehLhLch
0
1
)1( ehL
n+
++
])1()1()1(1[
21 np
hLhLhLch +++++++≤≤
+
LK
0
1
1
1
)1(
1)1(
1)1(
ehL
hL
hL
ch
n
n
p +
+
+
++
+
+
=
2004-12-13 41
若则且
,)(
00
yxy =
0
0
=e
整体截断误差与局部截断误差的关系（续）
hL
ehLhLhL =+++≤+< L
2
)(
2
1
110
得：
hLnn
ehL
)1(1
)1(0
++
≤+<
[]1
1
)1(
)1(
1
1
=
≤
+
+
+
+
hLnP
hLn
p
n
eh
L
c
hL
e
che
故当固定时，1+
=
n
xx
abxxhn
n
≤?=+
+ 01
)1(
[ ]
pLabp
n
hceh
L
c
e
1
)(
1
1 =?≤
+
所以证毕
2004-12-13 42
整体截断误差与局部截断误差的关系（续）
Remark1：该定理说明，数值方法的整体截断误差比局部截断误差低一阶。收敛的方法至少是一阶方法。
Remark2：在该定义的条件下，欧拉方法是一阶方法，欧拉预估-校正方法是二阶方法。当f(x,y)
关于y满足利普希茨条件时，r级龙格-库塔方法中的φ关于y也满足利普希茨条件，所以定理中的条件得到满足，解的收敛性得到了保证。
2004-12-13 43
二.相容性设单步法为，
),,(
1
hyxhyy
nnnn
φ+=
+
)],,([)(
1
hyxhyxy
nnnn
φ+?=
+
)),(,()()(
1
hxyxhxyxy
nnnn
φ=
+
11,
)(
++
=
nnhn
yxyR
由于,且为任意点，故该式相当于用近似方程
）（00
,
→→ hR
hn
n
x
)),(,(
)()(
hxyx
h
xyhxy
φ≈
+
))(,( xyxf
dx
dy
=来替代通过处求解近似方程而获得原方程的近似解。因此，必须要求当时，近似方程应逼近于原微分方程。
n
xx =
0→h
2004-12-13 44
相容性（续）
因此要使时，近似方程的极限状态成为原微分方程，需且只需极限
0→h
)(
)()(
lim
0
xy
h
xyhxy
h
′
=
+
→
由于
))(,()),(,(lim
0
xyxfhxyx
h
=
→
φ
成立。
)),(,( hxyxφ
由于假设为连续函数，因而上式可以表示为
))(,()0),(,( xyxfxyx =φ
定义：如果当时，近似方程能逼近微分方程，则称数值公式与原微分方程相容。
0→h
),,(
1
hyxhyy
nnnn
φ+=
+
相容的充要条件为
))(,()0),(,( xyxfxyx =φ
2004-12-13 45
相容性（续）
Remark：可以证明，若单步法的阶大于或等于1，则单步法与微分方程相容；反之，如果单步法与微分方程相容，且关于h满足Lipschitz条件，则单步法至少为一阶方法。事实上：
),,( hyxφ
⑴若单步法的阶大于或等于1，由
0)),(,()()(
1,
→=
+
hxyxhxyxyR
nnnnhn
φ
（当）知
0→h
)),(,(lim
)()(
lim
00
hxyx
h
xyhxy
hh
φ
→→
=
+
)0),(,())(,()(' xyxxyxfxy φ==
故有即单步法与微分方程相容。
2004-12-13 46
相容性（续）
⑵若单步法与微分方程相容，且关于
h满足Lipschitz条件，则
),,( hyxφ
),,()()(
111,
hyxhyxyyxyR
nnnnnnhn
φ=?=
+++
)),(,()()(
1
hxyxhxyxy
nnnn
φ=
+
)()()()(
2
nnn
xyhOxyhxy?+
′
+=
)]()0),(,([ hOxyxh
nn
+? φ
)())(,()(
2
hOxyxhfxyh
nnn
+?
′
=
)(
2
hO=
即与微分方程相容的单步法至少为一阶方法。
2004-12-13 47
相容性（续）
通过上面的讨论可以得到下面的定理：
定理：设增量函数在区域中连续，并对变量y满足利普希茨条件，
则单步法收敛的充要条件为相容性条件成立。
),,( hyxφ
}0,,),{(
0
hhRybxayxD ≤≤∈≤≤=
Remark：在满足定理的条件下，Euler方法，
Euler预估-校正格式，Runge-Kutta方法等都与原微分方程相容。
2004-12-13 48
三、稳定性定义1 *用一个数值方法求解微分方程初值问题时，对给定的步长h>0，若在计算时引入误差（也称扰动），
由此引起计算后面的时的误差按绝对值均不增加，则称这个数值方法是绝对稳定的。
n
y
n
δ
),2,1( L=
+
ky
kn
关于单步法收敛性的概念和收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的。但实际计算时通常会有舍入误差及其累积。数值求解微分方程的过程是一个递推公式，必须考虑误差积累是否得到控制。
设f(x,y)关于y满足Lipschitz条件，我们总是针对模型方程进行讨论，其中为复常数。
λ
yy λ=
′
0Re <）（λ
为保证微分方程本身的稳定性，这里假定。
2004-12-13 49
稳定性（续）
定义1设步长为h>0的单步法用于求解时，
yy λ=
′
中由引起的误差满足，则称
n
δ
m
δ
)( nm
nm
>< δδ
单步法对于所用步长h和复数是绝对稳定的。
λ
若在计算时有误差,但在计算后面的
)( nmy
m
>
n
δ
n
y
Remark1：在上面的定义中，可以取小于或等于关系符。取小于号是为了和线性多步法相一致。
Remark2：从上面的定义可知，单步法是否稳定，
与模型方程中的复数λ以及所用步长h有关。若对复平面上的某个区域G，当时，单步法绝对稳定，则称G为单步法的绝对稳定区域，G与实轴的交集为绝对稳定区间。
Gh∈λ
2004-12-13 50
稳定性（续）
①Euler显式公式将Euler显式公式用于模型方程，则有
yy λ=
′
nnnn
yhyhyy )1(
1
λλ +=+=
+
下面考察几个常用公式的稳定性。
设有误差，
n
y
n
δ nnn
yy δ+=
~
参与运算的量为，
误差为，则实际得到近似的量为
1+n
y
1+n
δ
由此引起的要求误差不增加，即，nn
δδ <
+1
1+n
y
111
~
+++
+=
nnn
yy δ
，即。与前面公式
))(1(
11 nnnn
yhy δλδ ++=+
++
nn
h δλδ )1(
1
+=
+
相减，有。
必须。
11
1
<+=
+
λ
δ
δ
h
n
n
是保证绝对稳定性对11 <+ λh
步长h加的限制。当λ为实数时，可以得到用hλ表示的绝对稳定的区间(-2,0)。
2004-12-13 51
稳定性（续）
②隐式Euler公式后退的Euler公式用到上，有，
yy λ=
′
11 ++
+=
nnn
yhyy λ
故
11 ++
+=
nnn
hδλδδ
即
nn
h
δ
λ
δ
=
+
1
1
1
令
1
1
1
1
<
=
+
λδ
δ
h
n
n
11 >? λh
得绝对稳定区域为。
可见，若取λ为实数，则对于任意的h都是绝对稳定的。
2004-12-13 52
稳定性（续）
③梯形公式梯形公式用到上，有yy λ=
′
)(
2
11 ++
++=
nnnn
yy
h
yy
λ
故)(
2
11 ++
++=
nnnn
h
δδ
λ
δδ
即
nn
h
h
δ
λ
λ
δ
+
=
+
2
1
2
1
1
设
iyxh +=λ
2004-12-13 53
稳定性（续）
iyx
iyx
++
=
2
2
2
1
2
1
λ
λ
h
h
+
=
n
n
δ
δ
1+
由于
2
1
22
22
}
)2(
)2(
{
yx
yx
+?
++
=
2
1
22
22
}
44
44
{
yxx
yxx
++?
+++
=
当x=Re( )<0,上式右端总小于1，故梯形法的绝对稳
hλ
hλ
定区域为Re( )<0,即左半平面。
2004-12-13 54
稳定性（续）
④四阶经典R-K方法四阶经典R-K方法用到上，有
yy λ=
′
n
yk λ=
1
)
2
1
()
2
1
(
2
12
λλλ hyhkyk
nn
+=+=
)]
2
1
(
2
1
[)
2
1
(
2
23
hhyhkyk
nn
λλλλλ ++=+=
)]
4
1
2
1
([)(
322
34
λλλλλλ hhhyhkyk
nn
+++=+=
)}
4
1
2
1
()
4
1
2
1
(2
43322322
λλλλλλλ hhhhh +++++++
)22(
6
43211
kkkk
h
yy
nn
++++=
+
)
2
1
(2{
6
2
λλλ hy
h
y
nn
+++=
})(
24
1
)(
6
1
)(
2
1
1{
432
λλλλ hhhhy
n
++++=
因此
2004-12-13 55
稳定性（续）
扰动满足
})(
24
1
)(
6
1
)(
2
1
1{
432
1
λλλλδδ hhhh
nn
++++=
+
令,得四阶经典R－K方法的绝对稳定区域为
1
1
<
+
n
n
δ
δ
1)(
24
1
)(
6
1
)(
2
1
1
432
<++++ λλλλ hhhh
Back
2004-12-13 56
§9.5线性多步法一、线性多步法的一般理论步的信息来预测，则可以期望获得较高的精度。
kn
y
+
线性多步法的基本思想：如果充分利用前面多
)(
n
xy
记的近似值为,并记，
则k步线性多步方法一般形式为
nhxxy
nn
+=
0
,),(
nnn
yxff =
)(
110
0
knknn
k
j
jnj
fffhfh
++
=
+
+++==
∑
ββββ L
knknnnn
k
j
jnj
yyyyyy
++++
=
+
+++++=
∑
αααααα L
3322110
0
jj
βα,
00
,,0 βαα ≠
k
其中为常数，不全为零。
2004-12-13 57
线性多步法的一般理论（续）
因为上式两边可同乘一常数，故可规定。
1=
k
α
则上式可写成
)()(
11011110 knknnknknnkn
fffhyyyy
++?+?++
+++++++?= βββααα LL
若则为隐式方法，若则为显式方法
0≠
k
β 0=
k
β
=
+=
∑
=
++
+++
1
0
)(
),(
k
j
jnjjnj
knknkkn
yfhg
gyxfhy
αβ
β
或写成上式称为线性k步法。
Remark：前面的R－K方法是增加一些非节点处的函数值的计算来提高单步法的精度的，这样使计算量增加了许多。而线性多步法每步只需要计算一个函数值。
2004-12-13 58
线性多步法的一般理论（续）
对于隐式情形（）的公式，由于f(x,y)一般是非线性函数，故难以求解到y
n+k
的显示表达式，
故常用迭代法求解：
0≠
k
β
gyxfhy
s
knknk
s
kn
+=
++
+
+
),(
)()1(
β
其中任意给出，s＝0,1,2,…,迭代到满足给定的精度要求为止。
)0(
kn
y
+
可以证明，当f(x,y)满足Lipschitz条件（或）
时，只要，迭代关系式就是收敛的。
L
y
f
≤
k
L
h
β
1
<
2004-12-13 59
线性多步法的一般理论（续）
定义：对于线性多步法定义处的局部截断误差为kn
x
+
)]()()([
110 knknn
xyxyxyh
++
′
++
′
+
′
βββ L
)()()(
110 knknn
xyxyxy
++
+++= ααα L
∑
=
+++
′
=
k
j
jnjjnjkn
xyhxyR
0
)]()([ βα
)(
1103322110 knknnknknnnn
fffhyyyyy
++++++
+++=+++++ βββααααα LL
kn
R
+按h展开的首项称主局部截断误差。
)1,,2,1,0(),(?==
++
kjxyy
jnjn
L
knknkn
yxyR
+++
= )(
kn
x
+
或在的假定下，定义处的局部截断误差为。
2004-12-13 60
线性多步法的一般理论（续）
若,则称线性多步法为p阶方法。
)(
1+
+
=
p
kn
hOR
若线性多步法为p阶方法，则
)()(
2)1(1
1
+++
++
+=
p
n
pp
pkn
hOxyhCR
设则有,1,,1,0),(?==
++
kjxyy
jnjn
L
knknknn
Rxyxyxyh
+++
+
′
++
′
+
′
= ))()()((
110
βββ L
∑
=
+++
+=
k
j
knjnjnj
Rxyxfh
0
))(,(β
)()()()(
110
0
knnn
k
j
jnj
xyxyxyxy
++
=
+
+++=
∑
Lααα
1+p
C其中称为主局部截断误差系数。
)(
)1(1
1 n
pp
p
xyhC
++
+
即主局部截断误差为
2004-12-13 61
线性多步法的一般理论（续）
将上式与有,1,,1,0),(?==
++
kjxyy
jnjn
L
相减，并注意
∑∑
=
++
=
+
=
k
j
jnjnj
k
j
jnj
yxfhy
00
),(βα
knknknknknk
Ryxfxyxfh
+++++
+?= )],())(,([β
knkn
yxy
++
)(
利用微分中值定理
))((
),(
knkn
knkn
yxy
y
xf
++
++
=
η
),())(,(
knknknkn
yxfxyxf
++++
其中介于与之间。
kn+
η kn
y
+ )(
kn
xy
+
))(](
),(
1[
knkn
knkn
kkn
yxy
y
xf
hR
++
++
+
=
η
β
故
2004-12-13 62
线性多步法的一般理论（续）
故在的假定下，若
1,,1,0),(?==
++
kjxyy
jnjn
L
0=
k
β
（显示公式），
knknkn
yxyR
+++
= )(
则：
若,且
（隐式公式）0≠
k
β
)()(
2)1(1
1
+++
++
+=
p
n
pP
pkn
hOxyhCR
即为p阶方法，则当y(x)充分可微时，
)]()(][
),(
1[
2)1(1
1
+++
+
++
++
+=
p
n
pp
p
knkn
k
hOxyhC
y
xf
h L
η
β
)]()([
),(
1
1
2)1(1
1
+++
+
++
+
=
p
n
pp
p
knkn
k
hOxyhC
y
xf
h
η
β
knkn
yxy
++
)(
)()(
2)1(1
1
+++
+
+=
p
n
pp
p
hOxyhC
2004-12-13 63
线性多步法的一般理论（续）
即的首项与的首项相同，因此两种局部截断误差的定义是相同的。
knkn
yxy
++
)(
kn
R
+
Remark1：可以利用此处的截断误差定义分析前面的单步隐式方法。对于Euler方法，其主局部截断误差为，而对于梯形方法，
其主局部截断误差为。
)(
2
1
2
n
xyh
′′
)(
12
1
3
n
xyh
′′′
Remark2：可以证明，显式线性多步法的整体截断误差比局部截断误差低一阶。
2004-12-13 64
二、基于数值积分的构造方法将方程两端从积分得))(,( xyxfy =
′
knjn
xx
+?
到
=
+=
∫
+
+
00
)(
))(,()()(
yxy
dttytfxyxy
kn
jn
x
x
jnkn
（*）
对，取等距插值节点，对应))(,( tytf
pnnn
xxx
,,,
1
L
))(,(,)),(,(
pnpnnn
xyxfxyxf
L
的函数值为，
构造p次Lagrange插值多项式：
))(,()()(
0
inin
p
i
ip
xyxftltL
=
∑
=
2004-12-13 65
基于数值积分的构造方法（续）
代入（*）式有
∑
=
+=
p
i
ininpijn
xyxfhxy
0
))(,()( β
∫
∑
+
=
+=
kn
jn
x
x
i
p
i
ininjn
dttlxyxfxy )())(,()(
0
∫
+
+≈
+
kn
jn
x
x
pjnkn
dttLxyxy )()()(
其中
，
ii
xxh?=
+1
∫
+
=
kn
jn
x
x
ipi
dttl
h
)(
1
β
),2,1,0(
)(
)( pi
xx
xt
tl
p
il
il
lnin
ln
i
L=
=
∏
≠
=
2004-12-13 66
基于数值积分的构造方法（续）
若作变量代换
。则],[],[,kjxxshxt
knjnn
+=
+?
hlslhxshxxshxxt
nnlnnln
)()( +=+=?+=?
hlilhxihxxx
nnlnin
)()()( +?==?
代入
∫
+
=
kn
jn
x
x
ipi
dttl
h
)(
1
β
则
∫
∏
≠
=
=
+?
+
=
k
j
p
il
l
pi
pids
li
ls
0
),,2,1,0( Lβ
，表示代替用),(),(
ininininin
yxffxyy
∑
=
+
+=
p
i
inpijnkn
fhyy
0
β
对k，j和p选择不同的值，即可以得到不同类型的具体公式。
则可得线性多步法显式公式
2004-12-13 67
基于数值积分的构造方法（续）
取可得到Adams显式公式
0,1 == jk
)(
1101 pnppnpnpnn
fffhyy
+
++++= βββ L
其中。最常用的是p＝3，
∫
∏
≠
=
+?
+
=
1
0
0
p
il
l
pi
ds
li
ls
β
]9375955[
24
3211+
+?+=
nnnnnn
ffff
h
yy
取k=0,j=1可得到Adams隐式公式
)(
*
1
*
1
*
01 pnppnpnpnn
fffhyy
++++= βββ L
再用n＋1代替n，得到
)(
1
**
11
*
01 +?++
++++=
pnppnpnpnn
fffhyy βββ L
其中
∫
∏
≠
=
+?
+
=
0
1
0
*
p
il
l
pi
ds
li
ls
β
2004-12-13 68
基于数值积分的构造方法（续）
最常用的是p＝3的情况：
]5199[
24
2111++
+?++=
nnnnnn
ffff
h
yy
取k＝1，j=1，得到Nystr?m显式公式：
=
+?
+
=
++++=
∫
∏
≠
=
+
1
1
0
11011
,,2,1,0,
)(
p
il
l
pi
pnppnpnpnn
pids
li
ls
fffhyy
L
L
β
βββ
取k=0，j=2，得到Milne隐式公式（用n＋1代替n）：
=
+?
+
=
++++=
∫
∏
≠
=
+?+?+
0
2
0
111011
,,2,1,0,
)(
p
il
l
pi
pnppnpnpnn
pids
li
ls
fffhyy
L
L
β
βββ
2004-12-13 69
基于数值积分的构造方法（续）
∏
=
+
+
=
p
i
in
p
p
xtyf
0
)1(
)!1(
1
)())(,( ξξ
)())(,( tLtytf
p
由
),max(),min(
npn
xtxt <<
ξ
于是线性多步法
==
+=
∫
∑
+
=
+
),2,1,0()(
1
0
pidttl
h
fhyy
kn
jn
x
x
ipi
p
i
inpijnkn
Lβ
β
的局部截断误差为：
2004-12-13 70
基于数值积分的构造方法（续）
knknkn
yxyR
+++
= )(
∫∫
+
+
+=
jn
jn
kn
n
x
x
p
x
x
jnjn
dttLydttytfxy )())(,()(
j－
∫
+
=
kn
jn
x
x
p
dttLtytf ))())(,((
dtxtyf
p
R
in
x
x
p
i
p
kn
kn
jn
)())(,(
)!1(
1
0
)1(
=
+
+
Π
+
=
∫
+
ξξ
即
∫
=
+
+
+Π
+
=
k
j
p
i
p
p
dsissysf
p
h
)()))((),((
)!1(
0
)1(
2
ξξ
2004-12-13 71
基于数值积分的构造方法（续）
对于Adams显式公式与隐式公式，由于)(
0
is
p
i
+Π
=
积分中值定理在[-1，0] (j=1，k=0，隐式公式)上恒负，故根据在[0，1]（j=0，k=1，显式公式）上恒正，
∫
+Π
+
=
=
+
+
+
1
0 0
)1(
2
,
)())(,(
)!1(
dsisyf
p
h
R
p
i
p
p
Ekn
ηη
∫
+Π
+
=
=
+
+
1
0 0
)2(
2
)()(
)!1(
dsisy
p
h
p
i
p
p
η
∫
=
+
+
+
+Π
+
=
0
1 0
)2(
2
,
)()(
)!1(
dsisy
p
h
R
p
i
p
p
Ikn
η
为某中间点，E表示Explicit，I表示Implicit。
η
2004-12-13 72
基于数值积分的构造方法（续）
当p=3时,
)()(
720
251
)(
720
251
6)5(5)5(5
,1
hOxyhyhR
nEn
+==
+
η
)()(
720
19
)(
720
19
6)5(5)5(5
,
hOxyhyhR
nIkn
+?=?=
+
η
由局部截断误差的分析可知，在f(x,y(x))具有p+1
阶连续偏导数的条件下，p+1步Adams显式方法与p步Adams隐式方法的局部截断误差与h
p+2
同阶，
即它们都是p＋1阶方法。特别地，当p＝3时，
Adams显、隐方法都是四阶的。
2004-12-13 73
三、基于Taylor 展开的构造方法
Taylor展开法更具一般性。下面构造四步方其为四阶方法，并求其局部截断误差的主项。
法，即k＝4。用Taylor展开法构造下述公式，使
33221104 ++++
+++=
nnnnn
yyyyy αααα
)(
443322110 ++++
+++++
nnnnn
fffffh βββββ
)3,2,1,(),( +++== nnnnixyy
ii
在上式右边近似认为
。)4,3,2,1,(),())(,( ++++=
′
== nnnnnixyxyxff
iiii
以及
)3,2,01()( =
+
ixy
in
)4,3,2,01()( =
′
+
ixy
in
n
x
然后将函数及在点按 的幂展开：
h
2004-12-13 74
基于Taylor 展开的构造方法（续）
])(
!5
)(
!4
)(
!3
)(
2
)()([)(
)5(
5
)4(
432
104
L+++
′′′
+
′′
+
′
++=
+ nnnnnnnn
xy
h
xy
h
xy
h
xy
h
xyhxyxyy
！
αα
])(
!4
)(
!3
)(
2
)()([)(
)5(
4
)4(
32
10
L+++
′′′
+
′′
+
′
+
′
+
nnnnnn
xy
h
xy
h
xy
h
xyhxyhxyh
！
ββ
])(
!5
)3(
)(
!4
)3(
)(
!3
)3(
)(
2
)3(
)(3)([
)5(
5
)4(
432
3
L
nnnnnn
xy
h
xy
h
xy
h
xy
h
xyhxy +++
′′′
+
′′
+
′
++
！
α
])(
!5
)2(
)(
!4
)2(
)(
!3
)2(
)(
2
)2(
)(2)([
)5(
5
)4(
432
2
L+++
′′′
+
′′
+
′
++
nnnnnn
xy
h
xy
h
xy
h
xy
h
xyhxy
！
α
])(
!4
)4(
)(
!3
)4(
)(
2
)4(
)(4)([
)5(
4
)4(
32
4
L+++
′′′
+
′′
+
′
+
nnnnn
xy
h
xy
h
xy
h
xyhxyh
！
β
])(
!4
)3(
)(
!3
)3(
)(
2
)3(
)(3)([
)5(
4
)4(
32
3
L+++
′′′
+
′′
+
′
+
nnnnn
xy
h
xy
h
xy
h
xyhxyh
！
β
])(
!4
)2(
)(
!3
)2(
)(
2
)2(
)(2)([
)5(
4
)4(
32
2
L+++
′′′
+
′′
+
′
+
nnnnn
xy
h
xy
h
xy
h
xyhxyh
！
β
2004-12-13 75
基于Taylor 展开的构造方法（续）
)()(
3210 n
xyαααα +++=
hxy
n
)()32(
43210321
′
++++++++ βββββααα
2
4321321
)()432
2
9
2
2
1
( hxy
n
′′
+++++++ ββββααα
3)3(
432
1
321
)()8
2
9
2
22
9
3
4
6
1
( hxy
n
βββ
β
ααα +++++++
(1)
4)4(
4321321
)()
3
32
2
9
3
4
6
1
8
27
3
2
24
1
( hxy
n
ββββααα +++++++
L++++++++
5)5(
4321321
)()
3
32
8
27
3
2
24
1
40
81
15
4
120
1
( hxy
n
ββββααα
2004-12-13 76
基于Taylor 展开的构造方法（续）
另一方面有
)()4(
!3
1
)()4(
!2
1
)('4)()(
32
4 nnnnn
xyhxyhxhyxyxy
′′′
+
′′
++=
+
（2）
L+++ )()4(
!5
1
)()4(
!4
1
)5(5)4(4
nn
xyhxyh
比较（1），（2）两式中h的同次幂，由于
)(xy
的任意性，可得：
2004-12-13 77
基于Taylor 展开的构造方法（续）
0
h
1
3210
=+++ αααα
432
43210321
=+++++++ βββββααα
1
h
8432
2
9
2
2
1
4321321
=++++++ ββββααα
2
h
3
32
8
2
9
2
22
9
3
4
6
1
432
1
321
=++++++ βββ
β
ααα
3
h
3
32
3
32
2
9
3
4
6
1
8
27
3
2
24
1
4321321
=++++++ ββββααα
4
h
15
128
3
32
8
27
3
2
24
1
40
81
15
4
120
1
4321321
=++++++ ββββααα
5
h
2004-12-13 78
基于Taylor 展开的构造方法（续）
如果要求前六个方程成立，得到的公式的局部截断误差为。也可以只要求前面的几个方程成立。
)(
6
hO
1、Simpson公式要求前四个方程成立。令,，
0
310
=== ααα
0
10
==ββ
.
3
1
,
3
4
,
3
1
,1
4322
==== βββα
]4[
3
23424 +++++
+++=
nnnnn
fff
h
yy
]4[
3
122 nnnnn
fff
h
yy +++=
+++
可得，
公式为即
2004-12-13 79
基于Taylor 展开的构造方法（续）
第五个方程恰好也成立。
将得
)1()2(?
44,4
)(
+++
=
nnhn
yxyR
)()()]
3
32
8
27
3
2
24
1
40
81
15
4
120
1
(
!5
4
[
65)5(
4321321
5
hOhxy
n
+++++++?= ββββααα
)()(
90
1
65)5(
hOhxy
n
+?=
2004-12-13 80
基于Taylor 展开的构造方法（续）
2.隐式的Adams公式若要求前五个方程成立，令0,0
0210
==== βααα
.
24
9
,
24
19
,
24
5
,
24
1
,1
43213
==?=== ββββα
得隐式的Adams公式为
]5199[
24
12334 nnnnnn
ffff
h
yy +?++=
+++++
将得：
)1()2(?
44,4
)(
+++
=
nnhn
yxyR
)()(
720
19
6)5(5
hOxyh
n
+
=
2004-12-13 81
基于Taylor 展开的构造方法（续）
3.显式的Adams公式要求前五个方程成立，令
0,0
4210
==== βααα
24
55
,
24
59
,
24
37
,
24
9
,1
32103
=?==?== ββββα
得显式Adams公式：
)9375955(
24
12334 nnnnnn
ffff
h
yy?+?+=
+++++
解得：
44,4
)(
+++
=
nnhn
yxyR
21321
5
3
2
24
1
40
81
15
4
120
1
(
!5
4
[ ββααα ++++?=
将得：
)1()2(?
)()(
720
251
65)5(
hOhxy
n
+=
)()()]
3
32
8
27
65)5(
43
hOhxy
n
+++ ββ
2004-12-13 82
基于Taylor 展开的构造方法（续）
4.Milne公式（四步四阶）
要求前四个方程成立，令
0,0
40210
===== ββααα
此时第五方程恰好成立。
得Milne公式：
)22(
3
4
1234 ++++
+?+=
nnnnn
fff
h
yy
44,4
)(
+++
=
nnhn
yxyR
)()()]
3
32
8
27
3
2
24
1
40
81
15
4
120
1
(
!5
4
[
65)5(
43
21321
5
hOhxy
n
+++
++++?=
ββ
ββααα
)()(
45
14
6)5(5
hOxyh
n
+=
解得
3
4
,
8
3
,1
2310
==== βββα
将得
)1()2(?
2004-12-13 83
基于Taylor 展开的构造方法（续）
5,Hamming公式要求前五个方程成立,令
0,0
1020
==== ββαα
可得Hamming公式：
解得
8
3
,
8
6
,
8
3
,
8
9
,
8
1
43232
==?=== βββαα
)2(
8
3
)9(
8
1
234134 ++++++
++?=
nnnnnn
fff
h
yyy
)2(
8
3
)9(
8
1
12323 +++++
++?=
nnnnnn
fff
h
yyy
将得：
)1()2(?
44,4
)(
+++
=
nnhn
yxyR
321
5
40
81
15
4
120
1
(
!5
4
[ ααα ++?=
)()(
40
1
65)5(
hOhxy
n
+?=)()()]
3
32
8
27
3
2
24
1
65)5(
4321
hOhxy
n
+++++ ββββ
2004-12-13 84
基于Taylor 展开的构造方法（续）
Remark1：用Taylor展开法比用数值积分法推导线性多步法更加灵活，推导出的方法也比积分法的更多。用数值积分法可以得到的方法，用Taylor
展开法都可以得到，但用Taylor展开法得到的方法，用数值积分法却不一定能够得到。
Back
Remark2：应用线性多步法求解初值问题时，开头几点处的数值要用别的方法先计算出来。一般选用与多步法同阶或更高阶的单步法，如Runge-Kutta法、
Taylor展开法等。
Remark3：对线性隐式多步方法，除开头几点的函数值需要单独计算外，还需要迭代求解或采用预估-校正方法求解。
2004-12-13 85
数值算例求解常微分方程初值问题
=
=
1)0(
2
y
y
x
y
dx
dy
在[0,1]上的解，取步长h＝0.1。
计算结果如下图所示：
2004-12-13 86
2004-12-13 87
