第八章 常微分方程的数值方法
( Numerical Methods for Ordinary
Differential Equations )
常微分方程分为
( 1)初值问题( 8.1节 )
( 2)边值问题( 8.2节 )
一阶常微分方程初值问题的一般形式是,
0
(,),
( 1 )
()
{ (,),}
y f x y a x b
y a y
D x y a x b c y d



8.1 初值问题的数值方法称 f(x,y)在区域 D上对 y满足
Lipschitz条件是指,
1 2 1 2
12
0,,
(,) (,),
[,],,[,]
L s t
f x y f x y L y y
x a b y y c d



利用 Picard逼近容易证明,
Th8.1.1 若 f(x,y)在区域 D上连续,且对 y满足 Lipschitz条件,则初值问题 (1)在 [a,b]上存在唯一的连续可微解 y.
利用 Gronwall不等式易证解连续依赖于初值条件,
()
1 2 1 2
T h 8,1,2 (,)
( ; ) ( ; ),
L x a
f x y D
y x s y x s e s s


设 在 上 连 续,且 对 y 满 足 Lipschitz
条 件,若 y(x;s) 是 初 值 问 题
y = f ( x,y ),a < x < b
y(a)=s
的 解,则 有定理 8.1.2的意义在于:若初值问题 (1),
(2)中的初始值有一微小扰动,则解的扰动也是微小的,也就是解连续依赖于初始条件.通常将具有这种特性的初值问题称为是 适定的.
(稳定的)
数值解和精确解用数值方法求解初值问题,不是求出它的解析解或其近似解析式,而是给出它的解在某些离散节点上的近似值用 y(x)表示问题 (1),(2)的准确解
y(x0),y(x1),y(xN) 表示解 y(x)在节点 x0,x1,…,
xN处的准确值
y0,y1,…,y N表示数值解,即问题 (1),(2)的解 y(x) 在相应节点处的近似值,
单步法和多步法单步法:在计算 yi+1 时只利用 y i
多步法:在计算 yi+1 时不仅利用 y i,还要利用
yi?1,yi?2,…,
k步法:在计算 yi+1 时要用到 yi,yi?1,…,yi?k+1
显式计算公式 可写成,yk+1=yk+hΦ f(xk,yk;h)
隐式格式,yk+1=yk+hΦ f( xk,yk,yk+1;h)
它每步求解 yk+1需要解一个隐式方程一,Euler方法
0 1 2 1
0
,
,,1,2,,.
NN
j
a x x x x x b
ba
x x j h h j N
N


8.1.2 单步法
Euler方法是一种最简单的单步法
00
1
( ),
(,),0,1,,1i i i i
y y x
y y h f x y i N?


