第九章 常微分方程的数值解法 主 要 内 容 §1、引言 §2、初值问题的数值解法--单步法 §3、龙格-库塔方法 §4、收敛性与稳定性 §5、初值问题的数值解法―多步法 §6、方程组和刚性方程 §7、习题和总结

§1、 引 言 主 要 内 容 研究的问题 数值解法的意义

找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系 1.什么是微分方程 ? 内部联系非常复杂 现实世界中大多数事物 其状态随着 时间、地点、条件 的不同而不同 找出其状态和状态变化规律之间的相互联系，也即一个或一些函数与他们的导数之间的关系 微分方程 此种关系的数学表达就为

2.数值求解微分方程的意义 如何建立数学模型已在建模课程中得到讨论，各类微分方程本身和他们的解所具有的特性 已在常微分方程及数学物理方程中得以解释， 本章专门讨论 如何利用数值方法求解微分方程（组）的问题。

3.什么是微分方程 (组)的解析解? 3.什么是微分方程(组)的解析解? 一个或一组具有所要求阶连续导数的解析函数，将 它代入微分方程（组），恰使其所有条件都得到满 足的解称为解析解(或古典解),称为真解或解。 寻找解析解的过程称为求解微分方程组。

4.什么是微分方程的数值解? 虽然求解微分方程有许多解析方法,但解析方法 只能够求解一些特殊类型的方程,从实际意义 上来讲 我们更关心的是某些 特定的自变量在某一个 定义范围内的一系列离散点上的近似值. 把这样一组近似解称为 微分方程在该范围内的 数值解 寻找数值解的过程称为数值求解微分方程。

在大量的实际方程中出现的函数起码的连续性都 无法保证，更何况要求阶的导数 很多微分方程 根本求不到 问题的解析解！ 求解数值解 重要手段。

5.常微分方程数值解法的特点 常微分方程的数值解法常用来求近似解 数值解法得到的近似 解(含误差)是一个 离散的函数表. 根据提供的算法 通过计算机 便捷地实现

6.基本知识 本章主要讨论一阶常微方程的初值问题 各种数值解法 其中f (x,y)是已知函数，(1.2)是定解条件也称为 初值条件。

常微分方程的理论指出： 当 f (x,y) 定义在区域 G=(a≤x≤b,｜y｜＜∞) (Lipschitz)条件 若存在正的常数 L 使： 则称 f (x,y) 对y 满足李普希兹条件，L 称为Lipschitz常数. 就可保证方程解的存在唯一性

若 f (x,y) 在区域 G连续，关于y 一阶常微分方程的初值问题的解存在且唯一. 我们以下的讨论，都在满足上述条件下进行. 满足李普希兹 条件 若 f (x,y) 在区域 G连续，关于y 一阶常微分方程的初值问题的解存在且唯一. 我们以下的讨论，都在满足上述条件下进行. 一阶常微分方程组常表述为： 方程组 初值条件

写成向量形式为 高阶常微分方程定解问题如二阶定解问题：

这些解法都可以写成向量形式 用于一阶常微分方程组的初值问题. 也就解决了高阶方程的定解问题.

简单的数值方法与基本概念 §2、初值问题的数值解法―单步法 1. 简单欧拉法（Euler） 2．后退的欧拉法 3．梯形法

弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路. 1. 简单的欧拉(Euler)方法 考虑模型： 欧拉方法 最简单而直观 实用方法 弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路. 在精度要求不高时 通过欧拉方法的讨论

2. 欧拉方法的导出 把区间［a,b］ 分为n个小区间 N等分 步长为 节点 要计算出解函数 y(x) 在一系列节点 处的近似值

对微分方程(1.1)两端从 进行积分

右端积分用左矩形数值求积公式 得

/* Euler’s polygonal arc method*/ x0 x1 或用向前差 商近似导数 依上述公式逐次计算可得： 亦称为欧拉折线法 /* Euler’s polygonal arc method*/ 例题 每步计算 只用到

