（ Numerical Methods for Ordinary Differential Equations ）

常微分方程分为 （ 1 ）初值问题 （ 2 ）边值问题 一、初值问题的数值解法 §1 引 言 一阶常微分方程初值问题的一般形式是 :

称 在区域 D 上对 满足 Lipschitz 条件 是指 :

二、初值问题解的存在性  考虑一阶常微分方程的初值问题 /* Initial-Value Problem */ : 则上述 IVP 存在唯一解。 只要 在 上连续, 且关于 y 满足 Lipschitz 条件， 即存在与 无关的常数 L 使 对任意定义在 上的 都成立，

要计算出解函数 y(x) 在一系列节点 a = x 0 < x 1 <…< x n = b 处的近似值 采用离散化方法。 称节点间距 为步长， 通常采用等距节点，即取 h i = h ( 常数 ) 。

三、初值问题的离散化方法 离散化方法的基本特点是依照某一递推公式， 值, 取 。 按节点从左至右的顺序依次求出 的近似 如果计算 ，只用到前一步的值 ，则称这类方 法为单步方法。 如果计算 需用到前 r 步的值,, 则称这类方法为 r 步方法。

§2 欧拉方法 /* Euler’s Method */ /* Euler’s Method */  欧拉公式 ( 单步显示公式）： 向前差商近似导数 记为 x0x0 x1x1

亦称为欧拉折线法 /* Euler’s polygonal arc method*/ 在假设 y i = y(x i ) ，即第 i 步计算是精确的前提 下，考虑的截断误差 R i = y(x i+1 )  y i+1 称为局部截断 误差 /* local truncation error */ 。 定义 若某算法的局部截断误差为 O(h p+1 ) ，则称该 算法有 p 阶精度。 定义

 欧拉法的局部截断误差： R i 的主项 /* leading term */ 欧拉法具有 1 阶精度。

例 1: 用欧拉公式求解初值问题 取步长 。 解 : 应用 Euler 公式于题给初值问题的具体形式为： 其中 。计算结果列于下表：

可用来检验近似解的准确程度。 进行计算，数值解已达到了一定的精度。 这个初值问题的准确解为 ， 从上表最后一列，我们看到取步长

 欧拉公式的改进：  隐式欧拉法 /* implicit Euler method */ 向后差商近似导数 x0x0 x1x1 )(,()( 1101 xyxfhyxy 

由于未知数 y i+1 同时出现在等式的两边，不能直接 得到，故称为隐式 /* implicit */ 欧拉公式，而前者 称为显式 /* explicit */ 欧拉公式。

一般先用显式计算一个初值，再迭代求解。  隐式欧拉法的局部截断误差： 即隐式欧拉公式具有 1 阶精度。

 梯形公式 /*trapezoid formula */ — 显、隐式两种算法的平均 注：的确有局部截断误差 ， 即梯形公式具有 2 阶精度，比欧拉方法有了进步。 但注意到该公式是隐式公式，计算时不得不用到 迭代法，其迭代收敛性与欧拉公式相似。

 中点欧拉公式 /* midpoint formula */ 中心差商近似导数 x0x0 x2x2 x1x1 假设, 则可以导出 即中点公式具有 2 阶精度。

方 法  显式欧拉 隐式欧拉 梯形公式 中点公式 简单精度低 稳定性最好精度低, 计算量大 精度提高计算量大 精度提高, 显式 多一个初值, 可能影响精度 Can’t you give me a formula with all the advantages yet without any of the disadvantages? Do you think it possible? Well, call me greedy… OK, let’s make it possible.

 改进欧拉法 /* modified Euler’s method */ Step 1 : 先用显式欧拉公式作预测，算出 Step 2 : 再将 代入隐式梯形公式的右边作校正，得到 1  i y

注 : 此法亦称为预测 - 校正法 /* predictor-corrector method */ 可以证明该算法具有 2 阶精度，同时可以看到它 是个单步递推格式，比隐式公式的迭代求解过程 简单。后面将看到，它的稳定性高于显式欧拉法。 改进的欧拉法

在实际计算时，可将欧拉法与梯形法则相结合， 计算公式为 应用改进欧拉法, 如果序列 收敛, 它的极限便满足方程

改进欧拉法的截断误差 因此，改进欧拉法公式具有 2 阶精度

例 2: 用改进 Euler 公式求解例 1 中的初值问题， 取步长 。 解：对此初值问题采用改进 Euler 公式， 其具体形式为 计算结果列于下表：

改进的 Euler 法 Euler 法

通过计算结果得比较可以看出，改进的 Euler 方法 的计算精度比 Euler 方法要高。

§3 龙格 - 库塔法 /* Runge-Kutta Method */ /* Runge-Kutta Method */ 建立高精度的单步递推格式。 单步递推法的基本思想是从 ( x i, y i ) 点出发， 欧拉法及其各种变形所能达到的最高精度为 2 阶。 以某一斜率沿直线达到 点。

 考察改进的欧拉法，可以将其改写为： 斜率 一定取 K 1 K 2 的平均值吗？ 步长一定是一个 h 吗？

将改进欧拉法推广为： ),( ),( ][ phKyphphxfK yxfK KKhyy ii ii ii     首先希望能确定系数 1 、 2 、 p ，使得到的算法格 式有 2 阶精度，即在 的前提假设下，使得 Step 1: 将 K 2 在 ( x i, y i ) 点作 Taylor 展开

Step 2: 将 K 2 代入第 1 式，得到 Step 3: 将 y i+1 与 y( x i+1 ) 在 x i 点的泰勒展开作比较

要求 ，则必须有： 这里有 个未知 数， 个方程。 3 2 存在无穷多个解。所有满足上式的格式统称为 2 阶 龙格 - 库塔格式。 注意到， 就是改进的欧拉法。

为获得更高的精度，应该如何进一步推广？ 其中 i ( i = 1, …, m ) ，  i ( i = 2, …, m ) 和  ij ( i = 2, …, m; j = 1, …, i  1 ) 均为待定系数，确定这些 系数的步骤与前面相似。

 考虑一阶常微分方程初值问题

将区域 [a,b] 进行分划：

若

则

n 级显式 Runge-Kutta 方法

二阶 Runge-Kutta 方法 取 n=2 记

由此得

另一方面 为使局部截断误差为, 应取

改进的 Euler 方法 取取

中点方法 取

二阶 Heun 方法 取

二级 Runge-Kutta 方法不超过二阶 记记 则

因此局部截断误差只能达到

三级 Runge-Kutta 方法 取 n=3

记

 又由于

因此要使局部截断误差为 ，必须

Kutta 方法 取

三阶 Heun 方法 取取

三级 Runge-Kutta 方法不超过三阶  完全类似于二级 Runge-Kutta 方法的分析 只能达到 三级 Runge-Kutta 方法的局部截断误差 将 和 都展开到 项, 易证

四级 R-K 方法 取 n=4

局部截断误 差为 O(h 5 )  最常用为四级 4 阶经典龙格 - 库塔法 /* Classical Runge-Kutta Method */ ：

附注：  二阶 Runge-Kutta 方法的局部截断误差 只能达到  五阶 Runge-Kutta 方法的局部截断误差 只能达到  四阶 Runge-Kutta 方法的局部截断误差 只能达到  三阶 Runge-Kutta 方法的局部截断误差 只能达到

注：  龙格 - 库塔法的主要运算在于计算 的值，即计 算 的值。 Butcher 于 1965 年给出了计算量与可达 到的最高精度阶数的关系： 753 可达到的最高精度 642 每步须算 K i 的个数  由于龙格 - 库塔法的导出基于泰勒展开，故精 太好的解，最好采用低阶算法而将步长 h 取小。 度主要受解函数的光滑性影响。对于光滑性不

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

解：该问题的精确解为 欧拉公式为 对任意固定的 x = x i = i h ，有 

 稳定性 /* Stability */ 例：考察初值问题 在区间 [0, 0.5] 上 的解. 分别用欧拉显、隐式格式和改进的欧拉格式 计算数值解。 精确解改进欧拉法 欧拉隐式欧拉显式 节点 x i    10 1    10   10   10   10   10      10   10   10   10   10  7 What is wrong ??! An Engineer complains: "Math theorems are so unstable that a small perturbation on the conditions will cause a crash on the conclusions!"

定义 若某算法在计算过程中任一步产生的误差在 以后的计算中都逐步衰减，则称该算法是绝对稳定 的 /*absolutely stable */ 。 一般分析时为简单起见，只考虑试验方程 /* test equation */ 常数，可以 是复数

我们称算法 A 比算法 B 稳定，就是指 A 的绝对稳定 区域比 B 的大。 当步长取为 h 时，将某算法应用于上式，并假设只在 初值产生误差 ，则若此误差以后逐步衰减, 就称该算法相对于 绝对稳定， 的全体构成 h hh  h 绝对稳定区域。

例：考察显式欧拉法 Re Img 由此可见，要保证初始误差  0 以后逐步衰减， 必须满足：

例：考察隐式欧拉法 可见绝对稳定区域为： 210Re Img 注：一般来说，隐式欧拉法的绝对稳定性比同阶 的显式法的好。

例：隐式龙格 - 库塔法 其中 2 阶方法 的绝对稳定 区域为 0 Re Img 无条件稳定

而显式 1~ 4 阶方法的绝对稳定区域为 k=1 k=2 k=3 k= Re Img