用差商近似导数
h
xyxyxy )()()( 01
0

x0 x1
),()()()( 000001 yxfhyxyhxyxy
)1,...,0(),(1 niyxfhyy iiii
继续这一过程,得到
1 0 0 0(,)y y h f x y
记为从而得到求解初值问题 (1),(2)的公式
Euler公式局部截断误差
1
1
(,,),
( ) [ ( ) (,( ),) ]
i i i i
i i i i
y y h x y h
y x y x h x y x h


i+1
对于数值方法局部截断误差定义为:
e
定义 在假设 yi = y(xi),即第 i 步计算是精确的前提下,
考虑的截断误差 Ri = y(xi+1)? yi+1 称为 局部截断误差
/* local truncation error */。
假定,yi = y(xi)”称为 局部化假定估计局部截断误差的主要方法是 Taylor展开法定义 若某算法的局部截断误差为 O(hp+1),则称该算法有 p阶精度 。
Euler方法的局部截断误差
2
2
3
11 2
3
1 2
1
( ) ( ) ( ) ( ) ( ),,
3!
( ) ( ) (,( ) ) ( ) ( )
h
i i i i i i i i
h
i i i i i
y x y x h y x y x h y x x
y x y x h f x y x y x O h



2
1 1 1
3
2
()
[ ( ) ( ) ( ) ( ) ] [ (,) ]
i i i
h
i i i i i i
e y x y
y x h y x y x O h y h f x y


)()( 32 2 hOxy ih 欧拉法具有 1 阶精度。
二,改进的 Euler方法
1
1
11
1
1 1 1
( ) ( ) (,( ) )
[ (,( ) ) (,( ) ) ]
2
(,),
[ (,) (,),
2
0,1,,1
i
i
x
ii
x
i i i i
i i i i
i i i i i i
y x y x f x y x d x
h
f x y x f x y x
y y h f x y
h
y y f x y f x y
iN







因为得到改进的 Euler方法的局部截断误差
2
1
3
(,)
( ) ( ) (,( ) ) [ (,( ) )
2
(,( ) ) (,( ) ) ] ( )
x y x y
i i i i x i i
y i i i i
y f x y
y f f y f f f
h
y x y x h f x y x f x y x
f x y x f x y x O h




1
2
2
3
[ (,) (,) (,)
2
(,) (,) ( ) ]
(,) [ (,) (,) (,) ]
2
( ),
i i i i i i x i i
i i y i i
i i i x i i i i y i i
h
y y f x y f x y h f x y
h f x y f x y O h
h
y h f x y f x y f x y f x y
Oh



对于改进的E u l e r 方法
8.1.2 一阶常微分方程初值问题的
Runge-Kutta方法考虑一阶常微分方程初值问题
0
(,),,
( ),
y f x y a x b
y a y


1
1
( ) ( ) (,( ) )
( ) (,( ) )
i
i
x
ii
x
i
y x y x f x y x d x
y x h f y


1
( ) (,( ) )
n
i j i j i j
j
y x h b f x c h y x c h

1
1
n
j
j
b

将区域 [a,b]进行分划:
0,,1,2,.,,,.jx a x a jh j N
ba
h
N

称为步长若则
1 1 1
2 2 2
(,( ) )
(,( ) )
......
(,( ) ),
ii
ii
n i n i n
k f x c h y x c h
k f x c h y x c h
k f x c h y x c h



(,( ) )
(,( ) ( ) )
(,( ) (,( ) ) )
j i j i j
i j i j i
i j i j i i
k f x c h y x c h
f x c h y x hc y x
f x c h y x hc f x y x



1
1
(,)
j
i j i jm m
m
f x c h y h a k
-1
1
j
acjm j
m

1
1
1
1
1
,
(,),
(,),2,3,.,,,
n
i i j j
j
ii
j
j i j i jm m
m
y y h b k
k f x y
k f x c h y h a k j n



n级显式 Runge-Kutta方法
1
1
j
jm j
m
ac

二级 Runge-Kutta方法取 n=2
1 1 1 2 2
1
2 2 2 1
( ),
(,),
(,),
ii
ii
ii
y y h b k b k
k f x y
k f x c h y c hk


(,),(,),(,),
(,),(,),(,),
i i x x i i y y i i
x x x x i i y y y y i i x y x y i i
f f x y f f x y f f x y
f f x y f f x y f f x y


记由此得
1
2 2 2 1
2
21
23
1 1 2 2 2 1
,
(,)
(,) ( ) ( ),
( ) ( ) ( )
ii
i i x y
i i x y
kf
k f x c h y c h k
f x y c h f k f O h
y y b b h f c b h f k f O h



另一方面
2
3
1
2
3
(,),,
( ) ( ) ( ) ( ) ( )
2
( ) ( )
2
i i x y
i i i i
i x y
y f x y f y f f f
h
y x y x hy x y x O h
h
y hf f f f O h



为使局部截断误差为,应取
12
22
1
1
2
bb
cb

3()Oh
改进的 Euler方法

2 1 2
1,1
2
b b c
1 1 2
1
21
( / 2 / 2),
(,),
(,),
0,1,...,1.
ii
ii
ii
y y h k k
k f x y
k f x h y hk
iN



中点方法取
2 1 2
11,0,
2
b b c
12
1
21
,
(,),
(,),
22
0,1,...,1.
ii
ii
ii
y y hk
k f x y
hh
k f x y k
iN



二阶 Heun方法取
2 1 2
3 1 2,,
4 4 3
b b c
1 1 2
1
21
13
( ),
44
(,),
22
(,),
33
0,1,..,,1,
ii
ii
ii
y y h k k
k f x y
hh
k f x y k
iN



二级 Runge-Kutta方法不超过二阶


2,2
x y x x x y y yF f f f G f f f f f
2 2 3
2 2 2
1 ()
2
k f c h F c h G O h
因此局部截断误差只能达到
2 2 3 4
1 1 2 2 2 2 2
1( ) ( ),
2ii
y y b b h f b c h F b c h G O h
23
4
1
23
4
( ) ( ) ( ) ( ) ( ) ( )
2 ! 3!
( ) ( ),
26
i i i i i
iy
hh
y x y x hy x y x y x O h
hh
y hf F G f F O h


3( ).Oh
三级 Runge-Kutta方法取 n=3 1 1 1 2 2 3 3
1
2 2 2 1
3 3 31 1 32 2
3 31 32
( ),
(,),
(,),
(,),
.
ii
ii
ii
ii
y y h b k b k b k
k f x y
k f x c h y c h k
k f x c h y a h k a h k
c a a





2,2
x y x x x y y yF f f f G f f f f f
2 2 2 1
2 2 3
22
(,)
1
()
2
iik f x c h y c hk
f c hF c h G O h


31 1 32 2
2 2 3
3 2 32 2 32
1
()
2
a k a k
c f c a hF c a h G O h

3 3 31 1 32 2
3 31 1 32 2
22
3 3 31 1 32 2
23
31 1 32 2
2 2 3
3 2 32 3
(,)
[ ( ) ]
1
[ 2 ( )
2
( ) ] ( )
1
( ) ( )
2
ii
xy
x x x y
yy
y
k f x c h y a h k a h k
f h c f a k a k f
h c f c a k a k f
a k a k f O h
f c h F h c a F f c G O h





又由于
2
1 1 2 3 2 2 3 3
3 2 2 4
3 2 32 2 2 3 3
( ) ( )
1
[ 2 ( ) ] ( ),
2
ii
y
y y b b b hf b c b c h F
h b c a Ff b c b c G O h


23
4
1( ) ( ) ( ),26i i y
hhy x y h f F G f F O h

因此要使局部截断误差为 O(h4),必须
1 2 3
2 2 3 3
22
2 2 3 3
3 2 3 2
1,
1 / 2,
1 / 3,
1 / 6,
b b b
b c b c
b c b c
b c a



Kutta方法取
1 2 3 2 3
3 2 3 1
1 / 6,2 / 3,1 / 6,1 / 2,1,
2,1
b b b c c
aa

1 1 2 3
1
21
3 1 2
1 2 1
( ),
6 3 6
(,),
11
(,),
22
(,2 ),
ii
ii
ii
ii
y y h k k k
k f x y
k f x h y h k
k f x h y h k h k



三阶 Heun方法

1 2 3 2
3 3 2 3 1
1 / 4,0,3 / 4,1 / 3,
2 / 3,2 / 3,0
b b b c
c a a


1 1 3
1
21
32
13
( ),
44
(,),
11
(,),
33
22
(,),
33
ii
ii
ii
ii
y y h k k
k f x y
k f x h y hk
k f x h y hk



三级 Runge-Kutta方法不超过三阶
完全类似于二级 Runge-Kutta方法的分析将 和 都展开到 项 易证 三级
Runge-Kutta方法的局部截断误差只能达到
4( ).Oh
1iy? 1()iyx? 4h
,
,
四级 R-K方法取 n=4
1 1 1 2 2 3 3 4 4
1
2 2 2 1
3 3 31 1 32 2
4 4 41 1 42 2 43 3
3 31 32 4 41 42 43
( ),
(,),
(,),
(,),
(,),
,.
ii
ii
ii
ii
ii
y y h b k b k b k b k
k f x y
k f x c h y c hk
k f x c h y a hk a hk
k f x c h y a hk a hk a hk
c a a c a a a





2 2 3
3 1 3 2 3 4 1 4 2 4 3 4
1 2 3 4 2 2 3 3 4 4
2 2 2 3 3 3
2 2 3 3 4 4 2 2 3 3 4 4
3 2 3 2 4 2 4 2 3 4 3
2 2 2
3 3 2 4 4 2 4 3
3 2 3 3 2 4 4 2 4 2 3 4 3
42
,,
1,1 / 2,
1 / 3,1 / 4,
( ) 1 / 6,
( ) 1 / 1 2,
( ) 1 / 8,
a a c a a a c
b b b b b c b c b c
b c b c b c b c b c b c
b c a b c a c a
b c a b c a c a
b c c a b c c a c a
b c a






3 2 4 3
1 / 2 4,a?
经典 R-K方法
1 1 2 3 4
1
21
32
43
( 2 2 ) / 6,
(,),
( / 2,/ 2),
( / 2,/ 2),
(,),
ii
ii
ii
ii
ii
y y h k k k k
k f x y
k f x h y hk
k f x h y hk
k f x h y hk




局部截断误差为 O(h5)
附注
二阶 Runge-Kutta方法的局部截断误差只能达到
三阶 Runge-Kutta方法的局部截断误差只能达到
四阶 Runge-Kutta方法的局部截断误差只能达到
五阶 Runge-Kutta方法的局部截断误差只能达到
3( );Oh
4( );Oh
5( );Oh
5( ).Oh