故也称Euler为单步法。 公式右端只含有已知项 所以又称为显格式的单步法。 3.欧拉公式有明显的几何意义 依此类推得到一折线

欧拉方法 就是用这条折线近似地代替曲线 也称欧拉折线法.

4.欧拉法的局部截断误差： 从上述几何意义上得知，由Euler法所得的 折线明显偏离了积分曲线，可见此方法 非常粗糙。 在假设第 i 步计算是精确的前提下，考虑 截断误差 称为局部截断误差 /* local truncation error */。

定义 　　　在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑的截断误差 Ri = y(xi+1)  yi+1 称为局部截断误差 /* local truncation error */。 定义 　　　若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。  欧拉法的局部截断误差： Ri 的主项 /* leading term */ 欧拉法具有 1 阶精度。

Taylor展开式,一元函数的Taylor展开式为: 如果单步差分公式的局部截断误差为O(hp+1), 则称该公式为p阶方法.这里p为非负整数. 显然,阶数越高,方法的精度越高. Taylor展开式,一元函数的Taylor展开式为: 若某算法的局部截断误差为O(hp+1)，则称该算法有p 阶精度。 Ri 的主项 /* leading term */

 隐式欧拉法 /* implicit Euler method */ 5. 欧拉公式的改进： x0 x1  隐式欧拉法 /* implicit Euler method */ 向后差商近似导数 )) ( , ) 1 x y f h +  ) 1 , ... ( - = + n i y x f h 由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。  隐式欧拉法的局部截断误差： 即隐式欧拉公式具有 1 阶精度。

6.梯形公式 /* trapezoid formula */ — 显、隐式两种算法的平均 注：的确有局部截断误差 ， 即梯形公式具有2 阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。 需要2个初值 y0和 y1来启动递推 过程，这样的算法称为双步法 /* double-step method */，而前面的三种算法都是单步法 /* single-step method */。  中点欧拉公式 /* midpoint formula */ x0 x2 x1 中心差商近似导数 假设 ，则可以导出 即中点公式具有 2 阶精度。

  方 法 显式欧拉 隐式欧拉 梯形公式 中点公式 简单 精度低 稳定性最好 精度低, 计算量大 精度提高 计算量大 精度提高, 显式 方 法   显式欧拉 隐式欧拉 梯形公式 中点公式 简单 精度低 稳定性最好 精度低, 计算量大 精度提高 计算量大 精度提高, 显式 多一个初值, 可能影响精度

 改进欧拉法 /* modified Euler’s method */ Step 1: 先用显式欧拉公式作预测，算出 ) , ( 1 i y x f h + = Step 2: 再将 代入隐式梯形公式的右边作校正，得到 1 + i y )] , ( ) [ 2 = x f h 注：此法亦称为预测-校正法 /* predictor-corrector method */。可以证明该算法具有 2 阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。

§3 龙格 - 库塔法 /* Runge-Kutta Method */ 建立高精度的单步递推格式。 单步递推法的基本思想是从 ( xi , yi ) 点出发，以某一斜 率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及其各种变形所 能达到的最高精度为2阶。  考察改进的欧拉法，可以将其改写为： 斜率 一定取K1 K2 的平均值吗？ 步长一定是一个h 吗？

首先希望能确定系数 1、2、p，使得到的算法格式有2阶精度，即在 的前提假设下，使得 将改进欧拉法推广为： ) , ( ] [ 1 2 phK y ph x f K h i + = l 首先希望能确定系数 1、2、p，使得到的算法格式有2阶精度，即在 的前提假设下，使得 Step 1: 将 K2 在 ( xi , yi ) 点作 Taylor 展开 Step 2: 将 K2 代入第1式，得到

存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。 Step 3: 将 yi+1 与 y( xi+1 ) 在 xi 点的泰勒展开作比较 要求 ，则必须有： 这里有 个未知数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库塔格式。 注意到， 就是改进的欧拉法。 Q: 为获得更高的精度，应该如何进一步推广？

 最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ： ) ... , ( ] [ 1 2 32 31 3 21 - + = m m m i hK y h x f K b a l 其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i1 ) 均为待定系数，确定这些系数的步骤与前面相似。  最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ：

 龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的值。Butcher 于1965年给出了计算量与可达到的最高精度阶数的关系： 注：  龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的值。Butcher 于1965年给出了计算量与可达到的最高精度阶数的关系： 7 5 3 可达到的最高精度 6 4 2 每步须算Ki 的个数  由于龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用低阶算法而将步长h 取小。 深入研究龙格- 库塔法请看此处！

 §4 收敛性与稳定性 /* Convergency and Stability */ 定义 　　　 若某算法对于任意固定的 x = xi = x0 + i h，当 h0 ( 同时 i  ) 时有 yi  y( xi )，则称该算法是收敛的。 例：就初值问题 考察欧拉显式格式的收敛性。 解：该问题的精确解为 欧拉公式为 对任意固定的 x = xi = i h ，有  

