第 6章
常微分方程数值解法
绪论
在工程和科学计算中，所建立的各
种常微分方程的初值或边值问题，除很
少几类的特殊方程能给出解析解，绝大
多数的方程是很难甚至不可能给出解析
解的，其主要原因在于积分工具的局限
性。因此，人们转向用数值方法去解常
微分方程，并获得相当大的成功，讨论
和研究常微分方程的数值解法是有重要
意义的。
6.1 初值问题的 Euler方法
,.,, )2,1,0()(
,)(
0,.,, ),2,1,0(
)(
),(
00
??
????
?
?
?
?
?
?
?
nyxy
yxxy
hnnhax
yxy
yxf
dx
dy
nn
nn
n
即上的近似值在点去获得解
值解是指通过某种方法假定为常数。该式的数
为步长，一般总记
问题设一阶常微分方程初值
初值问题的 Euler方法
。算出
出发，逐次公式，它可以从这就是显式的
）（
的近似值，则有表示以
于是该式可离散为：代替散化，并用
方法首先将微分算子离为实现这一目标，
...,,
` E u le r
)1(,.,,2,1,0),(
)(
))(,(
)()(
,
E u le r
321
0
1
0
yyy
y
nyxhfyy
xyy
xyxf
h
xyhxy
xx
nnnn
nn
nn
nn
n
???
?
??
?
初值问题的 Euler方法
。
才能得到步要解函数方程的不同在于，它每算一
方法，它与显式公式或向后这就是隐式的
）（
的近似值，则有表示以
于是该式可离散为：代替如果用
1
111
11
01
)2(
E u le r` E u le r
)2(,.,,2,1,0),(
)(
))(,(
)()(
,
?
???
??
?
???
?
??
n
nnnn
nn
nn
nn
n
y
nyxhfyy
xyy
xyxf
h
xyhxy
xx
初值问题的 Euler方法
)1(
11
)()1(
1
)(
11
)1(
1
)0(
1
111
,||
3
,.,,2,1,0)],(),([
2
),(
,.,,2,1,0)],(),([
2
?
??
?
?
??
?
?
?
???
???
?
?
?
?
?
????
??
????
k
nn
k
n
k
n
k
nnnnn
k
n
nnnn
n
nnnnnn
yyyy
kyxfyxf
h
yy
yxhfyy
y
nyxfyxf
h
yy
取时当
）（
）（
时常用以下迭代式：计算
称为梯形公式。
）（
平均值的结果，则得如果取以上两式的算术
?
初值问题的 Euler方法
收敛。
）产生的序列时，则由（即
满足常数。如果步长为条件，
满足对变量设函数定理
...)2,1,0}({3
2
,1
2
0L ip s c h it z
L ip s c h it zy),(1.1.6
)(
1
??
??
?
ky
L
h
hL
hL
yxf
k
n
初值问题的 Euler方法
。故有由假设知
）有）和（证明：由式（
1
)1(
1
1
)0(
11
1
)(
11
)(
1111
1)(k
1n1n
lim,0)
2
(lim:
||)
2
(
.,,,,,
||
2
|),(),(|
2
|yy|
32
?
?
?
??
?
??
??
?
??
????
?
??
??
??
??
???
n
k
n
k
k
k
nn
k
k
nn
k
nnnn
yy
hL
yy
hL
yy
hL
yxfyxf
h
初值问题的 Euler方法
。并取
即校正算法构成一类预估迭代一次
一般只由于迭代工作量较大计算对于
)c(
11
)(
11
)(
`
)(
1
1
)],(),([
2
),(
,,
,,)2(
??
???
?
?
?
?
?
?
?
?
???
??
?
nn
p
nnnnn
c
n
nnn
p
n
n
yy
yxfyxf
h
yy
yxhfyy
y
初值问题的 Euler方法
))],(,()([
2
,E u le r
,.,, )2,1,0(),(
),(
)(
2
1
11
12
1
211
nnnnnnnn
nn
nn
nn
yxhfyxfyxf
h
yy
nkyhxhfk
yxhfk
kkfyy
?????
?
?
?
?
?
?
?
????
?
???
??
?
亦可写成方法该式称为改进
上式还常写成
初值问题的 Euler方法
进行比较。精确解
并与法求解法和改进试分别用
设初值问题例
xy
y
y
x
y
dx
dy
21
,E u le rE u le r
1)0(
2
1.1.6
??
?
?
?
?
?
?
??
初值问题的 Euler方法
计算结果如下表所示：
法：改进的
法：
上结果，此时计算解：取
?
?
?
?
?
?
?
?
?
?
?
?
???
??
???
????
??
?
?
,.,, )2,1,0()
)1.0(2
(1.0
)
2
(1.0
)(
2
1
E u le r
,.,, )2,1,0()
2
(1.0
]1,0[,1.0
1
12
1
211
1
n
ky
x
kyk
y
x
yk
kkyy
n
y
x
yyyE u l e r
xh
n
n
n
n
n
n
nn
n
n
nnn
x Euler法 y 改进的 Euler法 y 精确解
0 1.000000 1.000000 1.000000
0.1 1.000000 1.095909 1.095445
0.2 1.191818 1.184097 1.183216
0.3 1.277438 1.266201 1.264911
0.4 1.358213 1.343360 1.341641
0.5 1.435133 1.416402 1.414214
0.6 1.508966 1.485956 1.483240
0.7 1.580338 1.552514 1.549193
0.8 1.649783 1.616475 1.612452
0.9 1.717779 1.678166 1.673320
1.0 1.784770 1.737867 1.732051
6.1.2 误差概述
有关，称为增量函数。与函数
为而隐式单步法一般形式
显式单步法一般形式为
),(
),,,,(
),,(
111
1
yxf
hyxyxhyy
hyxhyy
nnnnnn
nnnn
?
?
?
???
?
??
??
误差概述
上的局部截断误差。称为点则
由单步法计算出值计算没有误差，即
的步在点果第上的整体截断误差。如称为在点
则的值或隐式逐步计算，得
出发，由单步法显式从初值定义
11
~
11
1
~
1
11111
00
,)(
,),(
)(,
)(6, 1, 1
????
?
?
?????
??
?
??
?
nnnn
nnn
nn
nnnnn
xyxyT
yxyy
xnx
yxyeyx
yxy
误差概述
称 为为局部截断误差主))(,(则非零项
)())((
阶精度。如果则称该算法具有
误差为
截断如果 给果给定的算法的21..6定 义
1
21
1
1
1
?
??
?
?
?
??
?
p
nn
pp
nnn
p
n
hxyxg
hOhx,yxgT
p
)O(hT
误差概述
)()(
2
)O()(y
2
E u le r式对于
)O ( h)(
2
E u le r对于显式
则精度相对的越高。
越大，误差阶,一个算法,局部截断一般
4
3
1
3
2
1
3
2
1
hOxy
h
T
hx
h
T
xy
h
T
ｐ
nn
nn
nn
??????
?????
????
?
?
?
：对于梯形公式
：法隐
：法
说来
6.1.3 数值稳定性分析
值不大,则，则称该算法的舍入 误的舍入误差
。 如果算则称该算法是数值稳定
||||
，得
值方法又算,如果结果为
计算,则时有一舍入误差假 设
1
111
数
，也成绝对稳定
且
通过某种数
实际计算
nn
nnn
~
nnn
~
nn
ρρ
ρyy
ρyy
ρy
?
??
??
?
???
数值稳定性分析
? 定义 6.1.3 若某数值算法的绝对稳定性区
域包含 hλ平面上的左半平面 Re(hλ)＜ 0，
则称该方法是 A
? 隐式 Euler法是 A
6.2 Runge-Kutta方法
。这里仍假定
，使局部截断误差，，，适当选择参数
可设为方法启发，更一般算式受改进的
)(),()(
,.,, )2,1,0(
),(
),(
E u le r
3
111
21
12
1
22111
nnnnn
nn
nn
nn
xyyhOyxyT
n
kyhxhfk
yxhfk
kkyy
????
?
?
?
?
?
?
???
?
???
???
?
????
??
??
Runge-Kutta方法
)()),(),(),(()(
)(),(),(),(
:T a y lo r
32
3
1
2
2
hOyxfyxfyxfhxyh
hOyxfhkyxfhyxhfk
nnynnnnxn
nnynnxnn
???????
??????
??
??
展开式由二元函数
)()],(),(
),([)()()(
32
2
2211
hOhyxfyxf
yxfxyhxyy
nnynn
nnxnnn
???
??????
?
??
????
于是
Runge-Kutta方法
一。自由参数，即解答不唯
程，因此有一个由于四个参数，三个方
展式相比较得与
?
?
?
?
?
?
?
?
?
?
?
??
2
1
2
1
1
:T a y lo r
2
2
21
??
??
??
方法。这是改进的
此时算式为可得取
E u le r
),(
),(
)(
2
1
,1,
2
1
,
2
1
)1(
12
1
211
21
?
?
?
?
?
?
?
???
?
???
????
?
kyhxhfk
yxhfk
kkyy
nn
nn
nn
????
Runge-Kutta方法
.K-R
)
2
1
,
2
1
(
),(
,
2
1
,1,0)2(
12
1
21
21
方法这是二阶
此时算式为可得取
?
?
?
?
?
?
?
???
?
??
????
?
kyhxhfk
yxhfk
kyy
nn
nn
nn
????
方法。这也是二阶
（
又有算式可得取
K-R
)
3
2
,
3
2
(
),(
)3
4
1
,
3
2
,
4
3
,
4
1
)3(
12
1
211
21
?
?
?
?
?
?
?
???
?
???
????
?
kyhxhfk
yxhfk
kkyy
nn
nn
nn
????
6.2.2 四阶 Runge-Kutta方法
434241432313212
34324214143
23213133
12122
1
443322111
,.,, )2,1,0(
),(
),(
),(
),(
KR
?????????
????
???
??
????
??????
?
?
?
?
?
?
?
?
?
?
?????
????
???
?
?????
?
?
，，约定
分方程初值问题可设为方法类似，解一阶常微与二阶
n
kkkyhxhfk
kkyhxhfk
kyhxhfk
yxhfk
kkkkyy
nn
nn
nn
nn
nn
四阶 Runge-Kutta方法
稳定。但不是该公式结构整齐
算式求得经典四阶类似以上推导过程，可
A,
),(
)
2
1
,
2
1
(
)
2
1
,
2
1
(
),(
)22(
6
1
KR
33
23
12
1
43211
?
?
?
?
?
?
?
?
?
?
?
???
???
???
?
?????
?
?
kyhxhfk
kyhxhfk
kyhxhfk
yxhfk
kkkkyy
nn
nn
nn
nn
nn
6.2.3 R-K法的稳定性
?
?
?
?
?
?
?
?
?
????
???
??
?
?
n
n
n
n
yhhhhk
yhhhk
yhhk
yhk
y
x
y
))(
4
1
)(
2
1
)((
))(
4
1
)(
2
1
(
))(
2
1
(
d
d
,K-R
432
4
32
3
2
2
1
????
???
??
?
? 有应用于实验方程法为例以经典
R-K法的稳定性
nn
n
nn
hhh
h
y
hhh
h
kkkkyy
?
???
??
???
?
]
!4
)(
!3
)(
!2
)(
)(1[
:
]
!4
)(
!3
)(
!2
)(
)(1[
)22(
6
1
432
1
432
43211
?????
?????
?????
?
?
相应的误差为
于是
R-K法的稳定性
.
,,,
,078.2:,
1|
!4
)(
!3
)(
!2
)(
)(1|
:|,|||
432
1
稳定性必须很小才能保证算法
即限制很大步长的绝对值较大时当因此
可得稳定性区间为负实数时当
于是得稳定性区域为由稳定性要求
h
h
y
f
h
hhh
h
nn
?
?
?
???
?????
?
?
?
??
???
?
??
6.2.5 隐式 R-K法
.,
,,.,
A,,
E u le r
.
.,||,
,,K-R
1
1
又能保证数值稳定精度提高
既可以构造高阶隐式方法因此限制较小即对步长定的
稳它们均是但从稳定性分析知迭代求解是它们的缺点
需要多次法和梯形方法计算隐式格式的向后
否则数值不稳定
有较大的限制对步长较大时当但从稳定性分析知
使用方便可以直接计算出法的优点是从显式
h
y
h
y
f
yy
n
nn
?
?
?
?
隐式 R-K法
?
?
?
??
??
?
?
?
?
?
????
????
???
?
?
22212
12111
22212122
21211111
22111
),(
),(
KR
???
???
???
???
??
约定
方法形式为：设二阶隐式
kkyhxhfk
kkyhxhfk
kkyy
nn
nn
nn
隐式 R-K法
?
?
?
?
?
?
?
?
?
????
??
??
??
?
6
1
)()(
3
1
2
1
1
T a y lo r),(
T a y lo r),(
22221121221111
2
22
2
11
2211
21
1
21
??????????
????
????
??
的同次幂系数得：比较
展式的在再反复带入上式，并与
，展开求得的做二元函数在点，将
h
yxy
kyxkk
nnn
jnn
隐式 R-K法
?
?
?
?
?
?
?
???
???
???
??
?
)
3
2
,
3
2
(
),(
)3(
4
1
K-R
3
2
0)1(
22
211
211
21
kyhxhfk
kkyxhfk
kkyy
nn
nn
nn
方法，得二阶三级隐式，令 ??
隐式 R-K法
?
?
?
?
?
?
?
?
?
??????
??????
???
????
?
)
4
1
)
6
3
4
1
(,)
6
3
2
1
((
))
6
3
4
1
(
4
,)
6
3
2
1
((
)(
2
1
K-R
6
3
2
1
6
3
2
1
)2(
212
2
1
1
211
21
kkyhxhfk
k
k
yhxhfk
kkyy
nn
nn
nn
方法得二阶四级隐式，令 ??
6.3 线形多步法
? 单步法主要依据 yn的信息去计算 yn+1 。线性
多步法是想依据 yn,yn-1…, yn-r (r≥1)的信息
去计算 yn+1。
? 考虑到线性组合较为方便，因此，线性多步法
一般形式可设为
,,,0
).,.,,,1,0,1)(,(,
)...(...
1
10
0111101
否则称为隐式称为显式时当
其中
?
???
??
????????
?
???
??
?
?
?
??????
??
?
??
??????
rkyxff
fhy
fffhyyyy
knknkn
r
k
knh
r
k
knk
rnrnnrnrnnn
6.3.1 基于数值积分的方法
?
?
??
?
??
?
?
?
?
?
?
???
1
0
0
))(,()()(
))(,()()(
)(
),(
01
1
0
00
x
x
x
x
dttytfxyxy
xx
dttytfxyxy
yxy
bxayxf
dx
dy
得取
的等价积分方程形式为
初值问题
基于数值积分的方法
公式。这就是著名的
故有只能是近似值由于每一步得到的
略去误差项有形公式对右端定积分采用左矩
有对于一般区间
E u le r
,.,, )2,1,0(),(
,)(
))(,()()(
,)1(
))(,()()(
],[
1
1
1
1
1
1n
???
??
??
?
?
?
?
?
?
?
nyxhfyy
xy
xyxhfxyxy
dttytfxyxy
xx
nnnn
n
nnnn
x
x
nn
nn
n
有望提高精度。
则值多项式代替被积函数）如果对右端用高次插（
。这就是熟知的梯形公式
（
略去误差项有公式对右端定积分采用梯形
,3
,.,, )2,1,0()),(),(
2
,)2(
111
????
???
nyxfyxf
h
yy
nnnnnn
基于数值积分的方法
??
?
??
???
?
?
?
?
?
?
?
?
? ?
??
??
????
?
??
??
?
?
?
????
1
0
1
0
1
)1(2
11
0
11
rr11
)1,,.,,,1,0(
!
)1),,, (1(
)1(
),(),()(
)()()(
:N e w t o n
),,(),.,,,,(),,(1
rrjds
j
jsss
ds
j
s
b
xxfhbxT
xTfbhxyxy
fxfxfxr
j
j
nrn
rr
r
r
j
n
j
jnn
nnnnnn
??其中
后插值公式代入得到
由个点：取
基于数值积分的方法
算式时称为四阶显式
则），，，并逐次取式中的代替用
Ad a m s3
)(
7 2 0
2 5 1
)(
)9375955(
24
)(
8
3
)()51623(
12
)(
12
5
)()3(
2
)(
2
)(
,3210(),(
)5(5
1
3211
)4(
4
1211
)3(
3
111
''
2
11
j
?
?
?????
?????
????
???
??
????
???
??
?
j
yhxT
ffff
h
yy
y
h
xTfff
h
yy
y
h
xTff
h
yy
y
h
xThfyy
jfxyy
nnnnnn
nnnnn
nnnn
nnn
nnn
?
?
?
?
基于数值积分的方法
.A d a m s
)(
7 2 0
19
)5199(
24
:,
),(),.,,,,(),,(:1
)5(5
2111
1r1r11
算式称为四阶隐式
则可同上求得次插值多项式做
个点若取
?yh
ffff
h
yy
r
fxfxfxr
nnnnnn
nnnnnn
?
?????
?
????
??????
基于数值积分的方法
? Adams预估 — 校正法
? 预估
? 校正
? 并取
)9375955(24 321)( 1 ???? ????? nnnnnpn ffffhyy
]519),(9[24 21)( 11)( 1 ????? ????? nnnpnnncn fffyxfhyy
)( 11 cnn yy ?? ?
6.3.2 基于 Taylar展开式的方法
)...(...
,3),.,,,1,0,1(),(,
)...(...
33011331101
10
0111101
??????
???
??
?
?
?
??????
????????
????
??
????????
??
nnnnnnn
knknkn
r
k
knh
r
k
knk
rnrnnrnrnnn
fffhyyyy
rrkyxff
fhy
fffhyyyy
??????
??
??????
即设中取其中
式在线形多步法的一般形
)3,2,1,0,1()()(
!4
)(
)(
!3
)(
)(
!2
)(
)()()(
)()())(,(
)3,2,1,0()()(
!5
)(
)(
!4
)(
)(
!3
)(
)(
!2
)(
)()()(
)()(
,T a y lo r
)())(,()(4
5)5(
4
)4(
3
)3(
2
'''
''
6)5(
5
)4(
4
)3(
3
''
2
'
45
0
????
????
???
????
????
??
?
???
???
?
????
khOxy
kh
xy
kh
xy
kh
xykhxy
khxyxyxyxf
khOxy
kh
xy
kh
xy
kh
xy
kh
xykhxy
khxyxy
hhxx
xyxyxfxyp
n
nnnn
nknknkn
nn
nnnn
nkn
knknknkn
即项和处展到展式在用
及阶精度，将为使算法达到
基于 Taylar展开式的方法
?
?
?
?
?
?
?
?
?
???????
????????
???????
?????????
????
????
?
?
?
?
??
?
11 0 832448116
1271233278
1642294
132
1
,t a y lo r)(,
),())(,(),(
3211321
3211321
3211321
32101321
3210
11
1
???????
???????
???????
????????
????
可得到方程组
的系数展式比较的并与展开式则得到
，并令将其代入
j
nn
nnnnnnn
hxyy
xyxyxffxyyy
基于 Taylar展开式的方法
自由参数。
个方程，因此有四个个参数由于方程组是有
局部截断误差
59
)()(
!5
)4 0 5
8055-2 4 3321(
)(
6)5(
5
3
211321
111
hOxy
h
yxyT
n
nnn
??
??????
??
?
???
?
??????
算式和它的余式估计。阶显式这是
所以有
，，，，
可解得）取（
Ad a m s4
)()(
720
251
)9375955(
24
24
9
24
37
24
59
24
55
1
,0,01
6)5(5
1
3211
32100
3211
?
?
?
?
?
??
?????
???????
????
?
????
?
hOxyhT
ffff
h
yy
nn
nnnnnn
?????
????
算式和它的余式估计。阶隐式这是
所以有
，，，，
可解得）取（
Ad a m s4
)()(
7 2 0
19
)5199(
24
24
1
24
5
24
19
24
9
1
,0,02
6)5(5
1
2111
31010
3213
?
?
?
?
?
???
?????
??????
????
?
????
?
hOxyhT
ffff
h
yy
nn
nnnnnn
?????
????
?
?
?
?
?
??
????
?
?
?
?
?
??
????
?
????
?
????
)()(
45
14
)22(
3
4
M iln e4
)()(
90
1
)49(
3
算式si m p so n3
6)5(5
1
2131
6)5(5
1
1111
hOxyhT
fff
h
yy
hOxyhT
fff
h
yy
nn
nnnnn
nn
nnnnn
算式）（
）（
其他一些常见的算式：
?
?
?
?
?
?????
????
?
?
?
?
?
?
???
?????
?????
????
?
????
]2)),([
8
3
)9(
8
1
)22(
3
4
)6(
)()(
40
1
)2(
8
3
)9(
8
1
Ha m m in g5
1
)(
112
)(
1
213
)(
1
6)5(5
1
1121
nn
p
nnnn
C
n
nnnn
P
n
nn
nnnnnn
ffyxf
h
yyy
fff
h
yy
hOxyhT
fff
h
yyy
校正
预估
校正系统预估
算式）（
6.4 一阶常微分方程组数值解法
? 在许多实际问题中，常常出现高阶微分
方程和高阶微分方程组，通过引入新的
变量，总可化为一阶微分方程组。
? 由此可知，讨论一阶常微分方程组的数
值解法是很有意义的。
6.4.1 解一阶常微分方程组的 R-K方法
算式为四阶方法求解如果用经典的
时
设一阶方程组
,K-R
,,
),,(
d
d
),,(
d
d
2021010
212
2
211
1
?
?
?
?
?
?
?
?
?
???
?
?
yyyytt
yytf
t
y
yytf
t
y
一阶常微分方程组的 R-K方法
,.,, )2,1,0(
),,(),,,(
)
2
,
2
,
2
(),
2
,
2
,
2
(
)
2
,
2
,
2
(),
2
,
2
,
2
(
),,(),,(
)22(
6
1
)22(
6
1
232131221232131114
22
2
12
1223
22
2
12
1113
21
2
11
1222
21
2
11
1112
2122121111
24232221212
14131211111
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????????
????????
????????
??
?????
?????
?
?
n
kykyhthfkkykyhthfk
k
y
k
y
h
thfk
k
y
k
y
h
thfk
k
y
k
y
h
thfk
k
y
k
y
h
thfk
yythfkyythfk
kkkkyy
kkkkyy
nnnnnn
nnnnnn
nnnnnn
nnnnnn
nn
nn
，
。步求出
重复上述过程可逐，，，
的值；再以，时的，再计算出，
，，，，，代入上式，先计算
开始，将计算过程：从
,.,, )2,1,0(,
1
,,0
1211
2121111
211112414
2313122111202
1010
?
????
?
?
???
??
nyy
yyyyttn
yyttkk
kkkkkyy
yyttn
nn
一阶常微分方程组的 R-K方法
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
),.,,,2,1()(
),.,,,,(
d
d
.,,,,,
),.,,,,(
d
d
),.,,,,(
d
d
00
21
212
2
211
1
miyty
yyytf
t
y
yyytf
t
y
yyytf
t
y
m
ii
mm
m
m
m
值问题设一阶常微分方程组初
个方程的情形。到将两个方程的情形推广
一阶常微分方程组的 R-K方法
:K-R
,),.,,,,(
)(
),(
d
d
:
,),.,,,,(,),.,,,,(
T
020100
00
T
21
T
21
的向量形式是方法
于是经典的四阶，这里
向量形式为
则方程的若记
btayyy
t
t
t
fffyyy
m
mm
y
yy
yf
y
fy
???
?
?
?
?
?
?
?
??
一阶常微分方程组的 R-K方法
?
?
?
?
?
?
??
?
?
?
?
?
?
???
???
???
?
?????
?
,.,, )2,1,0
),(
)
2
,
2
(
)
2
,
2
(
),(
)22(
6
1
34
2
3
1
2
1
43211
n
hth
h
th
h
th
th
kyfk
k
yfk
k
yfk
yfk
kkkkyy
nn
nn
nn
nn
nn
（
,
),(
),(
),(
),(
.,,
1
12
11
2
1
1
2
1
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
mnnm
nn
nn
nn
mm
n
n
n
k
k
k
thf
thf
thf
th
y
y
y
y
y
y
yfky
??
，
这里向量
,
)
2
,
2
(
)
2
,
2
(
)
2
,
2
(
)
2
,
2
(
,
)
2
,
2
(
)
2
,
2
(
)
2
,
2
(
)
2
,
2
(
3
32
31
2
2
2
2
1
2
3
2
22
21
1
1
2
1
1
1
2
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
????
m
nnm
nn
nn
nn
m
nnm
nn
nn
nn
k
k
k
h
thf
h
thf
h
thf
h
t
k
k
k
h
thf
h
thf
h
thf
h
th
k
y
k
y
k
y
k
yfk
k
y
k
y
k
y
k
yfk
?
?
?
?
,.,, )2,1,0(),.,,2,1(
),.,,,,,,(
)
2
,.,,,
2
,
2
,
2
(
)
2
,.,,,
2
,
2
,
2
(
),.,,,,(
)2(
6
1
),(
),(
),(
),(
32321314
222
2
12
13
121
2
11
12
211
43211
4
42
41
3
32
31
34
??
?
?
?
?
?
?
?
?
?
?
?
?????
?????
?????
?
?????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
????
?
nmi
kykykyhthfk
k
y
k
y
k
y
h
thfk
k
y
k
y
k
y
h
thfk
yyythfk
kkkkyy
k
k
k
hthf
hthf
hthf
ht
mmnnnnii
m
mnnnnii
m
mnnnnii
mnnnnii
iiiiinin
m
nnm
nn
nn
nn
ky
ky
ky
kyfk
其分量形式是
??
一阶常微分方程组的 R-K算法;)2
);,.,,2,1(),.,,,,,(1
4,3,2,1)4(
);,.,,2,1(;;)3(;5.0;;;5.0;5.0)2(;),,.,,2,1(,,,)1(
21
54321
j
mii
iiii
ei
utt
mitfk
j
miyytt
huhuhuhuhu
tmiythm
??
??
?
????
?????
?
?
?
???
??
）
做对于
输入：
）。返回（停机
输出：
做对于
3e lset h e nif)6(
);,.,,2,1(,)5(;
3
2;1
,.,,,2,1)3
10
0
e
i
ij
ii
ijii
tt
miyt
ku
yy
ku
mi
?
?
??
??
?
?
??
一阶常微分方程组的 R-K方法
6.4.2 刚性方程组
?
?
?
???
?
???
??
m
j
j
t
j
jj
mm
mm
tect
mjm
ttttyyyR
tA
t
ry
r
y
y
y
j
1
T
21
T
21
)()(
:
),,.,,,2,1(A
))(),.,,,(),(()(),.,,,,(;A,
)(
d
d
?
?
????
?
则通解为
和相应的特征向量个互不相同的特征值有假设
。；其中
组设常系数线形微分方程
刚性方程组
即认为是刚性问题。
时当为微分方程组的刚性比称
称为慢瞬态解称为快瞬态解其中
为瞬态解称因此
且
即部的所有特征值都有非实假设
,1S.|
)(Re
)(Re
|S
,,,
).(0
|)(Re|.,,|)(Re||)(Re|
0)(Re,A
m
1
mm11
11
m21
m11
???
????
???
?
??
??
?
?
???
?
??
??
recrec
rectrec
tt
m
j
j
t
j
m
j
j
t
j
j
jj
刚性方程组
常收到较好的效果。
方法用隐式的才能保证其稳定性。应小时
只有当步长很法求解刚性问题用显式
,K-R,
,K-R
6.5 常微分方程边值问题的数值解法
? 设二阶线性常微分方程为
? 常见边界条件有三类,
],[),()()( baxxryxqyxpy ???????
1010
11
)()(,)()(
)(,)(
)(,)(
????
??
??
??????
????
??
bybyayay
byay
byay
6.5.1 差分方程的建立
。其中
。由数值微分公式
作等距分划对区间
),(,,
)(
12
)()(2)(
)(
)(
62
)()(
)(
),.,,2,1,0(:],[
11
)4(
2
2
11
)3(
2
11
??
??
??
?
?
??
???
?
?
??
?
?
???
jjjj
j
jjj
j
j
jj
j
j
xx
y
h
h
xyxyxy
xy
y
h
h
xyxy
xy
n
ab
h
njjhaxba
??
?
?
差分方程的建立
程由边界条件补出。
还缺的两个方个方程的这是求
得差分方程：代入
,1),.,,,2,1,0(
)1,.,,2,1(
2
2
],[),()()(
11
2
11
??
????
?
?
??
???????
????
nnjy
njryq
h
yy
p
h
yyy
baxxryxqyxpy
j
jjj
jj
j
jjj
差分方程的建立
由追赶法不难求解。其中
这时整理后有两个未知量的解
即已给出对于第一类边界条件：
,;
2
1
2;
2
1
,
,,
2
2
11
2
2
11
1
2
2
1
11
222
222
11
0
?
?
?
?
?
???
?????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
?
?
??
???
jjjj
jjjj
nn
n
n
n
nn
nnn
n
rhdp
h
c
qhbp
h
a
cd
d
d
ad
y
y
y
y
ba
cba
cba
cb
yy
?
?
??
?????
差分方程的建立
求解差分方程组联立，既可将上述两个方程分别与
具体是将边界条件也离散化，
式界条件，由数值微分公对于第二类，第三类边
?
?
?
?
?
??
??
??
???
?
?
?
?
?
?
??
?
???
??
??
10
12
100
210
1
12
1
210
2
34
2
43
2
34
2
43
??
??
?
?
n
nnn
nnn
y
h
yyy
y
h
yyy
h
yyy
h
yyy
算法
二阶常微分方程边值问题的差分算法
。）用追赶法解（
）（
）
）
）
做对
输入：
)1,.,,,2,1(5;;4;;
2
1;2;
2
13
);();();(2;1
1,.,,2,1)3(;)2(;,,,,)1(
111111
22
??
????
???????
???
?
??
?
?
???
njy
cddadd
rhdp
h
cqhbp
h
a
xrrxqqxpp
jhx
nj
n
ab
h
nba
j
nnn
jjjj
??
??
例题
? 例 6.5.1 用差分法求解边值问题
??
?
?
?
???
????????
2ln44)2(,2)1(
765
62 2
2
yy
xx
x
y
x
y
。此方程的解析解为解 xxxxxy ln21,2432 ????
例题
如下表：
果，由差分算法求解的结
，，
，，，
，，此时现取步长
2ln44
2),10,.,,,1,0(1.012
1765)(
6
)(
2
)(1.0
2
2
??
?????
??????
??
?
?jjxb
axxxr
x
xq
x
xph
j
x y y(x) e=y(x)-y
1.0 0.5 0.5 0
1.1 0.72798569 0.72637532 -1.610*10-3
1.2 1.0140390 1.0113430 -2.696*10-3
1.3 1.3678241 1.3644456 -3.378*10-3
…… …… …… ……
1.7 3.6896236 3.6865656 -3.058*10-3
1.8 4.5635316 4.5612288 -2.303*10-3
1.9 5.5854269 5.5841425 -1.284*10-3
2.0 6.7725887 6.7725887 0