数值计算方法-常微分方程(组).docx

科学计算一理论、方法
及其基于MATLAB的程序卖现与分析
微分方程（组」教值解法
§ 1常微分方程初值问题的数值解法
微分方程（组）是科学研究和工程应用中最常用的数学模型之一。
如揭示质点运动规律的Newton第二定律：
和刻n （19）
血）=%
例2 Euler方法与改进的Euler方法的比较
下图是当* = ：
open
§ Euler方法和改进的Euler方法的误差分析
由Taylor公式
y妃J = 乂九）+V（4）妃1 一九）+云矿'（4）妃1 一九）2 +… （迪）
=业）+—（4,旭））人+。倨）
说明Euler方法的截断误差是。伊），类似地，由
y 妃 i）=y（4）+f（4，y（4））打+»&贸+。（/） （20）
y^k ） = v 妃1） - •/妃1，布+1））"+§ V'妃1）层 + <?（/） （21）
以及
y"（4+i）= y"（4）+y"'（商+.•• =矿'妃1）=矿'（4）+。（力） （22）
让式（20）的两端减式（21）的两端，可得
y 妃1）=血）+g ［f^，））+y 妃i，y 妃i））］+:俨。3）+。便）
4 （23）
=血）+； ［f 妃血））+y 妃 1，y 妃1））］+。俯）
从上述推导Euler方法、改进的Euler方法的过程以及例1、例 2容易看出，改进的Euler方法Euler方法的精度高，其原因在于： 1在推导Euler方法时，我们是用待求解函数尸如在一点处的变化 率也，y（tk））代替y =如在区间妇京上的平均变彳匕率：
r+，= hf{t, y（F）） （24）
Jtk
2而在推导改进的Euler方法时，我们是用待求解函数y =市）在两点 处变化率的平均值；［儿,布））+ /■妃Q妃））］代替"业）在区间购扁上 的平均变化率；
显然，通常;［/■妃 y{tk）） + f（tk+l, ））］比 f（tk, y（tk））更接近于 y =业）在 区间此,匕J上的平均变化率。由此启发人们：适当地选取区间［今,知］上 函数y =业）若干点处的变化率，用它们加权平均值代替y =业）在区间 ［W+J上的平均变化率，近似解的精度应更高。
下面将要介绍的Runge- Kutta法就是基于上述想法得到的。
§ 2 Runge- Kutta 法
Runge- Kutta法是按选取区间L，知］上函数V =如变化率/'（3）的
个数的多少和截断误差的阶数O（hm）来区分的一系列方法，如
1 二阶的Runge- Kutta法（改进的Euler方法）
& = , y k)
（25）
（26）
< y顼心,％+咨)
%+1 = % +-(^1 + 匕)
2 三阶的Runge- Kutta法
'K[=低，%)
K2 = f\^k +5，为 +^K1 J
K3 = f*k+h，Vk +砍2)
h
>k+i = % + 工(K[ + 4K2 + K3)
L 6
四阶的Runge- Kutta法
1）古典形式
(27)
(28)
(29)
(30)
Y =/■妃 ％)
灼 捉J
匕=f(tk +h,yk+hK3)
%+i = % +*(Y + 2K? + 2K3 + 码)
I 6
2) Gill公式(具有减小舍入误差的优点)
妇=儿,光)
y =/"+?，% + 与1 hK} + 砍2
[22 2
( 皿 2 +而 )
K& = f h +h,yk - hK2 + —.一 hK3
\ J 7
y*+i = » +?(K1 +(2-扼)眩2 +(2 +zk +K4)
4 Runge—Kutta法的一般形式
n
yk+l = yk + "(4，%，") =处 + * atKi
i=l
其中，光，时称为增量函数(Increment Function)以及
K[ =f(tk,yk)
弓=f(tk+ pxh,yk+qnhK^
，K3 = + P2々，yk + QnhKi +q22hK2)
Kn = f^k + P 〃一也治 +q"K[ + ■■■ + qn_u_lhKii_l)
需要特别指出的是：在确定（29）、（30）中的参数a” p,.和心时,
应该使（29）的右端和适当阶的（如〃阶）Taylor展式一致，这样,
至少对于底阶的R-K法来说，参与加权的斜率个数"与方法的阶数是一 致的。
例如：一阶的R-K法即Euler法，局部的截断误差（Truncation Error）是。例），所以，当微分方程的阶是一次函数时是精确的；二阶 R-K法即改进的Euler法，局部的截断误