What is wrong ??!  2.稳定性 /* Stability */ 例：考察初值问题 在区间[0, 0.5]上的解。 例：考察初值问题 在区间[0, 0.5]上的解。 分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 0.0 0.1 0.2 0.3 0.4 0.5 精确解 改进欧拉法 欧拉隐式 欧拉显式 节点 xi 1.0000 2.0000 4.0000 8.0000 1.6000101 3.2000101 1.0000 2.5000101 6.2500102 1.5625102 3.9063103 9.7656104 1.0000 2.5000 6.2500 1.5626101 3.9063101 9.7656101 1.0000 4.9787102 2.4788103 1.2341104 6.1442106 3.0590107 What is wrong ??!

定义 　　　若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的 /*absolutely stable */。 常数，可以是复数 一般分析时为简单起见，只考虑试验方程 /* test equation */ 当步长取为 h 时，将某算法应用于上式，并假设只在初值产生误差 ，则若此误差以后逐步衰减，就称该算法相对于 绝对稳定， 的全体构成绝对稳定区域。我们称算法A 比算法B 稳定，就是指 A 的绝对稳定区域比 B 的大。 h l h =

由此可见，要保证初始误差0 以后逐步衰减， 必须满足： 例：考察显式欧拉法 - 1 2 Re Img 由此可见，要保证初始误差0 以后逐步衰减， 必须满足： 例：考察隐式欧拉法 2 1 Re Img 可见绝对稳定区域为： 注：一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。

例：隐式龙格-库塔法 其中2阶方法 的绝对稳定区域为 而显式 1~ 4 阶方法的绝对稳定区域为 无条件稳定 Img Re Img Re - 其中2阶方法 的绝对稳定区域为 Re Img 而显式 1~ 4 阶方法的绝对稳定区域为 k = 1 2 3 4 - Re Img 无条件稳定

例1 用欧拉方法，隐式欧拉方法和欧拉中点公式计算 的近似解，取步长h=0.1，并与精确值比较 解:欧拉方法的算式为： 欧拉隐式方法在本题可解出方程，不必迭代，公式为： 欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式

本题精确解为y=e-x,计算结果见表9-1 例2 用欧拉公式和梯形公式的预估校正法计算： 的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为: 解: 本例中欧拉公式为：

梯形公式只校正一次的格式为 结果列入表9.2

龙格－库塔法深入研究 1. Runge-Kutta 法的一般形式 2. 2阶Runge-Kutta 方法 4．R-K-Fehhlberg 方法 5. 隐式R-K方法 6. 变步长方法

1.Runge-Kutta 法的一般形式 Runge-Kutta 法是用区间上若干个点上的导数的线性组合得到平均斜率，以提高方法的阶。其一般形式为 ：

2. 2级2阶Runge-Kutta方法 (9.11) 式(9.11) 称L级p阶Runge-Kutta方法(简称R-K法)。 其中 它的局部截断误差是 (9.11) 式(9.11) 称L级p阶Runge-Kutta方法(简称R-K法)。 当L＝1就是欧拉法，当L＝2时为改进的欧拉法。 2. 2级2阶Runge-Kutta方法 令 L=2，则

