第八章 常微分方程数值解法
§ 1,引言
0 0
:
(,)
()
dy
f x y a x b
dx
y yx
?
? ? ??
?
? ?
?
一阶常微分方程的初值问题
''
,x y - 2 y = 4 x y =
f ( x,y ) = y ( 1 ) = - 3
(,)
()
2y
4
x
2y
4
x
d y 2 y
f x y 4
d x x
y 1 3
??
?
?
? ? ? ?
?
? ??
?
例 方程
令,且给出初值
就得到一阶常微分方程的初值问题：
(,)
L
(,) (,)
f x y
f x y f x y L y y? ? ?
只要函数 适当光滑连续，且关于y 满足李普
希兹( L i p s c h i t z ) 条件，即存在常数,使得
由常微分方程理论知，初值问题的解必存在且唯一。
微分方程的数值解：设方程问题的解 y(x)的存在区间
是 [a,b]，令 a= x0< x1<… < xn =b,其中 hk=xk+1-xk,如是等距
节点 h=(b-a)/n,h称为步长。
y(x)的解析表达式不容易得到或根本无法得到，我们用数
值方法求得 y(x)在每个节点 xk上 y(xk)的近似值，用 yk表示，
即 yk≈ y(xk)，这样 y0,y1,...,yn称为微分方程的数值解。
微分方程离散化常用方法
? ?
? ? ? ?
1
1
,
A
(,( ) )
n n
nn
nn
nn
yy
f x y x
yx
dy xx
xxdx
?
?
?
?
??
?
用差商代替微商
? ? ? ?
? ?
? ? ?,2,1,0,
,
,,
1
1
1
1
1
???
?
?
????
?
?
?
?
?
nhf
f
h
yyh
yxyy
yx
yy
xyxyxx
n
n
nn
n
n
nn
n
n
n
n
nn
代替，则：用
11
B,
:
(,) ( 0,
1,)
nn
nn
xx
xx
dy
dx f x y dx n
dx
??
???? L
数值积分
用数值积分方法离散化
1
11
1
,( ),( ),
(,) (,)
( 0,1,(,))
n
n
n n n n
x
nn
x
n n n n
y y y x y x
f x y dx h f x y
y y hf x ny
?
??
?
?
???
?
L
用 代替 对右端积分采用
取左端点的矩形公式
则有
2
/ //
2
//
C ( )
( ) ( ) ( ) ( )
2
( ) (,( ) ) ( )
2
n
n n n n
n n n n
y x Tay lor
y h y h
y h f y
x
h
yyx x x x
h
yx x x x
?
? ? ? ? ?
? ? ? ?
L
L
在 附近 的 展开：
1
1
( ) ( )
(,) 0,1,2,
nnn
nn n n
h y y
h f n
Tay lor
y xx
y y yx
?
?
?
? ? ? L
取 的线性部分，且 得 的近似值：
展开法不仅可得到求数值解的公式，且容易估计
截断误差。
0 0
(,)
()
dy
f x y a x b
dx
y yx
???
? ? ????
?
? ?
???
§ 2 尤拉（ Eular)方法
? ?
1
0 0
0 0
,0 1,(,)
()()
nn n n
E u le r
dy
h f n,f x y a x b
dx
yy
y y yx
yxyx
?
? ?
? ? ?? ? ?
??
???
??? ?
??
L
一,公式
计算公式：
? ? ? ?
? ?
? ?
? ? ? ?
0 0
0
0
0
0
,
00
00
1
,
,
f,,
y y x
f
yx
yx
dy
yx
dx
yyxx
?
?
，几何意义。
由 出发取曲线 的切线（存在！），则 斜率
由于 及 已知，必有切线方程。
? ?
? ?
),(f)x(xy yxxy
dx
dy
xy
y,x
000000
00
??????
：由点斜式写出切线方程
）（
：，可由切线算出，则为等步长
,yxyy
yxx
hf
hh
0001
101
??
??
? 2 1 0
1
1
，，，）（
：，点值在）（逐步计算出
,???
?
?
?
nhf
xyy
yxyy
x
nnnn
n
注意：这是“折线法”而非“切线法”
除第一个点是曲线切线外，其他点不是！
y
χ 0 χ 1 χ 2 χ 3 χ
2,Euler方法的误差估计
? ?
1 1
1
1
1)
~
~
( )
n n
n
n
n
n
n
y
y
E ule r
xT
y
y x
y
y
? ?
?
?
??
?
局部截断误差。
在一步中产生的误差而非累积误差：
其中 （精确解！）时
由 法求出的值，即 无误差！
是当
? ?
? ? xxyh
xxxxx
xx
nn
nnnnn
nn
yhfyhyy
Tay l o ry
1
//
2
1
1
2
))(,()()(
)(
?
?
?
???
?????
??
展开：点将 在
? ?
? ?? ?
? ? ? ?
1
1
2
//
1 11
1
~
,
~
( ),
~
2
nn n
n
n n n
n
n n n n
n
y hf y
y hf y
y
x
y
x x x
y
h
yx x xT
y
??
?
?
? ??
?
??
??
? ? ? ? ?则
? ?hhMT
yM
O
xyx
n
bxa
2
2
21
//
2
2
)(,)( m ax
??
?
?
??
充分光滑，则：令
3,总体方法误差 (1)
? ?? ? ? ? ? ? yxyxxx
n
n
n
nnn
yL,fy,f
L i p s c h i t z
yxf,
???
条件：充分光滑，或满足
），（自然假定一性由微分方程解的存在唯
关系。推出总体误差与步长的
邻步的总体误差关系递推方法：从任意两相
? ?
? ? ? ?
? ? ? ?? ?
11 1 11 1 1
1 1 1
1
11
n 1
~ ~ ~
~~
,
nn n
nn n nn n n
n n n
n n nn
nn
yn
yy
y h f y
yex
y y ye x x T
y y y
y x x x
yy
?? ? ?? ? ?
? ? ?
?
??
? ? ?
? ? ? ? ? ? ? ? ?
? ? ?
第 步的总体截断误差记为 则对 步：
以下估计 其中
总体方法误差 (2)
? ? ? ?? ? ? ?
? ? ? ? ? ?? ?
? ? ? ? ? ?
1
1
~
,,
y,,
1 ( 1 )
n n n nn n n
n
n n n nnn
nn n
y hf y hf
h f y f
L ipsc hitz hL y hL
y y yx x x x
y
yyx x x x
yxe
?
?
? ? ? ? ?
? ? ? ?
? ? ? ? ?由 条件
? ?
? ?
? ? ? ? ? ?
? ? ? ? ? ?
11
00 0
112
21
1 2 1
1 n
N 0
1 1 ( 1 )
1 11
nnn
N N NN N N
N
N N N
hL
y
hL hL hL
hL
ee T
yex
e e eT T T
hL hLT T T T
??
???
?
??
? ? ? ?
? ? ?
? ? ? ? ? ? ? ? ?
? ? ? ? ? ???
L
L
对一切 成立。
对取定,由, 则：
? ?
? ? ? ?
2
11
11
1
n
NnN
O
NhL
hT
hLe T T T?
?
?? ? ? ? ?
?L
由局部截断误差,则
? ? ? ? ? ?11 2
00
1 1
NN
NK
kk
kk Oh L h LhT??
?
??
??????
? ? ? ? ? ?1 2
11
1 N
O O h
hL
hL
h
?
? ? ?
??
? ? ? ? ?
? ? ? ?
? ?
0
00
0
0
1
[]
0
1l im l im 1
l im 1
N
hh
L
h
xx
N N
h
xxLx x h
NhL
hL hL
ehL
??
?
?
? ?
? ?
?
?
?
? ???
?
?
? ?
?
与步长 无关常数
? ?,0,0 ?? e NhhOE u l e r 速率收敛：方法以
4、微分方程数值解的稳定性
累会淹没真值？的近似计算值，误差积是其中
差：稳定性分析，对计算误
1
*
1
*
111
yy
yy
nn
nnn
??
???
???
? ?
'
n n k
n + k n
=,
0
1 2 h
y
yy
h
k
??
?
?
??
?
?
??L
定义：一种数值方法求解模型方程 其中 是复常
数。对给定的步长,若计算误差 在计算
，,时不产生增大的误差，即,称对 与
这种方法是绝对稳定的。
则称绝对稳定区域。
定的，的允许范围内是绝对稳，对 h?
Euler法的绝对稳定区域
*
n
*
n
*
n
nnn
/
yyy
yyy
y
h
h
E u l ery
???
???
??
?
?
计算值：
算法：的
1
1
??? ? nnn h??? 1 误差方程：
是绝对稳定区域当
从而
11
1
1
??
??
?
h
h
nn
?
???
稳定的。对稳定区域称如果整个左半平面是绝
越强。可选大些，方法适应性绝对稳定区域越大，
?A
h
I
m
( λ h )
-2 - 1 0
R
e
( λ h )
二、向后（后退的） Euler 方法
? ? ? ? ? ?,11 nnn yy hxxy x ?? ??用向后差商：
? ?
? ???
?
?
?
?
??
???
yx
yxyy
y
hf
nnnn
00
111
,
则隐式算法：
? ?
? ?
? ? ? ?
? ???
?
?
?
??
??
? ?
?
?
?
? ?,,,khf
,hf
E ul er
y,xyy
yxyy
n
k
nn
k
n
n
n
nn
210
1 1
1
1
0
1
法结合：与为避免解非线性方程，
? ? ? ? ? ?
? ?
? ?
? ?
? ? ? ?
11
111 1 1 1
1
11
,,
0 1
k k k k
nnn n n n
kk
nn
h f f
hL hL
y y y yxx
yy
??
??? ? ? ?
?
??
??
? ? ? ?
?
收敛条件
? ?
? ?
2
1
0 < h L < 1
n
O
Oh
hT ? ?局部截断误差,整体截断误差，
当 时：
? ?
? ?
? ? ? ?? ?
0
1
1
11
,
0,1,2,
1,
nn n n
kk
n n n
hf
h f k
n
y y yx
y y yx
?
?
??
? ??
?
?
? ??
?? ?
L
向后 Euler 方法收敛条件与截断误差
向后 Euler 法的稳定性
yyyy nnn hE u l e ry 11/ ?? ??? ??,法用向后对
??? ? 11 ?? ?? nnn h误差公式：
? ?? ?? ?hhhn
n
???
?
? ??
?
?
??
11 2
1
1 1
1
1 则
? ?? ? ?
?
??
?
? ?????
??
??????
2
2
2
)(211
2
1
2
1
11
hhRhh e
n1
n
( ) 0 1heR ?
?
?
???只要 则
仍受限制。要求
稳定的。但收敛法是因此向后
,10
hhL
AEu l e r
??
?
三、梯形公式
? ? ? ? ? ?
1
1,
n
n
nn
x
y y f x y d t
x
xx
?
? ?? ?由积分途径：
? ? ? ?
? ? ? ?? ?yxyxyy
xyxy
nnnnnn
nnnn
ff
h
yy
111
11
,
2
,
,
???
??
???
??
则得：
：积分用梯形公式，且令
? ?
? ?
? ?
? ?
? ?? ?? ?
?
?
?
?
?
????
??
?
??
?
?
?
?
?
,,,kff
h
hf
nE u l e r
yxyxyy
yxyy
k
nnnnn
k
n
nnnn
210,
2
210
11
1
1
0
1
,
,
，，，，对法结合，形成迭代算法同样与
(,) dy f x ydx ?（）
? ?2 0NE u le r eh ?梯形公式比 法的局部与总体误差均高一阶，,
但每次迭代均多算一次函数值—提高精度的计算代价。
? ?
? ?
? ?
? ?
? ?
? ?
0
1
1
111
0 1 2
,0 1 2
2
,
,
nn n n
kk
nnn n n n
n
hf
h
f f k,,,
y y yx
y y y yxx
?
?
???
?
?
??
?
?
??
? ? ? ? ???
???
L
L
由迭代算法，对,，，
? ? ? ? ? ?
? ?
? ?
? ?
? ? ? ?
111 1 1 1
1
11
,,
2
0 1
22
k k k k
nnn n n n
kk
nn
h
ff
hL hL
y y y yxx
yy
??
??? ? ? ?
?
??
??
? ? ? ?
?
收敛条件
梯形公式的收敛性
梯形公式的稳定性
? ?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
?
?
??
???
?
??
??
??
?
?
????
?
?
?
2
2
2
2
4
1
)(1
4
1
)(1
2
1
1
11
2
1
2
1
2
hhR
hhR
e
e
h
h
h
n
n
nnnn
-A 1 0)( 1 稳定的。梯形公式是，时当 ?? ?
?
?
?
n
n
e hR
? ? ? ?? ?111,2,nnn n n nh ffy y y yxx ???? ? ?
/
1 1
2nn n
hy
nyyy y y?? ? ?? ? ? ?对 用梯形公式,） （
四、改进的尤拉公式
梯形公式虽然提高了精度，但使算法复杂。而在实际计算中
只迭代一次，这样建立的预测 — 校正系统称作改进的尤拉公式。
1
11 1
1
1
1
(,) ;
[ (,) (,) ],
2
[ (,) (,(,) ) ],
2
(,) ;
(,) ;
( ) / 2.
n n nn
n n n n n n
n n n n n n n n
p n n n
c n n p
n p c
y y hf x y
h
y y f x y f x y
h
y y f x y f x h y hf x y
y y hf x y
y y hf x y
y y y
?
?? ?
?
?
?
??
? ? ?
? ? ? ? ?
? ??
?
???
?
??
?
预测
校正
这一计算公式也可表为
或表为下列平均化形式
'
1
1
1
2
( 0 1 );
( 0 ) 1,
0,1
2
( ).
2
( );
2
( );
1
( ).
2
n
n n n
n
n
p n n
n
n
c n p
p
n p c
x
y y x
y
y
h
x
y y h y
y
x
y y h y
y
x
y y h y
y
y y y
?
?
?
?
? ? ? ?
?
?
?
?
?
?
? ? ?
?
? ? ?
?
?
?
? ? ?
?
?
?
? ??
?
例：用尤拉公式和改进的尤拉公式解初值问题
解：取步长,
尤拉公式为,
改进的尤拉公式为：
计算结果略。
尤拉两步公式
'11
11
11
( ) ( )
( ),
2
2 (,),
nn
n
n n n n
nn
y x y x
yx
h
y y hf x y
yy
??
??
??
?
??
??
n- 1 n
上面介绍的改进尤拉公式，其预测公式（尤拉公式）
的精度差，与校正公式（梯形公式）不匹配。用中心差
商 替代 即
称其为尤拉两步公式（首先要确定y,y 两个值)。
用尤拉两步公式与梯形公式相匹配，得到下列
预测—校正系统：
预测
11 1
2 (,) ;
[ (,) (,) ],
2
nn
n n n n n n
hf x y
h
y y f x y f x y
?? ?
? ? ?校正
11
3
'"
1 1
( ),( ),
( ) ( ),
3
n n n n
nn n
y y x y y x
h
y x y y x
??
? ?
??
??
考虑上面给定预测－校正系统的局部截断误差：
对预测公式：假定
其局部截断误差,? ? ? ?
? ? ? ?
2
//
1
34
/// ( 4 )
1
( ) ( ) (,( ) )
2
3 ! 4 !
n n n n n n
n n n
y y h y h f y h yx x x x x x
hh yy
x x x??
?
?
? ? ? ? ?
? ? ? ?
''
111
2 3 4
' '' ''' ( 4 ) '
2 3 4
' '' ''' ( 4 )
2 (,) ( ) 2 ( ) ( ) 2 ( )
( ) ( ) ( ) ( ) ( ) 2 ( )
2 3! 4 !
( ) ( ) ( ) ( ) ( )
2 3! 4 !
n n n n n n nn
n n n n n
n n n n
y y hf x y y x hy x y x h hy x
h h h
y x hy x y x y x y hy x
h h h
y x hy x y x y x y
?
?
???
? ? ? ? ? ? ?
? ? ? ? ? ? ?
? ? ? ? ? ?
L
L
11
1 1
11 11
1 1 11
() 1
,
4()
:
4
( ) ( ) ;
5
1
( ) ( ),
5
nn
n n
nn nn
n n nn
y x y
y x y
y x y y y
y x y y y
??
? ?
?? ??
? ? ??
?
??
?
? ? ? ?
? ? ?
由上两式有
由此导出下列事后估计式
111
3
'"
11
( ),
( ) ( ),
12
nnn
n n n
y y y x
hy x y y x
???
??
?
? ? ?
对校正公式：假定 是准确的，即
其局部截断误差为：
'
1 1 1
2
' ' ' '' '''
23
' '' '''
[ (,) (,) ] ( ) [ ( )
22
( ) ] ( ) [ ( ) ( ) ( ) ( )
22
] ( ) ( ) ( ) ( )
24
n n n n n n n n
n n n n n n
n n n n
hh
y y f x y f x y y x y x
hh
y x h y x y x y x hy x y x
hh
y x hy x y x y x
? ? ?
? ? ? ? ?
? ? ? ? ? ? ?
? ? ? ? ? ?LL
'
11
11
'
1 1 1
''
11
,
,:
2 ;
4
( ) ;
5
(,) ;
( ) ;
2
n n n
n n n n
n n n
n n n n
p y hy
m p P c
m f x m
h
c y y m
??
??
? ? ?
??
??
? ? ?
?
? ? ?
利用估计出的误差作为计算结果的一种补偿 有可
能使精度得到改善此时的计算公式为
预测
改进
计算
校正
1 1 1 1
'
1 1 1
1
( ) ;
5
(,),
n n n n
n n n
y c p c
y f x y
? ? ? ?
? ? ?
? ? ?
?
改进
计算
§ 3.龙格 — 库塔方法
一,Runge-Kutta法的基本思想（ 1）
1
2
' " ( )
11
' '' ' '
( )
( ) ( ) ( ) ( ) ( )
2 ! !
( ) (,),( ) (,) (,) (,),
n
p
p
n n n n n n
xy
p T ay l or y x
hh
y y x y x hy x y x y x
P
y x f x y y x f x y f x y f x y
?
??
? ? ? ? ? ?
? ? ?
L
LL
若用 阶 多项式近似函数 有：
其中 。
但由于公式中各阶偏导数计算复杂，不实用。
' ( 0 )
( 0 ) ( 0 )
'' ( 1 )
( 1 ) ( 1 )
''' ( 2 )
( 2 ) ( 2 )
( ) ( 1 ); 2,3,
jj
jj
y f f
ff
y f f
xy
ff
y f f
xy
ff
y f f j
xy
??
?
??
??
? ? ?
??
??
? ? ?
??
??
? ? ? ?
??
L
一般地有
11
1
1 1 2
1
21
(,)
11
()
22
(,)
(,)
nn
nn
nn
nn
nn
Eule r Eule r
y y hK
Eule r
K f x y
y y h K K
Eule r
K f x y
K f x h y hK
?
?
???
?
??
?
? ? ?
?
?
?
?
?
? ? ? ?
?
如果将 公式与改进 公式写成下列形式：
公式
改进 公式
Runge-Kutta法的基本思想（ 2）
11
(,)
()
(,)
(,)
nn
f x y
y x y
f x y
f x y
??
以上两组公式都使用函数 在某些点上的
值的线性组合来计算 的近似值 。
Euler 公式：每步计算一次 的值，为一阶方法。
改进Euler公式：需计算两次 的值，二阶方法。
Runge-Kutta法的基本思想（ 3）
(,)
(,) ( )
n n n
f x y
x y Tay lo r y x x Tay lo r
于是可考虑用函数 在若干点上的函数值的
线性组合来构造近似公式，构造是要求近似公式在
处的 展开式与解 在 处的 展开
式的前面几项重合，从而使近似公式达到所需要的阶
数。即避免求偏导，又提高了方法的精度，此为RK 方
法的基本思想。
1
1
1
1
1
(,)
(,) ( 2,3,)
p
n n i i
i
nn
i
i n i n ij j
j
y y h c K
K f x y
K f x a h y h b K i p
?
?
?
?
?
???
?
?
??
?
? ? ? ? ?
??
?
? L
二、二阶龙格－库塔方法
1
1
1
1
1
RK
(,)
(,) ( 2,3,)
,
(,) ( )
p
n n i i
i
nn
i
i n i n ij j
j
i ij i
n n n
y y h c K
K f x y
K f x a h y h b K i p
a b c
x y T ay l or y x x T ay l or
?
?
?
?
?
??
?
?
?
??
?
? ? ? ? ?
?
?
?
? L
一般地,方法设近似公式为
其中, 都是参数，确定它们的原则是使近似公式
在 处的 展开式与 在 处的 展开式
的前面项尽可能多地重合。
1 1 1 2 2
1
2 2 21 1
()
2 (,)
(,)
nn
nn
nn
y y h c K c K
K f x y
K f x a h y hb K
? ? ? ??
?
?? ?
? ? ? ?
?
当p 时，近似公式为
1 1 2 2 21
'
1 2 2
'3
21
12
''
2 2 21
(,)
[ (,) (,(,) ) ]
{ (,) [ (,) (,)
(,) (,) ] } ( )
( ) (,)
[ (,) (,) (,
nn
n n n n n n n n
n n n n n x n n
y n n n n
n n n
x n n y n n n
x y T ay l or
y y h c f x y c f x a h y hb f x y
y h c f x y c f x y a hf x y
hb f x y f x y O h
y c c f x y h
c a f x y b f x y f x y
?
? ? ? ? ?
? ? ? ?
??
? ? ?
??
上式在 处的 展开式为
23
) ] ( )
n
h O h?
1
2
' " 3
1
2
' ' 3
()
( ) ( ) ( ) ( ) ( )
2
(,)
[ (,) (,) (,) ] ( )
2
nn
n n n n
n n n
x n n y n n n n
y x x Ta y lo r
h
y x y x h y x y x O h
y f x y h
h
f x y f x y f x y O h
?
?
? ? ? ?
??
? ? ?
在 处的 展开式为
12
22
2 21
1
1 / 2
1 / 2 ),
cc
ca
c b O
???
?
?? ?
? ?
?
3
有无穷多组解，每一组解得一
近似公式，局部截断误差均为
(h 这些方法统称二阶方法。
1 2 2 2 1
1 1 2
1
21
1
,1,
2
( ) / 2
(,)
(,)
nn
nn
nn
c c a b E u le r
y y h K K
K f x y
K f x h y h K
?
? ? ? ?
? ? ??
?
??
? ? ? ?
?
取 此为改进 公式。
近似公式为
1 2 2 2 1
12
1
21
1
0,1,,
2
(,)
( 2,2 )
nn
nn
nn
c c a b
y y h K
K f x y
K f x h y h K
?
? ? ? ?
???
?
??
?
? ? ??
取 此为常用的二阶公式，
称为中点公式。
／／
三、三阶龙格－库塔方法
1 1 2 3
1
21
3 1 2
3
( 4 )
6
(,)
(
,)
22
(,2 )
nn
nn
nn
nn
p
RK
RK
h
y y K K K
K f x y
hh
K f x y K
K f x h y hK hK
?
?
?
? ? ? ?
?
?
??
?
?
? ? ?
?
?
? ? ? ?
?
类似地，对,即三个点，通过更复杂的计算，
可导出三阶 公式。
常用的三阶 公式为：
四、四阶龙格－库塔方法
1 1 2 3 4
1
21
32
43
4
( 2 2 )
6
(,)
(,)
2
2
(,)
22
(,)
nn
nn
nn
nn
nn
p R K
RK
h
y y K K K K
K f x y
hh
K f x y K
hh
K f x y K
K f x h y hK
?
?
?
? ? ? ? ?
?
?
?
?
?
?
? ? ??
?
?
? ? ?
?
?
? ? ??
?
对,即四个点，可导出四阶 公式。
常用的四阶 公式为：
'
11
1 1 2 3 4
1
21
32
43
h= 0,2,x= 0 x= 1
2
( 0 1 ) ;
( 0) 1.
(
6
( 2 2 )
6
(,)
(,)
22
(,)
22
(,)
nn
nn
nn
nn
nn
nn
x
y y x
y
y
h
y y K
h
y y K K K K
K f x y
hh
K f x y K
hh
K f x y K
K f x h y hK
?
?
?
? ? ? ?
?
?
?
?
?
??
?
? ? ? ? ?
?
?
?
?
?
?
?? ? ?
?
?
?
? ? ?
?
?
? ? ??
?
例：设取步长 从 直到 用四阶龙格—库塔
方法求解初值问题
解：由经典的四阶龙格—库塔公式得
2 3 4
1
21
1
32
2
43
3
2 2 ) ;
2;
2;
2
2
2;
2
2
2( )
.
n
n
n
n
n
n
n
n
n
n
n
n
K K K
x
Ky
y
xhh
K y K
h
yK
xhh
K y K
h
yK
xh
K y hK
y hK
?
? ? ?
?
?
?
??
?
?
?
?
? ? ?
??
??
?
? ?
? ? ?
?
??
?
? ?
? ? ?
?
??
?
1 R K
RK
4 6 R K 5
2 R K
R K
RK
R K
两点说明：
）当p=1,2,3,4 时,公式的最高阶数恰好是p,
当p>4 时,公式的最高阶数不是p,如p=5 时仍
为,p= 时 公式的最高阶数为 。
） 方法的导出基于Tay lor 展开，故要求所求问
题的解具有较高的光滑度。
当解充分光滑时，四阶 方法确实优于改进
Euler 法。对一般实际问题，四阶 方法一般可达
到精度要求。
如果解的光滑性差，则用四阶 方法解的效果
不 如改进Eul er法 。
五、变步长的龙格 — 库塔方法
()
1
( ) 5
11
1
5
( 2 )
1
5
( 2 )
11
( 2 )
11
()
11
,
( ),
2
,,
2
( ) 2,
2
() 1
.
( ) 1 6
h
n
h
nn
nn
h
n
h
nn
h
nn
h
nn
h
y
y x y ch
h
xx
h
yc
h
y x y c
y x y
y x y
?
??
?
?
??
??
??
??
??
??
??
??
??
??
??
?
?
?
n
／
／
／
以经典四阶龙格—库塔公式为例。从节点x 出发，以 为
步长求一近似值
将步长折半，即取 为步长从 跨两步到,求一近似
值 每跨一步的截断误差是 因此有
由上两式
( 2 ) ( 2 ) ( )
1 1 1 1
1
( ) [ ].
15
h h h
n n n n
y x y y y
? ? ? ?
? ? ? ?
／／
R-K方法的绝对稳定区域
/
2
211
12
23
1
3
3 2 4 3
2
(,),(,)
22
(,)
(,)
1
,
22
11
,(,)
)
2 2 4
2
2
() hK
( ( ) ( )
n n n n
nn
n n n
n
nn
n
f x y y R K
h
h h h
hh
K f x y K f x y K
hh
K f x y K K f x
hh
K
h
h y h
y
y y y hhK K
y y h hhK K
?
? ? ?
??
?
??
? ? ? ?
? ? ? ? ?
? ? ?
? ? ? ?
? ? ? ? ?
? ? ? ?
? ? ? ?
??
? ? ? ? ?
?
?
?
??
代入 公式：将
234
4 3
11
h )
24
( ( ) ( ) ( )
nn
hhy y h h h
hKK
?? ? ? ?
??
? ? ? ? ? ?
??
??
2 3 411
2 3 4
2 )
6
1 1 1
1
2 6 2 4
( 2
( ) ( ) ( )
nn
n
h
h
y y K K K K
y h h h? ? ? ?
?
? ? ? ? ?
??
? ? ? ? ???
??
?
?
??
?
? ?????
?
)()()( 432
1 24
1
6
1
2
11 hhhh
nn
????? ?则
1
24
1
6
1
2
1
1 )()()(
432
????? hhhh ????
绝对稳定区域：
2
1
- 3 - 2 - 1
0
- 1
- 2
§ 4.线性多步法
1
1 - r
1
RK
,,
nn
n n n
n
yy
y y y
y
?
?
L
－
单步法在计算 时，只用到前一步的信息 。为提高精度，
需重新计算多个点处的函数值，如 方法，计算量较大。
如何通过较多地利用前面的已知信息，如,, 来
构造高精度的算法计算,这就是多步法的基本思想。
1 1 1
1
01
11
(,),,(,)
,(,) ( 1,,,)
00
T a y l o
n n n r n n n r n r
rr
n i n i i n i
ii
i i k k k
y y y f x y f x y
y y h f
f f x y k n n n r
??
??
??
? ? ? ? ? ?
? ? ?
? ? ?
??
? ? ? ?
?
??
LL
L
－－
多步法中最常用的是线性多步法，它的计算公式中只
出现,,, 及 的一次
项，其一般形式为
其中 均为常数,。
若 ＝,显式；,隐式。
构造线性多步公式常用 r 展开和数值积分方法。
线性多步公式的导出
nn
Ta y l o r
x Ta y l o r ) x Ta y l o r
,ii??
n+1
利用 展开导出的基本方法是：将线性多步公式
在 处进行 展开，然后与y ( x 在 处的 展
开式相比较，要求它们前面的项重合，由此确定参数 。
1
01
1 0 1 1 1 1 0 1 1
1 ( )
( )
rr
n i n i i n i
ii
n n n n n n
r y x
y y h f
y y y h f f f
??
? ? ? ? ?
? ? ?
? ? ?
? ? ? ? ?
?
? ? ?
? ? ? ? ?
??
以 为例：设初值问题的解 充分光滑，待定的
两步公式为
( ) ( )
'' ( )
'2
1
( ) ( 1,2,),( )
( ) ( ) ( ) ( )
2!
( ( ) )
kk
n n n
p
pnn
n n n n n
p
n
y y x k y x x T ay lor
yy
y x y y x x x x x x
p
O x x
?
??
? ? ? ? ? ? ? ?
??
L
L
记 则 在 处的 展开为
'
'' '''
' 2 3
1
( 4 ) ( 5 )
4 5 ( 6 )
' ' ''
1 1 1 1
( ),
( ) (,) ( ),
( )
2 ! 3!
( )
4 ! 5 !
(,) ( )
ii
i i i
nn
n n n n
nn
n n n n n n
n y y x
y x f x y i n
yy
y y x h y y h h h
yy
h h O h
f f x y y x y y h
?
? ? ? ?
?
??
? ? ? ? ? ?
? ? ?
? ? ? ?
假设前 步计算结果都是准确的，即
则有
'''
2
( 4 ) ( 5 )
3 4 ( 5 )
'
'''
' ' '' 2
1 1 1 1
( 4 ) (
3
2!
( )
3! 4 !
(,)
(,) ( )
2!
3!
n
nn
n n n n
n
n n n n n n
nn
y
h
yy
h h O h
f f x y y
y
f f x y y x y y h h
yy
h
? ? ? ?
?
? ? ?
??
? ? ? ? ?
??
5)
4 ( 5 )
()
4!
h O h?
' '' 21
1 0 1 1 1 0 1 1 1
''' 3 ( 4 ) 41 1 1 1 1 1
( 5 ) 5 61 1 1
( ) ( ) ( )
2
( ) ( )
6 2 2 2 4 6 6
( ) ( )
1 2 0 2 4 2 4
n n n n
nn
n
y y y h y h
y h y h
y h O h
?
? ? ? ? ? ? ? ?
? ? ? ? ? ?
? ? ?
? ? ?
??
?
? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ?
? ? ? ? ?
将以上各公式代入并整理，得
1
'' ( 5 )
' 2 5 6
1
()
( ) ( )
2 ! 5 !
p + 1
nn
nn
n n n
p y x x
T a y lo r
yy
y x y y h h h O h
?
?
? ? ? ? ? ?L
为使上式有 阶精度，只须使其与 在 处的
展开式
的前 项重合。
01
0 1 0 1
1 1 1
1 1 1
1 1 1
1
1
11
22
1 1 1 1
6 2 2 6
1 1 1 1
24 6 6 24
aa
a
a
a
a
? ? ?
??
??
??
?
?
?
?
???
?
? ? ? ? ?
?
?
? ? ?
?
?
? ?
?
? ? ? ?
?
?
?
? ? ?
??
5,5
,1ii P?? ?
个参数只须 个条件。由推导知，如果选取参数
，使其满足前 个方程（p = 1,2,3,4 ),则
近似公式为p 阶公式。
1 1 0 1
11
1
1,0,,0
2
( )
2
n n n n
h
y y f f
? ? ? ? ?
?
??
? ? ? ? ?
? ? ?
0
如 满足方程组前三个
方程，故公式
此为二阶公式。
0 1 1 1 0
1 1 1 1
5 ( 5 ) 6
1
14
0,1,,
33
( 4 )
3
1
( )
90
n n n n n
nn
h
y y f f f
R h y O h
? ? ? ? ?
?
? ? ? ?
?
? ? ? ? ?
? ? ? ?
? ? ?
又如：解上面方程组得
相应的线性二步四阶公式（S i m p s o n 公式) 为
其截断误差为
由此可知，线性二步公式至多是四阶公式。
1
0
1
01
1
1
1
2 3
()
1
( ) ( ) 1 ( 1,2,)
[ 1 ( ) ( 1 ) (
n n n
r
i
i
rr
kk
ii
ii
r
p
ni
i
r
x T ay l or y x x T ay l or
i k i k p
p
R i p
?
??
?
?
?
?
? ? ?
?
?
?
?
?
?
?
?
?
?
? ? ? ? ?
?
?
? ? ? ? ?
?
??
?
L
一般地，线性多步公式中有 个待定参数，如令其
右端在 处的 展开式与 在 处的 展开式
的前p+1项系 数对应相等，可得方程组
其解所对应的公式具有 阶精度，局部截断误差为
1
( 1 )
1
2
)]
( 1 ) !
( )
pr
pp
in
i
p
h
iy
p
Oh
?
?
?
??
?
?
?
?
?
显然，线性多步公式至多可达到2r +2 阶精度。
二、常用的线性多步公式
1 2 3 1
0
1
01
0 0 1 2 3
1 1 2 3
( A da m s)
r =3 0,
1
( ) ( ) 1 ( 1,2,3 4)
55 59 37 9
=1,,,
24 24 24 24
( 55 59 37 9
24
r
i
i
rr
kk
ii
ii
n n n n n n
i k i k
h
y y f f f f
? ? ? ?
?
??
? ? ? ? ?
?
?
?
? ? ?
? ? ? ?
? ? ? ?
?
?
?
?
?
?
? ? ? ? ?
?
?
? ? ? ? ? ?
? ? ? ? ?
?
??
（一）阿达姆斯 公式
取,并令 由方程组
，
可解得,
相应的线性多步公式为
1
)
=0 A da m s?
－
因,此式称为 显式公式，是四阶公式.
533
5 4 ( 5 ) 6
1
11
5 ( 5 ) 6
[ 1 ( ) 5 ( ) ] ( )
5!
251
()
720
n i i n
ii
n
h
R i i y O h
h y O h
??
?
? ? ?
? ? ? ? ? ?
??
??
局部截断误差为
1 2 3 3
0 1 0 1 2
1 1 1 2
5 ( 5 ) 6
1
0,
9 19 5 1
= 1,,,
24 24 24 24
( 9 19 5 )
24
A da m s
19
()
720
n n n n n n
nn
h
y y f f f f
R h y O h
? ? ? ?
? ? ? ? ?
?
? ? ? ?
?
? ? ? ?
? ? ? ? ?
? ? ? ? ?
? ? ?
如果令 由方程组可解得
，
相应的线性多步公式为
称其为四阶 隐式公式，其局部截断误差为
利用数值积分方法求线性多步公式
11
1
1 + 1 1
( ) ( ) (,( ) ) ( )
,,,,,,( )
( ) ( ),1
nn
nn
xx
nn
xx
n n n r n n n r
y x y x f x y x dx F x dx
x x x x x x F x
x F x r?
??
?
? ? ? ?
? ? ?
?
??
LL
r
基本思想是首先将初值问题化成等价的积分形式
用过节点 或 的 的r
次插值多项式 代替 求积分即得 阶的
线性多步公式。
1 2 3
3
3
0
1 2 3
3
0
3,,,,( )
( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( 0,1,2,3 )
( ) ( )
n n n n
i n i
i
n n n n
i
n i n i n j
j
ji
r x x x x F x
L x l x F x
x x x x x x x x
l x i
x x x x
? ? ?
?
?
? ? ?
? ? ?
?
?
?
?
? ? ? ?
??
??
?
?
例如 时，过节点 的三次
插值多项式为
其中
11
1
1
1
3
13
0
1 2 3
3
23
1
3
13
2
3
3
( ) ( ) ( ) [ ( ) ] ( )
( ) ( ) ( )
()
6
( ) ( ) ( )
()
2
( ) ( ) ( )
()
2
( ) (
()
nn
nn
n
n
n
n
n
n
xx
n n i n i
xx
i
x
n n n
n
x
x
n n n
n
x
x
n n n
n
x
n
n
y x y x L x d x l x d x F x
x x x x x x
F x d x
h
x x x x x x
F x d x
h
x x x x x x
F x d x
h
x x x
Fx
??
?
?
?
??
?
? ? ?
??
?
??
?
?
? ? ?
? ? ?
?
? ? ?
?
?
? ? ?
?
?
?
???
?
?
?
1
12
3
1 2 3
) ( )
6
[ 5 5 ( ) 5 9 ( ) 3 7 ( ) 9 ( ) ]
24
n
n
x
nn
x
n n n n
x x x
dx
h
h
F x F x F x F x
?
??
? ? ?
??
?
? ? ? ?
?
11
1 1 2 3
3
,( ),( ),(,)
( ) (,( )) (,1,2,3 ),
(55 5 9 3 7 9 )
24
[,]
n n n n k k k
k k k
n n n n n n
nn
y y y x y x f x y
F x f x y x k n n n n
h
y y f f f f
Ad a ms
x x Ad a ms
??
? ? ? ?
?
? ? ? ? ?
? ? ? ? ?
对上式用 代替 用 代替
则得
这就是四阶 显式公式。由于积分区间在插值
区间 外面，又称为四阶 外插公式。
11
1
( 4 ) ( 5 )33
1
00
31
( 5 ) 3
5 ( 5 )
1
0
( ) ( )
( ) ( )
4 ! 4 !
,),
( ) 25 1
( ) ( )
4 ! 72 0
nn
nn
n
n
xx
xx
n n j n j
xx
jj
nn
x
n n j
x
j
Fy
R x x dx x x dx
xx
y
R x x dx h y
??
?
?
?
??
?
? ? ?
??
??
??
?
? ? ? ?
?
? ? ?
????
??
由插值余项公式可得其局部截断误差为
由积分中值定理，存在 （ 使得
1 1 2
2
3
1
1 1 2
3
1
1
,,,( )
( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( 1,0,1,2)
( ) ( )
()
n n n n
i n i
i
n n n n
i
n i n i n j
j
ji
nn
x x x x F x
L x l x F x
x x x x x x x x
l x i
x x x x
F x A da m s
yy
??
?
??
? ? ?
? ? ?
??
?
?
?
? ? ? ?
? ? ?
??
?
?
?
＋
同样，如果过节点 的 三次插值多
项式为
其中
代替 求积分，即得四阶 隐式公式
1 1 2
5 ( 5 )
1 2 1
21
( 9 19 5 )
24
19
()
720
[,]
n n n n
n n n n
nn
h
f f f f
R h y x x
x x A da m s
A da m s
??
? ? ?
? ? ?
??
? ? ? ?
? ? ? ?其局部截断误差为
由于积分区间在插值区间 内，故 隐式
公式又称为 内插公式
1 2 1
3 0 1 2 3
1 3 1 2
5 ( 5 ) 6
1
( )
0,
8 4 8
1,,,,0,
3 3 3
4
( 2 2 )
3
14
( )
45
n n n n n
nn
M i l i ne
y y h f f f
M i l i ne
R h y O h
M i l i ne
? ? ? ?
? ? ? ? ?
?
? ? ? ?
?
? ? ? ?
? ? ? ? ? ?
? ? ? ?
??
0
（二）米尔尼 公式
取r =3,并令 由方程组可
解出
相应的线性多步式
称为 公式，其局部截断误差为
公式是四阶四 步显式公式。
12
1 2 1 1
5 ( 5 ) 6
1
( m in )
0,m i n
13
( 9 ) ( 2 )
88
1
( )
40
n n n n n n
nn
H am g
H am g
y y y h f f f
R h y O h
??
? ? ? ?
?
??
? ? ? ? ?
? ? ?
（三）哈明 公式
取r=2,并令 可得到 公式
其局部截断误差为
Hamming 公式是四阶三步隐式公式。
1
1
,
)
,
)
n
n
y
?
?
n+1
n+1
n+1
隐式法与显式法的比较
一般地，同阶的隐式法比显式法精确，而且数值
稳定性也好。但在隐式公式中，通常很难解出y 需
要用迭代法求解，这样又增加了计算量。
在实际计算中，很少单独用显式公式或隐式公式，
而是将它们联合使用：先用显式公式求出y(x 的预
测值，记作y 再用隐式公式对预测值进行校正，求
出y (x 的近似值 。
§ 5.预测 — 校正系统
用显式公式计算预测值，然后用隐式公式进行校正，
得到近似值 yn+1这样一组计算公式称为预测 — 校正系统。
一般采用同阶的隐式公式与显式公式。常用的预
测 — 校正系统有两种,
1 2 31
1 1 1 21
( 5 5 5 9 3 7 9 )
24
[ 9 (,) 1 9 5 ]
24
n n n n nn
n n n n n nn
A d a m s
h
y y f f f f
h
y y f x y f f f
? ? ??
? ? ? ??
?
? ? ? ? ?
??
?
? ? ? ? ? ?
??
预测—校正公式
预测
校正
3 1 21
1 2 1 11
l i m i n
4
( 2 2 )
3
13
( 9 ) [ (,) 2 ) ]
88
n n n nn
n n n n n nn
M i ne H am g
y y h f f f
y y y h f x y f f
? ? ??
? ? ? ??
?
?
? ? ? ?
?
?
?
? ? ? ? ? ?
??
预测—校正公式
( 1)
RK
3
说明：
以上两种预测—校正公式均为四阶公式，其起步值
通常用四阶 公式计算。
(2) 有时为提高精度，校正公式可迭代进行多次，但迭
代次数一般不超过 次。
用局部截断误差进一步修正预测－校正公式
5 ( 5 ) 6
1 1
5 ( 5 ) 6
11
5 ( 5 ) 6
1 1
5 ( 5 )
1 1
11 11
1
251
( ) ( )
720
19
( ) ( )
720
270
()
720
720
()
270
251
( ) ( )
270
()
nn n
n n n
nn n
nn n
nn nn
n
Adams
y x y h y O h
y x y h y O h
y y h y O h
h y y y
y x y y y
yx
? ?
??
? ?
? ?
?? ??
?
? ? ?
? ? ? ?
? ? ? ?
? ? ?
? ? ? ?
?
由 公式的局部截断误差公式
两式相减
11 1
19
()
270
nn n
y y y
?? ?
? ? ?
1 1 2 3
11
1 1 1 1 2
1 1 1 1
( 55 59 37 9 )
24
251
()
270
[ 9 (,) 19 5 ]
24
19
()
270
n n n n n n
n n n n
n n n n n n n
n n n n
A da m s
h
p y f f f f
m p c p
h
c y f x m f f f
y c c p
? ? ? ?
??
? ? ? ? ?
? ? ? ?
?
? ? ? ? ?
?
?
?
? ? ?
??
?
?
? ? ? ? ?
?
?
? ? ? ?
??
由上面就得到多环节的 预测－校正公式
预测
改进
校正
改进
1 3 1 2
11
1 2 1 1 1
1 1 1 1
l n m i n
4
( 2 2 )
3
112
()
121
13
( 9 ) [ (,) 2 ]
88
9
()
121
n n n n n
n n n n
n n n n n n n
n n n n
M i e H am g
p y h f f f
m p c p
c y y h f x m f f
y c c p
? ? ? ?
??
? ? ? ? ?
? ? ? ?
?
?
? ? ? ?
?
?
?
? ? ?
??
?
?
? ? ? ? ?
?
?
? ? ? ?
?
?
完全类似，可以导出多环节的
预测－校正公式
预测
改进
校正
改进
00
1 1 1
1
1 1 2 1 1
1,,(,),,( 2),,1
( 3 ) (,)
(,)
22
n n n
n n n
ba
a b f x y N y h x a n
N
f f x y
hK
K hf K hf x y
K
? ? ?
? ? ?
?
? ? ?
?
? ? ? ?
算法
（ ）输入 置
计算
2
3 1 1 4 1 1 3
1 1 2 3 4
00
(,) (,)
22
1
( 2 2 )
6
( 4) (,)
( 5 ) 3,1,3
1,0,0,6
n n n n
n n n
nn
hK
hf x y K hf x h y K
y y K K K K x a nh
xy
n n n
n n p c
? ? ? ?
?
? ? ? ? ? ?
? ? ? ? ? ? ?
? ? ?
? ? ? ?
输出
若 置 返回 ；
否则，置 转 。
3 3 3 3
0 3 2 1 0 0
3 1 3 2
1 1 1
3 3 0
( 6 )
(,)
4 1 1 2
( 2 2 ) ( )
3 1 2 1
13
( 9 ) [ (,) 2 ]
88
9
( ),)
121
( 7 ) 1,,,
( 0,1,2 ),,,,
j j j j j j
f f x y x x h
p y f f f m p c p
c y y h f x m f f
y c c p x y
n N n n x x y y f f
j x x y y p p
? ? ?
? ? ?
? ? ? ? ? ? ?
? ? ? ? ?
? ? ?
? ? ? ? ? ?
? ? ? ?
计算
输出(
若,置
0
,6cc ? 转；
否则停机。
§ 5,常微分方程组与高阶方程的数值解
'
12
0
0
(,,,,) ( 1,2,,)
( ) ( 1,2,,)
i i N
ii
y f x y y y i N
y x y i N
? ???
?
????
LL
L
一、一阶方程组
一阶方程组的初值问题
12
'
0 0 0 0
12 0
0
12
y f
y (,,) ;
y f (,y ) ;
y (,,) ;
y ( ) y ;
f (,,) ;
T
N
T
N
T
N
y y y
x
y y y
x
f f f
?
? ?
?? ?
??
?
L
L
L
若对 和 采用向量的记号
'
0
0
n+ 1 n 1 2 3 4
1 n 2 n 1
3 n 2 4 n 3
y f (,y ) ;
y ( ) y ;
y y ( k 2 k 2 k k )
6
k f (,y ) ; k f (,y k ) ;
22
k f (,y k ) ; k f (,y k ) ;
22
nn
nn
x
x
h
hh
xx
hh
x x h h
? ?
?
?
?
? ? ? ? ?
? ? ? ?
? ? ? ? ? ?
求解这一初值问题的四阶龙格—库塔公式为
其中
,1 1 2 3 4
1 1 2
2 1 1 1 2 2 1 1
3 1 1 2
( 2 2 )
6
( 1,2,,),
(,,,,) ;
(,,,,) ;
2 2 2 2
(,,
22
i n in i i i i
i i n n n Nn
i i n n n Nn N
i i n n
h
y y K K K K
iN
K f x y y y
h h h h
K f x y K y K y K
hh
K f x y K
?
? ? ? ? ?
?
?
? ? ? ? ?
? ? ?
L
L
L
或表达为
其中
2 2 2 2
4 1 1 3 2 2 3 3
,,) ;
22
(,,,,) ;
n Nn N
i i n n n Nn N
hh
y K y K
K f x h y hK y hK y hK
??
? ? ? ? ?
L
L
? ? ? ?
? ? ? ?
'
0 0
'
00
,,
,,
dy
y f x y z y
dx
dz
z g x y z z
dx
yx
x z
?
? ? ?
?
?
?
? ? ? ?
??
以两个方程组为例
? ? ? ?
? ? ? ?
?,2,1,0
,,,
,,,
,1
001
0
0
1
?
?
?
?
?
?
???
???
?
?
n
zhg
yhf
E ul e r
zxzyxzz
yxzyxyy
n
n
nnn
n
n
n
nn
法方程组的
2) 方程组的 R-k 法
? ?
? ?
1 2 3 41
1 1 2 3 4
2 2
6
2 2
6
nn
nn
h
h
yy K K K K
z z M M M M
?
?
? ? ? ? ?
? ? ? ? ?
? ? ? ?11,,,,,Mnnnnnnfg yyxxK z z??
11
2
11
2
,,
2 2 2
,,
2 2 2
nn n
nn n
h
f
h
g
h K h M
yxKz
h K h M
yxMz
??
? ? ? ???
??
??
? ? ? ???
??
22
3
22
3
,,
2 2 2
,,
2 2 2
nn n
nn n
h
f
h
g
h K h M
yxKz
h K h M
yxMz
??
? ? ? ???
??
??
? ? ? ???
??
? ?
? ?
4 33
4 33
,,
,,
0,1,2,
nn n
nn n
fh
gh
n
yx hK hMKz
yx hK hMMz
? ? ? ?
? ? ? ?
? L
二、化高阶方程为一阶方程组
( ) ' ( 1 )
' ' ( 1 ) ( 1 )
0 0 0 0 0 0
' '' ( 1 )
1 2 3
'
1 2 1 0
(,,,,);
( ),( ),,( ),
,,,,,
()
mm
mm
m
m
m
y f x y y y
y x y y x y y x y
y y y y y y y y
m
y y y x y
?
??
?
? ?
?
? ? ?
?
? ? ? ?
??
L
L
L
高阶微分方程的初值问题可通过变量代换化为一阶
微分方程组的初值问题。设有 阶常微分方程初值问题
引入新变量
则可将 阶方程化为如下一阶方程组：
0
''
2 3 2 0 0
' ( 2 )
1 1 0 0
' ( 1 )
1 2 0 0
()
()
(,,,,) ( )
m
m m m
m
m m m
y y y x y
y y y x y
y f x y y y y x y
?
??
?
?
?
??
?
?
?
?
??
?
? ??
?
L
L
"'
'
''
0 0 0 0
'
00
''
00
1 1 2 3 4
1 1 2 3 4
(,,) ;
( ),( ),
,( ) ;
(,,),( ),
( 2 2 ) ;
6
( 2 2 )
6
,1,2,3,4)
nn
nn
ii
y f x y y
zy
y x y y x y
y z y x y
z f x y z z x y
h
y y K K K K
h
z z M M M M
K M i
?
?
? ?
?
?
??
?
? ??
?
?
??
?
?
? ? ? ? ?
?
?
?
?
? ? ? ? ?
?
?
?
以二阶方程的初值问题为例：
引入新变量
对其用四阶龙格－库塔公式
（ 计算公式略.
? ?
? ?
" 2 '
'
'
'
'2
0 0 0
0
2
1
1
5 sin 2 2 0 1
( 0 ) 2,( 0 ) 3
( 0 ) 2
,
5 sin 2 2 ( 0 ) 3
,,,
0,2,3
,
2
,
3
,2
,
M
x
x
nn
n
nn
n
n
n
RK
y e x y y x
yy
y z y
zy
z e x y z z
x y z R K
zf
h
f
g
yxKz
z
xK
yx
? ? ? ? ? ?
?
? ? ? ?
?
? ? ? ?
??
?
? ? ? ? ?
?
? ? ? ?
?
?
?
?
?
??
?
?
?
例：用四阶 方法求解初值问题 ( 取h = 0, 1 )
解：令
因 由 公式得
1
11
1
2
0
1
M 3,1
2
1,6 2 3 8 2 2 4
,
22
,,
2 2 2
n
nn
n
h
z
h
g
hK hM
y z
hK hM
yxMz
??
??
??
??
??
? ? ? ?
??
? ? ? ?
??
??
LL 继续计算（略）。
§ 1,引言
0 0
:
(,)
()
dy
f x y a x b
dx
y yx
?
? ? ??
?
? ?
?
一阶常微分方程的初值问题
''
,x y - 2 y = 4 x y =
f ( x,y ) = y ( 1 ) = - 3
(,)
()
2y
4
x
2y
4
x
d y 2 y
f x y 4
d x x
y 1 3
??
?
?
? ? ? ?
?
? ??
?
例 方程
令,且给出初值
就得到一阶常微分方程的初值问题：
(,)
L
(,) (,)
f x y
f x y f x y L y y? ? ?
只要函数 适当光滑连续，且关于y 满足李普
希兹( L i p s c h i t z ) 条件，即存在常数,使得
由常微分方程理论知，初值问题的解必存在且唯一。
微分方程的数值解：设方程问题的解 y(x)的存在区间
是 [a,b]，令 a= x0< x1<… < xn =b,其中 hk=xk+1-xk,如是等距
节点 h=(b-a)/n,h称为步长。
y(x)的解析表达式不容易得到或根本无法得到，我们用数
值方法求得 y(x)在每个节点 xk上 y(xk)的近似值，用 yk表示，
即 yk≈ y(xk)，这样 y0,y1,...,yn称为微分方程的数值解。
微分方程离散化常用方法
? ?
? ? ? ?
1
1
,
A
(,( ) )
n n
nn
nn
nn
yy
f x y x
yx
dy xx
xxdx
?
?
?
?
??
?
用差商代替微商
? ? ? ?
? ?
? ? ?,2,1,0,
,
,,
1
1
1
1
1
???
?
?
????
?
?
?
?
?
nhf
f
h
yyh
yxyy
yx
yy
xyxyxx
n
n
nn
n
n
nn
n
n
n
n
nn
代替，则：用
11
B,
:
(,) ( 0,
1,)
nn
nn
xx
xx
dy
dx f x y dx n
dx
??
???? L
数值积分
用数值积分方法离散化
1
11
1
,( ),( ),
(,) (,)
( 0,1,(,))
n
n
n n n n
x
nn
x
n n n n
y y y x y x
f x y dx h f x y
y y hf x ny
?
??
?
?
???
?
L
用 代替 对右端积分采用
取左端点的矩形公式
则有
2
/ //
2
//
C ( )
( ) ( ) ( ) ( )
2
( ) (,( ) ) ( )
2
n
n n n n
n n n n
y x Tay lor
y h y h
y h f y
x
h
yyx x x x
h
yx x x x
?
? ? ? ? ?
? ? ? ?
L
L
在 附近 的 展开：
1
1
( ) ( )
(,) 0,1,2,
nnn
nn n n
h y y
h f n
Tay lor
y xx
y y yx
?
?
?
? ? ? L
取 的线性部分，且 得 的近似值：
展开法不仅可得到求数值解的公式，且容易估计
截断误差。
0 0
(,)
()
dy
f x y a x b
dx
y yx
???
? ? ????
?
? ?
???
§ 2 尤拉（ Eular)方法
? ?
1
0 0
0 0
,0 1,(,)
()()
nn n n
E u le r
dy
h f n,f x y a x b
dx
yy
y y yx
yxyx
?
? ?
? ? ?? ? ?
??
???
??? ?
??
L
一,公式
计算公式：
? ? ? ?
? ?
? ?
? ? ? ?
0 0
0
0
0
0
,
00
00
1
,
,
f,,
y y x
f
yx
yx
dy
yx
dx
yyxx
?
?
，几何意义。
由 出发取曲线 的切线（存在！），则 斜率
由于 及 已知，必有切线方程。
? ?
? ?
),(f)x(xy yxxy
dx
dy
xy
y,x
000000
00
??????
：由点斜式写出切线方程
）（
：，可由切线算出，则为等步长
,yxyy
yxx
hf
hh
0001
101
??
??
? 2 1 0
1
1
，，，）（
：，点值在）（逐步计算出
,???
?
?
?
nhf
xyy
yxyy
x
nnnn
n
注意：这是“折线法”而非“切线法”
除第一个点是曲线切线外，其他点不是！
y
χ 0 χ 1 χ 2 χ 3 χ
2,Euler方法的误差估计
? ?
1 1
1
1
1)
~
~
( )
n n
n
n
n
n
n
y
y
E ule r
xT
y
y x
y
y
? ?
?
?
??
?
局部截断误差。
在一步中产生的误差而非累积误差：
其中 （精确解！）时
由 法求出的值，即 无误差！
是当
? ?
? ? xxyh
xxxxx
xx
nn
nnnnn
nn
yhfyhyy
Tay l o ry
1
//
2
1
1
2
))(,()()(
)(
?
?
?
???
?????
??
展开：点将 在
? ?
? ?? ?
? ? ? ?
1
1
2
//
1 11
1
~
,
~
( ),
~
2
nn n
n
n n n
n
n n n n
n
y hf y
y hf y
y
x
y
x x x
y
h
yx x xT
y
??
?
?
? ??
?
??
??
? ? ? ? ?则
? ?hhMT
yM
O
xyx
n
bxa
2
2
21
//
2
2
)(,)( m ax
??
?
?
??
充分光滑，则：令
3,总体方法误差 (1)
? ?? ? ? ? ? ? yxyxxx
n
n
n
nnn
yL,fy,f
L i p s c h i t z
yxf,
???
条件：充分光滑，或满足
），（自然假定一性由微分方程解的存在唯
关系。推出总体误差与步长的
邻步的总体误差关系递推方法：从任意两相
? ?
? ? ? ?
? ? ? ?? ?
11 1 11 1 1
1 1 1
1
11
n 1
~ ~ ~
~~
,
nn n
nn n nn n n
n n n
n n nn
nn
yn
yy
y h f y
yex
y y ye x x T
y y y
y x x x
yy
?? ? ?? ? ?
? ? ?
?
??
? ? ?
? ? ? ? ? ? ? ? ?
? ? ?
第 步的总体截断误差记为 则对 步：
以下估计 其中
总体方法误差 (2)
? ? ? ?? ? ? ?
? ? ? ? ? ?? ?
? ? ? ? ? ?
1
1
~
,,
y,,
1 ( 1 )
n n n nn n n
n
n n n nnn
nn n
y hf y hf
h f y f
L ipsc hitz hL y hL
y y yx x x x
y
yyx x x x
yxe
?
?
? ? ? ? ?
? ? ? ?
? ? ? ? ?由 条件
? ?
? ?
? ? ? ? ? ?
? ? ? ? ? ?
11
00 0
112
21
1 2 1
1 n
N 0
1 1 ( 1 )
1 11
nnn
N N NN N N
N
N N N
hL
y
hL hL hL
hL
ee T
yex
e e eT T T
hL hLT T T T
??
???
?
??
? ? ? ?
? ? ?
? ? ? ? ? ? ? ? ?
? ? ? ? ? ???
L
L
对一切 成立。
对取定,由, 则：
? ?
? ? ? ?
2
11
11
1
n
NnN
O
NhL
hT
hLe T T T?
?
?? ? ? ? ?
?L
由局部截断误差,则
? ? ? ? ? ?11 2
00
1 1
NN
NK
kk
kk Oh L h LhT??
?
??
??????
? ? ? ? ? ?1 2
11
1 N
O O h
hL
hL
h
?
? ? ?
??
? ? ? ? ?
? ? ? ?
? ?
0
00
0
0
1
[]
0
1l im l im 1
l im 1
N
hh
L
h
xx
N N
h
xxLx x h
NhL
hL hL
ehL
??
?
?
? ?
? ?
?
?
?
? ???
?
?
? ?
?
与步长 无关常数
? ?,0,0 ?? e NhhOE u l e r 速率收敛：方法以
4、微分方程数值解的稳定性
累会淹没真值？的近似计算值，误差积是其中
差：稳定性分析，对计算误
1
*
1
*
111
yy
yy
nn
nnn
??
???
???
? ?
'
n n k
n + k n
=,
0
1 2 h
y
yy
h
k
??
?
?
??
?
?
??L
定义：一种数值方法求解模型方程 其中 是复常
数。对给定的步长,若计算误差 在计算
，,时不产生增大的误差，即,称对 与
这种方法是绝对稳定的。
则称绝对稳定区域。
定的，的允许范围内是绝对稳，对 h?
Euler法的绝对稳定区域
*
n
*
n
*
n
nnn
/
yyy
yyy
y
h
h
E u l ery
???
???
??
?
?
计算值：
算法：的
1
1
??? ? nnn h??? 1 误差方程：
是绝对稳定区域当
从而
11
1
1
??
??
?
h
h
nn
?
???
稳定的。对稳定区域称如果整个左半平面是绝
越强。可选大些，方法适应性绝对稳定区域越大，
?A
h
I
m
( λ h )
-2 - 1 0
R
e
( λ h )
二、向后（后退的） Euler 方法
? ? ? ? ? ?,11 nnn yy hxxy x ?? ??用向后差商：
? ?
? ???
?
?
?
?
??
???
yx
yxyy
y
hf
nnnn
00
111
,
则隐式算法：
? ?
? ?
? ? ? ?
? ???
?
?
?
??
??
? ?
?
?
?
? ?,,,khf
,hf
E ul er
y,xyy
yxyy
n
k
nn
k
n
n
n
nn
210
1 1
1
1
0
1
法结合：与为避免解非线性方程，
? ? ? ? ? ?
? ?
? ?
? ?
? ? ? ?
11
111 1 1 1
1
11
,,
0 1
k k k k
nnn n n n
kk
nn
h f f
hL hL
y y y yxx
yy
??
??? ? ? ?
?
??
??
? ? ? ?
?
收敛条件
? ?
? ?
2
1
0 < h L < 1
n
O
Oh
hT ? ?局部截断误差,整体截断误差，
当 时：
? ?
? ?
? ? ? ?? ?
0
1
1
11
,
0,1,2,
1,
nn n n
kk
n n n
hf
h f k
n
y y yx
y y yx
?
?
??
? ??
?
?
? ??
?? ?
L
向后 Euler 方法收敛条件与截断误差
向后 Euler 法的稳定性
yyyy nnn hE u l e ry 11/ ?? ??? ??,法用向后对
??? ? 11 ?? ?? nnn h误差公式：
? ?? ?? ?hhhn
n
???
?
? ??
?
?
??
11 2
1
1 1
1
1 则
? ?? ? ?
?
??
?
? ?????
??
??????
2
2
2
)(211
2
1
2
1
11
hhRhh e
n1
n
( ) 0 1heR ?
?
?
???只要 则
仍受限制。要求
稳定的。但收敛法是因此向后
,10
hhL
AEu l e r
??
?
三、梯形公式
? ? ? ? ? ?
1
1,
n
n
nn
x
y y f x y d t
x
xx
?
? ?? ?由积分途径：
? ? ? ?
? ? ? ?? ?yxyxyy
xyxy
nnnnnn
nnnn
ff
h
yy
111
11
,
2
,
,
???
??
???
??
则得：
：积分用梯形公式，且令
? ?
? ?
? ?
? ?
? ?? ?? ?
?
?
?
?
?
????
??
?
??
?
?
?
?
?
,,,kff
h
hf
nE u l e r
yxyxyy
yxyy
k
nnnnn
k
n
nnnn
210,
2
210
11
1
1
0
1
,
,
，，，，对法结合，形成迭代算法同样与
(,) dy f x ydx ?（）
? ?2 0NE u le r eh ?梯形公式比 法的局部与总体误差均高一阶，,
但每次迭代均多算一次函数值—提高精度的计算代价。
? ?
? ?
? ?
? ?
? ?
? ?
0
1
1
111
0 1 2
,0 1 2
2
,
,
nn n n
kk
nnn n n n
n
hf
h
f f k,,,
y y yx
y y y yxx
?
?
???
?
?
??
?
?
??
? ? ? ? ???
???
L
L
由迭代算法，对,，，
? ? ? ? ? ?
? ?
? ?
? ?
? ? ? ?
111 1 1 1
1
11
,,
2
0 1
22
k k k k
nnn n n n
kk
nn
h
ff
hL hL
y y y yxx
yy
??
??? ? ? ?
?
??
??
? ? ? ?
?
收敛条件
梯形公式的收敛性
梯形公式的稳定性
? ?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
?
?
??
???
?
??
??
??
?
?
????
?
?
?
2
2
2
2
4
1
)(1
4
1
)(1
2
1
1
11
2
1
2
1
2
hhR
hhR
e
e
h
h
h
n
n
nnnn
-A 1 0)( 1 稳定的。梯形公式是，时当 ?? ?
?
?
?
n
n
e hR
? ? ? ?? ?111,2,nnn n n nh ffy y y yxx ???? ? ?
/
1 1
2nn n
hy
nyyy y y?? ? ?? ? ? ?对 用梯形公式,） （
四、改进的尤拉公式
梯形公式虽然提高了精度，但使算法复杂。而在实际计算中
只迭代一次，这样建立的预测 — 校正系统称作改进的尤拉公式。
1
11 1
1
1
1
(,) ;
[ (,) (,) ],
2
[ (,) (,(,) ) ],
2
(,) ;
(,) ;
( ) / 2.
n n nn
n n n n n n
n n n n n n n n
p n n n
c n n p
n p c
y y hf x y
h
y y f x y f x y
h
y y f x y f x h y hf x y
y y hf x y
y y hf x y
y y y
?
?? ?
?
?
?
??
? ? ?
? ? ? ? ?
? ??
?
???
?
??
?
预测
校正
这一计算公式也可表为
或表为下列平均化形式
'
1
1
1
2
( 0 1 );
( 0 ) 1,
0,1
2
( ).
2
( );
2
( );
1
( ).
2
n
n n n
n
n
p n n
n
n
c n p
p
n p c
x
y y x
y
y
h
x
y y h y
y
x
y y h y
y
x
y y h y
y
y y y
?
?
?
?
? ? ? ?
?
?
?
?
?
?
? ? ?
?
? ? ?
?
?
?
? ? ?
?
?
?
? ??
?
例：用尤拉公式和改进的尤拉公式解初值问题
解：取步长,
尤拉公式为,
改进的尤拉公式为：
计算结果略。
尤拉两步公式
'11
11
11
( ) ( )
( ),
2
2 (,),
nn
n
n n n n
nn
y x y x
yx
h
y y hf x y
yy
??
??
??
?
??
??
n- 1 n
上面介绍的改进尤拉公式，其预测公式（尤拉公式）
的精度差，与校正公式（梯形公式）不匹配。用中心差
商 替代 即
称其为尤拉两步公式（首先要确定y,y 两个值)。
用尤拉两步公式与梯形公式相匹配，得到下列
预测—校正系统：
预测
11 1
2 (,) ;
[ (,) (,) ],
2
nn
n n n n n n
hf x y
h
y y f x y f x y
?? ?
? ? ?校正
11
3
'"
1 1
( ),( ),
( ) ( ),
3
n n n n
nn n
y y x y y x
h
y x y y x
??
? ?
??
??
考虑上面给定预测－校正系统的局部截断误差：
对预测公式：假定
其局部截断误差,? ? ? ?
? ? ? ?
2
//
1
34
/// ( 4 )
1
( ) ( ) (,( ) )
2
3 ! 4 !
n n n n n n
n n n
y y h y h f y h yx x x x x x
hh yy
x x x??
?
?
? ? ? ? ?
? ? ? ?
''
111
2 3 4
' '' ''' ( 4 ) '
2 3 4
' '' ''' ( 4 )
2 (,) ( ) 2 ( ) ( ) 2 ( )
( ) ( ) ( ) ( ) ( ) 2 ( )
2 3! 4 !
( ) ( ) ( ) ( ) ( )
2 3! 4 !
n n n n n n nn
n n n n n
n n n n
y y hf x y y x hy x y x h hy x
h h h
y x hy x y x y x y hy x
h h h
y x hy x y x y x y
?
?
???
? ? ? ? ? ? ?
? ? ? ? ? ? ?
? ? ? ? ? ?
L
L
11
1 1
11 11
1 1 11
() 1
,
4()
:
4
( ) ( ) ;
5
1
( ) ( ),
5
nn
n n
nn nn
n n nn
y x y
y x y
y x y y y
y x y y y
??
? ?
?? ??
? ? ??
?
??
?
? ? ? ?
? ? ?
由上两式有
由此导出下列事后估计式
111
3
'"
11
( ),
( ) ( ),
12
nnn
n n n
y y y x
hy x y y x
???
??
?
? ? ?
对校正公式：假定 是准确的，即
其局部截断误差为：
'
1 1 1
2
' ' ' '' '''
23
' '' '''
[ (,) (,) ] ( ) [ ( )
22
( ) ] ( ) [ ( ) ( ) ( ) ( )
22
] ( ) ( ) ( ) ( )
24
n n n n n n n n
n n n n n n
n n n n
hh
y y f x y f x y y x y x
hh
y x h y x y x y x hy x y x
hh
y x hy x y x y x
? ? ?
? ? ? ? ?
? ? ? ? ? ? ?
? ? ? ? ? ?LL
'
11
11
'
1 1 1
''
11
,
,:
2 ;
4
( ) ;
5
(,) ;
( ) ;
2
n n n
n n n n
n n n
n n n n
p y hy
m p P c
m f x m
h
c y y m
??
??
? ? ?
??
??
? ? ?
?
? ? ?
利用估计出的误差作为计算结果的一种补偿 有可
能使精度得到改善此时的计算公式为
预测
改进
计算
校正
1 1 1 1
'
1 1 1
1
( ) ;
5
(,),
n n n n
n n n
y c p c
y f x y
? ? ? ?
? ? ?
? ? ?
?
改进
计算
§ 3.龙格 — 库塔方法
一,Runge-Kutta法的基本思想（ 1）
1
2
' " ( )
11
' '' ' '
( )
( ) ( ) ( ) ( ) ( )
2 ! !
( ) (,),( ) (,) (,) (,),
n
p
p
n n n n n n
xy
p T ay l or y x
hh
y y x y x hy x y x y x
P
y x f x y y x f x y f x y f x y
?
??
? ? ? ? ? ?
? ? ?
L
LL
若用 阶 多项式近似函数 有：
其中 。
但由于公式中各阶偏导数计算复杂，不实用。
' ( 0 )
( 0 ) ( 0 )
'' ( 1 )
( 1 ) ( 1 )
''' ( 2 )
( 2 ) ( 2 )
( ) ( 1 ); 2,3,
jj
jj
y f f
ff
y f f
xy
ff
y f f
xy
ff
y f f j
xy
??
?
??
??
? ? ?
??
??
? ? ?
??
??
? ? ? ?
??
L
一般地有
11
1
1 1 2
1
21
(,)
11
()
22
(,)
(,)
nn
nn
nn
nn
nn
Eule r Eule r
y y hK
Eule r
K f x y
y y h K K
Eule r
K f x y
K f x h y hK
?
?
???
?
??
?
? ? ?
?
?
?
?
?
? ? ? ?
?
如果将 公式与改进 公式写成下列形式：
公式
改进 公式
Runge-Kutta法的基本思想（ 2）
11
(,)
()
(,)
(,)
nn
f x y
y x y
f x y
f x y
??
以上两组公式都使用函数 在某些点上的
值的线性组合来计算 的近似值 。
Euler 公式：每步计算一次 的值，为一阶方法。
改进Euler公式：需计算两次 的值，二阶方法。
Runge-Kutta法的基本思想（ 3）
(,)
(,) ( )
n n n
f x y
x y Tay lo r y x x Tay lo r
于是可考虑用函数 在若干点上的函数值的
线性组合来构造近似公式，构造是要求近似公式在
处的 展开式与解 在 处的 展开
式的前面几项重合，从而使近似公式达到所需要的阶
数。即避免求偏导，又提高了方法的精度，此为RK 方
法的基本思想。
1
1
1
1
1
(,)
(,) ( 2,3,)
p
n n i i
i
nn
i
i n i n ij j
j
y y h c K
K f x y
K f x a h y h b K i p
?
?
?
?
?
???
?
?
??
?
? ? ? ? ?
??
?
? L
二、二阶龙格－库塔方法
1
1
1
1
1
RK
(,)
(,) ( 2,3,)
,
(,) ( )
p
n n i i
i
nn
i
i n i n ij j
j
i ij i
n n n
y y h c K
K f x y
K f x a h y h b K i p
a b c
x y T ay l or y x x T ay l or
?
?
?
?
?
??
?
?
?
??
?
? ? ? ? ?
?
?
?
? L
一般地,方法设近似公式为
其中, 都是参数，确定它们的原则是使近似公式
在 处的 展开式与 在 处的 展开式
的前面项尽可能多地重合。
1 1 1 2 2
1
2 2 21 1
()
2 (,)
(,)
nn
nn
nn
y y h c K c K
K f x y
K f x a h y hb K
? ? ? ??
?
?? ?
? ? ? ?
?
当p 时，近似公式为
1 1 2 2 21
'
1 2 2
'3
21
12
''
2 2 21
(,)
[ (,) (,(,) ) ]
{ (,) [ (,) (,)
(,) (,) ] } ( )
( ) (,)
[ (,) (,) (,
nn
n n n n n n n n
n n n n n x n n
y n n n n
n n n
x n n y n n n
x y T ay l or
y y h c f x y c f x a h y hb f x y
y h c f x y c f x y a hf x y
hb f x y f x y O h
y c c f x y h
c a f x y b f x y f x y
?
? ? ? ? ?
? ? ? ?
??
? ? ?
??
上式在 处的 展开式为
23
) ] ( )
n
h O h?
1
2
' " 3
1
2
' ' 3
()
( ) ( ) ( ) ( ) ( )
2
(,)
[ (,) (,) (,) ] ( )
2
nn
n n n n
n n n
x n n y n n n n
y x x Ta y lo r
h
y x y x h y x y x O h
y f x y h
h
f x y f x y f x y O h
?
?
? ? ? ?
??
? ? ?
在 处的 展开式为
12
22
2 21
1
1 / 2
1 / 2 ),
cc
ca
c b O
???
?
?? ?
? ?
?
3
有无穷多组解，每一组解得一
近似公式，局部截断误差均为
(h 这些方法统称二阶方法。
1 2 2 2 1
1 1 2
1
21
1
,1,
2
( ) / 2
(,)
(,)
nn
nn
nn
c c a b E u le r
y y h K K
K f x y
K f x h y h K
?
? ? ? ?
? ? ??
?
??
? ? ? ?
?
取 此为改进 公式。
近似公式为
1 2 2 2 1
12
1
21
1
0,1,,
2
(,)
( 2,2 )
nn
nn
nn
c c a b
y y h K
K f x y
K f x h y h K
?
? ? ? ?
???
?
??
?
? ? ??
取 此为常用的二阶公式，
称为中点公式。
／／
三、三阶龙格－库塔方法
1 1 2 3
1
21
3 1 2
3
( 4 )
6
(,)
(
,)
22
(,2 )
nn
nn
nn
nn
p
RK
RK
h
y y K K K
K f x y
hh
K f x y K
K f x h y hK hK
?
?
?
? ? ? ?
?
?
??
?
?
? ? ?
?
?
? ? ? ?
?
类似地，对,即三个点，通过更复杂的计算，
可导出三阶 公式。
常用的三阶 公式为：
四、四阶龙格－库塔方法
1 1 2 3 4
1
21
32
43
4
( 2 2 )
6
(,)
(,)
2
2
(,)
22
(,)
nn
nn
nn
nn
nn
p R K
RK
h
y y K K K K
K f x y
hh
K f x y K
hh
K f x y K
K f x h y hK
?
?
?
? ? ? ? ?
?
?
?
?
?
?
? ? ??
?
?
? ? ?
?
?
? ? ??
?
对,即四个点，可导出四阶 公式。
常用的四阶 公式为：
'
11
1 1 2 3 4
1
21
32
43
h= 0,2,x= 0 x= 1
2
( 0 1 ) ;
( 0) 1.
(
6
( 2 2 )
6
(,)
(,)
22
(,)
22
(,)
nn
nn
nn
nn
nn
nn
x
y y x
y
y
h
y y K
h
y y K K K K
K f x y
hh
K f x y K
hh
K f x y K
K f x h y hK
?
?
?
? ? ? ?
?
?
?
?
?
??
?
? ? ? ? ?
?
?
?
?
?
?
?? ? ?
?
?
?
? ? ?
?
?
? ? ??
?
例：设取步长 从 直到 用四阶龙格—库塔
方法求解初值问题
解：由经典的四阶龙格—库塔公式得
2 3 4
1
21
1
32
2
43
3
2 2 ) ;
2;
2;
2
2
2;
2
2
2( )
.
n
n
n
n
n
n
n
n
n
n
n
n
K K K
x
Ky
y
xhh
K y K
h
yK
xhh
K y K
h
yK
xh
K y hK
y hK
?
? ? ?
?
?
?
??
?
?
?
?
? ? ?
??
??
?
? ?
? ? ?
?
??
?
? ?
? ? ?
?
??
?
1 R K
RK
4 6 R K 5
2 R K
R K
RK
R K
两点说明：
）当p=1,2,3,4 时,公式的最高阶数恰好是p,
当p>4 时,公式的最高阶数不是p,如p=5 时仍
为,p= 时 公式的最高阶数为 。
） 方法的导出基于Tay lor 展开，故要求所求问
题的解具有较高的光滑度。
当解充分光滑时，四阶 方法确实优于改进
Euler 法。对一般实际问题，四阶 方法一般可达
到精度要求。
如果解的光滑性差，则用四阶 方法解的效果
不 如改进Eul er法 。
五、变步长的龙格 — 库塔方法
()
1
( ) 5
11
1
5
( 2 )
1
5
( 2 )
11
( 2 )
11
()
11
,
( ),
2
,,
2
( ) 2,
2
() 1
.
( ) 1 6
h
n
h
nn
nn
h
n
h
nn
h
nn
h
nn
h
y
y x y ch
h
xx
h
yc
h
y x y c
y x y
y x y
?
??
?
?
??
??
??
??
??
??
??
??
??
??
??
?
?
?
n
／
／
／
以经典四阶龙格—库塔公式为例。从节点x 出发，以 为
步长求一近似值
将步长折半，即取 为步长从 跨两步到,求一近似
值 每跨一步的截断误差是 因此有
由上两式
( 2 ) ( 2 ) ( )
1 1 1 1
1
( ) [ ].
15
h h h
n n n n
y x y y y
? ? ? ?
? ? ? ?
／／
R-K方法的绝对稳定区域
/
2
211
12
23
1
3
3 2 4 3
2
(,),(,)
22
(,)
(,)
1
,
22
11
,(,)
)
2 2 4
2
2
() hK
( ( ) ( )
n n n n
nn
n n n
n
nn
n
f x y y R K
h
h h h
hh
K f x y K f x y K
hh
K f x y K K f x
hh
K
h
h y h
y
y y y hhK K
y y h hhK K
?
? ? ?
??
?
??
? ? ? ?
? ? ? ? ?
? ? ?
? ? ? ?
? ? ? ? ?
? ? ? ?
? ? ? ?
??
? ? ? ? ?
?
?
?
??
代入 公式：将
234
4 3
11
h )
24
( ( ) ( ) ( )
nn
hhy y h h h
hKK
?? ? ? ?
??
? ? ? ? ? ?
??
??
2 3 411
2 3 4
2 )
6
1 1 1
1
2 6 2 4
( 2
( ) ( ) ( )
nn
n
h
h
y y K K K K
y h h h? ? ? ?
?
? ? ? ? ?
??
? ? ? ? ???
??
?
?
??
?
? ?????
?
)()()( 432
1 24
1
6
1
2
11 hhhh
nn
????? ?则
1
24
1
6
1
2
1
1 )()()(
432
????? hhhh ????
绝对稳定区域：
2
1
- 3 - 2 - 1
0
- 1
- 2
§ 4.线性多步法
1
1 - r
1
RK
,,
nn
n n n
n
yy
y y y
y
?
?
L
－
单步法在计算 时，只用到前一步的信息 。为提高精度，
需重新计算多个点处的函数值，如 方法，计算量较大。
如何通过较多地利用前面的已知信息，如,, 来
构造高精度的算法计算,这就是多步法的基本思想。
1 1 1
1
01
11
(,),,(,)
,(,) ( 1,,,)
00
T a y l o
n n n r n n n r n r
rr
n i n i i n i
ii
i i k k k
y y y f x y f x y
y y h f
f f x y k n n n r
??
??
??
? ? ? ? ? ?
? ? ?
? ? ?
??
? ? ? ?
?
??
LL
L
－－
多步法中最常用的是线性多步法，它的计算公式中只
出现,,, 及 的一次
项，其一般形式为
其中 均为常数,。
若 ＝,显式；,隐式。
构造线性多步公式常用 r 展开和数值积分方法。
线性多步公式的导出
nn
Ta y l o r
x Ta y l o r ) x Ta y l o r
,ii??
n+1
利用 展开导出的基本方法是：将线性多步公式
在 处进行 展开，然后与y ( x 在 处的 展
开式相比较，要求它们前面的项重合，由此确定参数 。
1
01
1 0 1 1 1 1 0 1 1
1 ( )
( )
rr
n i n i i n i
ii
n n n n n n
r y x
y y h f
y y y h f f f
??
? ? ? ? ?
? ? ?
? ? ?
? ? ? ? ?
?
? ? ?
? ? ? ? ?
??
以 为例：设初值问题的解 充分光滑，待定的
两步公式为
( ) ( )
'' ( )
'2
1
( ) ( 1,2,),( )
( ) ( ) ( ) ( )
2!
( ( ) )
kk
n n n
p
pnn
n n n n n
p
n
y y x k y x x T ay lor
yy
y x y y x x x x x x
p
O x x
?
??
? ? ? ? ? ? ? ?
??
L
L
记 则 在 处的 展开为
'
'' '''
' 2 3
1
( 4 ) ( 5 )
4 5 ( 6 )
' ' ''
1 1 1 1
( ),
( ) (,) ( ),
( )
2 ! 3!
( )
4 ! 5 !
(,) ( )
ii
i i i
nn
n n n n
nn
n n n n n n
n y y x
y x f x y i n
yy
y y x h y y h h h
yy
h h O h
f f x y y x y y h
?
? ? ? ?
?
??
? ? ? ? ? ?
? ? ?
? ? ? ?
假设前 步计算结果都是准确的，即
则有
'''
2
( 4 ) ( 5 )
3 4 ( 5 )
'
'''
' ' '' 2
1 1 1 1
( 4 ) (
3
2!
( )
3! 4 !
(,)
(,) ( )
2!
3!
n
nn
n n n n
n
n n n n n n
nn
y
h
yy
h h O h
f f x y y
y
f f x y y x y y h h
yy
h
? ? ? ?
?
? ? ?
??
? ? ? ? ?
??
5)
4 ( 5 )
()
4!
h O h?
' '' 21
1 0 1 1 1 0 1 1 1
''' 3 ( 4 ) 41 1 1 1 1 1
( 5 ) 5 61 1 1
( ) ( ) ( )
2
( ) ( )
6 2 2 2 4 6 6
( ) ( )
1 2 0 2 4 2 4
n n n n
nn
n
y y y h y h
y h y h
y h O h
?
? ? ? ? ? ? ? ?
? ? ? ? ? ?
? ? ?
? ? ?
??
?
? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ?
? ? ? ? ?
将以上各公式代入并整理，得
1
'' ( 5 )
' 2 5 6
1
()
( ) ( )
2 ! 5 !
p + 1
nn
nn
n n n
p y x x
T a y lo r
yy
y x y y h h h O h
?
?
? ? ? ? ? ?L
为使上式有 阶精度，只须使其与 在 处的
展开式
的前 项重合。
01
0 1 0 1
1 1 1
1 1 1
1 1 1
1
1
11
22
1 1 1 1
6 2 2 6
1 1 1 1
24 6 6 24
aa
a
a
a
a
? ? ?
??
??
??
?
?
?
?
???
?
? ? ? ? ?
?
?
? ? ?
?
?
? ?
?
? ? ? ?
?
?
?
? ? ?
??
5,5
,1ii P?? ?
个参数只须 个条件。由推导知，如果选取参数
，使其满足前 个方程（p = 1,2,3,4 ),则
近似公式为p 阶公式。
1 1 0 1
11
1
1,0,,0
2
( )
2
n n n n
h
y y f f
? ? ? ? ?
?
??
? ? ? ? ?
? ? ?
0
如 满足方程组前三个
方程，故公式
此为二阶公式。
0 1 1 1 0
1 1 1 1
5 ( 5 ) 6
1
14
0,1,,
33
( 4 )
3
1
( )
90
n n n n n
nn
h
y y f f f
R h y O h
? ? ? ? ?
?
? ? ? ?
?
? ? ? ? ?
? ? ? ?
? ? ?
又如：解上面方程组得
相应的线性二步四阶公式（S i m p s o n 公式) 为
其截断误差为
由此可知，线性二步公式至多是四阶公式。
1
0
1
01
1
1
1
2 3
()
1
( ) ( ) 1 ( 1,2,)
[ 1 ( ) ( 1 ) (
n n n
r
i
i
rr
kk
ii
ii
r
p
ni
i
r
x T ay l or y x x T ay l or
i k i k p
p
R i p
?
??
?
?
?
?
? ? ?
?
?
?
?
?
?
?
?
?
?
? ? ? ? ?
?
?
? ? ? ? ?
?
??
?
L
一般地，线性多步公式中有 个待定参数，如令其
右端在 处的 展开式与 在 处的 展开式
的前p+1项系 数对应相等，可得方程组
其解所对应的公式具有 阶精度，局部截断误差为
1
( 1 )
1
2
)]
( 1 ) !
( )
pr
pp
in
i
p
h
iy
p
Oh
?
?
?
??
?
?
?
?
?
显然，线性多步公式至多可达到2r +2 阶精度。
二、常用的线性多步公式
1 2 3 1
0
1
01
0 0 1 2 3
1 1 2 3
( A da m s)
r =3 0,
1
( ) ( ) 1 ( 1,2,3 4)
55 59 37 9
=1,,,
24 24 24 24
( 55 59 37 9
24
r
i
i
rr
kk
ii
ii
n n n n n n
i k i k
h
y y f f f f
? ? ? ?
?
??
? ? ? ? ?
?
?
?
? ? ?
? ? ? ?
? ? ? ?
?
?
?
?
?
?
? ? ? ? ?
?
?
? ? ? ? ? ?
? ? ? ? ?
?
??
（一）阿达姆斯 公式
取,并令 由方程组
，
可解得,
相应的线性多步公式为
1
)
=0 A da m s?
－
因,此式称为 显式公式，是四阶公式.
533
5 4 ( 5 ) 6
1
11
5 ( 5 ) 6
[ 1 ( ) 5 ( ) ] ( )
5!
251
()
720
n i i n
ii
n
h
R i i y O h
h y O h
??
?
? ? ?
? ? ? ? ? ?
??
??
局部截断误差为
1 2 3 3
0 1 0 1 2
1 1 1 2
5 ( 5 ) 6
1
0,
9 19 5 1
= 1,,,
24 24 24 24
( 9 19 5 )
24
A da m s
19
()
720
n n n n n n
nn
h
y y f f f f
R h y O h
? ? ? ?
? ? ? ? ?
?
? ? ? ?
?
? ? ? ?
? ? ? ? ?
? ? ? ? ?
? ? ?
如果令 由方程组可解得
，
相应的线性多步公式为
称其为四阶 隐式公式，其局部截断误差为
利用数值积分方法求线性多步公式
11
1
1 + 1 1
( ) ( ) (,( ) ) ( )
,,,,,,( )
( ) ( ),1
nn
nn
xx
nn
xx
n n n r n n n r
y x y x f x y x dx F x dx
x x x x x x F x
x F x r?
??
?
? ? ? ?
? ? ?
?
??
LL
r
基本思想是首先将初值问题化成等价的积分形式
用过节点 或 的 的r
次插值多项式 代替 求积分即得 阶的
线性多步公式。
1 2 3
3
3
0
1 2 3
3
0
3,,,,( )
( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( 0,1,2,3 )
( ) ( )
n n n n
i n i
i
n n n n
i
n i n i n j
j
ji
r x x x x F x
L x l x F x
x x x x x x x x
l x i
x x x x
? ? ?
?
?
? ? ?
? ? ?
?
?
?
?
? ? ? ?
??
??
?
?
例如 时，过节点 的三次
插值多项式为
其中
11
1
1
1
3
13
0
1 2 3
3
23
1
3
13
2
3
3
( ) ( ) ( ) [ ( ) ] ( )
( ) ( ) ( )
()
6
( ) ( ) ( )
()
2
( ) ( ) ( )
()
2
( ) (
()
nn
nn
n
n
n
n
n
n
xx
n n i n i
xx
i
x
n n n
n
x
x
n n n
n
x
x
n n n
n
x
n
n
y x y x L x d x l x d x F x
x x x x x x
F x d x
h
x x x x x x
F x d x
h
x x x x x x
F x d x
h
x x x
Fx
??
?
?
?
??
?
? ? ?
??
?
??
?
?
? ? ?
? ? ?
?
? ? ?
?
?
? ? ?
?
?
?
???
?
?
?
1
12
3
1 2 3
) ( )
6
[ 5 5 ( ) 5 9 ( ) 3 7 ( ) 9 ( ) ]
24
n
n
x
nn
x
n n n n
x x x
dx
h
h
F x F x F x F x
?
??
? ? ?
??
?
? ? ? ?
?
11
1 1 2 3
3
,( ),( ),(,)
( ) (,( )) (,1,2,3 ),
(55 5 9 3 7 9 )
24
[,]
n n n n k k k
k k k
n n n n n n
nn
y y y x y x f x y
F x f x y x k n n n n
h
y y f f f f
Ad a ms
x x Ad a ms
??
? ? ? ?
?
? ? ? ? ?
? ? ? ? ?
对上式用 代替 用 代替
则得
这就是四阶 显式公式。由于积分区间在插值
区间 外面，又称为四阶 外插公式。
11
1
( 4 ) ( 5 )33
1
00
31
( 5 ) 3
5 ( 5 )
1
0
( ) ( )
( ) ( )
4 ! 4 !
,),
( ) 25 1
( ) ( )
4 ! 72 0
nn
nn
n
n
xx
xx
n n j n j
xx
jj
nn
x
n n j
x
j
Fy
R x x dx x x dx
xx
y
R x x dx h y
??
?
?
?
??
?
? ? ?
??
??
??
?
? ? ? ?
?
? ? ?
????
??
由插值余项公式可得其局部截断误差为
由积分中值定理，存在 （ 使得
1 1 2
2
3
1
1 1 2
3
1
1
,,,( )
( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( 1,0,1,2)
( ) ( )
()
n n n n
i n i
i
n n n n
i
n i n i n j
j
ji
nn
x x x x F x
L x l x F x
x x x x x x x x
l x i
x x x x
F x A da m s
yy
??
?
??
? ? ?
? ? ?
??
?
?
?
? ? ? ?
? ? ?
??
?
?
?
＋
同样，如果过节点 的 三次插值多
项式为
其中
代替 求积分，即得四阶 隐式公式
1 1 2
5 ( 5 )
1 2 1
21
( 9 19 5 )
24
19
()
720
[,]
n n n n
n n n n
nn
h
f f f f
R h y x x
x x A da m s
A da m s
??
? ? ?
? ? ?
??
? ? ? ?
? ? ? ?其局部截断误差为
由于积分区间在插值区间 内，故 隐式
公式又称为 内插公式
1 2 1
3 0 1 2 3
1 3 1 2
5 ( 5 ) 6
1
( )
0,
8 4 8
1,,,,0,
3 3 3
4
( 2 2 )
3
14
( )
45
n n n n n
nn
M i l i ne
y y h f f f
M i l i ne
R h y O h
M i l i ne
? ? ? ?
? ? ? ? ?
?
? ? ? ?
?
? ? ? ?
? ? ? ? ? ?
? ? ? ?
??
0
（二）米尔尼 公式
取r =3,并令 由方程组可
解出
相应的线性多步式
称为 公式，其局部截断误差为
公式是四阶四 步显式公式。
12
1 2 1 1
5 ( 5 ) 6
1
( m in )
0,m i n
13
( 9 ) ( 2 )
88
1
( )
40
n n n n n n
nn
H am g
H am g
y y y h f f f
R h y O h
??
? ? ? ?
?
??
? ? ? ? ?
? ? ?
（三）哈明 公式
取r=2,并令 可得到 公式
其局部截断误差为
Hamming 公式是四阶三步隐式公式。
1
1
,
)
,
)
n
n
y
?
?
n+1
n+1
n+1
隐式法与显式法的比较
一般地，同阶的隐式法比显式法精确，而且数值
稳定性也好。但在隐式公式中，通常很难解出y 需
要用迭代法求解，这样又增加了计算量。
在实际计算中，很少单独用显式公式或隐式公式，
而是将它们联合使用：先用显式公式求出y(x 的预
测值，记作y 再用隐式公式对预测值进行校正，求
出y (x 的近似值 。
§ 5.预测 — 校正系统
用显式公式计算预测值，然后用隐式公式进行校正，
得到近似值 yn+1这样一组计算公式称为预测 — 校正系统。
一般采用同阶的隐式公式与显式公式。常用的预
测 — 校正系统有两种,
1 2 31
1 1 1 21
( 5 5 5 9 3 7 9 )
24
[ 9 (,) 1 9 5 ]
24
n n n n nn
n n n n n nn
A d a m s
h
y y f f f f
h
y y f x y f f f
? ? ??
? ? ? ??
?
? ? ? ? ?
??
?
? ? ? ? ? ?
??
预测—校正公式
预测
校正
3 1 21
1 2 1 11
l i m i n
4
( 2 2 )
3
13
( 9 ) [ (,) 2 ) ]
88
n n n nn
n n n n n nn
M i ne H am g
y y h f f f
y y y h f x y f f
? ? ??
? ? ? ??
?
?
? ? ? ?
?
?
?
? ? ? ? ? ?
??
预测—校正公式
( 1)
RK
3
说明：
以上两种预测—校正公式均为四阶公式，其起步值
通常用四阶 公式计算。
(2) 有时为提高精度，校正公式可迭代进行多次，但迭
代次数一般不超过 次。
用局部截断误差进一步修正预测－校正公式
5 ( 5 ) 6
1 1
5 ( 5 ) 6
11
5 ( 5 ) 6
1 1
5 ( 5 )
1 1
11 11
1
251
( ) ( )
720
19
( ) ( )
720
270
()
720
720
()
270
251
( ) ( )
270
()
nn n
n n n
nn n
nn n
nn nn
n
Adams
y x y h y O h
y x y h y O h
y y h y O h
h y y y
y x y y y
yx
? ?
??
? ?
? ?
?? ??
?
? ? ?
? ? ? ?
? ? ? ?
? ? ?
? ? ? ?
?
由 公式的局部截断误差公式
两式相减
11 1
19
()
270
nn n
y y y
?? ?
? ? ?
1 1 2 3
11
1 1 1 1 2
1 1 1 1
( 55 59 37 9 )
24
251
()
270
[ 9 (,) 19 5 ]
24
19
()
270
n n n n n n
n n n n
n n n n n n n
n n n n
A da m s
h
p y f f f f
m p c p
h
c y f x m f f f
y c c p
? ? ? ?
??
? ? ? ? ?
? ? ? ?
?
? ? ? ? ?
?
?
?
? ? ?
??
?
?
? ? ? ? ?
?
?
? ? ? ?
??
由上面就得到多环节的 预测－校正公式
预测
改进
校正
改进
1 3 1 2
11
1 2 1 1 1
1 1 1 1
l n m i n
4
( 2 2 )
3
112
()
121
13
( 9 ) [ (,) 2 ]
88
9
()
121
n n n n n
n n n n
n n n n n n n
n n n n
M i e H am g
p y h f f f
m p c p
c y y h f x m f f
y c c p
? ? ? ?
??
? ? ? ? ?
? ? ? ?
?
?
? ? ? ?
?
?
?
? ? ?
??
?
?
? ? ? ? ?
?
?
? ? ? ?
?
?
完全类似，可以导出多环节的
预测－校正公式
预测
改进
校正
改进
00
1 1 1
1
1 1 2 1 1
1,,(,),,( 2),,1
( 3 ) (,)
(,)
22
n n n
n n n
ba
a b f x y N y h x a n
N
f f x y
hK
K hf K hf x y
K
? ? ?
? ? ?
?
? ? ?
?
? ? ? ?
算法
（ ）输入 置
计算
2
3 1 1 4 1 1 3
1 1 2 3 4
00
(,) (,)
22
1
( 2 2 )
6
( 4) (,)
( 5 ) 3,1,3
1,0,0,6
n n n n
n n n
nn
hK
hf x y K hf x h y K
y y K K K K x a nh
xy
n n n
n n p c
? ? ? ?
?
? ? ? ? ? ?
? ? ? ? ? ? ?
? ? ?
? ? ? ?
输出
若 置 返回 ；
否则，置 转 。
3 3 3 3
0 3 2 1 0 0
3 1 3 2
1 1 1
3 3 0
( 6 )
(,)
4 1 1 2
( 2 2 ) ( )
3 1 2 1
13
( 9 ) [ (,) 2 ]
88
9
( ),)
121
( 7 ) 1,,,
( 0,1,2 ),,,,
j j j j j j
f f x y x x h
p y f f f m p c p
c y y h f x m f f
y c c p x y
n N n n x x y y f f
j x x y y p p
? ? ?
? ? ?
? ? ? ? ? ? ?
? ? ? ? ?
? ? ?
? ? ? ? ? ?
? ? ? ?
计算
输出(
若,置
0
,6cc ? 转；
否则停机。
§ 5,常微分方程组与高阶方程的数值解
'
12
0
0
(,,,,) ( 1,2,,)
( ) ( 1,2,,)
i i N
ii
y f x y y y i N
y x y i N
? ???
?
????
LL
L
一、一阶方程组
一阶方程组的初值问题
12
'
0 0 0 0
12 0
0
12
y f
y (,,) ;
y f (,y ) ;
y (,,) ;
y ( ) y ;
f (,,) ;
T
N
T
N
T
N
y y y
x
y y y
x
f f f
?
? ?
?? ?
??
?
L
L
L
若对 和 采用向量的记号
'
0
0
n+ 1 n 1 2 3 4
1 n 2 n 1
3 n 2 4 n 3
y f (,y ) ;
y ( ) y ;
y y ( k 2 k 2 k k )
6
k f (,y ) ; k f (,y k ) ;
22
k f (,y k ) ; k f (,y k ) ;
22
nn
nn
x
x
h
hh
xx
hh
x x h h
? ?
?
?
?
? ? ? ? ?
? ? ? ?
? ? ? ? ? ?
求解这一初值问题的四阶龙格—库塔公式为
其中
,1 1 2 3 4
1 1 2
2 1 1 1 2 2 1 1
3 1 1 2
( 2 2 )
6
( 1,2,,),
(,,,,) ;
(,,,,) ;
2 2 2 2
(,,
22
i n in i i i i
i i n n n Nn
i i n n n Nn N
i i n n
h
y y K K K K
iN
K f x y y y
h h h h
K f x y K y K y K
hh
K f x y K
?
? ? ? ? ?
?
?
? ? ? ? ?
? ? ?
L
L
L
或表达为
其中
2 2 2 2
4 1 1 3 2 2 3 3
,,) ;
22
(,,,,) ;
n Nn N
i i n n n Nn N
hh
y K y K
K f x h y hK y hK y hK
??
? ? ? ? ?
L
L
? ? ? ?
? ? ? ?
'
0 0
'
00
,,
,,
dy
y f x y z y
dx
dz
z g x y z z
dx
yx
x z
?
? ? ?
?
?
?
? ? ? ?
??
以两个方程组为例
? ? ? ?
? ? ? ?
?,2,1,0
,,,
,,,
,1
001
0
0
1
?
?
?
?
?
?
???
???
?
?
n
zhg
yhf
E ul e r
zxzyxzz
yxzyxyy
n
n
nnn
n
n
n
nn
法方程组的
2) 方程组的 R-k 法
? ?
? ?
1 2 3 41
1 1 2 3 4
2 2
6
2 2
6
nn
nn
h
h
yy K K K K
z z M M M M
?
?
? ? ? ? ?
? ? ? ? ?
? ? ? ?11,,,,,Mnnnnnnfg yyxxK z z??
11
2
11
2
,,
2 2 2
,,
2 2 2
nn n
nn n
h
f
h
g
h K h M
yxKz
h K h M
yxMz
??
? ? ? ???
??
??
? ? ? ???
??
22
3
22
3
,,
2 2 2
,,
2 2 2
nn n
nn n
h
f
h
g
h K h M
yxKz
h K h M
yxMz
??
? ? ? ???
??
??
? ? ? ???
??
? ?
? ?
4 33
4 33
,,
,,
0,1,2,
nn n
nn n
fh
gh
n
yx hK hMKz
yx hK hMMz
? ? ? ?
? ? ? ?
? L
二、化高阶方程为一阶方程组
( ) ' ( 1 )
' ' ( 1 ) ( 1 )
0 0 0 0 0 0
' '' ( 1 )
1 2 3
'
1 2 1 0
(,,,,);
( ),( ),,( ),
,,,,,
()
mm
mm
m
m
m
y f x y y y
y x y y x y y x y
y y y y y y y y
m
y y y x y
?
??
?
? ?
?
? ? ?
?
? ? ? ?
??
L
L
L
高阶微分方程的初值问题可通过变量代换化为一阶
微分方程组的初值问题。设有 阶常微分方程初值问题
引入新变量
则可将 阶方程化为如下一阶方程组：
0
''
2 3 2 0 0
' ( 2 )
1 1 0 0
' ( 1 )
1 2 0 0
()
()
(,,,,) ( )
m
m m m
m
m m m
y y y x y
y y y x y
y f x y y y y x y
?
??
?
?
?
??
?
?
?
?
??
?
? ??
?
L
L
"'
'
''
0 0 0 0
'
00
''
00
1 1 2 3 4
1 1 2 3 4
(,,) ;
( ),( ),
,( ) ;
(,,),( ),
( 2 2 ) ;
6
( 2 2 )
6
,1,2,3,4)
nn
nn
ii
y f x y y
zy
y x y y x y
y z y x y
z f x y z z x y
h
y y K K K K
h
z z M M M M
K M i
?
?
? ?
?
?
??
?
? ??
?
?
??
?
?
? ? ? ? ?
?
?
?
?
? ? ? ? ?
?
?
?
以二阶方程的初值问题为例：
引入新变量
对其用四阶龙格－库塔公式
（ 计算公式略.
? ?
? ?
" 2 '
'
'
'
'2
0 0 0
0
2
1
1
5 sin 2 2 0 1
( 0 ) 2,( 0 ) 3
( 0 ) 2
,
5 sin 2 2 ( 0 ) 3
,,,
0,2,3
,
2
,
3
,2
,
M
x
x
nn
n
nn
n
n
n
RK
y e x y y x
yy
y z y
zy
z e x y z z
x y z R K
zf
h
f
g
yxKz
z
xK
yx
? ? ? ? ? ?
?
? ? ? ?
?
? ? ? ?
??
?
? ? ? ? ?
?
? ? ? ?
?
?
?
?
?
??
?
?
?
例：用四阶 方法求解初值问题 ( 取h = 0, 1 )
解：令
因 由 公式得
1
11
1
2
0
1
M 3,1
2
1,6 2 3 8 2 2 4
,
22
,,
2 2 2
n
nn
n
h
z
h
g
hK hM
y z
hK hM
yxMz
??
??
??
??
??
? ? ? ?
??
? ? ? ?
??
??
LL 继续计算（略）。
