数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
第 8章 常微分方程
实际中，很多问题的数学模型都是微分方程。我们可以研究它们的一些
性质。但是，只有极少数特殊的方程有解析解。对于绝大部分的微分方程是
没有解析解的。
常微分方程作为微分方程的基本类型之一，在自然界与工程界有很广泛
的应用。很多问题的数学表述都可以归结为常微分方程的定解问题。很多偏
微分方程问题，也可以化为常微分方程问题来近似求解。
本章讨论常微分方程的数值解法
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
对于一个常微分方程：
],[,),(' baxyxfdxdyy ???
通常会有无穷个解。如：
Raaxyxdxdy ??????,)s i n ( )c o s (
因此，我们要加入一个限定条件。通常会在端点出给出，如下面的初值问题：
??
?
?
?
?
??
0)(
],[,),(
yay
baxyxf
dx
dy
为了使解存在唯一，一般，要加限制条件在 f上，要求 f对 y满足 Lipschitz条件：
2121 ),(),( yyLyxfyxf ???
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
常微分方程的解是一个 函数,但是，计算机没有办法对函数进行运算。
因此，常微分方程的数值解并不是求函数的近似，而是求解函数在某些节
点的近似值。
例：我们对区间做等距分割：,( ) /
ix h i h b a m? ? ?
设解函数在节点的近似为 }{
iy
由数值微分公式，我们有
i
i
xx
xx
yxfdxdy ?
?
? ),(
，则：
),(1 iiii yxfh yy ???
向前差商公式
),( 1 iiii yxfhyy ???
可以看到，给出初值，就可以用上式求出所有的 }{
iy
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
基本步骤如下：
③ 解差分方程，求出格点函数
① 对区间作分割,bxxxa nI ?????? ?10:
求 在 上的近似值 。)(xy
ix iy }{ iy
称为分割 I? 上的格点函数
② 由微分方程出发，建立求格点函数的差分方程。这个方程应该满足：
A、解存在唯一； B、稳定，收敛； C、相容
数值方法，主要研究步骤②，即如何建立差分方程，并研究差分方程的性质。
这种方法,称为 数值离散方法 。求的是在一系列离散点列上，求未知函数 y在这些
点上的值的近似。
我们的目的，就是求这个格点函数
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
为了考察数值方法提供的数值解，是否有实用价值，需要知道如下几个结论：
① 步长充分小时，所得到的数值解能否逼近问题得真解；即收敛性问题
② 误差估计
③ 产生得舍入误差，在以后得各步计算中，是否会无限制扩大；稳定性问题
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
8.1 Euler公式
做等距分割，利用数值微分代替导数项，建立差分方程。
im abx iI ???,
1、向前差商公式
)(''2)(')()( 1 nnnn yhxyh xyxy ?????
)(''2))(,()()( 1 nnnnn yhxyxfh xyxy ?????
)(''2))(,()()(
2
1 nnnnn y
hxyxhfxyxy ????
?
所以，可以构造差分方程
),(1 nnnn yxhfyy ???
称为局部截断误差。
显然，这个误差在逐
步计算过程中会传播，
积累。因此还要估计
这种积累
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
定义 在假设 yi = y(xi),即第 i 步计算是精确的前提下, 考虑
的截断误差 Ri = y(xi+1) ? yi+1 称为 局部截断误差 /* local
truncation error */。
定义 若某算法的局部截断误差为 O(hp+1),则称该算法有 p
阶精度 。
记为
2、收敛性
)(''2))(,()()(
2
1 nnnnn y
hxyxhfxyxy ????
?
),(1 nnnn yxhfyy ???
1?nhT
考察局部误差的传播和积累
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
111 )( ??? ?? nnn yxye
1),())(,()( ?????? nnnnnnn Thyxfxyxfhyxy
1)( ????? nnnn ThyxyhLe
jjn TThTehL m a x,)1( ????
? ? hThTehLhL n ????? ? 1)1()1(
? ?hThLehL n 1)1()1( 12 ????? ?
? ? ? ?hThLhTehLhL n 1)1()1()1( 22 ??????? ?
? ?hThLhLehL n 1)1()1()1( 223 ??????? ?
??
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
? ?hThLhLehL nn 1)1()1()1( 01 ???????? ? ?
hThLhLehL
n
n
)1(1
)1(1)1( 1
0
1
??
????? ??
hThLhLehL
n
n
1
0
1 )1()1(
?
? ????
?????? ??? ? LTehL n 01)1(
??????? ? LTe hLn )1(
)( 1 hOe n ?? ?
)(
)1(
00
hOT
ex
e
nxn
?
??
?
称为整体截断误差是 1阶方法
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
3、稳定性－误差在以后各步的计算中不会无限制扩大。是格式对舍入误差的抑止作用
我们考虑一种简单情况，即仅初值有误差，而其他计算步骤无误差。
设 }{
iz
是初值有误差后的计算值，则
),(
),(
1
1
nnnn
nnnn
zxhfzz
yxhfyy
??
??
?
?
所以，我们有：
),(),(111 nnnnnnnn zxfyxfhezye ????? ???
)1( hLezyhLe nnnn ?????
hLnn eehLe )1(010 )1( ?? ???? ?
可以看出，向前差商公式关于初值是稳定的。当初始误差充分小，以后各步的误差
也充分小
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
4、向后差商公式
)(''2)(')()( 11 nnnn yhxyh xyxy ???? ??
)(''2))(,()()( 111 nnnnn yhxyxfh xyxy ???? ???
)(''2))(,()()(
2
111 nnnnn y
hxyxhfxyxy ????
???
),( 111 ??? ?? nnnn yxhfyy
是隐格式，要迭代求解
)0(
1
)(
11
)1(
1 ),(
?
??
?
? ??
n
k
nnn
k
n
y
yxhfyy
可以由向前差商公式求出
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
5、中心差商公式
)(''2)(')()( 111 nnnn yhxyh xyxy ???? ???
),( 1111 ???? ?? nnnn yxhfyy
是多步,2阶格式，该格式不稳定
6、梯形法－基于数值积分的公式
对微分方程
],[,),(' baxyxfdxdyy ???
做积分，则：
? ? ?1 ),(nnxx yxfdxdy
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
)( 21 hOe n ??
类似，可以算出其误差估计式：
2阶的方法
所以，有格式为：
)],(),([2 111 ??? ??? nnnnnn yxfyxfhyy
是个隐式的方法，要用迭代法求解
1
1 ( ) ( ) (,( ) )
n
n
x
nn xy x y x f x y x d x
?
?? ? ? ?
???
?
???
? ?????
??? )('''12))](,())(,([2)()(
2
111 ?f
hxyxfxyxfhxyxy
nnnnnn
局部截断误差
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
8.2 Runge－ Kutta法
由 Taylor展开
)()!1()(!)(')()( )1(
1
)(
1 n
k
k
n
k
k
nnn yk
hxy
k
hxhyxyxy ???
? ?????? ?
记为
1?nhT
),()(' yxfxy ?
'),(),()('' yyxfyxfxy yx ???
??)(''' xy
所以，可以构造格式
? ? ???????? ),(),(),(!2),( 21 nnnnynnxnnnn yxfyxfyxfhyxhfyy
这种格式使用到了各阶偏导数，使用不便。
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
从另一个角度看，
11 )),(,,()()( ?? ??? nnnnn hTfxyxhhFxyxy
取 (x,y)及其附近的点做线性组合，表示 F，问题就好办了。当然，要求此时的展开精
度相同。这种方法称为 Runge－ Kutta法
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
在 (x,y)处展开，
? ?)(),(),(),(),(),(
),,,(
2
21221 hOyxfyxhfbyxhfayxfcyxfc
fyxhF
yx ?????
比较
? ? 1
1
),(),(),(
!2
),(
)()(
?
?
??
?
?
??
? ????
?
nnnnnynnxnn
nn
hTyxfyxfyxf
h
yxfh
xyxy
以 2阶为例，设
)),(,(),(),,,( 21221 yxhfbyhaxfcyxfcfyxhF ????
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
有：
?
?
?
?
?
?
?
??
2/1
2/1
1
212
22
21
bc
ac
cc
?
?
?
?
?
?
?
??
1
1
2/1
21
2
21
b
a
cc
1、改进的 Euler公式
?
?
?
?
?
???
?
???
?
?
),(
),(
][2/
12
1
211
hKyhxfK
yxfK
KKhyy
nn
nn
nn
?
?
?
?
?
?
?
??
3/2
3/2
4/3,4/1
21
2
21
b
a
cc
2,Heun公式
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
一般的 Runge－ Kutta法构造
?
?
?
?
?
?
?
?
?
???
???
?
????
?
?
?
1
1
12122
1
2211
,(
),(
),(
),,,(
m
i
imimm
mm
KbhyhaxfK
hKbyhaxfK
yxfK
KcKcKcfyxhF
?
?
常见的为 3阶,4阶公式
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
8.3 线性多步法
用 若干 节点处的 y及 y’ 值的 线性组合 来近似
y(xn+1)。
)...(..,110111101 knknnnknknnn ffffhyyyy ??????? ????????? bbbbaaa
其通式可写为：
),( nnn yxff ?
当 b?1?0 时，为 隐式公
式 ; b?1=0 则为 显式公式 。
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
?基于数值积分的构造法
将 在 上积分，得到),( yxfy ?? ],[
1?? npn xx
? ???? ?? 1 ))(,()()( 1 i pnxxpnn dxxyxfxyxy
只要 近似地算出右边的积分,则可通
过 近似 y(xn+1) 。而 选用不同近似式 Ik，可得到不
同的计算公式 。
? ??? 1 ))(,(n pnxxk dxxyxfI
kpnn Iyy ?? ?? 1
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
? ?? ?
?
1 )(
)!1(
)()2(n
pn
x
x q
q
dxxqy ??
? ?? 1 )('n pnxx dxxy
若积分 用节点
qnnn xxx ??,,,1 ?
作为积分点，则有
1110 )](')(')('[)('
1
??? ??????
?
? nqnqnn
x
x hTxyaxyaxyahdxxy
n
pn
?
))(,()(' nnn xyxfxy ?
1
0
1 ))(,()()( ?
?
???? ??? ? n
q
j
jnjnjpnn hTxyxfahxyxy
积分系数
? ??? 1 )(n pnxx jj dxxlha
这是显格式,q+1阶 r+1步格式。 r=max{p,q}
11,,,??? qnnn xxx ?
为积分节点，可以构造 r+1步 q+1阶隐格式
局部截断误差
同样，若以
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
例：建立 p=1,q=2的显格式
p=1，
q=2，显格式，
? ?? 11 )('nnxx dxxy积分区间为
积分节点为
21,,?? nnn xxx
hdxxxxx xxxxha n
n
x
x nnnn
nn
3
7
))((
))((1
1 21
21
0 ???
??? ? ?
? ??
??
所以
hdxxxxx xxxxha n
n
x
x nnnn
nn
3
2
))((
))((1
1 211
2
1 ????
??? ? ?
? ???
?
hdxxxxx xxxxha n
n
x
x nnnn
nn
3
1
))((
))((1
1 122
1
2 ???
??? ? ?
? ???
?
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
)(31))()(()!3( )( )4(421
)4(
1
1
1
?? yhdxxxxxxxyT n
n
x
x nnnn
????? ? ?
?
???
例：建立 p=2,q=2的隐格式
p=2，
q=2，隐格式，
? ?? 12 )('nnxx dxxy积分区间为
积分节点为
11,,?? nnn xxx
hdxxxxx xxxxha n
n
x
x nnnn
nn
4
3
))((
))((1
2 111
1
0 ???
??? ? ?
? ???
?
所以
0))(( ))((1
2 11
11
1 ???
??? ? ?
? ??
??n
n
x
x nnnn
nn dx
xxxx
xxxxha
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
)(83))()(()!3( )( )4(411
)4(
1
1
1
?? yhdxxxxxxxyT n
n
x
x nnnn
?????? ? ?
?
???
hdxxxxx xxxxha n
n
x
x nnnn
nn
4
9
))((
))((1
2 111
1
0 ???
??? ? ?
? ???
?
它的截断误差较 显格式 小，通常也具有更好的稳定性。
? Adams公式 －－ p=0 时候的多步法
参见书
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
§ 8.4 方程组和高阶方程的数值解法
bxa
ay
ay
yyxf
dx
dy
yyxf
dx
dy
mm
mm
m
m
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)(
),,,(
),,,(
11
1
11
1
?
?
?
?
?
?
写成向量的形式：
??
?
?
?
?
?
?)(
),(
aY
YxF
dx
dY
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
各种方法都可以直接运用过来。
bxa
zaz
yay
zyxg
dx
dz
zyxf
dx
dy
??
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)(
),,(
),,(
0
0
Euler公式
?
?
?
??
??
?
?
),,(
),,(
1
1
nnnnn
nnnnn
zyxhgzz
zyxhfyy
以两个方程的方程组为例
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
Runge-Kutta公式
1
1 2 3 4
1
( 2 2 )6nn
nn
yy h K K K K
zz
?
?
? ? ? ?? ? ? ? ?
? ? ? ?? ? ? ?
)2,2( 12 KhYhxFK ???
1
( 1 ) ( 2 )
11
2
( 1 ) ( 2 )
11
(,,)
(,,)
(,,)
2 2 2
(,,)
2 2 2
n n n
n n n
n n n
n n n
f x y z
K
g x y z
h h h
f x y K z K
K
h h h
g x y K z K
??
? ??
??
??
? ? ?
??
?
??
? ? ?
??
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS ( 1 ) ( 1 )
22
3
( 1 ) ( 1 )
22
( 1 ) ( 1 )
33
4 ( 1 ) ( 1 )
33
(,,)
2 2 2
(,,)
2 2 2
(,,)
(,,)
n n n
n n n
n n n
n n n
h h h
f x y K z K
K
h h h
g x y K z K
f x h y h K z h K
K
g x h y h K z h K
??
? ? ?
??
?
??
? ? ?
??
??? ? ?
? ??
? ? ???
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
0.0 5 ( 1 ) 0.0 02
20
0.0 9 ( 1 ) 0.1 5
15
( 0) 0.1 93
( 0) 0.0 83
du u
u uv
dt
dv v
v uv
dt
u
v
?
? ? ?
?
?
?
? ? ??
?
??
?
??
1、
(,,) 0, 0 5 ( 1 / 2 0 ) 0, 0 0 2
(,,) 0, 0 9 ( 1 / 1 5 ) 0, 1 5
f u v t u u u vF
g u v t v v u v
??? ? ? ???
? ? ? ???? ? ? ?
2、确定方法，然后求解 （ 0.20276 0.0881157）
（ 0.213007 0.0934037）
（ 0.223763 0.0988499）
（ 0.235052 0.104437）
（ 0.246902 0.110146）
4阶 Runge-Kutta法,h=1
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
高阶方程
bxa
ay
ay
ay
yyyxf
dx
yd
m
m
m
m
m
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)('
)(
),,',,(
)1(
2
1
)1(
1
?
?
?
?
?
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
则有：
bxa
ay
ay
yyxf
dx
dy
y
dx
dy
mm
m
m
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)(
),,,(
11
1
2
1
?
?
?
?
?
令
?
?
?
?
??
?
?
?
?
?
?
?
m
m
y
dx
dy
y
dx
dy
yy
1
2
1
1
?
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
例,考察初值问题 在区间 [0,0.5]上的解。
分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。
??
?
?
???
1)0(
)(30)(
y
xyxy
0.0
0.1
0.2
0.3
0.4
0.5
精确解改进欧拉法欧拉隐式欧拉显式节点 xi xey 30??
1.0000
?2.0000
4.0000
?8.0000
1.6000?101
?3.2000?101
1.0000
2.5000?10?1
6.2500?10?2
1.5625?10?2
3.9063?10?3
9.7656?10?4
1.0000
2.5000
6.2500
1.5626?101
3.9063?101
9.7656?101
1.0000
4.9787?10?2
2.4788?10?3
1.2341?10?4
6.1442?10?6
3.0590?10?7
What is wrong!
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
§ 8.5 差分方程的绝对稳定性
对于一般的差分方程
??
?
??
?
? ?
k
j
jnjnj
k
j
jnj yxfbhya
00
),(
由初始误差产生了差分解的误差，实际上是同一差分方程，取不同初值所得到的 2组
差分解之间的差。这个差不仅于差分方程本身有关，而且与微分方程本身有关。如果
微分方程本身是不稳定，那就没理由要求这 2组解充分接近。因此，差分方程的稳定性
概念是建立在微分方程稳定的基础上的。把这个典型微分方程规定为：
仍然考虑最简单的模型，即只有初值产生误差，看看这个误差的传播。
??
?
?
?
?
??
0)(
)0( R e
yay
y
dx
dy ??
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
差分方程运用到如上的微分方程后，可以得到
0)R e(,
00
?? ??
?
?
?
? ??
k
j
jnj
k
j
jnj ybhya
对于给定的初始误差
110,,,?keee ?
，误差方程具有一样的形式
0)R e(,
00
?? ??
?
?
?
? ??
k
j
jnj
k
j
jnj ebhea
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
定义：差分方程称为绝对稳定的，若差分方程作用到微分方程
)0( R e ?? ?? ydxdy
时，对任意的初值，总存在左半复平面上的一个区域，当 在这个区域时，差分
方程的解趋于 0。这个区域称为稳定区域 h?
例：向后 Euler公式的稳定性
11 ?? ?? nnn hyyy ?
误差方程：
11 ?? ?? nnn heee ?
he
e
n
n
??
??
1
11 1? 210 Re
Img
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
考察隐式欧拉法
11 ?? ?? iii yhyy ?
1
1
1iiyyh??
??? ??
???
1
10
1
1
i
i h???
?
?
??? ??
???
可见绝对稳定区域为：
| 1 | 1h???
210 Re
Img
注,一般来说，隐式欧拉法的绝对稳定性比同阶的显式
法的好。
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
3阶 Runge－ Kutta
? ?? ?
2
1
1
1
3211
1
)
2
1
1(
]4[
6
hhyK
hyK
yK
KKK
h
yy
n
n
n
nn
???
??
?
???
??
?
????
?
? ? ? ? ?
?
?
??
? ????
?
32
1 6
1
2
11 hhhyy
nn ???
显式 1~ 4 阶方法的绝对稳定区域为
k=1
k=2
k=3
k=4
-1-2-3
-
-
-
1
2
3
Re
Img
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
第 8章 常微分方程
实际中，很多问题的数学模型都是微分方程。我们可以研究它们的一些
性质。但是，只有极少数特殊的方程有解析解。对于绝大部分的微分方程是
没有解析解的。
常微分方程作为微分方程的基本类型之一，在自然界与工程界有很广泛
的应用。很多问题的数学表述都可以归结为常微分方程的定解问题。很多偏
微分方程问题，也可以化为常微分方程问题来近似求解。
本章讨论常微分方程的数值解法
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
对于一个常微分方程：
],[,),(' baxyxfdxdyy ???
通常会有无穷个解。如：
Raaxyxdxdy ??????,)s i n ( )c o s (
因此，我们要加入一个限定条件。通常会在端点出给出，如下面的初值问题：
??
?
?
?
?
??
0)(
],[,),(
yay
baxyxf
dx
dy
为了使解存在唯一，一般，要加限制条件在 f上，要求 f对 y满足 Lipschitz条件：
2121 ),(),( yyLyxfyxf ???
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
常微分方程的解是一个 函数,但是，计算机没有办法对函数进行运算。
因此，常微分方程的数值解并不是求函数的近似，而是求解函数在某些节
点的近似值。
例：我们对区间做等距分割：,( ) /
ix h i h b a m? ? ?
设解函数在节点的近似为 }{
iy
由数值微分公式，我们有
i
i
xx
xx
yxfdxdy ?
?
? ),(
，则：
),(1 iiii yxfh yy ???
向前差商公式
),( 1 iiii yxfhyy ???
可以看到，给出初值，就可以用上式求出所有的 }{
iy
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
基本步骤如下：
③ 解差分方程，求出格点函数
① 对区间作分割,bxxxa nI ?????? ?10:
求 在 上的近似值 。)(xy
ix iy }{ iy
称为分割 I? 上的格点函数
② 由微分方程出发，建立求格点函数的差分方程。这个方程应该满足：
A、解存在唯一； B、稳定，收敛； C、相容
数值方法，主要研究步骤②，即如何建立差分方程，并研究差分方程的性质。
这种方法,称为 数值离散方法 。求的是在一系列离散点列上，求未知函数 y在这些
点上的值的近似。
我们的目的，就是求这个格点函数
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
为了考察数值方法提供的数值解，是否有实用价值，需要知道如下几个结论：
① 步长充分小时，所得到的数值解能否逼近问题得真解；即收敛性问题
② 误差估计
③ 产生得舍入误差，在以后得各步计算中，是否会无限制扩大；稳定性问题
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
8.1 Euler公式
做等距分割，利用数值微分代替导数项，建立差分方程。
im abx iI ???,
1、向前差商公式
)(''2)(')()( 1 nnnn yhxyh xyxy ?????
)(''2))(,()()( 1 nnnnn yhxyxfh xyxy ?????
)(''2))(,()()(
2
1 nnnnn y
hxyxhfxyxy ????
?
所以，可以构造差分方程
),(1 nnnn yxhfyy ???
称为局部截断误差。
显然，这个误差在逐
步计算过程中会传播，
积累。因此还要估计
这种积累
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
定义 在假设 yi = y(xi),即第 i 步计算是精确的前提下, 考虑
的截断误差 Ri = y(xi+1) ? yi+1 称为 局部截断误差 /* local
truncation error */。
定义 若某算法的局部截断误差为 O(hp+1),则称该算法有 p
阶精度 。
记为
2、收敛性
)(''2))(,()()(
2
1 nnnnn y
hxyxhfxyxy ????
?
),(1 nnnn yxhfyy ???
1?nhT
考察局部误差的传播和积累
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
111 )( ??? ?? nnn yxye
1),())(,()( ?????? nnnnnnn Thyxfxyxfhyxy
1)( ????? nnnn ThyxyhLe
jjn TThTehL m a x,)1( ????
? ? hThTehLhL n ????? ? 1)1()1(
? ?hThLehL n 1)1()1( 12 ????? ?
? ? ? ?hThLhTehLhL n 1)1()1()1( 22 ??????? ?
? ?hThLhLehL n 1)1()1()1( 223 ??????? ?
??
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
? ?hThLhLehL nn 1)1()1()1( 01 ???????? ? ?
hThLhLehL
n
n
)1(1
)1(1)1( 1
0
1
??
????? ??
hThLhLehL
n
n
1
0
1 )1()1(
?
? ????
?????? ??? ? LTehL n 01)1(
??????? ? LTe hLn )1(
)( 1 hOe n ?? ?
)(
)1(
00
hOT
ex
e
nxn
?
??
?
称为整体截断误差是 1阶方法
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
3、稳定性－误差在以后各步的计算中不会无限制扩大。是格式对舍入误差的抑止作用
我们考虑一种简单情况，即仅初值有误差，而其他计算步骤无误差。
设 }{
iz
是初值有误差后的计算值，则
),(
),(
1
1
nnnn
nnnn
zxhfzz
yxhfyy
??
??
?
?
所以，我们有：
),(),(111 nnnnnnnn zxfyxfhezye ????? ???
)1( hLezyhLe nnnn ?????
hLnn eehLe )1(010 )1( ?? ???? ?
可以看出，向前差商公式关于初值是稳定的。当初始误差充分小，以后各步的误差
也充分小
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
4、向后差商公式
)(''2)(')()( 11 nnnn yhxyh xyxy ???? ??
)(''2))(,()()( 111 nnnnn yhxyxfh xyxy ???? ???
)(''2))(,()()(
2
111 nnnnn y
hxyxhfxyxy ????
???
),( 111 ??? ?? nnnn yxhfyy
是隐格式，要迭代求解
)0(
1
)(
11
)1(
1 ),(
?
??
?
? ??
n
k
nnn
k
n
y
yxhfyy
可以由向前差商公式求出
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
5、中心差商公式
)(''2)(')()( 111 nnnn yhxyh xyxy ???? ???
),( 1111 ???? ?? nnnn yxhfyy
是多步,2阶格式，该格式不稳定
6、梯形法－基于数值积分的公式
对微分方程
],[,),(' baxyxfdxdyy ???
做积分，则：
? ? ?1 ),(nnxx yxfdxdy
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
)( 21 hOe n ??
类似，可以算出其误差估计式：
2阶的方法
所以，有格式为：
)],(),([2 111 ??? ??? nnnnnn yxfyxfhyy
是个隐式的方法，要用迭代法求解
1
1 ( ) ( ) (,( ) )
n
n
x
nn xy x y x f x y x d x
?
?? ? ? ?
???
?
???
? ?????
??? )('''12))](,())(,([2)()(
2
111 ?f
hxyxfxyxfhxyxy
nnnnnn
局部截断误差
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
8.2 Runge－ Kutta法
由 Taylor展开
)()!1()(!)(')()( )1(
1
)(
1 n
k
k
n
k
k
nnn yk
hxy
k
hxhyxyxy ???
? ?????? ?
记为
1?nhT
),()(' yxfxy ?
'),(),()('' yyxfyxfxy yx ???
??)(''' xy
所以，可以构造格式
? ? ???????? ),(),(),(!2),( 21 nnnnynnxnnnn yxfyxfyxfhyxhfyy
这种格式使用到了各阶偏导数，使用不便。
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
从另一个角度看，
11 )),(,,()()( ?? ??? nnnnn hTfxyxhhFxyxy
取 (x,y)及其附近的点做线性组合，表示 F，问题就好办了。当然，要求此时的展开精
度相同。这种方法称为 Runge－ Kutta法
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
在 (x,y)处展开，
? ?)(),(),(),(),(),(
),,,(
2
21221 hOyxfyxhfbyxhfayxfcyxfc
fyxhF
yx ?????
比较
? ? 1
1
),(),(),(
!2
),(
)()(
?
?
??
?
?
??
? ????
?
nnnnnynnxnn
nn
hTyxfyxfyxf
h
yxfh
xyxy
以 2阶为例，设
)),(,(),(),,,( 21221 yxhfbyhaxfcyxfcfyxhF ????
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
有：
?
?
?
?
?
?
?
??
2/1
2/1
1
212
22
21
bc
ac
cc
?
?
?
?
?
?
?
??
1
1
2/1
21
2
21
b
a
cc
1、改进的 Euler公式
?
?
?
?
?
???
?
???
?
?
),(
),(
][2/
12
1
211
hKyhxfK
yxfK
KKhyy
nn
nn
nn
?
?
?
?
?
?
?
??
3/2
3/2
4/3,4/1
21
2
21
b
a
cc
2,Heun公式
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
一般的 Runge－ Kutta法构造
?
?
?
?
?
?
?
?
?
???
???
?
????
?
?
?
1
1
12122
1
2211
,(
),(
),(
),,,(
m
i
imimm
mm
KbhyhaxfK
hKbyhaxfK
yxfK
KcKcKcfyxhF
?
?
常见的为 3阶,4阶公式
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
8.3 线性多步法
用 若干 节点处的 y及 y’ 值的 线性组合 来近似
y(xn+1)。
)...(..,110111101 knknnnknknnn ffffhyyyy ??????? ????????? bbbbaaa
其通式可写为：
),( nnn yxff ?
当 b?1?0 时，为 隐式公
式 ; b?1=0 则为 显式公式 。
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
?基于数值积分的构造法
将 在 上积分，得到),( yxfy ?? ],[
1?? npn xx
? ???? ?? 1 ))(,()()( 1 i pnxxpnn dxxyxfxyxy
只要 近似地算出右边的积分,则可通
过 近似 y(xn+1) 。而 选用不同近似式 Ik，可得到不
同的计算公式 。
? ??? 1 ))(,(n pnxxk dxxyxfI
kpnn Iyy ?? ?? 1
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
? ?? ?
?
1 )(
)!1(
)()2(n
pn
x
x q
q
dxxqy ??
? ?? 1 )('n pnxx dxxy
若积分 用节点
qnnn xxx ??,,,1 ?
作为积分点，则有
1110 )](')(')('[)('
1
??? ??????
?
? nqnqnn
x
x hTxyaxyaxyahdxxy
n
pn
?
))(,()(' nnn xyxfxy ?
1
0
1 ))(,()()( ?
?
???? ??? ? n
q
j
jnjnjpnn hTxyxfahxyxy
积分系数
? ??? 1 )(n pnxx jj dxxlha
这是显格式,q+1阶 r+1步格式。 r=max{p,q}
11,,,??? qnnn xxx ?
为积分节点，可以构造 r+1步 q+1阶隐格式
局部截断误差
同样，若以
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
例：建立 p=1,q=2的显格式
p=1，
q=2，显格式，
? ?? 11 )('nnxx dxxy积分区间为
积分节点为
21,,?? nnn xxx
hdxxxxx xxxxha n
n
x
x nnnn
nn
3
7
))((
))((1
1 21
21
0 ???
??? ? ?
? ??
??
所以
hdxxxxx xxxxha n
n
x
x nnnn
nn
3
2
))((
))((1
1 211
2
1 ????
??? ? ?
? ???
?
hdxxxxx xxxxha n
n
x
x nnnn
nn
3
1
))((
))((1
1 122
1
2 ???
??? ? ?
? ???
?
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
)(31))()(()!3( )( )4(421
)4(
1
1
1
?? yhdxxxxxxxyT n
n
x
x nnnn
????? ? ?
?
???
例：建立 p=2,q=2的隐格式
p=2，
q=2，隐格式，
? ?? 12 )('nnxx dxxy积分区间为
积分节点为
11,,?? nnn xxx
hdxxxxx xxxxha n
n
x
x nnnn
nn
4
3
))((
))((1
2 111
1
0 ???
??? ? ?
? ???
?
所以
0))(( ))((1
2 11
11
1 ???
??? ? ?
? ??
??n
n
x
x nnnn
nn dx
xxxx
xxxxha
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
)(83))()(()!3( )( )4(411
)4(
1
1
1
?? yhdxxxxxxxyT n
n
x
x nnnn
?????? ? ?
?
???
hdxxxxx xxxxha n
n
x
x nnnn
nn
4
9
))((
))((1
2 111
1
0 ???
??? ? ?
? ???
?
它的截断误差较 显格式 小，通常也具有更好的稳定性。
? Adams公式 －－ p=0 时候的多步法
参见书
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
§ 8.4 方程组和高阶方程的数值解法
bxa
ay
ay
yyxf
dx
dy
yyxf
dx
dy
mm
mm
m
m
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)(
),,,(
),,,(
11
1
11
1
?
?
?
?
?
?
写成向量的形式：
??
?
?
?
?
?
?)(
),(
aY
YxF
dx
dY
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
各种方法都可以直接运用过来。
bxa
zaz
yay
zyxg
dx
dz
zyxf
dx
dy
??
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)(
),,(
),,(
0
0
Euler公式
?
?
?
??
??
?
?
),,(
),,(
1
1
nnnnn
nnnnn
zyxhgzz
zyxhfyy
以两个方程的方程组为例
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
Runge-Kutta公式
1
1 2 3 4
1
( 2 2 )6nn
nn
yy h K K K K
zz
?
?
? ? ? ?? ? ? ? ?
? ? ? ?? ? ? ?
)2,2( 12 KhYhxFK ???
1
( 1 ) ( 2 )
11
2
( 1 ) ( 2 )
11
(,,)
(,,)
(,,)
2 2 2
(,,)
2 2 2
n n n
n n n
n n n
n n n
f x y z
K
g x y z
h h h
f x y K z K
K
h h h
g x y K z K
??
? ??
??
??
? ? ?
??
?
??
? ? ?
??
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS ( 1 ) ( 1 )
22
3
( 1 ) ( 1 )
22
( 1 ) ( 1 )
33
4 ( 1 ) ( 1 )
33
(,,)
2 2 2
(,,)
2 2 2
(,,)
(,,)
n n n
n n n
n n n
n n n
h h h
f x y K z K
K
h h h
g x y K z K
f x h y h K z h K
K
g x h y h K z h K
??
? ? ?
??
?
??
? ? ?
??
??? ? ?
? ??
? ? ???
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
0.0 5 ( 1 ) 0.0 02
20
0.0 9 ( 1 ) 0.1 5
15
( 0) 0.1 93
( 0) 0.0 83
du u
u uv
dt
dv v
v uv
dt
u
v
?
? ? ?
?
?
?
? ? ??
?
??
?
??
1、
(,,) 0, 0 5 ( 1 / 2 0 ) 0, 0 0 2
(,,) 0, 0 9 ( 1 / 1 5 ) 0, 1 5
f u v t u u u vF
g u v t v v u v
??? ? ? ???
? ? ? ???? ? ? ?
2、确定方法，然后求解 （ 0.20276 0.0881157）
（ 0.213007 0.0934037）
（ 0.223763 0.0988499）
（ 0.235052 0.104437）
（ 0.246902 0.110146）
4阶 Runge-Kutta法,h=1
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
高阶方程
bxa
ay
ay
ay
yyyxf
dx
yd
m
m
m
m
m
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)('
)(
),,',,(
)1(
2
1
)1(
1
?
?
?
?
?
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
则有：
bxa
ay
ay
yyxf
dx
dy
y
dx
dy
mm
m
m
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
,
)(
)(
),,,(
11
1
2
1
?
?
?
?
?
令
?
?
?
?
??
?
?
?
?
?
?
?
m
m
y
dx
dy
y
dx
dy
yy
1
2
1
1
?
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
例,考察初值问题 在区间 [0,0.5]上的解。
分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。
??
?
?
???
1)0(
)(30)(
y
xyxy
0.0
0.1
0.2
0.3
0.4
0.5
精确解改进欧拉法欧拉隐式欧拉显式节点 xi xey 30??
1.0000
?2.0000
4.0000
?8.0000
1.6000?101
?3.2000?101
1.0000
2.5000?10?1
6.2500?10?2
1.5625?10?2
3.9063?10?3
9.7656?10?4
1.0000
2.5000
6.2500
1.5626?101
3.9063?101
9.7656?101
1.0000
4.9787?10?2
2.4788?10?3
1.2341?10?4
6.1442?10?6
3.0590?10?7
What is wrong!
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
§ 8.5 差分方程的绝对稳定性
对于一般的差分方程
??
?
??
?
? ?
k
j
jnjnj
k
j
jnj yxfbhya
00
),(
由初始误差产生了差分解的误差，实际上是同一差分方程，取不同初值所得到的 2组
差分解之间的差。这个差不仅于差分方程本身有关，而且与微分方程本身有关。如果
微分方程本身是不稳定，那就没理由要求这 2组解充分接近。因此，差分方程的稳定性
概念是建立在微分方程稳定的基础上的。把这个典型微分方程规定为：
仍然考虑最简单的模型，即只有初值产生误差，看看这个误差的传播。
??
?
?
?
?
??
0)(
)0( R e
yay
y
dx
dy ??
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
差分方程运用到如上的微分方程后，可以得到
0)R e(,
00
?? ??
?
?
?
? ??
k
j
jnj
k
j
jnj ybhya
对于给定的初始误差
110,,,?keee ?
，误差方程具有一样的形式
0)R e(,
00
?? ??
?
?
?
? ??
k
j
jnj
k
j
jnj ebhea
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
定义：差分方程称为绝对稳定的，若差分方程作用到微分方程
)0( R e ?? ?? ydxdy
时，对任意的初值，总存在左半复平面上的一个区域，当 在这个区域时，差分
方程的解趋于 0。这个区域称为稳定区域 h?
例：向后 Euler公式的稳定性
11 ?? ?? nnn hyyy ?
误差方程：
11 ?? ?? nnn heee ?
he
e
n
n
??
??
1
11 1? 210 Re
Img
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
考察隐式欧拉法
11 ?? ?? iii yhyy ?
1
1
1iiyyh??
??? ??
???
1
10
1
1
i
i h???
?
?
??? ??
???
可见绝对稳定区域为：
| 1 | 1h???
210 Re
Img
注,一般来说，隐式欧拉法的绝对稳定性比同阶的显式
法的好。
数 学 系
University of Science and Technology of China
DEPARTMENT OF MATHEMATICS
3阶 Runge－ Kutta
? ?? ?
2
1
1
1
3211
1
)
2
1
1(
]4[
6
hhyK
hyK
yK
KKK
h
yy
n
n
n
nn
???
??
?
???
??
?
????
?
? ? ? ? ?
?
?
??
? ????
?
32
1 6
1
2
11 hhhyy
nn ???
显式 1~ 4 阶方法的绝对稳定区域为
k=1
k=2
k=3
k=4
-1-2-3
-
-
-
1
2
3
Re
Img