其局部截断误差是： 这是3个未知量的两个方程，其中有一个自由参数，方程组有无穷多解，从而得到一族2级2阶R-K方法 。

3. 经典Runge-Kutta方法 我们可以构造出一族3级3阶，一族4级4阶和一族5级4 阶等R-K方法。最常用的4级4阶是如下的经典R-K方法：

4．R-K-Fehhlberg 方法 R-K-Fehhlberg方法是在R-K方法的基础上引进误差和步长 控制的办法。即利用5阶R-K法

注：R-K-Fehhlberg比4阶R-K方法具有更大的优越性，他是计算稳定，高精度的方法，他的显著优点是，每一步仅需计算f的6个值；若用4阶R-K-L与5阶R-K-L在一起使用，则每步需要计算f的10个值。推荐使用！

5. 隐式R-K方法 类似于显式R-K公式，稍加改变，就得到隐式R-K方法。 它与显式R-K公式的区别在于：显式公式中对系数求和的上限是i-1，从而构成的矩阵是一个严格下三角阵。而在隐式公式中对系数求和的上限是L，从而构成的矩阵是方阵，需要用迭代法求出Ki，。推导隐式公式的思路和方法与显式R-K法类似。通常，同级的隐式公式获得比显式公式更高的阶。

通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式 R-K法有: 1级2阶中点公式 : 2级2阶梯形公式： 2级4阶R-K公式：

6.变步长方法 在单步法中每一积分步步长实际上是相互独立的，步长的 选择具有了灵活性。根据合理地选择每一积分步的步长， 既保证精度的要求，又可以减少计算量，从而减小舍入误 差。其方便的控制手段是基于误差的事后估计式。 对于给定的精度 ,如果 > ,反复将步长折半进行计算，直至 < 为止,这时取最终得到的作为结果； 如果 < ，将反复将步长加倍，直到 > 为止，这时再将步长折半一次，就得到所要的结果。 这种通过加倍或折半处理步长的计算方法称为变步长方法。 注：推荐使用精度好计算量低的变步长方法。 用四阶显式R-K方法做变步长方法是实践中较好的方法！

例 分别用改进的欧拉格式和四阶龙格—库塔格式解初值问题（取步长h=0.2）：

表 节点 改进欧拉法 四阶龙格—库塔法 准确解 0 1 1 1 0.2 1.186667 1.183229 1.183216 0.4 1.348312 1.341667 1.341641 0.6 1.493704 1.483281 1.483240 0.8 1.627861 1.612514 1.612452 1 1.754205 1.732142 1.732051 （注：已指出过准确解 ）

单步法的相容性、收敛性和稳定性 1.单步法的相容性 2.单步法的收敛性 3.单步法的稳定性 4.相容性和收敛性的关系 5.相容性和方法阶的关系 6.稳定性和收敛性 7.绝对稳定性和绝对稳定域

单步法的相容性 定义一：对于（9.1.1）常微分方程初值问题 单步法的形式可以变表示为 （9.2.19） 其中 h为步长 若对求解区间中任一固定的x，当 时皆成立 则称由（9.2.19）确定的单步法与微分方程初值问题是相容的 注意到上式左边极限为由（9.1.1）知它应等于从而由相容性 定义得 称相容性条件。

单步法的收敛性 定义二： 设 y(x) 是（9.1.1）的解， 是单步法（9.2.19）产生的数值解，对于每一个固定的 , ， 当 即 。若成立 ， 则称该方法是收敛的。

单步法的稳定性 定义三： 若一个数值方法在基点 处的值有 的扰动，在 此后各基点 （m>n）处的值产生的偏差均不超 若一个数值方法在基点 处的值有 的扰动，在 此后各基点 （m>n）处的值产生的偏差均不超 过 ，则称该方法是稳定的。 单步法的稳定性有以下定理 定理二： 若单步法的增量函数对变量y满足 Lipschtiz 条件， 则单步方法是稳定的。

