1、毕业论文 开题报告 信息与计算科学 常微分方程初值问题的 Runge-Kutta 解法 一、选题的背景、意义 常微分方程是指只有一个自变量的微分方程。而且 常微分方程在很多学科领域内有着重要的作用， 如 自动控制、各种电子学装置的设计、弹道的计算、飞机和导弹飞行的稳定性的研究、化学反应过程稳定性的研究等等，这些问题都可以化为求微分方程的解，或者化为研究解的性质的问题。 这些问题都包含某个变量关于另一个变量的变化。大多数这样的问题需要求解一个初值问题，即求解满足给定初值条件的微分方程 1 。 现代科学、技术、工程中的大量数学模型都可以用常微分方程的初值问题来描述，所以研究常微分方程初值问题的一些
2、性质、解法很有必要。 1.1 历史背景 常微分方程在微积分概念出现后即已出现。从莱布尼兹专门研究用变量变换解决一阶微分方程的求解问题的“求通解”时代，到 1841 年刘维尔的里卡蒂方程不存在一般初等解和柯西的初值问题的提出，常微分方程转向“求定解”时代。再到 19 世纪末天体力学的太阳系稳定性问题眼界需要常微分方程解的大范围性态，从而发展到了“求所有解问题”。最后到了 20 世纪六七十年代，常微分方程因为计算 机的发展进入“求特殊解”阶段，发现了具有新性质的特殊的解和方程，如混沌（解）、奇异吸引子及孤立子等等 2 。 常微分方程初值问题的数值求解是求近似解的一种经典方法。 其中常微分初值问题的
3、数值求解问题的方法又包括单步法和线性多步法 3 。单步法主要有欧拉法 84 和 Runge-Kutta法 105 。 Runge-Kutta 法是最重要也最便于应用的单步法。 1.2 国内外研究现状 常微分方程在科学和工程中常用于建立数学模型，通常它们没有解析解，而需要数值方法来近似求解。在科学的计算机化进程中，常微分方程也在不断地精进中，更在不同的领域得到了前所未有的发展和应用。求解常微分方程的初值问题的数值方法有 单步法和多步法。其中单步法中的 Runge-Kutta 方法是目前工程上求常微分方程数值解的基本方法。 现在的 Runge-Kutta 方法 运用非常广泛，国内外都对其 做的 研
4、究 也 非常多，以 前 Runge- Kutta 方法 仅仅是运用在一些简单的方程和实际问题中，现在科学家都希望将其能够 运用在更广泛和更复杂的情况下，且能构造 更优化的 Runge-Kutta 方程格式来 满足各种不同目的需要。 1.3 发展趋势 Runge-Kutta 方法是目前工程上求常微分方程数值解的基本方法，它的基本思想是在每个区间多取几点，将它们的斜率加权平均作为平均斜率，就有可能构造出精度更高的计算公式。但是这并不表示 Runge-Kutta 方法是最好的，它也存在着缺点 95 ，例如：它计算函数值 f 的次数较多。尤其是在级数高于 4 级时，这一缺点特别 明显：计算函数值 f
5、的工作量增加较快，但是精度即阶数提高较慢。还有各阶 Runge-Kutta 方法的绝对稳定区域不够大，对于坏条件的初值问题不尽适用等。但是它作为一类便于应用的单步方法，其实用价值还是很高的，特别是在现如今计算机发展飞快的时代。一定会有更先进的方法来使得 Runge-Kutta方法得到优化，从而得到更好也更为实用的数值方法。 1.4 计算前提 对于一阶的常微分方程初值问题 yxfdxdy , ， Tx0 00 yy （ 1） 的数值方法，其中 f 是 x 和 y 的已知函数， 0y 为给定的初值。为使上面的解存在、唯一且连续依赖初 值 0y ，则 f 必须关于 y 满足 Lipschitz 条件
6、：即存在常数 L，使得对于任意 Tx ,0 ，均成立不等式 2121 , yyLyxfyxf 单步显式公式的一般形式为 hyxhyy nnnn ,1 ， 00 yy （ 2） 称 hyx , 为增量函数。一般来说微分方程问题（ 1）的精确解 nxy 不满足式（ 2），即一般地有 hxyxhxyxy nnnn ,1 定义 1 称 hxyxhxyxyR nnnnn ,11 为单步显式公式（ 2）在 1nx 处的截断误差。 一个求解公式的局部截断误差刻画了其逼近微分方程的准确程度。根据上述定义可以直接求的各式的局部截断误差。例如欧拉公式的局部截断误差为 11 ,n n n n nR y x y x
7、hf x y x nnnnn xyhxyyhxyhxy 221 nyh 2211 nnn xx 定义 2 如果一个求解公式的局部截断误差为 1 ph ，则称该求解公式是 p 阶的，或称具有 p 阶精度。 二、研究的基本内容与拟解决的主要问题 2.1 基本内容 （ 1）欧拉法 84 单步法中最简单的数值解法是欧拉法。欧拉法的基本思想是由 ny 推导出下一个节点1ny ，其推导公式为 nnnn yxhfyy ,1 ， ,1,0 . 其局部截断误差为 21 hxR nh 。欧拉法的几何意义是用一条过 00,yx 的折线来近似替代过 00,yx 的积分曲线。 虽然这个方法的计算很简单，但精确度较低，所
8、以有了精确度较高的梯形方法。推导公式为 111 ,2 nnnnnn yxfyxfhyy， ,1,0n ，它的局部截断误差为 32 hxR nh 。此为隐型格式，迭代求解所 需计算量大。 结合上述两种方法的优点，从而得到改进欧拉方法。其计算公式为 nnnnnnnn yxhfyxfyxfhyy ,2 11 ， ,1,0n 。 （ 2）线性多步法 常用的多步法主要有 Adams 法和 Milne 法，其中 Adams 法又包括显式 Adams 法和隐式Adams 法。显式 Adams 法的线性 1k 步公式为 kj jnkjnn fhyy 01 其中 kjm mjkj jm 1 ， kj ,2,1,
9、0 。 且隐式 Adams 公式为 km nmmnn fhyy0 1*1 其中 01* 1 dsmsmm ， km ,1,0 2.2 主要内容 （ 1） Runge-Kutta 法的构造 105 用 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 xR
11、xyphxyhxyhxyxy !2 21 （ 5）（其中 11!1 pph yphR， 1 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） 再
12、将 Taylor 展式（ 3）改成 321 !2 hfffhhfxyxy nnynxnnn （ 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 这是改进欧拉公式
13、。 因为二级的 RK 方法最多是二阶的，若想得到更高阶的方法，必须增加级数，即取较大的 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 方法，其计算公式为 43
14、211 2261 kkkkhyy nn nn yxfk ,1 12 21,21 hkyhxfk nn 23 21,21 hkyhxfk nn 34 , hkyhxfk nn （ 9） 其局部截断误差为 5h 。 由于 RK 方法是通过对函数的 Taylor 展开来构造的，因此高阶 RK 方法要求微分方程的解具有高阶导数。如果解的光滑性较差，即使用高阶方法也不可能得到高精度的数值解，此时 应当用低阶方法。或者说从这些上述公式的讨论中可知，级数 n 表示计算函数值 f 的次数，当 4,3,2n 时， n 级的 RK 方法的阶数 nnp ，即计算函数值 f 的次数越多， RK 方法的精度即阶数越高。
15、但是当 4n 时，情况就大不一样了， J.C.Butcher 指出：当 7,6,5n时， 1nnp ；当 9,8n 时， 2nnp ；而且 ， 810p 。由此可见四级以上的RK 方法，计算函数值 f 的工作量增加较快，而阶数即精度提高较慢，故而，在实际问题中比较常用的是四级 RK 方法。 （ 2） Runge-Kutta 法的误差分析 98 I、相容性 当 phhxyL ; ， 2 时，称数值方法是相容的。 由于 hxyxh xyhxyhxyL ,; 2hhxyxyhxy hyxfwkwhxyx ni ini ii 11 , 于是 hyxfwxyhxyL ni i 1 ,; hyxfwxy
16、ni i 1 , hyxfwxy ni i ,1 1是微分方程的解当 hwni i 11故有结论： R-K 法相容 11 ni iw。即 p 阶的 Runge-Kutta 方法是相容的。 II、收敛性 定义 3 如果一个单步法满足：对 Tx,0 中任意一个固定的点 a 成立 ayynanhx h 00lim则称该方法是收敛的。 由于 anhx 0 ，当 0h 时， n ，也即当步长 h 减小时 ，从 0x 计算到点 a 的步数就增大，故有下面的定理。 定理 1 若公式（ 2）是 p 阶方法，且增量函数 hyx , 关于 y 满足 Lipschitz 条件，则该方法收敛。 对于 Runge-Ku
17、tta 方法，当函数 yxf , 关于 y 满足 Lipschitz 条件，可用数学归纳法证明所有的 ik 关于 y 满足 Lipschitz 条件，从而 增量函数 hyx , 关于 y 满足 Lipschitz 条件，所以 Runge-Kutta 方法收敛。而且对于这类满足 Lipschitz 条件的单步法而言，相容条件是收敛条件的充分必要条件。 III、稳定性 定义 4 如果一个单步法满足：存在正常数 c 和 0h ，使得用任意步长 0,0hh 以及任意两个初值 0y 和 0g 计算得到的两组数值 nyyy , 21 和 nggg , 21 成立 001m a x gycgy iini ，
18、 则称该方法是稳定的，其中 hTn。 定理 2 若 增量函数 hyx , 关于 y 满足 Lipschitz 条件，则公式（ 2）是稳定的。 由上可知，欧拉法，改进欧拉法和 Runge-Kutta 方法都是稳定的。 IV、绝对稳定性 对于模型问题： ydxdy ， Tx0 00 yy （ 10） 其中 是常数（可以是复数）。 定义 5 对初值问题（ 10），记 hh ， 若一个数值方法在某一步由计算产生的误差在以后各步中逐步减少，则称该方法关于 h 是绝对稳定的。复平面上所有这样的 h 组成的区域称为该方法的绝对稳定区域。 由此可知，当 1!0 niiih满足时， n 级 n 阶的 Runge
19、-Kutta 方法 (n=1,2,3,4)是绝对稳定的。 三、研究的方法与技术路线、研究难点，预期达到的目标 3.1 研究的方法和路线 本文中的讨论了常微分方程初值问题的数值方法的一些基本解法，具体描述构造Runge-Kutta 方法的过程，给出了这个 Runge-Kutta 方法的相容性分析，收敛性分析和稳定性分析。从 Runge-Kutta 方法的构造可以看出它的优势，阶数和级数的关系，及其收敛趋势。从对 Runge-Kutta 方法的误差分析，可以得到它的优缺点。 3.2 研究的难点 常微分方程的初值和边值是不同的，这里只考虑了初值问题的数值解法。论文的难点在于一些计算公式的计算量小但是
20、精确度低或者精确度高但是计算量增大很多。其次是常微分方程的初值问题的数值求解，还有其他的解法，如线性多步法等， Runge-Kutta 方法也有它的一些不足之处。 另外，怎样更好地对 Runge-Kutta 方法进行优化 ，使其在不过分增加计算工作量的前提 下使得计算简便和精确度得到提高，且计算的各点在总体上协调也是个难点。 3.3 预期目标 通过对 Runge-Kutta 方法的研究，得到它的相容性、 收敛性 和稳定分析， 通过程序（ C语言、 FORTRAN 语言或 MATLAB 语言）编程，给出数值例子，验证理论分析。求 Runge-Kutta方法的数值解在于对 Runge-Kutta
21、方法的计算中，阶数和级数已知时，求它的构造公式。再由这构造公式推导出确切的数值解。 四、论文详细工作进度和安排 第 7 学期第 9 周（ 2010 年 11 月 5 号）至第 7 学期第 19 周（ 2011 年 1 月 10 号） 完成毕业论文文献检索、文献综述、外文文献翻译及开题报告。 第 7 学期第 19 周（ 2011 年 1 月 10 号）至第 8 学期第 3 周（ 2011 年 3 月 11 号） 完成毕业论文的数据收集、论文初稿。 第 8 学期第 3 周（ 2011 年 3 月 11 号）至第 8 学期第 11 周（ 2011 年 5 月 3 号） 1、进入实习单位进行毕业实习，
22、对论文进行修改； 2、第 11 周（ 2011 年 5 月 3 日）前必须返校，完成毕业实习返校，并递交毕业实习报告，进一步完善毕业论文； 第 8 学期第 14 周（ 2011 年 5 月 23 号 2011 年 5 月 28 号）完成第一轮毕业论文答辩； 第 8 学期第 15 周（ 2011 年 5 月 28 号 2011 年 6 月 3 号） 第一轮毕业论文答辩未通过的学生完成第二轮毕业论文答辩，并随机抽取部分完成较好地毕业论文进行校级答辩 五、主要参考文献： 1 美 Richard L.Burden， J.Douglas Faires数值分析 M北京：高等教育出版社， 2005 2 王高雄，朱思铭，王寿松常微分方程 M北京：高等教育出版社， 2006 3 孙维夫，姜云义常微分方程初值问题 RK 法和多步法 J烟台职业学院学报， 2010，（ 2）： 84-88 4 余徳浩，汤华中微分方程数值解法 M 北京 ：科学出版社， 2003
