1、毕业论文文献综述 信息与计算科学 常微分方程初值问题的 Runge-Kutta 解法 一、前言部分 常微分方程在很多学科领域内有着重要的作用， 如 自动控制、各种电子学装置的设计、弹道的计算、飞机和导弹飞行的稳定性的研究、化学反应过程稳定性的研究等等，这些问题都可以化为求微分方程的解，或者化为研究解的性质的问题。 这些问题都包含某个变量关于另一个变量的变化。大多数这样的问题需要求解一个初值问题，即求解满足给定初值条件的微分方程 1 。我们 更多的是使用逼近原问题的解的方法来逼近原方程的解。因逼近方法给出更精确的结果和实际的误差信息，所以更经常被使用。一些典型的常微分初值问题的数值求解问题的方法
2、有：单步法和线性多步法。 在求解区间 ba, 上取定节点 bxxxxxa NN 1210 令 nnn xxx 1 ，称为积分网格的步长。用 Nyyy , 10 表示初值问题精确解 xy 在节点 Nxxx , 10 上函数值 Nxyxyxyxy , 210 的近似值。对给定的数值积分方法，ny 中的各个 ny 市按某一个递推算法确定的。一个递推算法，如果在用它计算 1ny 时只用到已经求出的诸值 nyyy , 10 中的 ny ，而无须使用其余值 110 , nyyy 中的任何一个，则称此算法为单步方法 2 。相反的则是多步法。单步法主要有欧拉法和 Runge-Kutta 法，多步法主要有 Ad
3、ams 法和 Milne 法等 3 。 本文综述 常微分初值问题 初值问题的数值解法及其误差估计（ 相容性、稳定性和收敛性分析 ），重点介绍了 Runge-Kutta 法 。 二、主题部分 2.1 常微分方程的初值问题概述 114 常微分方程在微积分概念出现后即已出现。从莱布尼兹专门研究用 变量变换解决一阶微分方程的求解问题的“求通解”时代，到 1841 年刘维尔的里卡蒂方程不存在一般初等解和柯西的初值问题的提出，常微分方程转向“求定解”时代。再到 19 世纪末天体力学的太阳系稳定性问题眼界需要常微分方程解的大范围性态，从而发展到了“求所有解问题”。最后到了 20 世纪六七十年代，常微分方程因
4、为计算机的发展进入“求特殊解”阶段，发现了具有新性质的特殊的解和方程，如混沌（解）、奇异吸引子及孤立子等等。 在数学分析（微积分）中研究了变量的各种函数及函数的微分与积分。如函数未知，但知道变量与函数的代数关系式，便 组成代数方程，通过求解代数方程解出未知函数。同样，如果知道自变量、未知函数及函数的导数（或微分）组成的关系式，得到的便是微分方程，通过求解微分方程求出未知函数。自变量只有一个的微分方程称为常微分方程。常微分方程是数学分析或基础数学的一个组成部分，在整个数学大厦中占据着重要位置。 在常微分方程中，我们知道只有少数简单的常微分方程才能用初等积分法求出它们的精确解析解，多数情形只能用近
5、似解法求其近似解。近似解法可以分成两类：解析方法和数值方法。其中前者寻求解的近似表达式，后者则是计算微分方程解在求解区域中一些离散 点上的近似值。 对于一阶的常微分方程初值问题 yxfdxdy , ， Tx0 00 yy （ 1） 的数值方法，其中 f 是 x 和 y 的已知函数， 0y 为给定的初值。为使上面的解存在、唯一且连续依赖初值 0y ，则 f 必须关于 y 满足 Lipschitz 条件：即存在常数 L，使得对于任意 Tx ,0 ，均成立不等式 2121 , yyLyxfyxf 单步显式公式的一般形式为 hyxhyy nnnn ,1 ， 00 yy （ 2） 称 hyx , 为增量
6、函数。一般来说微分方程问题（ 1）的精确解 nxy 不满足式（ 2），即一般地有 hxyxhxyxy nnnn ,1 定义 1 称 hxyxhxyxyR nnnnn ,11 为单步显式公式（ 2）在 1nx 处的截断误差。 一个求解公式的局部截断误差刻画了其逼近微分方程的准确程度。根据上述定义可以直接求的各式的局部截断误差。例如欧拉公式的局部截断误差为 11 ,n n n n nR y x y x hf x y x nnnnn xyhxyyhxyhxy 221 nyh 2211 nnn xx 定义 2 如果一个求解公式的局部截断误差为 1 ph ，则称该求解公式是 p 阶的，或称具有 p 阶精
7、度。 2.2 常微分方程的初值问题的数值方法 下面将先简单的介绍一下除 Runge-Kutta 法外的一些数值方法。 2.2.1 单步法中的欧拉法 116 单步法中最简单的数值解法是欧拉法。欧拉法的基 本思想是由 ny 推导出下一个节点1ny ，其推导公式为 nnnn yxhfyy ,1 ， ,1,0 . 其局部截断误差为 21 hxR nh 。欧拉法的几何意义是用一条过 00,yx 的折线来近似替代过 00,yx 的积分曲线。 虽然这个方法的计算很简单，但精确度较低，所以有了精确度较高的梯形方 法。推导公式为 111 ,2 nnnnnn yxfyxfhyy， ,1,0n ，它的局部截断误差为
8、 32 hxR nh 。此为隐型格式，迭代求解所需计算量大。 结合上述两种方法的优点，从而得到改进欧拉方法。其计算公式为 nnnnnnnn yxhfyxfyxfhyy ,2 11 ， ,1,0n 。 2.2.2 线性多步法 10,3 常用的多步法主要有 Adams 法和 Milne 法，其中 Adams 法又包括显式 Adams 法和隐式Adams 法。显式 Adams 法的线性 1k 步公式为 kj jnkjnn fhyy 01 其中 kjm mjkj jm 1 ， kj ,2,1,0 。 且隐式 Adams 公式为 km nmmnn fhyy0 1*1 其中 01* 1 dsmsmm ，
9、km ,1,0 2.3 Runge-Kutta 法 2.3.1 Runge-Kutta 法的构造 157 用 Taylor 级数法的公式 hyxhxy nnnn ,1 ， ,1,0n 可以构造高阶单步法，但增量函数 hyx , 是用 f 的各阶导数表示的，通常不易计算 。由于函数在一点的导数值可以用该点附近若干点的函数值近似表示，因此可以把 Taylor 级数法中的增量函数 改为f 在一些点上函数值的组合，然后用 Taylor 展开确定待定的系数，使方法达到一定的阶。这就是 Runge-Kutta 法的基本思想。或者是在每个区间多取几点，将它们的斜率加权平均作为平均斜率，就有可能构造出精度更高
10、的计算公式。通过构造高阶单步法来提高精度，而较高的精度意味着计算结果更加精确，误差会随着 h 的减小而迅速减小。 现在假设 Ni iinn kwhyy 11（ 3） 其中 nn yxfk ,1 ， 1, iinini khyhxfk ， ni ,3,2 （ 4） 此类递推算法为显示 n 级 Runge-Kutta 方法。上述算法公式中的系数 iw ， i 和 i 希望如此确定，使得（ 3）式的 Taylor 展开式中所有 h 幂次不超过 p 的那些项的系数与 Taylor 展开式 nhnppnnnn xRxyphxyhxyhxyxy !2 21 （ 5）（其中 11!1 pph yphR， 1
11、 nn xx ）中的相应项的系数相等。 当 n=1 时，则可知 1w =1，此时式（ 3）为欧拉公式。 当 n=2 时，确定系数 iw ， i 和 i 使得 nn xyhxy 与 Ni iikwh 1两者的 Taylor展开中含 rh 项的系数 ， pr ,2,1 ，对尽可能大的 p 相等。 利用二元的 Taylor 展开公式 nnnnnn yxfykxhyxfkhyhxf , nn yxfykxh ,!222 可将（ 3）式表示为 32222211 hfffwhfwwhyy nnynxnnn （ 6） 再将 Taylor 展式（ 3）改成 321 !2 hfffhhfxyxy nnynxnn
12、n （ 7） 比较（ 6）和（ 7）式中 2,1rhr 项的系数，得到 121 ww ， 2122 w， 2122 w（ 8） 因为有 4 个待定系数，可是就只有 3 个方程数目。因此具有无穷多个数组 2221 , ww 满足条件（ 8），即存在无穷多个二级 Runge-Kutta 方法，它们的截断误差 3hxR nh 。 若取 212w，便得到 211 2 kkhyy nn nn yxfk ,1 12 , hkyhxfk nn 或 nnnnnnnn yxhfyxfyxfhyy ,2 11 这是改进欧拉公式。 因为二级的 RK 方法最多是二阶的，若想得到更高阶的方法，必须增加级数，即取较大的
13、n。当 n 较大时，计算也就会更加的复杂。 当 n=3 时，和 n=2 的推导过程相相似，当 43,0,41321 www时，得到 311 34 kkhyy nn nn yxfk ,1 12 31,31 hkyhxfk nn 23 32,32 hkyhxfk nn或当 61,32,61321 www时，有 3211 46 kkkhyy nn nn yxfk ,1 12 21,21 hkyhxfk nn 213 2, hkhkyhxfk nn 同理， n=4 时，构造经典的四阶 RK 方法，其计算公式为 43211 2261 kkkkhyy nn nn yxfk ,1 12 21,21 hkyh
14、xfk nn 23 21,21 hkyhxfk nn 34 , hkyhxfk nn （ 9） 其局部截断误差为 5h 。 由于 RK 方法是通过对函数的 Taylor 展开来构造的，因此高阶 RK 方法要求微分方程的解具有高阶导数。如果解的光滑性较差，即 使用高阶方法也不可能得到高精度的数值解，此时 应当用低阶方法。或者说从这些上述公式的讨论中可知，级数 n 表示计算函数值 f 的次数，当 4,3,2n 时， n 级的 RK 方法的阶数 nnp ，即计算函数值 f 的次数越多， RK 方法的精度即阶数越高。但是当 4n 时，情况就大不一样了， J.C.Butcher 指出：当 7,6,5n时
15、， 1nnp ；当 9,8n 时， 2nnp ；而且 ， 810p 。由此可见四级以上的RK 方法，计算函数值 f 的工作量增加较快，而阶数即精度提高较慢，故而，在实际问题中比较常用的是四级 RK 方法。 2.3.2 Runge-Kutta 法的误差分析 1311 I、相容性 当 phhxyL ; ， 2 时，称数值方法是相容的。 由于 hxyxh xyhxyhxyL ,; 2hhxyxyhxy hyxfwkwhxyx ni ini ii 11 , 于是 hyxfwxyhxyL ni i 1 ,; hyxfwxy ni i 1 , hyxfwxy ni i ,1 1是微分方程的解当 hwni
16、i 11故有结论： R-K 法相容 11 ni iw。即 p 阶的 Runge-Kutta 方法是相容的。 II、收敛性 由于截断误差的存在，计数值 ny 只是精确解 nxy 的一个近似。故我们可以相当这样的问题：差 nn yxy 是不是随着步长 h 的减小而减小。这就是收敛性问题。 定义 3 如果一个单步法满足：对 Tx,0 中任意一个固定的点 a 成立 ayynanhx h 00lim则称该方法是收敛的。 由于 anhx 0 ，当 0h 时， n ，也即当步长 h 减小时，从 0x 计算到点 a 的步数就增大，故有下面的定理。 定理 1 若公式（ 2）是 p 阶方法，且增量函数 hyx ,
17、 关于 y 满足 Lipschitz 条件，则该方法收敛。 对于 Runge-Kutta 方法，当函数 yxf , 关于 y 满足 Lipschitz 条件，可用数学归纳法证明所有的 ik 关于 y 满足 Lipschitz 条件，从而 增量函数 hyx , 关于 y 满足 Lipschitz 条件，所以 Runge-Kutta 方法收敛。而且对于这类满足 Lipschitz 条件的单步法而言，相容条件是收敛条件的充分必 要条件。 III、稳定性 下面讨论稳定性问题，更准确地说就是数值解关于初值的连续依赖性问题。 定义 4 如果一个单步法满足：存在正常数 c 和 0h ，使得用任意步长 0,0
18、hh 以及任意两个初值 0y 和 0g 计算得到的两组数值 nyyy , 21 和 nggg , 21 成立 001m a x gycgy iini ， 则称该方法是稳定的，其中 hTn。 定理 2 若 增量函数 hyx , 关于 y 满足 Lipschitz 条件，则公式（ 2）是稳定的。 由上可知，欧拉法，改进欧拉法和 Runge-Kutta 方法都是稳定的。但是因为上面的稳定性定义中没有考虑计算过程中的舍入误差，使得最后得到 的结果可能会很差，所以又有了任何一步产生的误差在以后过程中足部减少的绝对稳定性。它不仅与函数 f 有关，还和 步长 h有关。 对于模型问题： ydxdy ， Tx0
19、 00 yy （ 10） 其中 是常数（可以是复数）。 定义 5 对初值问题（ 10），记 hh ， 若一个数值方法在某一步由计算产生的误差在以后各步中逐步减少，则称该方法关于 h 是绝对稳定的。复平面上所有这样的 h 组成的区域称为该方法的绝对稳定区域。 由此可知，当 1!0 niiih满足时， n 级 n 阶的 Runge-Kutta 方法 (n=1,2,3,4)是绝对稳定的。 如果一个方法的绝对稳定区域比较大，则可说它的绝对稳定性比较好，且隐式方法的绝对稳定性比显式方法好；对于 Runge-Kutta 方法，阶愈高绝对稳定性愈好 11 。 三、总结部分 常微分方程在科学和工程中常用于建立
20、数学模型，通常它们没有解析解，而需要数值方法来近似求解。在科学的计算机化进程中，常微分方程也在不断地精进中，更在不同的领域得 到了前所未有的发展和应用。求解常微分方程的初值问题的数值方法有 单步法和多步法。其中单步法中的 Runge-Kutta 方法是目前工程上求常微分方程数值解的基本方法。但是这并不表示 Runge-Kutta 方法是最好的，它也存在着缺点 1387 ，例如：它计算函数值 f 的次数较多。尤其是在级数高于 4 级时，这一缺点特别明显：计算函数值 f 的工作量增加较快，但是精度即阶数提高较慢。还有各阶 Runge-Kutta 方法的绝对稳定区域不够大，对于坏条件的初值问题不尽适
21、用等。但是它作为一类便于应用的单步方法，其实用价值还是很高的，特别是在现如今计算机发展飞快的时代。一定会有更先进的方法来使得 Runge-Kutta 方法得到完善，从而得到更好也更为实用的数值方法。 四、参考文献 1 美 Richard L.Burden， J.Douglas Faires 数值分析 M北京：高等教育出版社， 2005 2 黄明游 ， 刘播 ， 徐涛 数值计算方法 M 北京 ： 科学出版社 ， 2006 3 孙维夫，姜云义常微分 方程初值问题 RK 法和多步法 J烟台职业学院学报， 2010，（ 2）： 84-88 4 中国大百科全书数学北京：中国大百科全书出版社， 1988
22、5 王高雄，朱思铭，王寿松常微分方程 M北京：高等教育出版社， 2006 6 余徳浩，汤华中微分方程数值解法 M北京：科学出版社， 2003 7 胡健伟，汤怀民微分方程数值方法 M北京：科学出版社， 1999 8 李荣华，冯果忱微分方程数值解法 M北京：高等教育出版社， 1997 9 Arieh Iserles 微分方程数值分析基础教程 M北京：清华大学出版社， 2005 10 李忠杰常微分方程初值问题 RK 法和多步法 J黑龙江科技信息，科教文化，（ 1）： 199 11 李瑞遐，何志庆微分方程数值方法 M上海：华东理工大学出版社， 2005 12 刘利斌，王伟四阶 Runge-Kutta 算法的优化分析 J成都大学学报（自然科学版）， 2007，（ 1）：19-21 13 林群微分方程数值解法基础教程 M北京：科学出版社， 2003 10 马明书龙格 -库塔法的级数和阶数 J河南师范大学河南机专学报， 1994，（ 2）： 5-7 15 吴怡，张向辉，熊超西三阶 Runge-Kutta 最优算法 J河北秦皇岛职业技术学院基础部， 2005，（ 18）： 57-61