相容性和收敛性的关系 定理一： 若单步法的增量函数对变量y满足Lipschtiz 条件， 即存在与 h , x 无关的常数 L， 对区域D= 任意两点 （x，y1）,(x,y2)成立，则单步法收敛的充分必要 条件是相容性条件成立。（读者自证）

相容性和方法阶的关系 若单步法是p阶方法则成立 若单步法满足相容性条件，得 所以 =0 也就是说单步法的阶数一定要是正数。由于我们考虑 所以 =0 也就是说单步法的阶数一定要是正数。由于我们考虑 的单步法皆为正整数，p至少为一。因此我们考虑的 单步法都满足相容性条件。

稳定性和收敛性的关系 若单步法的增量函数满足定理二的条件即单 步法是稳定的则单步法收敛的充分必要条件 是 相容性条件成立。

绝对稳定性和绝对稳定域 稳定性问题是一个比较复杂的问题。为了简化讨论一 般仅对试验方程 进行考察。这里假定 般仅对试验方程 进行考察。这里假定 Re<0,即试验方程本身是稳定的。 定义三： 若求解微分方程的数值方法对试验方程和给定 的步长h，在计算时引入舍入误差 ，这个 误差在计算后继的 ，k＝1,2,…所引起的 误差按绝对值均不增加，就称该方法是对这个 步长h和复数是绝对稳定的。保证数值方法绝 对稳定性的h>0和的允许范围，称为该方法的 绝对稳定域。

绝对稳定性和绝对稳定域2 将Euler方法应用到试验方程得 误差方程是 要求误差不增长则

则Euler 方法的绝对稳定域是以 为半径的圆的内部。为了保证稳定性步长有所控制。 假如当 时h应满足 ，当 时， h 应满 足 等等。 同样我们可以将试验方程用到其它各种单步法当中，求出其绝对稳定域。在实际应用中必须在这个范围内，否则误差传播相当大，即不稳定。

§5 初值问题的数值解法 ―多步法 1.Adams方法 2.米尔尼方法、汉明方法及辛普森方法 3.预测校正方法 §5 初值问题的数值解法 ―多步法 1.Adams方法 2.米尔尼方法、汉明方法及辛普森方法 3.预测校正方法 4.多步法的相容性、稳定性和收敛性

1.Adams方法 考虑型如 的k步法， 称为阿当姆斯（Adams）方法 拿其中一种为例，推导其公式。对上面所说的法一般形 式若取 ，有方程组 同样解之，得到3步4阶隐式Adams公式和它的余项。

其它请读者自证。我们仅将结论列于下表。 Adams显式公式 k p 公式 1 2 3 4

注：在这些方法中，4阶的Adams预测校正方法具有方法 简洁、使用方便、易排序、高精度等优点。尤其当函数 f比较复杂时更能显出它的优越性。 k p 公式 1 2 3 4 5 注：在这些方法中，4阶的Adams预测校正方法具有方法 简洁、使用方便、易排序、高精度等优点。尤其当函数 f比较复杂时更能显出它的优越性。

2.米尔尼方法、汉明方法及 辛普森方法 同理得到5个待定参数方程组。解之得 ， ， ， ， 。 同理得到5个待定参数方程组。解之得 ， ， ， ， 。 构成著名的Miline 4步4阶显式公式和它的余项。 ，

米尔尼方法、汉明方法及 辛普森方法 同理得到5个参数方程组。求解后就构成著名的3步4阶隐式Hamming公式和它的余项。 若取，也得到5个参数方程组。求解后就构成Simpson隐式公式和其余项。

3. 预测校正方法 不论单步法或多步法，隐式公式比显式公式稳定性好，但在实际使用隐式公式时，都会遇到两个问题：一个是隐式公式如何能方便地进行计算；另一个是实际计算步长取多大。如隐式梯形公式，每往前推进一步，不必进行多次迭代，而是采用一阶显式Euler公式预测，二阶隐式梯形公式校正一次，构成显式改进Euler公式，能达到与梯形公式同阶的精度，即二阶精度。

