1、 本科毕业论文（设计） （ 20 届） 常微分方程初值问题的 Runge-Kutta 解法 所在学院 专业班级 信息与计算科学 学生姓名 学号 指导教师 职称 完成日期 年 月 I 摘要： 常微分方程是指只有一个自变量的微分方程。 常微分方程在很多学科领域内有着重要的作用， 如 自 动控制、各种电子学装置的设计 等 。 本文首先介绍了常微分方程以及初值问题的除 Runge-Kutta 解法外的数值解一些理论；接着，我们对常微分方程初值问题提出 Runge-Kutta 解法的构造公式；然后，我们再对这个公式进行误差分析，包括相容性、稳定性和收敛性分析；最后，通过一些数值实例，我们验证了Runge
2、-Kutta法的优缺点。 关键词： 常微分方程； Runge-Kutta解法；误差分析 II The Runge-Kutta Method for Initial Value Problems of Ordinary Differential Equations Abstract: The ordinary differential equations have only one independent variable in the equation. And the ordinary differential equations play an important role in many
3、 areas, such as automatic control, electronics equipment design and so on. In this paper, we firstly introduce some theories of the initial value problems for ordinary differential equations, except the Runge-Kutta method; Next, we construct some numerical formulations of the Runge-Kutta Method to s
4、olve the initial value problems for ordinary differential equations; Then, we get the error analysis of these formulations, including consistency, stability and convergence analysis; Finally, some numerical examples are provided to confirm the advantages and disadvantages of the Runge-Kutta method.
5、Keywords: ordinary differential equations; Runge-Kutta method; error estimates III 目 录 1 绪 论 . 1 1.1 问题的背景 . 1 1.2 常微分方程初值问题简介 . 1 2 常微分方程初值问题的一些数值解法 . 3 2.1 单步法中的欧拉解法 . 3 2.2 线性多步法 . 3 3 RUNGE-KUTTA解法 . 4 3.1 RUNGE-KUTTA解法的构造 . 4 3.2 误差分析 . 8 3.2.1 相容性 . 8 3.2.2 收敛性 . 8 3.2.3 稳定性 . 10 4 数值算例 . 14 5
6、 结 论 . 16 致 谢 . 17 参考文献 . 18 附录 . 19 1 1 绪 论 1.1 问题的 背景 常微分方程在微积分概念出现后即已出现。从莱布尼兹专门研究用变量变换解决一阶微分方程的求解问题的“求通解”时代，到 1841 年刘维尔的里卡蒂方程不存在一般初等解和柯西的初值问题的提出，常微分方程转向“求定解”时代。再到 19 世纪末天体力学的太阳系稳定 性问题眼界需要常微分方程解的大范围性态，从而发展到了“求所有解问题”。最后到了 20 世纪六七十年代，常微分方程因为计算机的发展进入“求特殊解”阶段，发现了具有新性质的特殊的解和方程，如混沌（解）、奇异吸引子及孤立子等等 1。 常微分
7、方程是指只有一个自变量的微分方程。而且 常微分方程在很多学科领域内有着重要的作用， 如 自动控制、各种电子学装置的设计、弹道的计算、飞机和导弹飞行的稳定性的研究、化学反应过程稳定性的研究等等 。 常微分方程初值问题的数值求解是求近似解的一种经典方法。 1.2 常微分方程初值问题简介 2-9 大多数学科领域内的问题需要求解一个初值问题，即求解满足给定初值条件的微分方程。我们更多的是使用逼近原问题的解的方法来逼近原方程的解。因逼近方法给出更精确的结果和实际的误差信息，所以更经常被使用。 现代科学、技术、工程中的大量数学模型都可以用常微分方程的初值问题来描述，所以研究常微分方程初值问题的一些性质、解
8、法很有必要。 一些典型的常微分初值问题的数值求解问题的方法有：单步法和线性多步法。 在求解区间 ba, 上取定节点 bxxxxxa NN 1210 令 nnn xxx 1 ，称为积分网格的步长。用 Nyyy , 10 表示初值问题精确解 xy 在节点Nxxx , 10 上函数值 Nxyxyxyxy , 210 的近似值。对给定的数值积分方法， ny 中的各个 ny 市按某一个递推算法确定的。一个递推算法，如果在用它计算 1ny 时只用到已经求出的诸值 nyyy , 10 中的 ny ，而无须使用其余值 110 , nyyy 中的任何一个，则称此算法为单步方法。相反的则是多步法。 对于一阶的常微
9、分方程初值问题 yxfdxdy , ， Tx0 00 yy （ 1） 2 的数值方法，其中 f 是 x 和 y 的已知函数， 0y 为给定的初值。为使上面的解存在、唯一且连续依赖初值 0y ，则 f 必须关于 y 满足 Lipschitz 条件：即存在常数 L，使得对于任意 Tx ,0 ，均成立不等式 2121 , yyLyxfyxf 单步显式公式的一般形式为 hyxhyy nnnn ,1 ， 00 yy （ 2） 称 hyx , 为增量函数。一般来说微分方程问题（ 1）的精确解 nxy 不满足式（ 2），即一般地有 hxyxhxyxy nnnn ,1 定义 1 称 hxyxhxyxyR nn
10、nnn ,11 为单步显式公式（ 2）在 1nx 处的截断误差。 一个求解公式的局部截断误差刻画了其逼近微分方程的准确程度。根据上述定义 可以直接求得各式的局部截断误差。 定义 2 如果一个求解公式的局部截断误差为 1 ph ，则称该求解公式是 p 阶的，或称具有 p阶精度。 在多步法中也有几个基本的理论问题的定义，定理。 定义 3 对于初值问题（ 1）的多步法 kj jnjkj jnj fhy 00 ， ,1,0n （ 3） 说 是相容的，如果它至少是一阶的。 定理 1 为使 k 步法（ 3）相容，必需且只需 11,01 3 2 常微分方程初值问题的一些数值解法 下面将先简单的介绍一下 除
11、Runge-Kutta法外的一些数值方法。 2.1 单步法中的欧拉解法 4-9 单步法中最简单的数值解法是欧拉法。欧拉法的基本思想是由 ny 推导出下一个节点 1ny ，其推导公式为 nnnn yxhfyy ,1 ， ,1,0 . 其局部截断误差为 21 hxR nh 。欧拉法的几何意义是用一条过 00,yx 的折线来近似替代过 00,yx 的积分曲线。 虽然这个方法的计算很简单，但精确度较低，所以有了精确度较高的梯形方法。推导公式为 111 ,2 nnnnnn yxfyxfhyy ， ,1,0n ，它的局部截断误差为 32 hxR nh 。此为隐型格式，迭代求解所需计算量大。 结 合 上 述
12、 两 种 方 法 的 优 点 ， 从 而 得 到 改 进 欧 拉 方 法 。 其 计 算 公 式 为 nnnnnnnn yxhfyxfyxfhyy ,2 11 ， ,1,0n 。 2.2 线性多步法 810 常用的多步法主要有 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 ， km ,1,0 4 3 Run
13、ge-Kutta 解法 3.1 Runge-Kutta 解法的构造 5-15 用 Taylor 级数法的公式 hyxhxy nnnn ,1 ， ,1,0n 可以构造高阶单步法，但增量函数 hyx , 是用 f 的各阶导数表示的，通常不易计算。由于函数在一点的导数值可以用该点附近若干点的函数值近似表示，因此可以把 Taylor 级数法中的增量函数 改为 f 在一些点上函数值的组合，然后用 Taylor 展开确定待定的系数，使方法达到一定的阶。这就是 Runge-Kutta 法的基本思想。或者是在每个区间多取几点，将它们的斜率加权平均作为平均斜率，就有可能构造出精度更高的计算公式。通过构造高阶单步
14、法来提高精度，而较高的精度意味着计算结果更加精确，误差会随着 h 的减小而迅速减小。 现在假设 Ni iinn kwhyy 11（ 4） 其中 nn yxfk ,1 ， 1, iinini khyhxfk ， ni ,3,2 （ 5） 此类递推算法为显示 N 级 Runge-Kutta 方法。上述算法公式中的系数 iw ， i 和 i 希望如此确定，使得（ 4）式的 Taylor 展开式中所有 h 幂次不超过 p 的那些项的系数与 Taylor 展开式 nhnppnnnn xRxyphxyhxyhxyxy !2 21 （ 6） （其中 11!1 pph yphR， 1 nn xx ）中的相应项
15、的系数相等。 当 N=1 时，则可知 1w =1，此时式（ 4）为欧拉公式。 当 N=2 时，确定系数 iw ， i 和 i 使得 nn xyhxy 与 Ni iikwh 1两者的 Taylor 展开中含 rh 项的系数， pr ,2,1 ，对尽可能大的 p 相等。 利用二元的 Taylor 展开公式 nnnnnn yxfykxhyxfkhyhxf , 5 nn yxfykxh ,!222 可将（ 4）式表示为 32222211 hfffwhfwwhyy nnynxnnn （ 7） 再将 Taylor 展式（ 4）改成 321 !2 hfffhhfxyxy nnynxnnn （ 8） 比较（
16、7）和（ 8）式中 2,1rhr 项的系数，得到 121 ww ， 2122 w， 2122 w（ 9） 因为有 4 个待定系数，可是就只有 3 个方程数目。因此具有无穷多个数组 2221 , ww 满足条件（ 9），即存在无穷多个二级 Runge-Kutta 方法，它们的截断误差 3hxR nh 。 若取 212w，便得到 211 2 kkhyy nn nn yxfk ,1 12 , hkyhxfk nn 或 nnnnnnnn yxhfyxfyxfhyy ,2 11 这是改进欧拉公式。 因为二级的 RK 方法最多是二阶的，若想得到更高阶的方法，必须增加级数，即取较大的 N。当 N 较大时，计
17、算也就会更加的复杂。 当 N=3 时，和 N=2 的推导过程相相似，当 43,0,41321 www时，得到 311 34 kkhyy nn nn yxfk ,1 12 31,31 hkyhxfk nn6 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 hkyhxfk nn 23
18、21,21 hkyhxfk nn 34 , hkyhxfk nn （ 10） 其局部截断误差为 5h 。 但是，以上的方法还是不够的。在涉及初值问题数值解的实际问题中，总要先设定一个精度允许值，数值解与真解的偏离不能超过这个允许值。因此，步长还受到这个允许值的控制，常常会发生这样的情形： 仅在求解的某一部分需要小的步长，而在其余部分则用较大的步长就可以了。 所以，各种用于算法中自动地调整步长的方法（称为自适应方法）就得到了发展。就如在四阶的 Runge-Kutta 方法中，为了从 x 到 hx ，我们可以用步长为 h 做一步，也可以用步长为 2h 做两步，如果没有截断误差，它们求出的 hxy 的值应是相同的 。但事实上它们不会相同，所以其差异可看成是对局部截断误差的一个估计：如果这个差异不超过预先给定的允许值，那么步长 h 就是满意的；如果它超过允许值，就将 h 减半；如果它远小于允许值，就将步长加倍。 虽然这样的过程很容易计算，但要浪费大量的时间，因而不宜提倡。
