Transcript 偏微分方程

第九章
常微分方程数值解法
微分方程
 常微分方程( ODEs 未知函数是一元函数)
dv
c
 g v
dt
m
 偏微分方程( PDEs
未知函数是多元函数)
u
u
 2u
u
 2
t
x
x
 2  2

0
2
y
x
常微分方程
 一阶方程：如
dv
c
 g v
dt
m
 二阶方程 ：如
d2x
dx
m 2  c  kx  0
dt
dt
y' ( x )  2 x的解,
y( x )  x 2  c
方程加上适当的定解条件,才能确定方程的解.
定解条件 :
(1)初始值条件
( 2)边值条件
同一个微分方程,
具有不同的初始
条件
考虑一阶常微分方程初值问题
(Initial Value Problems )的数值解
 dy
a xb
  f ( x , y( x ))
 dx
问题的解y( x )是[a , b]上的函数
 y(a )  y0
数值解是求y( x )在一些离散点上的近似值
即给定节点: a  x0  x1  x 2    x N  b,
求上述问题的解y( x )在x n ( n  1,2,  N )处
的近似值yn .
假定 : (1) f ( x , y )连续, 且关于y满足Lipschitz条件
由此保证初值问题的解存在唯一.
| f ( x , y )  f ( x , y ) | L | y  y |
(2)hn  xn1  xn (n  0,1,, N  1)为常量h.
a  x  b 求y ( k  1,2, N )
 y  f ( x , y )
n

 y ( a )  y0
后边节点值 前边节点值 斜率  步长
单步法一般形式: yn1  yn  n h
 yn 1  yn  hf ( x n , yn )

 y0  y ( a )
 yn1  yn  hf ( x n1 , yn1 )

 y0  y ( a )
h

 y n  1  y n   f ( x n , y n )  f ( x n  1 , y n  1 )

2
 y0  y(a )
 y  f ( x , y )
微分方程离散化的方法 
 y ( a )  y0
(1) 用差商近似导数
y( xn )  f ( xn , y( xn ))
a xb
y( x n 1 )  y( x n )
y( xn ) 
h
y( x n  1 )  y( x n )
 f ( xn , y( xn ))
h
y( xn1 )  y( xn )  hy( xn )  y( xn )  hf ( xn , y( xn ))
 yn 1  yn  hf ( x n , yn )

 y0  y ( a )
差分方程初值问题
Euler方法
 y  f ( x , y )

 y ( a )  y0
a xb
若用向后差商近似导数，即
y
(
x
)

y
(
x
)
n

1
n
y( xn1 )  f ( xn1 , y( xn1 )) y( xn1 ) 
h
y( x n  1 )  y( x n )
 f ( xn1 , y( xn1 ))
h
y( xn1 )  y( xn )  hf ( xn1 , y( xn1 ))
 yn1  yn  hf ( x n1 , yn1 )

 y0  y ( a )
向后Euler方法