对于线性多步公式，一般地，不难验证，如果 预测公式是阶或阶精度，校正公式为阶精度， 则用预测公式提供初值，校正公式迭代一次的 效果也能达到阶精度，再迭代下去，效果就不 明显了。预测-校正技术即保证了计算精度，又 使隐式计算显式化，克服了隐式公式要反复迭 代的困难。至于实际计算步长的选取，我们对 预测-校正公式使用外推原理，得到误差估计式 ，用来调整计算步长，使达到控制误差要求。 对于线性多步方法常用的预测-校正公式有Miline- Hamming方法和4阶Adams显隐式预测-校正公式

修正Miline-Hamming公式 将Miline 公式和Hamming公式结合，构成预测-校正公式

利用外推原理，即加上后消去局部截断误差的主项。使 说明经过外推后的算法其精度提高一阶。忽略误差项，上式可改写为

由于 和是 在计算过程中获得的数据，称为Miline 公式和Hamming公式的事后误差估计式。我们可以用它们来调节计算步长的大小，即选择一个合适的的步长，使 其中的是要求达到的计算精度。

又可得到 Miline公式和Hamming公式的修正公式，它们分别是 从而构成如下的修正Hamming预测-校正公式，简称修正Hamming公式：

预测： 修正： 校正： 在应用这套公式时，先由同阶单步法提供初值 ， ， ， 。计算 时可取。

Adams预测-校正公式 将4步4阶显式Adams公式作为预测公式和3步4阶式Adams公式作为校正公式，构成Adams预测-校正公式。 它们的局部截断误差分别是

利用外推原理，将上两式作线性组合消去局部截断误差的主项。使计算精度至少提高一阶，同时得到两个修正公式，将它们和上两式构成了如下修正Adams公式： 预测： 修正： 校正:

, 由同阶单步法提供初值 , ， ， 。计算 时，可取 。 修正： 同理，在计算时，调节计算步长 使 , 由同阶单步法提供初值 , ， ， 。计算 时，可取 。 理论分析得出它们的绝对稳定区域，修正Hamming法： ; 4阶Adams预测-校正法： 其中 我们也可以在教学演示上看出修正的预测校正格式的误差非常的小。

4.多步法的相容性、稳定性和收敛性 多步法的相容性条件 k步法的一般形式为 其中 由对单步法的讨论可知,当 时,方法阶数至少为1。即对 y＝1,y=x 应精确成立。令 y ＝1，所以令y=x得

所以 我们称为线性多步法的相容性条件。

多步法的收敛问题 k步法需要k个出发值，而初值问题只提供了一个初值，在使用k步法时尚需要其它方法作补充k-1个出发值，今假定它们 应收敛于共同的极限y（a）即 在讨论多步法收敛性时我们总假定（9.3.12）成立 定义四：

＝x＝a+nh. 的解 收敛于初值问题的解 y(x)。这里 定义五：称多项式 为k步法的特征多项式。如果特征多项式的零点皆不大于1，且等 于1的零点是单重的，称根条件成立。称多项式 为第二特征多项式。 显然根条件可以表示 定理三：k步法收敛的充要条件为： （1）相容性条件成立。 （2）特征多项式的零点皆不大于1，且等于1的零点是单重的。 （称2为）特征根条件。

多步法的稳定性 关于线性多步法成立以下定理： 定理四：若函数f(x,y)对变量y满足Lipschtiz 条件在与h,x无关的常数L，对区域D= 任意两点（x，y1）,(x,y2)成立 k步法的相容性、收敛性、稳定性有以下关系 对于常微分方程右端函数f(x,y)，在相容性条件成立情况下，k步法的收敛性和稳定性是等价的。 事实上在相容性条件成立时，收敛性和稳定性的充要条件都是特征根条件成立。

多步法的绝对稳定性和绝对稳定域 将线性多步法的公式应用到试验方程 进行考察。这里假定Re <0,即试验方程本身是稳定的。 记 得 记 得 是常系数差分方程，其特征方程为 记它的 k个特征根为 并设它们是互异的。显然根与 有关，不妨记为 注意到当 时 这时由特征方程得 由线性多步法的相容性条件得 是一个根。不妨设, 差分方程的解为 其中系数由线性多步法的出发值确定。

另一方面，y(0)=1的试验方程的精确解为 ， 设多步法截断误差为 ，由此可得 ， 我们称 为主根，其它根都为增根。 定义五：线性多步法的绝对稳定区域 对给定的 ，如果特征方程的特征根 皆按模小于1，则线性多步法关于u是绝对稳定的。使得 成立的 构成绝对稳定区域。 注：从误差角度来看绝对稳定区域的方法是一个理想的方法。这样，绝对稳定区域越大，方法适用性越广，因而越优越。

§6 方程组和刚性方程 9.1 一阶方程组 9.2 化高阶方程为一阶方程组 9.3 刚性方程组

9.1 一阶方程组 前面我们研究了单个方程y’=f的数值解法，只要把y和f理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形。 9.1 一阶方程组 前面我们研究了单个方程y’=f的数值解法，只要把y和f理解为向量，那么，所提供的各种计算公式即可应用到一阶方程组的情形。 考虑一阶方程组的出值问题 若采用向量的记号，记

则常微分方程组的出值问题可以表示为 前几节我们主要讨论了常微分方程的出值问题的数值解法。只要将 y 和f 改写为向量，那么前面提供的各种计算公式即可适用于一阶常微分方程组的出值问题。 Runge-Kutta方法 对于方程组 四级四阶显示Runge- Kutta公式

若写成分量形式就是 i=1,2,…,m 为了帮助理解这一公式的计算过程，我们再考虑两个方程的特殊情形：

这时四阶龙格－库塔公式具有形式

其中

这是一步法，利用节点 上的值 ， ，由上式顺序计 算 ，然后代入 即可求得 节点上的 。

9.2化高阶方程为一阶方程组 关于高阶微分方程的初值问题，原则上总可以归结为一阶方程组来求解，例如，考察下列m 阶微分方程 初始条件为 只要引进新的变量 即可将m 阶方程化为如下的一阶方程组

初始条件则相应的化为

9.3 刚性方程组 在求解方程组 时，经常出现解的分量数量级差别很大的情形，这给数值求解带来很大困难，这种问题称为刚性问题。 9.3 刚性方程组 在求解方程组 时，经常出现解的分量数量级差别很大的情形，这给数值求解带来很大困难，这种问题称为刚性问题。 若线性系统 ，其中， 中A的特征值 满足条件<0(j=1,…,N),

且 则称 系 统为刚性方程，称s为刚性比

§7 习题与总结 1数值例题 2数值练习 3总结

1、数值例题 我们已经学习了很多数值算法，他们的效果到底如何呢？ 下面我们来分析一道例题，看看那些方法，就这个问题，最能接近真实值 求初值问题 的数值解，取h=0.1并同精确解 比较 （1）用欧拉法来计算这个初值问题 （2）用各阶的Runge—Kutta方法来计算这个初值问题 （3）用四阶的Adams定步长，Adams定步长预测校正方 法来 计算这个初值问题。然后将数值结果列成表格：

自变量 解析解 欧拉法 后退欧拉法 梯形法 改进的欧拉法 0.1 1.0048 1.0000 1.0091 1.0050 0.2 1.0187 1.0100 1.0264 1.0186 1.0190 0.3 1.0408 1.0290 1.0513 1.0406 1.0412 0.4 1.0703 1.0561 1.0830 1.0701 1.0708 0.5 1.1065 1.0905 1.1209 1.1063 1.1071 0.6 1.1488 1.1314 1.1645 1.1485 1.1494 0.7 1.1966 1.1783 1.2132 1.1963 1.1972 0.8 1.2493 1.2305 1.2665 1.2490 1.2500 0.9 1.3066 1.2874 1.3241 1.3063 1.3072 1.0 1.3679 1.3488 1.3855 1.3676 1.3685