（2）用数值积分方法
x n 1
xn
y( x n 1 )  y( x n )  
x n 1
xn
y( x )dx  
x n 1
xn
f ( x , y( x ))dx
f ( x , y( x ))dx
对积分用点xn 处的矩形公式得
y( xn1 )  y( xn )  f ( xn , y( xn )( xn1  xn )
yn1  yn  hf ( xn , yn )
若对积分用梯形公式，则得
h
y( xn1 )  y( xn )   f ( xn , y( xn ))  f ( xn1 , y( xn1 ))
2
h

 y n  1  y n   f ( x n , y n )  f ( x n  1 , y n  1 )

2
 y0  y(a )
梯形公式
（3）用Taylor多项式近似
y( xn1 )  y( xn  h)  y( xn )  hy( xn )  
yn1  yn  hf ( xn , yn )
Euler方法
§1
Euler方法
1.Euler方法
以差分方程初值问题
 yn1  yn  hf ( xn , yn )

 y0  y ( a )
n  0,1,, N  1
的解作为微分方程初值问题的数值解，即
y( xn )  yn , 称为Euler方法.
单步法一般形式: yn1  yn  n h
用f ' ( xn )代替n
dy
 y   f ( x , y );
dx
y(x 0 )  y 0
用分段折线逼近曲线
y0
x0
h
x1
h
x2
h
x3
例 : 用欧拉方法求下列方程的数值解
 dy

 x y,
 dx

 y(0)  1,
0 x1
分别取h  0.5
h  0.25
解：Euler公式为
yn1  yn  hf ( xn , yn )  yn  hxn yn
n  0,1,
当h=0.5时
 y(0.5)  y1  y0  hx0 y0  1  (0.5)(0 1 )  1.0

 y(1.0)  y2  y1  hx1 y1  1  (0.5)( 0.5 1 )  1.25
yn1  yn  hxn yn
(n  0,1,)
当h=0.25时
y(0.25)  y1  y0  0.25 x0 y0  1  (0.25)(0 1 )  1
y(0.5)  y2  y1  0.25 x1 y1
 1  (0.25)( 0.25 1 )  1.0625
y(0.75)  y3  y2  0.25 x2 y2
 1.0625  (0.25)( 0.5 1.0625 )  1.191347
y(1.0)  y4  y3  0.25 x3 y3  y(0.75)  hf (0.75, 1.191347)
 1.191347  (0.25)( 0.75 1.191347)  1.39600
dy
 y  x y ;
dx
y(0 )  1
解析解为: y  (1  x 2 / 4)2
h = 0.25
h = 0.5
1
0
0.25
0.5
0.75
1.0
1.2
Euler方法的误差估计
局部截断误差：Rn1  y( xn1 )  yn1，其中yn  y( xn ).
整体截断误差: en1  y( xn1 )  yn1 .
dy
 y  f ( x , y );
dx
y
x0
h
y ( x 0 )  y0
x1
h
x2
h
x3
定义 : 若某种数值方法的局部截断误差
Rn1  O( h P 1 ), 则称该方法是 p 阶方法.
若Rn1   ( xn , y( xn ))h P 1  O( h P  2 ),
则称 ( xn , y( xn ))h P 1为局部截断误差的主项.
对Euler方法，局部截断误差
1 2
y( xn1 )  y( xn  h)  y( xn )  hy' ( xn )  h y( )
2
xn    xn1
yn1  yn  hf ( xn , yn )  yn  hy' ( xn )
 Rn1  y( xn1 )  yn1
1 2
 h y( )
2
xn    xn 1
Euler方法为1阶方法，局部截断误差： Rn 1  1 h2 y( )
2
Euler方法的整体截断误差
记：yn1  y( xn )  hf ( xn , y( xn ))
hM
L( b  a )
e n  1  y( x n  1 )  y n  1 

e
 1  O( h)
2L
若f ( x, y )充分光滑, 使得y( x )  C a, b, 则 M  0, 使得 | y( x ) | M
( 2)
向后Euler方法 :
 yn1  yn  hf ( xn1 , yn1 ) n  0,1,

 y0  y ( a )
这是隐式公式, 一般需用迭代法求解
 yn( 0)1  yn  hf ( xn , yn )
 ( k 1)
 yn1  yn  hf ( xn , yn( k)1 ) k  0,1,
局部截断误差:
Rn1  y( xn1 )  yn1
1 2
  h y( xn )  O( h3 )
2
§2
改进Euler方法
2.1 梯形公式
h

 y n  1  y n   f ( x n , y n )  f ( x n  1 , y n  1 )

2
 y0  a
h3
截 断 误 差: Rn1  y( xn1 )  yn1   y( )
12
xn    xn1
梯形公式为二阶方法, 且为隐式格式
 yn( 0)1  yn  hf ( xn , yn )
求解公式 
 ( k 1)
h
(k )
y

y

f
(
x
,
y
)

f
(
x
,
y
(k  0,1,2,)
 n1
n
n n
n1 n1 )
2



h

 y n  1  y n   f ( x n , y n )  f ( x n  1 , y n  1 )

2
 y0  a
求截断误差
h2
h3
h4 ( 4 )
y( xn1 )  y( xn )  hy( xn )  y( xn )  y( xn )  y (1 ) xn  1  xn1
2
6
24
若yn  y( xn ),
h
yn1  y( xn )  ( f ( xn , y( xn ))  f ( xn1 , yn1 ))
2
h
 y( xn )  ( f ( xn , y( xn ))  f ( xn1 , y( xn1 )))
2
h
 y( xn )  ( y( xn )  y( xn1 ))
2
h
yn1  y( xn )  ( y( xn )  y( xn1 ))
2
h2
h3
h4 ( 4 )
 y( xn )  hy( xn )  y( xn )  y( xn )  y ( 2 ) xn   2  xn1
2
4
12
h2
h3
h4 ( 4 )
y( xn1 )  y( xn )  hy( xn )  y( xn )  y( xn )  y (1 )
2
6
6
xn  1  xn1
h3
 Rn1  y( xn1 )  yn1   y( xn )  O( h4 )
12
或用梯形公式的误差
Rn1  y( xn )  
x n 1
xn
h
h3
f ( xn , y( xn ))  y( xn )   f ( xn , yn )  f ( xn1 , yn1 )   y( )
2
12
xn    xn1
2.2 改进Euler法
 yn1  yn  hf ( xn , yn )


h
 yn1  yn  2  f ( xn , yn )  f ( xn1 , yn1 )
预测
校正
称为Euler公式与梯形公式的预测—校正系统。
实际计算时，常改写成以下形式

 y p  yn  hf ( xn , yn )

 yq  yn  hf ( xn1 , y p )

 y n  1  1 ( y p  yq )

2
yi01  yi  hf ( xi , yi )
predictor


1
f ( xi , yi )  f ( xi 1 , yi01 )
2


h
yi 1  yi  f ( xi , yi )  f ( xi 1 , yi01 )
2
corrector
 y'  y
例 : ( p.265 13)写出初值问题
改进的Euler法
 y ( 0)  1
近似解的表达式, 并证明当h  0时, 收敛于真解.
解 : 改进的Euler法计算公式
 yn1  yn  hf ( xn , yn )


h
 f ( xn , yn )  f ( xn1 , yn1 )
y

y

n
 n1
2
预测
校正
 yn1  yn  hyn

2
h
h

yn1  yn  ( yn  yn  hyn )  yn  hyn 
yn

2
2

h2
h2 n
yn  yn1 (1  h  )  y0 (1  h  )
2
2
2
x
h
令x  nh
2
h
h
lim yn  lim (1  h  )  lim (1  h  )
h 0
h 0
h 0
2
2
2 h h2

2h
2 h h2
2x
 ex

例 : 用改进的Euler方法计算初值问题  y p  yn  hf ( xn , yn )

2
 yq  yn  hf ( xn1 , y p )
 y'   y  y sin x ,


1

 y(1)  1,
 yn1  2 ( y p  yq )
的解在1.2和1.4处的近似值, 取步长h  0.2,
小数点后至少保留5位.
解 : 改进的Euler方法
 yn1  yn  hf ( xn , yn )


h
 yn1  yn  2  f ( xn , yn )  f ( xn1 , yn1 )
预测
校正
f ( x , y )   y  y 2 sin x
y p  yn  h(  yn  yn2 sin x n )
y p  yn  h(  y p  y 2p sin x n1 )
yn1
1
 ( y p  yq )
2
y(0)  y0  1, h  0.2
y p  yn  h(  yn  yn2 sin x n )
y p  yn  h(  y p  y 2p sin x n1 )
1
y n  1  ( y p  yq )
2
n  0时, y p  0.631706
n  0,1
y p  0.799272
y(1.2)  y1  0.71549
n  1时, y p  0.476964
y(1.4)  y2  0.52611
y p  0.575259
y(0)  y0  1
§3
Runge—Kutta法
3.1 : 基本思想
以函数f ( x , y )在若干点上的函数值的线性组合构造公式,
组合系数与点的选取要使近似公式在( x n , yn )处的Taylor
展开与解y( x )在x n 在的Taylor展开式的前p  1项重合,
从而得到p阶公式.
单步法一般形式:
Runge  Kutta :
yn1  yn   ( xn , yn , h) h
  c1 K1  c2 K 2   cn K n
1. RK方法的构造
RK公式的一般形式为
yn1  yn   ( xn , yn , h)h
 ( xn , yn , h)  c1 K 1  c2 K 2    cn K n
 K1

 K2

 K3


 K n
 f ( xn , yn )
 f ( xn  a2 h, yn  b21 K 1h)
 f ( xn  a3 h, yn  b31 K 1h  b32 K 2 h)

 f ( xn  an h, yn  bn ,1 K 1h  bn , 2 K 2 h    bn ,n1 K n1h)
其中 ai , bij , ci 为待定系数.
以p  2为例, RK公式的几何意义:
 y n  1  y n  h( c1 K 1  c 2 K 2 )

 K 1  f ( x n , yn )
 K  f ( x  a h, y  hb K )
n
2
n
21
1
 2
K2
K1
c1 K1 + c2 K2
加权平均斜率
xn
xn+a2h
xn+1
= xn+h
Taylor展开公式:
h2
h3
y( x  h)  y( x )  hy( x ) 
y( x ) 
y( x )  
2
6
f ( x  h, y  k )  f ( x , y )  h
f
f
( x, y)  k ( x, y)
x
y
2
2
2

f

f

f
2
2
h
( x , y )  2hk
( x, y )  k
( x, y)  
2
2
x
xy
y
复合函数求导链式法则:
y( x )  f ( x , y( x ))
y( x )  f ( x , y( x )) 
y( x )  f ( x , y ) 
f f dy f f



f  f x  ff y
x y dx x y
  f f dy    f f dy  dy

  



x  x y dx  y  x y dx  dx
以p  2为例, RK公式的推导:
 y n  1  y n  h( c1 K 1  c 2 K 2 )

 K 1  f ( x n , yn )
 K  f ( x  a h, y  hb K )
n
2
n
21
1
 2
yn1  yn  hc1 f ( xn , yn )  c2 f  xn  a2 h, yn  hb21 f ( xn , yn )
c1 f ( xn , yn )

 yn  h
2 


 c2 f ( xn , yn )  a2 hf x ( xn , yn )  b21hf ( xn , yn ) f y ( xn , yn )  O(h ) 




 yn  (c1  c2 )hf ( xn , yn )  c2h2 a2 f x ( xn , yn )  b21 f ( xn , yn ) f y ( xn , yn )  O(h3 )


yn1  yn  (c1  c2 )hf ( xn , yn )  c2h2 a2 f x ( xn , yn )  b21 f ( xn , yn ) f y ( xn , yn )  O(h3 )
h2
y( xn1 )  y( xn )  hy( xn )  y( xn )  O( h3 ) 由yn  y( xn )得
2
h2
 yn  hf ( xn , yn ) 
f x ( xn , yn )  f y ( xn , yn ) f ( xn , yn )  O(h3 )
2

为使y( xn1 )  yn1  O( h3 ),


c1  c2  1

1

c 2 a 2 
2

1

c 2 b21 

2


c1  c2  1

1

c
a

 2 2
2

1

c 2 b21 

2

1
1
特别地, 令c2  , 则c1  , a2  b21  1
2
2
h
yn1  yn   f ( xn , yn )  f ( xn  h, yn  hf ( xn , yn ))
2
这就是改进Euler公式，故其为二阶方法。
h

 yn1  yn  2 ( K 1  K 2 )

 K 1  f ( xn , yn )
 K  f ( x  h, y  hK )
n
n
1
 2


 y p  yn  hf ( x n , yn )

 yq  yn  hf ( x n1 , y p )

 y n  1  1 ( y p  yq )

2
由

c1  c 2  1

1

c 2 a 2 
2

1

c 2 b21 

2

1
令c2  1, 则c1  0, a2  b21  ,
2
则得如下中点公式
K2

 yn1  yn  hK 2

 K 1  f ( x n , yn )

h
h
 K 2  f ( xn  , yn  K 1 )
2
2

K1
xn
xn+1/2
xn+1
三阶Runge  Kutta公式的一般形式：
 K 1  f ( xn , yn )

 K 2  f ( x n  a 2 h, yn  b21 K 1 h)

 K 3  f ( x n  a 3 h, yn  b31 K 1 h  b32 K 2 h)
 yn1  yn  (c1 K 1  c2 K 2  c3 K 3 )h
  c1 K1  c2 K 2  c3 K 3
1
y n  1  y n  ( 2 K 1  3 K 2  3 K 3 )h y n  1 
8


 K 1  f ( x n , yn )
 K1 


2
2


 K 2  f ( x n  h, yn  K 1h)
K2 
3
3


2
2


 K 3  f ( x n  3 h, yn  3 K 2 h)
 K 3 
1
y n  ( 2 K 1  3 K 2  4 K 3 )h
9
f ( x n , yn )
1
1
f ( xn  h, yn  K 1h)
2
2
3
3
f ( xn  h, yn  K 2 h)
4
4
经典三阶RK公式（右端类似Simpson公式）
h

 yn1  yn  6 ( K 1  4 K 2  K 3 )

 K 1  f ( xn , yn )

K  f ( x  h , y  h K )
n
n
1
 2
2
2
 K  f ( x  h, y  hK  2hK )
n
n
1
2
 3
经典四阶RK公式
h
yn1  yn  ( K 1  2 K 2  2 K 3  K 4 )
6
 K 1  f ( xn , yn )

 K 2  f ( xn  h , yn  h K 1 )

2
2

K  f ( x  h , y  h K )
n
n
2
 3
2
2
 K  f ( x  h, y  hK )
n
n
3
 4
经典四阶RK公式的几何意义
K2
K4
K3
K1
xn

1
  ( K1  2K 2  2K 3  K 4 )
6
xn + h/2
xn + h
说明： 当 p  4时, RK公式的最高阶数为p;
当 p  4时, RK公式的最高阶数小于p.
4阶Runge-Kutta方法
dy
  y  sin x ; y(0 )  1
dx
Euler
Euler_mod
midpoint
RK4
n = 10
Euler
Euler_mod
midpoint
RK4
n = 20
例 : 证明Runge  Kutta公式
h

yn1  yn  ( k 2  k 3 )

2

k1  f ( x n , y n )


k 2  f ( x n  th, yn  thk1 )

 k 3  f ( x n  (1  t )h, yn  (1  t )hk1 )
对任意参数t至少是二阶公式.
证明: 精确解y( xn1 )和数值解yn1展开为f ( xn , yn )及其各阶偏导的形式
y( xn1 )在xn 处Taylor展开
h2
yn  y( x n )
y( xn1 )  y( xn )  hy( xn )  y( xn )  O( h3 )
2
h2
 yn  hf ( xn , yn )   f1( xn , yn )  f 2( xn , yn ) f ( xn , yn )  O( h3 )
2
h

y

y

(k2  k3 )
n
 n1
2
k  f ( x , y )
n
n
 1
 k 2  f ( x n  th, yn  thk1 )

 k 3  f ( x n  (1  t )h, yn  (1  t )hk1 )
k2  f ( xn  th, yn  thk1 )
 f ( xn , yn )  thf '1 ( xn , yn )  thk1 f '2 ( xn , yn )  O( h2 )
 f ( xn , yn )  thf1 ' ( xn , yn )  thf ( xn , yn ) f 2 ' ( xn , yn )  O( h2 )
k3  f ( xn  (1  t )h, yn  (1  t )hk1 )
 f ( xn , yn )  (1  t )hf1 ' ( xn , yn )  (1  t )hk1 f 2 ' ( xn , yn )  O( h2 )
 f ( xn , yn )  (1  t )hf1 ' ( xn , yn )  (1  t )hf ( xn , yn ) f 2 ' ( xn , yn )  O( h2 )
k2  f ( xn , yn )  thf1 ' ( xn , yn )  thf ( xn , yn ) f 2 ' ( xn , yn )  O( h2 )
k3  f ( xn , yn )  (1  t )hf1 ' ( xn , yn )  (1  t )hf ( xn , yn ) f 2 ' ( xn , yn )  O( h2 )
于是, 得到yn1在( xn , yn )处的展开式:
h
yn1  yn  ( k 2  k 3 )
2
h
 yn  2 f ( xn , yn )  hf1 ' ( xn , yn )  hf 2 ' ( xn , yn ) f ( xn , yn )  O( h2 )
2



h
yn1  yn  2 f ( xn , yn )  hf1 ' ( xn , yn )  hf 2 ' ( xn , yn ) f ( xn , yn )  O( h2 )
2

h2
y( xn1 )  yn  hf ( xn , yn )   f1( xn , yn )  f 2( xn , yn ) f ( xn , yn )  O( h3 )
2
y( x n  1 )  y n  1  O( h 3 )
对任意t上述Runge  Kutta方法至少是二阶公式.
注 : 事实上，可以证明局部截断误差h3的系数不为零，
是二阶公式
§4
线性多步法
单步法的一般形式：yn1  yn  h ( xn , yn , yn1 , h)
显式为: yn1  yn  h ( xn , yn , h)
其中称为增量函数.
r
r
i 0
i  1
多步法一般形式: yn1    i yn i  h   i f n i
其中f n i  f ( xn i , yn i ),  i ,  i 为待定系数.
若 1  0 为显式公式,  1  0 为隐式公式.
显式线性三步法的几何意义 :
2
2
i 0
i 0
yn1    i yn i  h  i f n i
fn
f n 1
f n 2
 0 f n   1 f n 1   2 f n  2
 0 y n   1 y n 1   2 y n  2
xn-2
xn-1
xn
xn+1 = xn+h
4.1 线性多步公式的导出
r
r
i 0
i  1
多步法一般形式: yn1    i yn i  h   i f n i
利用Taylor展开
p
h2
h
yn i  y( xn  ih)  yn  ( i )hyn  ( i )2 h2 yn    ( i ) p yn( p )  O(h p1 )
2
p!
h p 1 ( p )
f n i  y( xn i )  yn  ( i )hyn   
y n  O( h p )
( p  1)!
h p 1 ( p )
 f ( xn1 , y( xn1 ))  y( xn1 )  yn  hyn   
y n  O( h p )
( p  1)!
约定
f n1
h2
h p ( p)
y( xn1 )  yn  hyn 
yn   
y n  O( h p  1 )
2
p!
r
r
i 0
i  1
yn1    i yn i  h   i f n i
 yn1
r
 r

   i yn  h  (  i ) i    i  yn
i 0
i  1
 i 1

r
r

h2  r
2

(  i )  i  2  (  i ) i  yn  


2  i 1
i  1

r
 ( p)
hp  r
p
p 1
p 1

(

i
)


p
(

i
)

y

O
(
h
)


i
i n

p!  i 1
i  1

 r
  i  1
 i 0
 r
r
 (  i ) k   k (  i ) k 1   1


i
i

i  1
 i 1
()
( k  1,2,  , p)
p阶公式的局部截断误差:
r
 ( p 1)
h p 1  r
p 1
p
p 2
Rn1 
1

(

i
)


(
p

1
)
(

i
)

y

O
(
h
)


i
i n

( p  1)!  i 1
i 1

结论：r  1步公式至多可以达到2r  2阶精度.
例 : 讨论二步公式
h
y n  1  y n 1  ( f n  1  4 f n  f n 1 )
3
是几阶公式, 并求出它的局部截断误差.
h
y n  1  y n 1  ( f n  1  4 f n  f n 1 )
3
解 : 精确解y( xn1 )和数值解yn1展成y( xn )及其各阶导数形式
yn1  y( xn  h)
h2
h3
h4 ( 4 ) h5 ( 5 )
 yn  hyn  yn  yn  yn  yn  O( h6 )
2
3！
4!
5!
f n1  f ( xn1 , y( xn1 ))  y( xn1 )  y( xn  h)
h2
h 3 ( 4 ) h4 ( 5 )
 yn  hyn  yn  yn  yn  O( h5 )
2!
3！
4!
f n1  f ( xn1 , y( xn1 ))  y( xn1 )  y( xn  h)
h2
h 3 ( 4 ) h4 ( 5 )
 yn  hyn  yn  yn  yn  O( h5 )
2!
3！
4!
h
y n  1  y n 1  ( f n  1  4 f n  f n 1 )
3 1


1 1

2



yn1  yn  hyn   1  (1  4  1)  h yn   (1  1)
3


2 3

1 1 1 
 1 1 1 1 
3
4 (4)  1



 h y n    (  )  h y n   (  )
 3! 3 2 2 
 4! 3 3! 3! 
 1 1 1 1 
 h y    (  )  O ( h 6 )
 5! 3 4! 4! 
h2
h3
h4 ( 4 )
7 5 (5)
 yn  hyn 
yn 
yn 
yn 
h y n  O ( h6 )
2
6
4!
360
5
(5)
n
y( xn1 )  y( xn  h)  yn  hy
(1)
n
y( x n1 )  yn 1
h 2 ( 2 ) h 3 ( 3 ) h4 ( 4 ) h5 ( 5 )

yn 
yn 
yn 
yn  O( h6 )
2
6
4!
5!
1
7
( 
)h5 yn( 5 )  O( h6 )
5! 360
h5 ( 5 )

y n  O ( h6 )
90
4.2 常用的线性多步公式
1. Adams公式
h
yn1  yn  (55 f n  59 f n1  37 f n 2  9 f n 3 )
24
Rn1
251 5 ( 5 )

h y n  O ( h6 )
720
四阶Adams显式公式
h
yn1  yn  (9 f n1  19 f n  5 f n1  f n 2 )
24
19 5 ( 5 )
Rn1  
h y n  O ( h6 )
720
四阶Adams隐式公式
（二）Milne公式
4
yn1  yn 3  h(2 f n  f n1  2 f n 2 )
3
Rn1
14 5 ( 5 )

h y n  O ( h6 )
45
（三）Hamming公式
1
3
yn1  (9 yn  yn 2 )  h( f n1  2 f n  f n1 )
8
8
1 5 (5)
Rn1   h yn  O( h6 )
40
说明
1.线性多步公式不能自启动，
一般需用同阶单步法求得初值后
再用线性多步公式计算；
2. 一般地，同阶隐式公式比显式公式精确，
但隐式公式计算复杂，需用迭代法求解。
4.3 预测—校正系统
Adams预测—校正公式：
h

y

y

 n1 n 24 (55 f n  59 f n1  37 f n 2  9 f n 3 )

 y  y  h 9 f ( x , y )  19 f  5 f  f 
n1 n1
n
n 1
n 2
 n1 n 24
Milne—Hamming预测—校正公式：
4

 yn1  yn 3  3 h( 2 f n  f n1  2 f n 2 )

 y  1 (9 y  y )  3 h f ( x , y )  2 f  f 
n1
n1
n
n 1
 n1 8 n n 2 8
§5
相容性、收敛性与稳定性
5.1 相容性与收敛性
原方程：
 y  f ( x , y )

 y ( x 0 )  y0
()
 yn1  yn  h ( xn , yn , h)
单步法离散方程： 
 y ( x 0 )  y0
()
问题：
1. 当 h  0 时，差分方程是否无限逼近微分方程？
2. x  (a , b], 当 h  0时，问题（* *）的解
是否无限逼近问题（
*）的解？
相容性:
如果增量函数  ( x, y, h) 关于 h 连续且满足条件
 ( x , y,0)  f ( x, y )
则称单步法与问题（*）相容，
也称问题（**）与（*）相容。
结论：若 ( x , y , h)足够光滑(连续，对h偏导数存在有界)
则单步法为 p 阶 ( p  1)的充要条件问题 (*) 相容。
必要性 :
若 yn1  yn  h ( xn , yn , h) 为p阶方法
则 y( xn1 )  y( xn ) - h ( xn , y( xn ), h)  O( h p 1 )
也即 y( x  h)  y( x ) - h ( x, y( x ), h)  O( h p1 )
y( x  h)  y( x )
  ( x , y( x ), h)  O( h p )
h
两边取极限,由的连续性得
 ( x , y( x ),0)  lim  ( x , y( x ), h)  y( x )  f ( x , y( x ))
h 0
收敛性 : 如果某种数值方法对任意x  a , b,
lim yn  y( x )
x  a  nh
h 0
则称该数值方法是收敛的.
定理9.1 : 设单步法具有p阶精度, 其增量函数
 ( x , y , h)关于y满足Lipschitz条件, 问题(*)的
初值是精确的, 则单步法的整体截断误差为
en 1  y( x n 1 )  yn 1  O( h p )
定理9.2 : 设增量函数 ( x , y , h)在区域S 上连续,
S:
a  x  b,
   y  ,
0  h  h0
且关于y满足Lipschitz条件, 则单步法收敛的充要
条件是相容性条件成立.
参考李立康等
收敛性
《微分方程数值解》
定理9.1
相容性
结论
p阶方法( p  1)
结论 : 若f ( x , y )在区域D : a  x  b,  y  
连续, 且关于y满足L条件, 则Euler, 改进Euler
和各阶RK方法都是收敛的.
以上方法均为( p  1)方法,
且增量函数满足定理9.1的条件.
注 对形式简单的方程，可以由差分方程解的表达式
取极限导出收敛性。
x0
 y  ay
例如对初值问题： 
 y ( 0)  1
用Euler法得近似解表达式
yn  yn1  hayn1  (1  ha ) yn1  (1  ha )
对 x  nh， 当 h  0 时有
x
h
lim yn  lim (1  ha )  lim (1  ha )  e
n
h 0
h 0
ax
h 0
容易验证 y  e ax 是初值问题的解。
n
5.2 稳定性
 y( x )  30 y( x )
例 : 考虑初值问题
在区间[0,0.5]上的解.
 y ( 0)  1
分别用欧拉显式, 欧拉隐式和改进欧拉方法求解.
节点
欧拉显式
xi
0.0
1.0000
0.1
2.0000
0.2
4.0000
0.3
8.0000
0.4
1.6000101
0.5
3.2000101
欧拉隐式
改进欧拉法
精确解y  e 30 x
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
定义 : 若算法在计算过程中任一步产生的误差
在以后的计算中都逐步衰减, 则称该算法
是绝对稳定的.
一般分析时, 为了简单, 只考虑试验方程
y   y
为常数, 可以是复数
当步长取为h时, 将某算法应用于上式(试验方程),
假设yn 有误差 , 则若此误差以后逐步衰减, 则称
该算法相对于h  h绝对稳定, h 的全体构成稳定区域.
若算法A的稳定区域比算法B的稳定区域大,
则称算法A比算法B稳定.
特别地, 若绝对稳定区域包含左半平面, 则称
该方法是A  稳定的.
讨论各种方法的稳定的
.
Im( h)
1.显式Euler方法 :
用Euler方法求解试验方程: y   y
-2
-1
yn1  yn  hyn  (1  h) yn
若yn有舍入误差 , 则
yn1  (1  h)( yn   )  yn1  (1  h)
 n1  (1  h)
 n k  (1  h)k 
故Euler法的绝对稳定区域为： 1  h  1
0
Re(h)
2.隐式Euler方法 :
用隐式Euler方法求解试验方程: y   y
Im( h)
yn1  yn  hyn1
yn1
1

yn
1  h
故隐式Euler法的绝对稳定区域为：
0
1
2
1  h  1
注：一般来说，隐式欧拉法的绝对稳定性比同阶
的显式法的好。
Re(h)
3.改进的Euler方法 :
h
yn 1  yn  [ f ( x n , yn )  f ( x n 1 , yn 1 )]
2
h
 yn  [ f ( x n , yn )  f ( x n 1 , yn  hf ( x n , yn ))]
2
将改进Euler法用于试验方程，则有
2 h 2
h
) yn
yn1  yn  yn   ( yn  hyn )  (1  h 
2
2
2 h 2
 n1  (1  h 
)
2
故改进Euler法的绝对稳定区域为：
1  h 
2 h 2
2
1
4.梯形公式:
梯形公式用于模型方程则为
yn 1
h
 yn  (yn  yn 1 )
2
故其绝对稳定区域为
h
h
1
 1
2
2
1
yn1 
h
2 1
h
1
2
Re(h)  0
因此梯形公式是A―稳定的。
1
h
2 y
h n
1
2
Im( h)
0
Re(h)
对四阶RK方法，有
h
yn1  yn  ( K 1  2 K 2  2 K 3  K 4 )
6
h
2 h
K 1  yn K 2   ( yn  K 1 )  ( 
) yn
2
2
h
2 h 3 h 2
K 3   ( yn  K 2 )  ( 

) yn
2
2
4
3 2
4 3

h

h
2
K 4   ( yn  hK 3 )  (   h 

) yn
2
4
3 2
3 2
4 3
h

h

h

h
2
2
2
yn1  yn  (  2   h  2   h 
    h

) yn
6
2
2
4

(h)2 (h)3 (h)4 
 1  h 


 yn
2
6
24 

其绝对稳定区域为
(h)2 (h)3 (h)4
1  h 


1
2
6
24
(1)当为实数时, Euler方法, 改进的Euler方法,
梯形公式, 四阶RK法的绝对稳定区间分别为
 2  h  0
   h  0
(h)2
 2  h 
0
2
 2.78  h  0
f
( 2)对一般常微分方程y  f ( x , y ),  
y
( 3)只有当步长h的选取使得h在绝对稳定区域内时,
计算结果才能稳定; 否则, 数值计算是不稳定的.
例：有解初值问题
 y  f ( x , y )

 y ( x 0 )  y0
h
h
的龙格  库塔公式: yn1  yn  hf ( xn  , yn  f ( xn , yn ))
2
2
(1)证明该公式是二阶公式；
( 2)求该公式的绝对稳定的实数区间；
( 3)若f ( x , y )  y , y(0)  1, h  0.2求y(0.2)的近似值
解：
(2)将上述公式用于模型方程
 y  y

 y ( x 0 )  y0
h
( h )2
yn1  yn  h ( yn  yn )  yn (1  h 
)
2
2
( h )2
[2，
0]
稳定性区域：
| 1  h 
| 1 得实数区间：
2
( 3)1.22
§6一阶微分方程组和高阶微分方程的数值解法
6.1 一阶微分方程组的数值解法
一阶微分方程组的初值问题的一般形式为：
 dy1
 dx  f1 ( x , y1 , y2 ,  , ym )



 dym  f ( x , y , y ,  , y )
m
1
2
m
 dx

初值 : y1 ( x0 )  y10 , y2 ( x0 )  y20 , ym ( x0 )  ym 0
6.1.1 : 用显式Euler公式解微分方程组初值问题 :
以两个方程为例
 dy1
 dx  f1 ( x , y1 , y2 )
 dy
 1  f1 ( x , y1 , y2 )
 dx
初值y1 ( x0 )  y10 , y2 ( x0 )  y20

 y1,n1  y1,n  hf1 ( x n , y1,n , y2 ,n )


 y2 ,n1  y2 ,n  hf 2 ( x n , y1,n , y2 ,n )
6.1.2 : 用经典4阶RK方法解微分方程组初值问题 :
以两个方程为例
 k1,1 

 k1, 2 

 k 2 ,1 

k 
 2 , 2

 k 3 ,1 

k 
 3 , 2
 k4 ,1 

 k4 , 2 
f1 ( x n , y1,n , y2 ,n )
f 2 ( x n , y1,n , y2 ,n )
h
h
h
, y1, n  k1,1 , y2, n  k1,2 )
2
2
2
h
h
h
f 2 ( x n  , y1, n  k1,1 , y2, n  k1,2 )
2
2
2
h
h
h
f1 ( x n  , y1, n  k 2,1 , y2, n  k 2,2 )
2
2
2
h
h
h
f 2 ( x n  , y1, n  k 2,1 , y2, n  k 2,2 )
2
2
2
f1 ( x n  h, y1,n  k 3,1h, y2, n  k 3,2h)
f1 ( xn 
f 2 ( x n  h, y1,n  k 3,1h, y2, n  k 3,2h)
1

y

y

( k1 , 1  2 k 2 , 1  2 k 3 , 1  k 4 , 1 ) h
1, n
 1,n1
6

1
y

y

( k1 , 2  2 k 2 , 2  2 k 3 , 2  k 4 , 2 ) h
2,n
 2 ,n1
6
6.2高阶微分方程的数值解法
 y ( n )  f ( x , y, y, y,  , y ( n1) )

 y( x0 )   0 , y( x0 )   1 ,  , y ( n1) ( x0 )   n1
化为一阶微分方程组:
 y1  y
y1 ( x0 )  α0
 y1  y2 ,

 

y2 ( x0 )  α1
 y2  y
 y2  y3 ,


令  y3  y   y3  y4 ,
y3 ( x0 )  α2
 
 


 yn  y ( n1)
 yn  f ( x,y1 ,y2 , ,yn ) , yn ( x0 )  αn1