自变量 解析解 4阶显RK 4阶Adams显法 预测校正法 0.1 1.0048 0.2 1.0187 0.3 1.0408 0.4 1.0703 0.5 1.1065 0.6 1.1488 0.7 1.1966 0.8 1.2493 0.9 1.3066 1.0 1.3679

2、数值练习 2. 利用4阶Runge-Kutta算法求初值问题 计算实习题： 1.使用改进的Euler算法求初值问题 的数值解，取步长h=0.1 ,并同精确解 比较 2. 利用4阶Runge-Kutta算法求初值问题

4.使用4阶定步长Adams预测－校正算法求初值问题 的数值解，取步长：(1)h=0.1;(2)h=0.01.求数值解及精确解 3.使用4阶定步长Adams预测－校正算法求初值问题 的数值解，取h=0.1,并同精确解 比较 4.使用4阶定步长Adams预测－校正算法求初值问题 的数值解，使用h=0.5

5．综合题 请用已学的数值方法来计算下面初值问题，看看那种精度最好， 的数值解。

总 结 本章研究的是常微分方程初值问题的数值解法，我们从单步法、多步法两大类来介绍经典算法，并且对相容性、收敛性、稳定性这些理论进行系统的、全面介绍，为同学们学习偏微分方程数值解法作了较好的铺垫。 梯形法、R-K方法、变步长法以及Adams方法，都是一些高精度的方法，但是它们也有缺点，四阶R-K方法必须要求f(x,y)有较好的光滑性，若光滑性差，还不如欧拉法或改进欧拉公式。 它的另一个缺点是计算量较大，需要好较多的机器时间（每一步需要计算4次f(x,y)），相比之下汉明方法节省计算量，但它是四步方法不是自开始的需要借助四阶R-K方法提供初始值。 对数值方法的分析还涉及到局部截断误差，整体误差，相容性、收敛性、稳定性等概念，他们涉及到步长的选取，如步长不合适，舍入误差恶性增长，结果会完全错误。大家在实际应用时一定要注意这些问题。

刚性方程组具有很重要的应用价值，它的理论和解法很多，并且还在进一步发展，我们仅是简单介绍一下，详细介绍请查看参考书目。

抽取客观事物的本质特征及其内部联系，用适当的数学语言和工具所建立的一个数学结构 什么是数学模型 抽取客观事物的本质特征及其内部联系，用适当的数学语言和工具所建立的一个数学结构 例如著名的Malthus人口模型

又如弱肉强食模型

例1 用欧拉方法，隐式欧拉方法和欧拉中点公式计算 的近似解，取步长h=0.1，并与精确值比较 解:欧拉方法的算式为： 欧拉隐式方法在本题可解出方程，不必迭代，公式为： 欧拉中点公式是两步法，第一步y1用欧拉公式，以后用公式

利用 Euler方法求初值问题 例1 的数值解.此问题的精确解是 解 此时的 Euler公式为 分别取步长h=0.2 ,0.1 ,0.05,计算结果如下

h xn yn y(xn) y(xn)-yn h=0.2 0.00 0.40 0.80 1.20 1.60 2.00 0.00000 0.37631 0.54228 0.52709 0.46632 0.40682 0.34483 0.48780 0.49180 0.44944 0.40000 -0.03148 -0.05448 -0.03529 -0.01689 -0.00682 h=0.1 0.36085 0.51371 0.50961 0.45872 0.40419 -0.01603 -0.02590 -0.01781 -0.00928 -0.00419 h=0.05 0.35287 0.50049 0.50073 0.45425 0.40227 -0.00804 -0.01268 -0.00892 -0.00481 -0.00227

例2 .用欧拉公式和梯形公式的预估校正法计算： 的数值解，取h=0.1，梯形公式只迭代一次，并与精确值比较.方程的解析解为: 解: 本例中欧拉公式为：