含有未知函数及其导数或微分的等式, 常称为微分方程. 科学研究和工程技术中的许多问题都可通过微分方程来建立数学模型进行分析. 一般而言, 除了少数特殊类型的微分方程能用解析方法求得精确解外, 多数情况找不到解的解析表达式. 因此, 采用适当的数值方法来求解微分方程已成为微分方程求解的主要手段.
本章主要研究两类常微分问题：一类是一阶常微分方程的初值问题
$$ \left\{\begin{array}{l l} \boldsymbol{y}^{\prime} = \boldsymbol{f} (x, \boldsymbol{y}), & x \in [ a, b ], \\ \boldsymbol{y} (\boldsymbol{a}) = \boldsymbol{y}_{0}, \end{array} \right.\tag{9.1} $$其中 $\boldsymbol{y}(x)$ 是定义在 $[a,b]$ 上的 m 维函数向量, $\boldsymbol{f}(x,\boldsymbol{y})$ 是定义在 $m+1$ 维区域
$$ G = \{(x, \boldsymbol{y}) \mid x \in [ a, b ], \boldsymbol{y} \in \mathbb{R}^{m} \}\tag{9.2} $$上的 $m$ 维已知函数向量. 由常微分方程理论知: 如果函数 $\pmb{f}(x, \pmb{y})$ 在区域 $G$ 中连续, 且关于 $y$ 满足利普希茨 (Lipschitz) 条件, 即对所有 $x \in [a, b]$ 及 $\pmb{y} \in \mathbb{R}^m$ , $\pmb{z} \in \mathbb{R}^m$ , 总存在常数 $L > 0$ , 使得
$$ \left\| \boldsymbol{f} (x, \boldsymbol{y}) - \boldsymbol{f} (x, \boldsymbol{z}) \right\| \leqslant L \| \boldsymbol{y} - \boldsymbol{z} \|,\tag{9.3} $$那么初值问题方程 (9.1) 存在唯一解, 且解连续依赖于初始条件和右端项.
另一类是两阶常微分方程边值问题
$$ \left\{\begin{array}{l l} - y^{\prime \prime} + q (x) y = f (x), & x \in [ a, b ], \\ y (a) = \alpha , y (b) = \beta , \end{array} \right.\tag{9.4} $$其中 $q(x)$ 和 $f(x)$ 在区间 $[a, b]$ 上连续, $q(x) > 0$ . 这里假设上述边值问题存在唯一解, 且解连续依赖于边界条件和右端项.
无论是初值问题还是边值问题, 其解 $\boldsymbol{y} = \boldsymbol{y}(x)$ 都是区间 $[a, b]$ 上关于变量 x 的函数或函数向量. 记 $a = x_{0} < x_{1} < \cdots < x_{N-1} < x_{N} = b$ 为求解区域中的一系列节点. 数值解就是要计算精确解 $\boldsymbol{y}(x)$ 在这些节点 $x_{n}$ 处的近似值 $y_{n}$ . 为简单起见, 本章假设点列满足
$$ h_{n} = x_{n + 1} - x_{n} = h = \frac{b - a}{N},\tag{9.5} $$其中 h 称为步长. 本章主要研究 m = 1 时的初值问题, 对于边值问题和 m > 1 的情况下的初值问题将分别在最后两节作简单介绍.
9.1 欧拉公式及其改进#
9.1.1 欧拉公式#
欧拉公式是求解初值问题 (9.1) 的一种简单而古老的数值方法. 其基本思路是把方程 (9.1) 中的导数项 $y'$ 用差商逼近, 从而把一个微分方程转化成为一个代数方程, 以便求解.
在节点 $x_{n}$ 处，(9.1)式可写成
$$ y^{\prime} (x_{n}) = f (x_{n}, y (x_{n})).\tag{9.6} $$由于
$$ y^{\prime} (x_{n}) \approx \frac{y (x_{n + 1}) - y (x_{n})}{h},\tag{9.7} $$故
$$ y (x_{n + 1}) \approx y (x_{n}) + h f (x_{n}, y (x_{n})).\tag{9.8} $$因为真解 $y(x)$ 是未知的, 所以将 $y(x_{i})$ 的近似值 $y_{i}, (i = n, n + 1)$ 代入 (9.8) 式. 这就得到欧拉公式
$$ y_{n + 1} = y_{n} + h f (x_{n}, y_{n}).\tag{9.9} $$欧拉公式的几何意义是十分清楚的．方程 $y' = f(x, y)$ 满足初始条件 $y(x_{0}) = y_{0}$ 的解 $y = y(x)$ 是 xOy 平面上过点 $P_{0}(x_{0}, y_{0})$ 的一条特殊积分曲线．欧拉方法就是用从 $P_{0}$ 出发的折线 $P_{0}P_{1}\cdots P_{N}$ 来作为积分曲线 $y = y(x)$ 的近似解，如图 9-1 所示．所以欧拉方法又称为欧拉折线法．
图 9-1 欧拉法的几何意义
以下是欧拉方法的算法描述:
算法 9.1.1(欧拉方法)
(1) 给定整数 N;
(2) $h = \frac{b - a}{N}$ ;
(3) $y(1) = y_0$ ;
(4) 对 $k = 2, 3, \cdots, N + 1$ 做:
$$ x = a + (k - 1) h; $$$$ y (k) = y (k - 1) + h f (x, y (k - 1)); $$(5) 返回 y.
基于欧拉方法的 Matlab 程序如下:
$$ \begin{array}{l} \text{for i = 1:n,} \\ \quad \mathrm{y(i + 1) = y(i) + h*feval(f,x(i),y(i));} \\ \text{end} \end{array} $$例 9.1.2 用欧拉方法求解初值问题
$$ \left\{\begin{array}{l} y^{\prime} = y - \frac{2 x}{y}, \\ y (0) = 1, \end{array} \right. $$在区间 $[0,1]$ 上，步长 $h = 0.1$ 的数值解.
解：由欧拉公式 (9.9) 得到具体计算公式
$$ y_{n + 1} = y_{n} + h \left(y_{n} - \frac{2 x_{n}}{y_{n}}\right).\tag{9.10} $$令 $y_{0}=y(x_{0})$ ，数值计算结果见表 9-1.
表 9-1 欧拉法的计算结果
| $x_{k}$ | $y_{k}$ | $y(x_{k})$ | $|y_{k}-y(x_{k})|$ |
| 0.0 | 1.000 0 | 1.000 0 | 0.000 0 |
| 0.1 | 1.100 0 | 1.095 4 | 0.004 6 |
| 0.2 | 1.191 8 | 1.183 2 | 0.008 6 |
| 0.3 | 1.277 4 | 1.264 9 | 0.012 5 |
| 0.4 | 1.358 2 | 1.341 6 | 0.016 6 |
| 0.5 | 1.435 1 | 1.414 2 | 0.020 9 |
| 0.6 | 1.509 0 | 1.483 2 | 0.025 7 |
| 0.7 | 1.580 3 | 1.549 2 | 0.031 1 |
| 0.8 | 1.649 8 | 1.612 5 | 0.037 3 |
| 0.9 | 1.717 8 | 1.673 3 | 0.044 5 |
| 1.0 | 1.784 8 | 1.732 1 | 0.052 7 |
该初值问题的解析解为 $y = \sqrt{1 + 2x}$ . 为便于比较, 表 9-1 中同时列出了解析解在各节点处的值. 两者相比较, 显然欧拉方法给出的数值解误差较大, 一般只有两位有效数字.
9.1.2 数值积分与多步法#
为了改进欧拉方法, 我们从另一角度来考察初值问题. 事实上, (9.1) 式可改写为
$$ \mathrm{d} y (x) = f (x, y (x)) \mathrm{d} x.\tag{9.11} $$两边关于 $x$ 在区间 $[x_{n}, x_{n+1}]$ 上积分, 得
$$ y (x_{n + 1}) = y (x_{n}) + \int_{x_{n}}^{x_{n + 1}} f (x, y (x)) \mathrm{d} x.\tag{9.12} $$要得到 $y(x_{n+1})$ 的值, 就必须计算右端的积分
$$ \int_{x_{n}}^{x_{n + 1}} f (x, y (x)) \mathrm{d} x.\tag{9.13} $$但积分式中含有未知函数 $y(x)$ ，无法直接计算，因此考虑采用数值积分方法计算积分式 (9.13) 的近似值。这种将初值问题方程 (9.1) 离散化的方法称为初值问题的数值积分方法。
选用不同的数值积分公式, 就会得到不同的计算公式.
例如, 用最简单的左矩形公式
$$ \int_{x_{n}}^{x_{n + 1}} f (x, y (x)) \mathrm{d} x \approx h f (x_{n}, y (x_{n})).\tag{9.14} $$将 (9.14) 式代入 (9.12) 式, 可得
$$ y (x_{n + 1}) \approx y (x_{n}) + h f (x_{n}, y (x_{n})).\tag{9.15} $$如用 $y_{n}$ 代替上式右端的 $y(x_{n})$ ，并将算出的右端值作为 $y(x_{n+1})$ 的近似值 $y_{n+1}$ ，这样建立起来的计算公式就是前述的欧拉公式.
类似地, 利用右矩形公式, 可得到后退欧拉公式
$$ y_{n + 1} = y_{n} + h f (x_{n + 1}, y_{n + 1}).\tag{9.16} $$矩形求积公式的代数精度较低. 为了提高计算精度, 我们用梯形求积公式在 $[x_n, x_{n+1}]$ 上计算右端的积分, 即
$$ \int_{x_{n}}^{x_{n + 1}} f (x, y (x)) \mathrm{d} x \approx \frac{h}{2} (f (x_{n}, y (x_{n})) + f (x_{n + 1}, y (x_{n + 1}))).\tag{9.17} $$于是得到计算初值问题的梯形公式
$$ y_{n + 1} = y_{n} + \frac{h}{2} (f (x_{n}, y_{n}) + f (x_{n + 1}, y_{n + 1})).\tag{9.18} $$由于梯形求积公式比矩形求积公式的代数精度高, 因此一般情况下梯形公式 (9.18) 的精度比欧拉公式 (9.9) 高.
以下是梯形方法的算法描述.
算法 9.1.3(梯形方法)
(1) 给定整数 N;
(2) $h = \frac{b - a}{N}$ ;
(3) $y(1) = y_0$ ;
(4) 对 $k = 2, 3, \cdots, N + 1$ 做:
$$ x = a + (k - 1) h; $$$$ \text{求解方程} y (k) = y (k - 1) + 0.5 h (f (x, y (k - 1)) + f (x + h, y (k)); $$(5) 返回 y.
例 9.1.4 用梯形公式 (9.18) 求解例 9.1.2 中的初值问题.
解：对于例 9.1.2 用公式 (9.18) 的具体形式为
$$ y_{n + 1} = y_{n} + \frac{h}{2} \left(\left(y_{n} - \frac{2 x_{n}}{y_{n}}\right) + \left(y_{n + 1} - \frac{2 x_{n + 1}}{y_{n + 1}}\right)\right).\tag{9.19} $$这是一个非线性方程, 可利用非线性方程的各种求解方法进行计算. 这里采用简单的迭代法
$$ y_{n + 1}^{(t + 1)} = y_{n} + \frac{h}{2} \left(\left(y_{n} - \frac{2 x_{n}}{y_{n}}\right) + \left(y_{n + 1}^{(t)} - \frac{2 x_{n + 1}}{y_{n + 1}^{(t)}}\right)\right).\tag{9.20} $$由压缩映像定理可以证明：当 h 充分小时，该迭代法收敛于唯一不动点.
取步长 h = 0.1, 初始值 $y_{0} = y(0)$ , 计算结果见表 9-2. 同例 9.1.2 中欧拉方法计算的结果比较, 梯形公式明显地改善了精度.
表 9-2 梯形公式的计算结果
| $x_{k}$ | $y_{k}$ | $y(x_{k})$ | $|y_{k}-y(x_{k})|$ |
| 0.0 | 1.000 0 | 1.000 0 | 0.000 0 |
| 0.1 | 1.095 7 | 1.095 4 | 0.000 2 |
| 0.2 | 1.183 6 | 1.183 2 | 0.000 4 |
| 0.3 | 1.265 4 | 1.264 9 | 0.000 5 |
| 0.4 | 1.342 3 | 1.341 6 | 0.000 7 |
| 0.5 | 1.415 1 | 1.414 2 | 0.000 8 |
| 0.6 | 1.484 3 | 1.483 2 | 0.001 0 |
| 0.7 | 1.550 4 | 1.549 2 | 0.001 2 |
| 0.8 | 1.613 9 | 1.612 5 | 0.001 5 |
| 0.9 | 1.675 1 | 1.673 3 | 0.001 8 |
| 1.0 | 1.734 1 | 1.732 1 | 0.002 1 |
利用高阶数值积分, 还可构造其他离散格式. 例如,
亚当斯－巴什福思（Adams-Bashforth）公式
$$ y_{n + 1} = y_{n} + \frac{h}{24} (55 f (x_{n}, y_{n}) - 59 f (x_{n - 1}, y_{n - 1}) + 37 f (x_{n - 2}, y_{n - 2}) - 9 f (x_{n - 3}, y_{n - 3})).\tag{9.21} $$亚当斯－莫尔顿（Adams-Moulton）公式
$$ y_{n + 1} = y_{n} + \frac{h}{24} (9 f (x_{n + 1}, y_{n + 1}) + 19 f (x_{n}, y_{n}) - 5 f (x_{n - 1}, y_{n - 1}) + f (x_{n - 2}, y_{n - 2})).\tag{9.22} $$欧拉公式与梯形公式仅含有 $y_{n}$ 和 $y_{n+1}$ ，只要知道 $y_{n}$ 即可求出 $y_{n+1}$ ，这类公式称为单步法。而用亚当斯-巴什福思公式与亚当斯-莫尔顿公式计算 $y_{n+1}$ 时，除了要知道 $y_{n}$ 外，还需知道 $y_{n-1}, y_{n-2}$ 等，这类公式称为多步法。
欧拉公式和亚当斯-巴什福思公式右端不含有 $y_{n + 1}$ ，把右端的已知量算出来即可求出 $y_{n + 1}$ ，这类方法称为显式格式.梯形公式与亚当斯-莫尔顿公式的右端隐含 $y_{n + 1}$ ，求 $y_{n + 1}$ 时通常需要解方程,此类公式称为隐式格式.特别是当 $f(x,y)$ 关于 $y$ 为非线性时，就需要解非线性方程
以上这些不同的离散格式在计算精度上的表现是不同的. 下面介绍局部截断误差的概念. 局部截断误差可用于表征求解初值问题的数值方法的计算精度.
一般地, 常微分方程初值问题的数值解满足形如
$$ y_{n + 1} = y_{n} + h g (y_{n + 1}, y_{n}, \dots , y_{n - r})\tag{9.23} $$的等式, 其中 $y_{n}, \cdots, y_{n-r}$ 为 y 在 $r+1$ 个节点 $x_{n}, \cdots, x_{n-r}$ 处的数值解.
假设 (9.23) 中等号右端的数值解都是准确的, 则利用该式得到的数值解 $y_{n+1}$ 与真解 $y(x_{n+1})$ 之间的差异
$$ \varepsilon_{n + 1} = y (x_{n + 1}) - y (x_{n}) - h g (y (x_{n + 1}), y (x_{n}), \dots , y (x_{n - r})),\tag{9.24} $$称为数值方法的局部截断误差. 而真解与近似解之间的差异 $e_n = y(x_n) - y_n$ 称为数值方法的整体截断误差. 在第 §9.3 节, 我们将看到局部截断误差与整体截断误差有着密切的联系, 并在一定程度上刻画了数值方法的计算精度.
例 9.1.6 欧拉公式是一阶方法.
证：设 $y(x)$ 在 $[x_0, x_N]$ 上充分光滑， $M = \max_{x_0 \leqslant x \leqslant x_N} |y''(x)|$ ，由泰勒展开可得
$$ y (x_{n + 1}) = y (x_{n}) + h y^{\prime} (x_{n}) + \frac{h^{2}}{2} y^{\prime \prime} (\xi), \quad x_{n} < \xi < x_{n + 1}.\tag{9.25} $$假设式 (9.9) 右端的 $y_{n}=y(x_{n})$ ，则
$$ \begin{array}{c} | \varepsilon_{n + 1} | = | y (x_{n + 1}) - y (x_{n}) - h f (x_{n}, y (x_{n})) | \=| y (x_{n + 1}) - y (x_{n}) - h y^{\prime} (x_{n}) | \=frac{h^{2}}{2} | y^{\prime \prime} (\xi) | \leqslant \frac{M}{2} h^{2} = O (h^{2}) \end{array}\tag{9.26} $$这表明欧拉公式的局部截断误差是 $h^{2}$ 的同阶无穷小, 即欧拉公式是一阶方法.
类似地, 可以证明梯形公式是二阶方法, 而亚当斯-巴什福思公式与亚当斯-莫尔顿公式均为四阶方法.
9.1.3 预估校正格式#
由例 9.1.2 和例 9.1.4 可以看出, 梯形方法相对于欧拉方法虽然精度提高了, 但是每一步都要用迭代法解方程. 每迭代一次, 都要重新计算函数 $f(x)$ 的值, 所以计算量也增加很多. 实际计算时, 一种有效简化计算的方法是: 当 h 较小时, 先用显格式计算一个合适的预估值, 然后利用隐格式迭代一二次得到校正值. 这就是所谓的预估校正方法.
例如, 先用欧拉公式求得一个初始的近似值 $\overline{y_{n+1}}$ , 称之为预估值. 预估值 $\overline{y_{n+1}}$ 的精度可能很差. 然后以预估值作为初始迭代值代入梯形公式 (9.18) 迭代一次得 $y_{n+1}$ (与例 9.1.4 中简单迭代类似), 计算结果称为校正值, 并以之作为初值问题在节点 $x_{n+1}$ 处的数值解. 实际计算公式如下
$$ \left\{\begin{array}{l l} {\text{预估}} & {\overline{{{y_{n + 1}}}} = y_{n} + h f (x_{n}, y_{n}),} \\ {\text{校正}} & {y_{n + 1} = y_{n} + \frac{h}{2} (f (x_{n}, y_{n}) + f (x_{n + 1}, \overline{{{y_{n + 1}}}})),} \end{array} \right.\tag{9.27} $$这就是一个预估校正公式. (9.27) 式通常称为改进的欧拉公式. 该公式亦可表示为
$$ y_{n + 1} = y_{n} + \frac{h}{2} \left(f (x_{n}, y_{n}) + f (x_{n} + h, y_{n} + h f (x_{n}, y_{n}))\right)\tag{9.28} $$或等价的平均化形式
$$ \left\{\begin{array}{l} y_{p} = y_{n} + h f (x_{n}, y_{n}), \\ y_{c} = y_{n} + h f (x_{n} + h, y_{p}), \\ y_{n + 1} = \frac{1}{2} (y_{p} + y_{c}). \end{array} \right.\tag{9.29} $$以下是改进欧拉方法的算法描述.
算法 9.1.7(改进欧拉方法)
(1) 给定整数 N;
(2) $h = \frac{b - a}{N}$ ;
(3) $y(1) = y_0$ ;
(4) 对 $k = 2, 3, \cdots, N + 1$ 做:
$$ \begin{array}{l} {x = a + (k - 1) h;} \\ {y_{p} = y (k - 1) + h f (x, y (k - 1));} \\ {y_{c} = y (k - 1) + h f (x + h, y_{p});} \\ {y (k) = 0.5 (y_{c} + y_{p});} \end{array} $$(5) 返回 $y$ .
基于改进欧拉方法的 Matlab 程序如下:
$$ \begin{array}{l} \text{function [x,y] = odeIEuler(f,y0,a,b,n)} \\ \quad \mathrm{y(1) = y0;} \\ \quad \mathrm{h = (b - a) / n;} \\ \quad \mathrm{x=a:h:b;} \\ \quad \text{for i = 1:n,} \\ \quad \quad \mathrm{yp = y(i) + h*feval(f,x(i),y(i));} \\ \quad \quad \mathrm{yc = y(i) + h*feval(f,x(i + 1),yp);} \\ \quad \quad \mathrm{y(i + 1) = 0.5*(yp + yc);} \\ \quad \text{end} \end{array} $$例 9.1.8 用改进的欧拉公式 (9.27) 求解例 9.1.2 中的初值问题.
解：对于例 9.1.2, 改进的欧拉公式的具体形式为
$$ \left\{\begin{array}{l} \overline{{y_{n + 1}}} = y_{n} + h \left(y_{n} - \frac{2 x_{n}}{y_{n}}\right), \\ y_{n + 1} = y_{n} + \frac{h}{2} \left(\left(y_{n} - \frac{2 x_{n}}{y_{n}}\right) + \left(\overline{{y_{n + 1}}} - \frac{2 x_{n + 1}}{\overline{{y_{n + 1}}}}\right)\right). \end{array} \right.\tag{9.30} $$取步长 h = 0.1, 初始值 $y_{0} = y(0)$ , 计算结果见表 9-3.
表 9-3 改进欧拉法的计算结果
| $x_{k}$ | $y_{k}$ | $y(x_{k})$ | $|y_{k}-y(x_{k})|$ |
| 0.0 | 1.000 0 | 1.000 0 | 0.000 0 |
| 0.1 | 1.095 9 | 1.095 4 | 0.000 5 |
| 0.2 | 1.184 1 | 1.183 2 | 0.000 9 |
| 0.3 | 1.266 2 | 1.264 9 | 0.001 3 |
| 0.4 | 1.343 4 | 1.341 6 | 0.001 7 |
| 0.5 | 1.416 4 | 1.414 2 | 0.002 2 |
| 0.6 | 1.486 0 | 1.483 2 | 0.002 7 |
| 0.7 | 1.552 5 | 1.549 2 | 0.003 3 |
| 0.8 | 1.616 5 | 1.612 5 | 0.004 0 |
| 0.9 | 1.678 2 | 1.673 3 | 0.004 8 |
| 1.0 | 1.737 9 | 1.732 1 | 0.005 8 |
同例 9.1.2 和例 9.1.4 相比, 改进的欧拉方法的计算精度明显要优于欧拉方法, 但略差于梯形公式. 另一方面, 由于改进的欧拉公式是显格式, 不用求解非线性方程, 所以计算量要远小于梯形公式. 因此, 改进的欧拉公式可视为欧拉公式和梯形公式在计算精度与计算量方面的折中方案.
一般地, 较为简单的预估校正格式都包含两个计算公式, 一个是显式公式, 作为预估公式; 另一个是隐式公式, 作为校正公式, 当然也可以构造包含多个计算公式的预估校正格式. 构造预估校正格式时, 应该注意阶数的匹配, 例如在式 (9.27) 中, 校正公式具有二阶精度, 而预估公式仅具有一阶精度. 因为提供的预估精度较差, 且仅经一次校正, 校正值的精度不会太高.
下面的预估校正格式
$$ \left\{\begin{array}{l l} {\text{预估}} & {\overline{{{y_{n + 1}}}} = y_{n - 1} + 2 h f (x_{n}, y_{n}),} \\ {\text{校正}} & {y_{n + 1} = y_{n} + \frac{h}{2} (f (x_{n}, y_{n}) + f (x_{n + 1}, \overline{{{y_{n + 1}}}}))} \end{array} \right.\tag{9.31} $$要好一些. 因为式 (9.31) 中的预估公式也具有二阶精度, 与校正公式是一致的. 预估值精度提高了, 校正值精度也会更高一些. 式 (9.31) 中的预估公式可用数值积分的中矩形公式推得.
此外, 利用亚当斯-巴什福思公式和亚当斯-莫尔顿公式, 可构造常用的四阶亚当斯预估校正系统
$$ \left\{\begin{array}{l} \overline{{y_{n + 1}}} = y_{n} + \frac{h}{24} (55 f (x_{n}, y_{n}) - 59 f (x_{n - 1}, y_{n - 1}) + 37 f (x_{n - 2}, y_{n - 2}) - 9 f (x_{n - 3}, y_{n - 3})), \\ y_{n + 1} = y_{n} + \frac{h}{24} (9 f (x_{n + 1}, \overline{{y_{n + 1}}}) + 19 f (x_{n}, y_{n}) - 5 f (x_{n - 1}, y_{n - 1}) + f (x_{n - 2}, y_{n - 2})). \end{array} \right.\tag{9.32} $$9.2 龙格－库塔公式#
下面我们介绍一类应用较广泛的高精度的显式单步法——龙格-库塔 (Runge-Kutta) 方法, 简称 R-K 方法. 为了导出龙格-库塔公式, 我们先对前面介绍的欧拉公式和预估校正公式作进一步的分析.
欧拉公式 (9.9) 可改写成
$$ \left\{\begin{array}{c} y_{n + 1} = y_{n} + k_{1}, \\ k_{1} = h f (x_{n}, y_{n}). \end{array} \right.\tag{9.33} $$用它计算 $y_{n+1}$ ，需计算一次 $f(x, y)$ 的值。若设 $y_n = y(x_n)$ ，则 $y_{n+1}$ 的表达式与 $y(x_{n+1})$ 的泰勒展开式的前两项完全相同，即局部截断误差为 $O(h^2)$ 。
预估校正公式 (9.31) 可改写为
$$ \left\{\begin{array}{c} y_{n + 1} = y_{n} + \frac{1}{2} k_{1} + \frac{1}{2} k_{2}, \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{2} = h f (x_{n} + h, y_{n} + k_{1}). \end{array} \right.\tag{9.34} $$用它计算 $y_{n+1}$ ，需计算两次 $f(x, y)$ 的值。若设 $y_n = y(x_n)$ ，则 $y_{n+1}$ 的展开式与 $y(x_{n+1})$ 的泰勒展开式的前三项完全相同，即局部截断误差为 $O(h^3)$ 。
上述两组公式在形式上有一个共同点：都是用 $f(x, y)$ 在某些点上的值的线性组合得出 $y(x_{n+1})$ 的近似值 $y_{n+1}$ ，而且可以看出，增加计算 $f(x, y)$ 的次数，可提高局部截断误差的阶.
龙格-库塔公式的基本思想就是设法计算 $f(x,y)$ 在某些点上的函数值，然后对这些函数值作线性组合，构造近似计算公式，再把近似公式和解的泰勒展开式相比较，使前面的若干项吻合，从而获得达到一定精度的数值计算公式.
具体构造时, 先引进若干参数, 例如一般的显式龙格-库塔公式的形式为
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + \sum_{i = 1}^{r} \omega_{i} k_{i} \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{i} = h f \left(x_{n} + \alpha_{i} h, y_{n} + \sum_{j = 1}^{i - 1} \beta_{ij} k_{j}\right), \quad i = 2, 3, \dots , r. \end{array} \right.\tag{9.35} $$其中参数 $\omega_{i}, \alpha_{i}, \beta_{ij}$ 是与初值问题方程 (9.1) 右端 $f(x, y)$ 和步长 h 均无关的常数. 计算公式 (9.35) 称为 r 段的龙格-库塔公式. 特别地, 若式 (9.35) 与 $y(x_{n+1})$ 的泰勒展开式的前 $p+1$ 项完全一致, 即局部截断误差达到 $O(h^{p+1})$ , 则称公式 (9.35) 为 p 阶 r 段龙格-库塔公式.
下面我们以二段龙格-库塔方法为例, 介绍格式的推导过程. 从龙格-库塔公式的一般形式 (9.35) 出发, 可知二段龙格-库塔格式为
$$ \left\{\begin{array}{c} y_{n + 1} = y_{n} + \omega_{1} k_{1} + \omega_{2} k_{2}, \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{2} = h f (x_{n} + \alpha_{2} h, y_{n} + \beta_{21} k_{1}). \end{array} \right.\tag{9.36} $$由式 (9.36) 得
$$ y_{n + 1} = y_{n} + \omega_{1} h f (x_{n}, y_{n}) + \omega_{2} h f (x_{n} + \alpha_{2} h, y_{n} + \beta_{21} h f (x_{n}, y_{n})).\tag{9.37} $$记 $f(x_{n},y_{n}) = f_{n}$ ，利用泰勒展开得
$$ f \left(x_{n} + \alpha_{2} h, y_{n} + \beta_{21} h f_{n}\right) = f_{n} + \alpha_{2} h f_{x} + \beta_{21} h f_{y} f_{n} + O \left(h^{2}\right),\tag{9.38} $$其中 $f_{x}$ 和 $f_{y}$ 分别表示 $\frac{\partial f}{\partial x}$ 和 $\frac{\partial f}{\partial y}$ 在 $(x_{n},y_{n})$ 处的值.
假设 $y_{n}=y(x_{n})$ ，并将上式代入式 (9.37)，化简得
$$ \begin{array}{c} y_{n + 1} = y (x_{n}) + \omega_{1} h f_{n} + \omega_{2} (h f_{n} + h^{2} (\alpha_{2} f_{x} + \beta_{21} f_{y} f_{n})) + O (h^{3}) \=y (x_{n}) + (\omega_{1} f_{n} + \omega_{2} f_{n}) h + \omega_{2} (\alpha_{2} f_{x} + \beta_{21} f_{y} f_{n}) h^{2} + O (h^{3}). \end{array}\tag{9.39} $$将 $y(x_{n + 1})$ 在点 $x = x_{n}$ 处的泰勒展开
$$ \begin{array}{r} y (x_{n + 1}) = y (x_{n} + h) = y (x_{n}) + y^{\prime} (x_{n}) h + \frac{1}{2 !} y^{\prime \prime} (x_{n}) h^{2} + O (h^{3}) \=y (x_{n}) + f_{n} h + \frac{1}{2} (f_{x} + f_{y} f_{n}) h^{2} + O (h^{3}). \end{array}\tag{9.40} $$逐项进行比较式 (9.39) 和式 (9.40), 为使 $y(x_{n+1}) - y_{n+1} = O(h^3)$ , 则有
$$ \left\{\begin{array}{c} \omega_{1} + \omega_{2} = 1, \\ \omega_{2} \alpha_{2} = \frac{1}{2}, \\ \omega_{2} \beta_{21} = \frac{1}{2}. \end{array} \right.\tag{9.41} $$上述方程组 (9.41) 由关于四个未知元 $\alpha_{2}, \beta_{21}, \omega_{1}, \omega_{2}$ 的三个方程组成. 事实上, 方程 (9.41) 有无穷多组解, 且每一组解都构成一个二阶二段的龙格-库塔公式.
令 $\alpha_{2} = 1$ ，从式(9.41)解得 $\omega_{1} = \omega_{2} = \frac{1}{2},\beta_{21} = 1,$ 则有
$$ \left\{\begin{array}{c} y_{n + 1} = y_{n} + \frac{1}{2} k_{1} + \frac{1}{2} k_{2}, \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{2} = h f (x_{n} + h, y_{n} + k_{1}). \end{array} \right.\tag{9.42} $$这就是前面的预估校正公式.
若令 $\omega_{1} = 0$ ，从式(9.41)解得 $\omega_{2} = 1,\alpha_{2} = \frac{1}{2},\beta_{21} = \frac{1}{2}$ ，则得到另一个二阶二段龙格-库塔公式
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + k_{2}, \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{2} = h f \left(x_{n} + \frac{1}{2} h, y_{n} + \frac{1}{2} k_{1}\right). \end{array} \right.\tag{9.43} $$该格式称为中点公式.
显然, 还有许多其他可能的取法, 相应的计算公式的局部截断误差都是 $O(h^{3})$ . 所以龙格-库塔公式是一类公式, 每确定一组特殊的系数, 就得到一个特殊的龙格-库塔公式.
值得注意的是：对于上述二阶龙格－库塔公式，若计算函数值的次数不增加，则不能通过选择其他参数，来得到更高阶的龙格－库塔公式．这是因为在 $y(x_{n+1})$ 的展开式中的 $f_{x}f_{y} + f_{y}^{2}f$ 是不能通过选择 $\alpha_{2}, \beta_{21}, \omega_{1}, \omega_{2}$ 来消掉的，所以不论四个参数如何选择都不能使局部截断误差的阶达到 $O(h^{4})$ ．这说明在计算二次函数值的情况下，局部截断误差的阶最高是 $O(h^{3})$ ，要进一步提高方法的阶，就必须继续增加计算函数值的次数．
高阶龙格-库塔公式的推导方法与上面讨论的二阶龙格-库塔公式推导类似, 只是随着阶数的增高, 推导的工作量也随着增大. 这里仅列举几个常用的公式.
(1) 三阶三段龙格-库塔公式
$$ \left\{\begin{array}{c} y_{n + 1} = y_{n} + \frac{1}{6} (k_{1} + 4 k_{2} + k_{3}), \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{2} = h f \left(x_{n} + \frac{1}{2} h, y_{n} + \frac{1}{2} k_{1}\right), \\ k_{3} = h f (x_{n} + h, y_{n} - k_{1} + 2 k_{2}). \end{array} \right.\tag{9.44} $$(2) 三阶三段 Heun 公式
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + \frac{1}{4} (k_{1} + 3 k_{3}), \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{2} = h f \left(x_{n} + \frac{1}{3} h, y_{n} + \frac{1}{3} k_{1}\right), \\ k_{3} = h f \left(x_{n} + \frac{2}{3} h, y_{n} + \frac{2}{3} k_{2}\right). \end{array} \right.\tag{9.45} $$(3) 标准四阶四段龙格-库塔公式
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + \frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \\ k_{1} = h f (x_{n}, y_{n}), \\ k_{2} = h f \left(x_{n} + \frac{1}{2} h, y_{n} + \frac{1}{2} k_{1}\right), \\ k_{3} = h f \left(x_{n} + \frac{1}{2} h, y_{n} + \frac{1}{2} k_{2}\right), \\ k_{4} = h f (x_{n} + h, y_{n} + k_{3}). \end{array} \right.\tag{9.46} $$(4) 四阶四段 Gill 公式
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + \frac{1}{6} (k_{1} + (2 - \sqrt{2}) k_{2} + (2 + \sqrt{2}) k_{3} + k_{4}), \\ \quad k_{1} = h f (x_{n}, y_{n}), \\ \quad k_{2} = h f \left(x_{n} + \frac{1}{2} h, y_{n} + \frac{1}{2} k_{1}\right), \\ \quad k_{3} = h f \left(x_{n} + \frac{1}{2} h, y_{n} + \frac{\sqrt{2} - 1}{2} k_{1} + \frac{2 - \sqrt{2}}{2} k_{2}\right), \\ \quad k_{4} = h f \left(x_{n} + h, y_{n} - \frac{\sqrt{2}}{2} k_{2} + \frac{2 + \sqrt{2}}{2} k_{3}\right). \end{array} \right.\tag{9.47} $$在实际应用中, 最常用的是标准四阶四段龙格-库塔公式.
从理论上说, 我们可构造任意高阶的计算方法, 但方法的阶数与计算函数值的次数之间的关系并非等量增加的, 它们之间的关系如表 9-4.
表 9-4 龙格－库塔公式阶数和次数的关系
| 计算函数值次数 r | 1 | 2 | 3 | 4 | 5 | 6 | 7 | r ≥ 8 |
| 阶数 p | 1 | 2 | 3 | 4 | 4 | 5 | 6 | r - 2 |
从表 9-4 可以看出, 四阶及四阶以下的龙格-库塔公式, 每一步需调用函数 $f(x,y)$ 的次数与阶数一致, 如二阶公式需调用二次 $f(x,y)$ , 而四阶公式需调用四次 $f(x,y)$ , 但对于更高阶的情形则不然, 需要调用 $f(x,y)$ 的次数, 远远大于方法的阶数. 由于计算量较大, 所以很少使用更高阶的龙格-库塔公式. 事实上, 对于大量的实际问题, 四阶的龙格-库塔公式已可满足对精度的要求.
以下是标准四阶四段龙格-库塔公式的算法描述.
算法 9.2.1(标准四阶四段龙格-库塔方法)
(1) 给定整数 N;
(2) $h = \frac{b - a}{N}$ ;
(3) $y(1) = y_0$ ;
(4) 对 $k = 2, 3, \cdots, N + 1$ 做:
$$ \begin{array}{r l} & {x = a + (k - 1) h;} \\ & {k_{1} = h f (x, y (k - 1));} \\ & {k_{2} = h f (x + 0.5 h, y (k - 1) + 0.5 k_{1});} \\ & {k_{3} = h f (x + 0.5 h, y (k - 1) + 0.5 k_{2});} \\ & {k_{4} = h f (x + h, y (k - 1) + k_{3});} \\ & {y (k) = y (k - 1) + (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}) / 6;} \end{array} $$(5) 返回 $y$ .
例 9.2.2 用标准四阶龙格-库塔方法 (9.46) 求解例 9.1.2 中的初值问题.
解：已知 $f(x,y) = y - \frac{2x}{y}, y(0) = 1$ ，并取 $h = 0.1$ . 将 $f, h$ 代入标准的四阶四段龙格-库塔格式 (9.46)，得具体计算公式为
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + \frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), \\ \quad k_{1} = 0.1 \left(y_{n} - \frac{2 x_{n}}{y_{n}}\right), \\ \quad k_{2} = 0.1 \left(y_{n} + \frac{k_{1}}{2} - \frac{4 (x_{n} + 0.5 h)}{2 y_{n} + k_{1}}\right), \\ \quad k_{3} = 0.1 \left(y_{n} + \frac{k_{2}}{2} - \frac{4 (x_{n} + 0.5 h)}{2 y_{n} + k_{2}}\right), \\ \quad k_{4} = 0.1 \left(y_{n} + k_{3} - \frac{2 (x_{n} + h)}{y_{n} + k_{3}}\right). \end{array} \right.\tag{9.48} $$计算结果见表 9-5.
表 9-5 标准四阶龙格－库塔方法的计算结果
| $x_k$ | $y_k$ | $y(x_k)$ | $|y_k - y(x_k)|$ |
| 0.0 | 1.000 000 00 | 1.000 000 00 | 0.000 0e-5 |
| 0.1 | 1.095 445 53 | 1.095 445 12 | 0.041 7e-5 |
| 0.2 | 1.183 216 74 | 1.183 215 96 | 0.078 9e-5 |
| 0.3 | 1.264 912 23 | 1.264 911 06 | 0.116 4e-5 |
| 0.4 | 1.341 642 54 | 1.341 640 79 | 0.156 7e-5 |
| 0.5 | 1.414 215 58 | 1.414 213 56 | 0.201 6e-5 |
| 0.6 | 1.483 242 22 | 1.483 239 69 | 0.252 5e-5 |
| 0.7 | 1.591 964 52 | 1.549 193 34 | 0.311 4e-5 |
| 0.8 | 1.612 455 35 | 1.612 451 55 | 0.380 0e-5 |
| 0.9 | 1.673 324 66 | 1.673 320 05 | 0.460 6e-5 |
| 1.0 | 1.732 056 37 | 1.732 050 81 | 0.555 8e-5 |
与例 9.1.2、例 9.1.4、例 9.1.8 的结果相比较, 可以看出, 四阶龙格-库塔公式的计算精
度比欧拉方法、梯形公式、改进的欧拉公式都高.
龙格-库塔公式在求解范围较大而精度要求较高时是比较好的方法, 它与欧拉公式和改进的欧拉公式一样都是显格式且可以自启动. 为了达到同样的精度, 它所用的步长可比欧拉公式和改进的欧拉公式所用的步长大得多. 值得指出的是, 由于龙格-库塔公式的推导基于泰勒展开方法, 因而它要求所求的解具有较好的光滑性质. 反之, 如果解的光滑性差, 那么, 使用四阶龙格-库塔公式求得的数值解, 其精度可能反而不如改进的欧拉公式. 在实际计算时, 我们应当针对问题的具体特点选择合适的算法.
9.3 收敛性与稳定性#
收敛性与稳定性从不同角度描述了微分方程数值解法的有效性．只有既收敛又稳定的方法,才能提供比较可靠的计算结果.以下将分别讨论单步法的收敛性与稳定性.
9.3.1 单步法的收敛性#
一般来说, 显式单步法总可写成如下形式:
$$ y_{n + 1} = y_{n} + h g (x_{n}, y_{n}, h),\tag{9.49} $$其中 $g(x_{n}, y_{n}, h)$ 称为增量函数, 它依赖于 f 且仅是关于 $x_{n}, y_{n}$ 和 h 的连续函数.
不同的单步法, 对应于不同的增量函数. 例如欧拉公式的增量函数为
$$ g (x, y, h) = f (x, y),\tag{9.50} $$改进的欧拉公式的增量函数为
$$ g (x, y, h) = \frac{1}{2} (f (x, y) + f (x + h, y + h f (x, y))).\tag{9.51} $$讨论数值方法的收敛性时总是假设计算过程中没有舍入误差, 因此只与方法的截断误差有关. 下面的定理给出了一个判断 p 阶方法收敛性的方法, 并指明了局部截断误差与收敛性的关系.
若 p 阶单步法 $y_{n+1}=y_{n}+hg(x_{n},y_{n},h)$ 中的增量函数 $g(x,y,h)$ 在区域 $a\leqslant x\leqslant b,$ $-\infty 其中 L > 0 为正常数, 则由该单步法得到的 $y_{n+1}$ 满足如下误差估计:
证：令 $\varepsilon_0 = y(x_0) - y_0 = 0$ ，则有
$$ \begin{array}{r l} & {e_{n + 1} = y (x_{n + 1}) - y_{n + 1}} \\ & {\qquad = \varepsilon_{n + 1} + y (x_{n}) + h g (x_{n}, y (x_{n}), h) - y_{n} - h g (x_{n}, y_{n}, h)} \\ & {\qquad = \varepsilon_{n + 1} + e_{n} + h (g (x_{n}, y (x_{n}), h) - g (x_{n}, y_{n}, h)),} \end{array}\tag{9.54} $$其中 $\varepsilon_{n+1}$ 表示局部截断误差, 即 $\varepsilon_{n+1}=y(x_{n+1})-(y(x_{n})+hg(x_{n},y(x_{n}),h))$ .
利用三角不等式, 可得 p 阶单步法的整体截断误差满足如下估计式:
$$ \begin{array}{l} | e_{n + 1} | \leqslant | \varepsilon_{n + 1} | + | e_{n} | + h L | e_{n} | \\ \quad \leqslant (1 + h L) | e_{n} | + C h^{p + 1} \\ \quad \leqslant (1 + h L)^{2} | e_{n - 1} | + ((1 + h L) + 1) C h^{p + 1} \\ \quad \leqslant (1 + h L)^{n + 1} | e_{0} | + ((1 + h L)^{n} + (1 + h L)^{n - 1} + \dots + (1 + h L) + 1) C h^{p + 1} \\ \quad = \frac{(1 + h L)^{n + 1} - 1}{hL} C h^{p + 1}. \end{array}\tag{9.55} $$注意到
$$ (1 + h L)^{n + 1} \leqslant (1 + h L)^{\frac{(b - a) L}{hL}} = \left((1 + h L)^{\frac{1}{hL}}\right)^{(b - a) L} \leqslant \mathrm{e}^{(b - a) L}.\tag{9.56} $$故 $|e_{n + 1}|\leqslant \frac{Ch^p}{L} (\mathrm{e}^{(b - a)L} - 1).$ 定理得证.
由定理 9.3.2 可知, 当 $p \geqslant 1$ 时, 单步法是收敛的, 且当局部截断误差 $\varepsilon_{n} = O(h^{p+1})$ 时, 其整体截断误差 $e_{n} = O(h^{p})$ . 由于整体截断误差比局部截断误差低一阶, 因此在构造高精度的计算方法时, 只要设法提高方法的局部截断误差.
利用定理 9.3.2 可以证明前面讲过的各单步法都是收敛的. 例如, 对于欧拉方法, 其增量函数为 $g(x,y,h)=f(x,y)$ , 所以当初值问题的右端函数 $f(x,y)$ 满足利普希茨条件时, 它是收敛的. 对于标准四阶四段龙格-库塔公式, 若 $f(x,y)$ 满足利普希茨条件, 则龙格-库塔公式的增量函数可写成
$$ g (x, y, h) = \frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}).\tag{9.57} $$设利普希茨常数为 L. 因为 $k_{1}=f(x,y)$ ，所以
$$ \begin{array}{l} | k_{1} (x, y_{1}, h) - k_{1} (x, y_{2}, h) | \leqslant L | y_{1} - y_{2} |, \\ | k_{2} (x, y_{1}, h) - k_{2} (x, y_{2}, h) | \leqslant L \left| y_{1} - y_{2} + \frac{h}{2} k_{1} (x, y_{1}, h) - \frac{h}{2} k_{1} (x, y_{2}, h) \right| \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \leqslant L \left(1 + \frac{1}{2} h L\right) | y_{1} - y_{2} |, \\ | k_{3} (x, y_{1}, h) - k_{3} (x, y_{2}, h) | \leqslant L \left(1 + \frac{1}{2} h L + \frac{1}{4} (h L)^{2}\right) | y_{1} - y_{2} |, \\ | k_{4} (x, y_{1}, h) - k_{4} (x, y_{2}, h) | \leqslant L \left(1 + h L + \frac{1}{2} (h L)^{2} + \frac{1}{4} (h L)^{3}\right) | y_{1} - y_{2} |. \end{array}\tag{9.58} $$故
$$ \begin{array}{c} | g (x, y_{1}, h) - g (x, y_{2}, h) | \leqslant L \left(1 + \frac{1}{2} h L + \frac{1}{6} (h L)^{2} + \frac{1}{24} (h L)^{3}\right) | y_{1} - y_{2} | \\ \leqslant L_{g} | y_{1} - y_{2} |. \end{array}\tag{9.59} $$由定理 9.3.2 知, 标准四阶四段龙格-库塔方法收敛.
9.3.2 单步法的稳定性#
前面讨论数值解的收敛性, 总是认为数值方法的计算是准确的. 但实际计算过程中总会有误差, 例如计算机的舍入误差. 一般地, 这种误差伴随着计算的过程的继续会不断积累并传递下去. 所谓的数值解的稳定性问题, 就是指误差的积累是否受到控制的问题. 粗略地说, 如果计算结果对初始数据的误差及计算过程中的误差不敏感, 也就是舍入误差不增长, 就称相应的计算方法是稳定的, 否则就称为不稳定的.
数值稳定性有各种各样的定义, 下面我们只简单地介绍 “绝对稳定性” 的概念. 设某数值方法在节点 $x_{n}$ 处的数值解应为 $y_{n}$ ，但实际计算时因为有误差，算得的值为 $\tilde{y}_{n}$ . 记 $\delta_{n} = \tilde{y}_{n} - y_{n}$ ，并称之为第 n 步数值解的扰动.
若一种数值方法在节点 $x_{n}$ 处的数值解 $y_{n}$ 伴随有大小为 $\delta_{n}$ 的扰动, 由 $\delta_{n}$ 的传播而导致以后各节点上的值 $y_{t}(t > n)$ 产生的扰动值 $\delta_{t}(t > n)$ , 其绝对值均不超过 $|\delta_{n}|$ , 即
$$ | \delta_{t} | \leqslant | \delta_{n} |, t = n + 1, n + 2, \dots ,\tag{9.60} $$则称该数值方法是绝对稳定的.
数值稳定性的分析是相当复杂的, 且总跟微分方程的右端 $f(x,y)$ 有关, 所以一般性地讨论常微分方程初值问题的稳定性很困难. 事实上, 我们不可能也不必对每个不同的右端函数进行讨论. 一个约定俗成的方法是对如下的“模型方程”
$$ y^{\prime} = \lambda y, \quad (\mathrm{Re} \lambda < 0)\tag{9.61} $$来讨论数值方法的绝对稳定性. 这是因为若一种数值方法对如此简单的问题都不是绝对稳定的话, 那么就难以用此数值方法来求解一般的微分方程初值问题. 显然, 若某数值方法对“模型方程”是绝对稳定的, 则对一般的微分方程, 该数值方法不一定也是绝对稳定的.
在把某数值方法应用于模型方程时, 为了保证数值稳定, 步长要受到一定的限制. 步长被允许的范围就称为该方法的绝对稳定区域. 一般来说, 绝对稳定区域越大, 算法的稳定性越好.
首先讨论欧拉公式的稳定性. 将欧拉公式用于模型方程可得
$$ y_{n + 1} = (1 + \lambda h) y_{n} = (1 + z) y_{n},\tag{9.62} $$其中, $z = \lambda h$ . 设在计算 $y_{n}$ 时有扰动 $\delta_{n}$ , 由 $\delta_{n}$ 的传播造成 $y_{n+1}$ 产生扰动值 $\delta_{n+1}$ , 则
$$ \delta_{n + 1} = (1 + z) \delta_{n},\tag{9.63} $$故要欧拉公式稳定, 只要
$$ \vert 1 + z \vert \leqslant 1.\tag{9.64} $$于是, 欧拉方法的绝对稳定区域是复平面上以 $z = -1$ 为圆心、1 为半径的圆盘.
其次考虑后退欧拉公式的稳定性. 将后退欧拉公式用于模型方程可得
$$ y_{n + 1} = y_{n} + \lambda h y_{n + 1},\tag{9.65} $$即
$$ y_{n + 1} = \frac{1}{1 - z} y_{n}.\tag{9.66} $$于是有
$$ \delta_{n + 1} = \frac{1}{1 - z} \delta_{n}.\tag{9.67} $$为使后退欧拉公式稳定, 只要
$$ \frac{1}{| 1 - z |} \leqslant 1,\tag{9.68} $$即 $|z - 1| \geqslant 1$ . 故后退欧拉公式的绝对稳定区域是复平面上以 $z = 1$ 为圆心、1 为半径的圆外区域. 显然, 它比欧拉方法的绝对稳定区域要大很多.
用类似的方法分析标准四阶四段龙格-库塔公式, 可得其绝对稳定区域为
$$ \left| 1 + z + \frac{z^{2}}{2} + \frac{z^{3}}{6} + \frac{z^{4}}{24} \right| \leqslant 1.\tag{9.69} $$若 $\lambda < 0$ 为实数, 则由上式可近似得到 $-2.78 \leqslant z < 0$ . 另外, 当 $\lambda$ 是负实数时, 亚当斯-巴什福思公式的绝对稳定区域为 $-0.55 < z < 0$ , 亚当斯-莫尔顿公式的绝对稳定区域为 $-3 < z < 0$ . 一般而言, 隐式格式的稳定性要比显式格式的稳定性要好.
9.4 微分方程组和刚性问题#
本节将主要讨论有关常微分方程组的数值解法. 事实上, 前三节在未知函数个数 $m = 1$ 情况下得到的大部分结论都可平行地推广到 $m > 1$ 的情况.
例如, 考虑如下一阶常微分方程组
$$ \binom{u^{\prime}}{v^{\prime}} = \left( \begin{array}{c c} 0 & 1 \-x^{2} & - x \end{array} \right) \binom{u}{v} + \binom{0}{x + 1}, \quad x \in (0, 10),\tag{9.70} $$初值条件为
$$ \binom{u (0)}{v (0)} = \binom{a}{b}.\tag{9.71} $$将区间 [0,10] 进行 $N$ 等分, 记 $h = \frac{10}{N}$ . 将欧拉公式 (9.9) 应用于该问题, 则有
$$ \binom{u_{n + 1}}{v_{n + 1}} = \binom{u_{n}}{v_{n}} + h \left(\left( \begin{array}{c c} 0 & 1 \-x_{n}^{2} & - x_{n} \end{array} \right) \binom{u_{n}}{v_{n}} + \binom{0}{x_{n} + 1}\right)\tag{9.72} $$且
$$ \binom{u_{0}}{v_{0}} = \binom{a}{b}.\tag{9.73} $$下面讨论应用欧拉方法求解 $(9.70) \sim (9.71)$ 的绝对稳定性. 简单计算可知
$$ \binom{\delta_{u}^{n + 1}}{\delta_{v}^{n + 1}} = \left( \begin{array}{c c} 0 & 1 \-x^{2} & - x \end{array} \right) \binom{\delta_{u}^{n}}{\delta_{v}^{n}},\tag{9.74} $$其中 $(\delta_u^{(n + 1)},\delta_v^{(n + 1)})^{\mathrm{T}}$ 和 $(\delta_u^{(n)},\delta_v^{(n)})^{\mathrm{T}}$ 分别表示第 $n$ 步和第 $n + 1$ 步数值解的扰动.由于矩阵
$\left( \begin{array}{cc}0 & 1\-x^{2} & -x \end{array} \right)$ 的特征值为 $\lambda_{1,2} = \mathrm{e}^{\pm \frac{2\pi}{3}\mathrm{i}}x,$ 故可以证明当 $|h\lambda_i| < 1$ (或 $h < 0.1$ )时，
$$ \left\| \binom{\delta_{u}^{n + 1}}{\delta_{v}^{n + 1}} \right\| \leqslant \left\| \binom{\delta_{u}^{n}}{\delta_{v}^{n}} \right\|,\tag{9.75} $$即该方法是绝对稳定的.
类似地, 前文介绍的梯形公式、预估校正格式和龙格-库塔公式均可应用于一阶常微分方程组的求解. 只是在进行理论分析时, 需要将绝对值替换为向量范数.
对于高阶常微分方程初值问题
$$ \left\{\begin{array}{l} y^{(n)} = f (x, y (x), y^{\prime} (x), \dots , y^{(n - 1)} (x)), \qquad x \in [ a, b ], \\ y (a) = \tilde{y}_{0}, y^{\prime} (a) = \tilde{y}_{0}^{\prime}, \dots , y^{(n - 1)} (a) = \tilde{y}_{0}^{(n - 1)}, \end{array} \right.\tag{9.76} $$若记
$$ \pmb{u} = (y (x), y^{\prime} (x), \dots , y^{(n - 1)})^{\mathrm{T}}, \quad \pmb{y_{0}} = (\tilde{y}_{0}, \tilde{y}_{0}^{\prime}, \dots , \tilde{y}_{0}^{(n - 1)})^{\mathrm{T}},\tag{9.77} $$则可将高阶微分方程化为一阶微分方程组, 从而利用前面介绍的方法近似求解.
例如, 考虑初值问题
$$ \left\{\begin{array}{l l} y^{\prime \prime} = - x y^{\prime} - x^{2} y + x + 1, & x \in [ 0, 10 ], \\ y (0) = a, \quad y^{\prime} (0) = b. \end{array} \right.\tag{9.78} $$令 y = u 和 $y' = v$ ，则该初值问题转化成一阶常微分方程组 (9.70)～(9.71).
在本节的最后, 我们简单介绍一类在化学反应、电子网络和自动控制等领域经常会遇到病态方程组——刚性方程组.
常微分方程初值问题 (9.1) 称为刚性方程组, 若对 $x \in [a, b]$ , 成立下列各式:
(1) $\operatorname{Re}(\mu_{j}(t)) < 0, j = 1, \cdots, m;$
$$ s = \frac{\max_{1 \leqslant j \leqslant m} | \mathrm{Re} (\mu_{j}) |}{\min_{1 \leqslant j \leqslant m} | \mathrm{Re} (\mu_{j}) |} \gg 1. $$其中 $\mu_{j}(t)$ 是雅可比矩阵 $\left(\frac{\partial f_{i}}{\partial y_{j}}\right)_{m\times m}$ 的特征值, s 称为刚性比.
例如, 考虑常微分方程组
$$ \binom{u^{\prime}}{v^{\prime}} = \left( \begin{array}{c c} 9 & 24 \-24 & - 51 \end{array} \right) \binom{u}{v} + \binom{5 \cos x - \frac{1}{3} \sin x}{- 9 \cos x + \frac{1}{3} \sin x}, \quad x \in (0, 1),\tag{9.79} $$初值条件为
$$ \binom{u (0)}{v (0)} = \left( \begin{array}{c} \frac{4}{3} \\ \frac{2}{3} \end{array} \right).\tag{9.80} $$这就是一个刚性方程组, 其雅可比矩阵为 $\left(\begin{array}{cc}9 & 24 \-24 & -51\end{array}\right)$ ，刚性比 $s = \frac{|-39|}{|-3|} = 13$ 。该问题存在唯一解
$$ u = 2 \mathrm{e}^{- 3 x} - \mathrm{e}^{- 39 x} + \frac{1}{3} \cos x, \qquad v = - \mathrm{e}^{- 3 x} + 2 \mathrm{e}^{- 39 x} - \frac{1}{3} \cos x. $$求解刚性问题的困难之处在于：为保证算法的稳定性，必须将步长限制在较小的范围内。因此，若需要计算到某一个较长的区间，则需要迭代非常多的时间步。这导致计算量大，并且由于舍入误差的累计，结果也很不准确。
应用标准四阶龙格-库塔公式求解问题 $(9.79)\sim (9.80)$ .由于绝对稳定性的限制，步长 $h$ 须满足 $|\mu_i h|\leqslant 2.78$ ，其中 $\mu_{i}$ 为雅可比矩阵的特征值.这要求 $h < \frac{2.78}{39}\approx 0.0713,$ 下面的数值结果验证了这一分析.取 $h = 0.05$ 和 $h = 0.1$ ，用标准四阶龙格-库塔公式求解该初值问题的计算结果见表9-6.由表9-6可以看出，若步长稍大并且不满足稳定性条件，则得到完全错误的数值结果.
表 9-6 龙格－库塔方法的计算结果
| h=0.05 | h=0.1 | |||||
| x | u(x) | v(x) | un | vn | un | vn |
| 0.1 | 1.793 061 | -1.032 001 | 1.712 219 | -0.870 315 2 | -2.645 169 | 7.844 527 |
| 0.2 | 1.423 901 | -0.874 680 9 | 1.414 070 | -0.855 014 8 | -18.451 58 | 38.876 31 |
| 0.3 | 1.131 575 | -0.724 998 4 | 1.130 523 | 0.722 891 0 | -87.472 21 | 176.482 8 |
| 0.4 | 0.909 408 6 | -0.608 214 1 | 0.909 276 3 | -0.607 947 5 | -934.072 2 | 789.354 0 |
| 0.5 | 0.738 787 7 | -0.515 657 5 | 9.738 750 6 | -0.515 581 0 | -1760.016 | 3 520.00 |
| 0.6 | 0.605 709 4 | -0.440 410 8 | 0.605 683 3 | -0.440 355 8 | -7848.550 | 15 697.84 |
| 0.7 | 0.499 860 3 | -0.377 403 8 | 0.499 836 1 | -0.377 354 0 | -34989.63 | 69 979.87 |
| 0.8 | 0.413 671 4 | -0.322 953 5 | 0.413 649 0 | -0.322 907 8 | -155 979.4 | 311 959.5 |
| 0.9 | 0.341 614 3 | -0.274 408 8 | 0.341 593 9 | -0.274 367 3 | -695 332.0 | 139 066 4 |
| 1.0 | 0.279 674 8 | -0.229 887 7 | 0.279 656 8 | -0.229 851 1 | -309 967 1 | 619 935 2 |
对于刚性问题, 一个最自然的想法就是扩大数值方法的绝对稳定区域. 显然, 若某数值方法是 A- 稳定的, 则应用该方法时步长可随意选取, 不再受稳定性限制. 遗憾的是, 早已证明, 显式多步法和显式龙格-库塔法不可能是 A- 稳定的, A- 稳定的隐式线性多步法的阶不超过 2, 而梯形公式是二阶隐式线性多步法中精度最高的一个. 因此, 为构造更实用的求解刚性问题的数值方法, 需要放松稳定性的限制.
实际计算时, 常采用隐式或半隐式的龙格-库塔公式求解刚性方程组. 下面列出几个常用的计算格式, 可以证明这些格式都是 A-稳定的.
(1) 一段二阶隐式龙格-库塔方法:
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + h k_{1}, \\ k_{1} = f \left(x_{n} + \frac{h}{2}, y_{n} + \frac{h}{2} k_{1}\right). \end{array} \right.\tag{9.81} $$(2) 二段二阶隐式龙格-库塔方法:
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2}), \\ k_{1} = f (x_{n}, y_{n}), \\ k_{2} = f \left(x_{n} + h, y_{n} + \frac{h}{2} (k_{1} + k_{2})\right). \end{array} \right.\tag{9.82} $$(3) 二段四阶隐式龙格-库塔方法:
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + \frac{h}{2} (k_{1} + k_{2}), \\ k_{1} = f \left(x_{n} + \left(\frac{1}{2} + \frac{\sqrt{3}}{6}\right) h, y_{n} + \frac{h}{4} \left(k_{1} + \left(1 + \frac{2 \sqrt{3}}{3} k_{2}\right)\right), \right. \\ k_{2} = f \left(x_{n} + \left(\frac{1}{2} - \frac{\sqrt{3}}{6}\right) h, y_{n} + \frac{h}{4} \left(\left(1 - \frac{2 \sqrt{3}}{3}\right) k_{1} + k_{2}\right)\right). \end{array} \right.\tag{9.83} $$(4) 半隐式龙格-库塔方法:
$$ \left\{\begin{array}{l} y_{n + 1} = y_{n} + k_{2}, \\ \quad k_{1} = f (x_{n}, y_{n}) + \left(1 - \frac{\sqrt{2}}{2}\right) h^{2} \frac{\partial f}{\partial x} (x_{n}, y_{n}) + \left(1 - \frac{\sqrt{2}}{2}\right) h \frac{\partial f}{\partial y} (x_{n}, y_{n}) k_{1}, \\ \quad k_{2} = h f \left(x_{n} + \frac{\sqrt{2} - 1}{2} h, y_{n} + \frac{\sqrt{2} - 1}{2} k_{1}\right) + \left(1 - \frac{\sqrt{2}}{2}\right) h^{2} \frac{\partial f}{\partial x} (x_{n}, y_{n}) \\ \quad + \left(1 - \frac{\sqrt{2}}{2}\right) h \frac{\partial f}{\partial y} (x_{n}, y_{n}) k_{2}. \end{array} \right.\tag{9.84} $$9.5 有限差分法#
本节主要讨论求解微分方程边值问题的有限差分法. 下面以两点边值问题 (9.4) 为例, 简单介绍该方法的基本思想.
通常, 利用有限差分法离散微分方程包含两步: 第一步是将求解区域进行网格剖分; 第二步是将微分方程在节点处进行离散化. 建立差分格式的离散化方法有多种, 这里仅介绍以差商代替微商的方法.
对于内部节点 $x_{n}(n=1,2,\cdots,N-1)$ ，由泰勒展开公式得
$$ \frac{y (x_{n + 1}) - 2 y (x_{n}) + y (x_{n - 1})}{h^{2}} = \left[ \frac{\mathrm{d}^{2} y}{\mathrm{d} x^{2}} \right]_{n} + \frac{h^{2}}{12} \left[ \frac{\mathrm{d}^{4} y}{\mathrm{d} x^{4}} \right]_{n} + O (h^{3}).\tag{9.85} $$其中 $[ ]_n$ 表示方括号内的函数在点 $x_{n}$ 取值. 于是在 $x_{n}$ 点处可将方程 (9.4) 写成
$$ \;-\; \frac{y (x_{n + 1}) - 2 y (x_{n}) + y (x_{n - 1})}{h^{2}} + q (x_{n}) y (x_{n}) = f (x_{n}) + R_{n} (y),\tag{9.86} $$其中
$$ R_{n} (y) = - \frac{h^{2}}{12} \left[ \frac{\mathrm{d}^{4} y}{\mathrm{d} x^{4}} \right]_{n} + O (h^{3}).\tag{9.87} $$显然, 当 $h$ 充分小时, $R_{n}(y)$ 是的二阶无穷小量. 若舍去 $R_{n}(y)$ , 则得到求解方程的差分方程:
$$ \;-\; \frac{y_{n + 1} - 2 y_{n} + y_{n - 1}}{h^{2}} + q_{n} y_{n} = f_{n},\tag{9.88} $$式中 $q_{n}=q(x_{n})$ , $f_{n}=f(x_{n})$ . 称 $R_{n}(y)$ 为差分方程 (9.88) 的截断误差.
对于边界节点, 由边界条件知 $y_{0} = \alpha, y_{N} = \beta$ .
于是可得到关于 $y_{n}$ 的线性代数方程组：
$$ \left\{\begin{array}{l} - \frac{y_{n + 1} - 2 y_{n} + y_{n - 1}}{h^{2}} + q_{n} y_{n} = f_{n}, \quad n = 1, \dots , N - 1, \\ y_{0} = \alpha , \quad y_{N} = \beta . \end{array} \right.\tag{9.89} $$若记方程组的未知向量 $\boldsymbol{y}_{h} = (y_{1}, y_{2}, \cdots, y_{N-1})^{\mathrm{T}}$ ，右端向量 $\boldsymbol{g} = \left( f_{1} + \frac{\alpha}{h^{2}}, f_{2}, \cdots, f_{N-1} + \frac{\beta}{h^{2}} \right)^{\mathrm{T}}$ ，系数矩阵
$$ \boldsymbol{H} = \left( \begin{array}{c c c c c} \frac{2}{h^{2}} + q_{1} & - \frac{1}{h^{2}} & 0 & \dots & 0 \-frac{1}{h^{2}} & \frac{2}{h^{2}} + q_{2} & - \frac{1}{h^{2}} & \dots & 0 \\ 0 & - \frac{1}{h^{2}} & \frac{2}{h^{2}} + q_{3} & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & \frac{2}{h^{2}} + q_{N - 1} \end{array} \right),\tag{9.90} $$则有
$$ \boldsymbol{H} \boldsymbol{y}_{h} = \boldsymbol{g}.\tag{9.91} $$易知, H 为对称正定矩阵, 故该方程组有唯一解. 此外, 由于矩阵 H 为三对角阵, 所以可以用稀疏数组的方式存储并用追赶法求解.
下面是求解两点边值问题的有限差分法的算法描述.
算法 9.5.1(有限差分法)
(1) 给定整数 N;
(2) $h = \frac{b - a}{N}$ ;
(3) 生成三对角阵 H 和向量 g;
(4) 求解方程组 $H y_{h} = g;$
(5) 返回向量 $y_{h}$ .
例 9.5.2 用有限差分法算法 9.5.1 求解常微分方程边值问题
$$ \left\{\begin{array}{l l} - y^{\prime \prime} + y = (\pi^{2} + 1) \sin (\pi x), \\ y (0) = 0, \quad y (1) = 0. \end{array} \right.\tag{9.92} $$解：该问题的真解为 $y(x)=\sin\pi x$ 。取 N=2,4,8,16, $\cdots$ ，则有限差分法计算结果见表 9-7，其中 $\|y_{h}-y\|_{\infty}=\max_{1\leqslant n\leqslant N}|y_{h}(x_{n})-y(x_{n})|$
表 9-7 不同 N 情况下有限差分解的误差
| N | 2 | 4 | 8 | 16 | 32 | 64 | 128 | 256 |
| $\|y_h - y\|_\infty$ | 0.207 7 | 0.047 9 | 0.011 7 | 0.002 9 | 7.295 9e-4 | 1.823 4e-4 | 4.558 2e-5 | 1.139 5e-5 |
从表 9-7 可以看出：随着等分节点个数的增加，有限差分方法得到的数值解与真解之间的误差不断减小。事实上，为控制计算误差在可接受的范围之内，可采用更多的网格节点。一般地，若有限差分法格式满足一定条件，则当步长 h 趋向于零时，由有限差分法得到的数值解收敛于真解。
有限差分方法还可以使用数值积分的方法进行离散. 更多内容参见偏微分方程数值解法的有关教材.
评 注#
本章主要介绍了几种求解常微分方程初边值问题的数值方法. 对于初边值问题, 离散常微分方程的常用方法主要有两种: 一种是基于泰勒展开, 另一种是基于数值积分. 这两种方法中, 利用泰勒展开来构造数值计算格式比较简单灵活, 但通常需要在真解比较光滑的前提下才能进行, 而数值积分方法对真解的光滑性要求较低. 此外, 这两种方法也可用于偏微分方程的数值离散, 如利用泰勒展开可构造有限差分法, 利用数值积分可构造有限体积法等.
单个方程的数值解法, 稍加改动都可平行地用于一阶方程组的情形. 构造常微分方程组的离散格式时, 需要将相应的单个方程计算公式中的实 (复) 函数替换成实 (复) 向量函数; 理论分析时, 需要将函数的绝对值换成向量范数. 高阶微分方程等价于一阶微分方程组, 故也可应用本章所介绍的数值方法进行计算.
对于刚性微分方程组, 一般的显格式计算效果很差, 故需要考虑 A-稳定、 $\mathrm{A}(\alpha)$ -稳定和刚性稳定的数值方法. 本章给出了几种常用的隐式或半隐式龙格-库塔公式, 用于求解此类问题.
对于初值问题, Matlab 软件提供了一系列的函数, 如 ode45 和 ode23 等, 其使用格式如下:
$$
[ t, y ] = \text{odesolver} (\text{odefun}, \text{tspan}, y 0) $$
其中 odefun 为函数名字符串, 表示 (9.1) 中的 $f(x,y)$ ; 向量 tspan=[a,b] 表示初始和结束时间; y0 表示初始条件. 返回值 t 和 y 分别表示区间 $[a,b]$ 上的节点以及处近似解. 可以直接用命令 plot(t,y) 画出近似解的折线图.
下面以例 9.1.2 中的初值问题为例介绍该函数的使用方法.
(1) 编写函数文件 f.m.
function r = f(x, y)
$$ r = y - 2 * x. / y; $$(2) 在命令窗口 (或脚本文件) 中, 执行如下代码.
$$
[ t, y ] = \text{ode45} (^{\prime} f^{\prime}, [ 0, 1 ], 1) $$
注意, odesolver 表示函数名, 实际使用时用 ode45 等替换. 这些函数也适用于求解常微分方程组.
习题九#
- 取 $h = 0.1$ , 分别用欧拉公式、梯形公式和改进的欧拉公式在 $0 \leqslant x \leqslant 1$ 上求解初值问题
- 取 h = 0.2, 用四阶龙格-库塔方法在 $0 \leqslant x \leqslant 1$ 上求解初值问题
- 证明对于任意参数 $\alpha$ , 下列格式是二阶的:
- 分析用中点公式求解初值问题
时绝对稳定性对步长的限制.
- 导出用梯形公式求解
的计算公式
$$ y_{n} = \left(\frac{2 + h}{2 - h}\right)^{n}. $$取 $h = \frac{1}{4}$ ，计算 $y(1)$ 的近似值.
- 取 $h = 0.5$ , 用有限差分方法在 $0 \leqslant x \leqslant 1$ 上求解边值问题
数值实验九#
- 分别用欧拉方法、梯形公式、改进的欧拉方法以及标准四阶龙格-库塔方法求解下列常微分方程初值问题
比较四种方法的计算精度, 并体会显式格式与隐式格式的区别.
- 用标准四阶龙格-库塔方法, 对 $x \geqslant 0$ 时的标准正态分布函数
产生一张在 $[0,5]$ 之间 80 个等距节点处的函数表.
- 考虑刚性问题
该问题的真解为 $y(x)=x-\mathrm{e}^{-5x}$ . 分别取步长 h=0.2 和 0.25, 用标准四阶龙格-库塔公式和梯形公式求解该初值问题, 并对计算结果进行分析.
- 尝试用不同方法求解下面的初值问题
初值条件为
$$ \binom{u (0)}{v (0)} = \left( \begin{array}{c} \frac{1}{3} \\ \frac{1}{3} \end{array} \right),\tag{9.94} $$比较各种方法的计算结果和计算时间, 并分析它们的精度. (该问题的精确解为 $u = \frac{2}{3}x + \frac{2}{3}e^{-x} - \frac{1}{3}e^{-100x}$ , $v = -\frac{1}{3}x - \frac{1}{3}e^{-x} + \frac{2}{3}e^{-100x}$ .)
- 分别取步长 $h = \frac{1}{2}, \frac{1}{4}, \cdots, \frac{1}{256}$ ，用有限差分方法求解边值问题
参考文献#
[1] Atkinson K.E., An Introduction to Numerical Analysis, John Wiley, New York, 1989.
[2] Richard L., Burden J., Douglas Faires, Numerical Analysis (Seventh Edition), 高等教育出版社, 2001.
[3] David, P.J. and Rabinowitz P., Methods of Numerical Integration, New York Academic Press, New York, 1975.
[4] Michael T. Heath, Scientific Computing: An Introductory Survey (Second Edition), 清华大学出版社, 2001.
[5] 胡健伟, 汤怀民, 微分方程数值方法 [M], 科学出版社, 1999.
[6] 贾现正, 王彦博, 程晋, 散乱数据的数值微分及其误差估计, 高等学校计算数学学报, vol 25, No.1, 2003.
[7] Kalos M. H., Whitlock P.A., Monte Carlo methods, John Wiley and Sons, New York, 1986.
[8] 李立康, 於崇华, 朱政华, 微分方程数值解 [M], 复旦大学出版社, 1999.
[9] Quarteroni A., Sacco R., Saleri F., Numerical Mathematics, Springer, Science and Business Media, Inc, Berlin, 2000.
[10] 同济大学计算数学教研室, 数值分析基础 [M], 同济大学出版社, 1998.
[11] 王仁宏. 数值逼近 [M]. 北京：高等教育出版社, 2001.
[12] 肖云茹. 概率统计计算方法 [M]. 南开大学出版社, 1994.
[13] 张平文, 李铁军. 数值分析 [M]. 北京大学出版社, 2007.
[14] 赵访熊, 李庆扬, 富利叶变换滤波在地震勘探处理中应用, 清华大学学报, 1978, 18(4).
[15] Jie Shen, Tao Tang, Spectral and High-Order Methods with Applications, Science Press, Beijing, 2006.
埃尔米特插值, 66 埃尔米特多项式, 85 艾特肯加速方法, 187
贝塞尔不等式, 103 边界条件, 71 变步长方法, 118 病态, 155 不动点, 183 不可约, 160 步长, 218
残差, 102 插值多项式, 56 插值基函数, 58 插值余项, 60 差分, 65 差商, 62 乘幂法, 208 传值, 21
代数精度, 111 单步法, 222 递归函数, 22 杜利脱尔分解, 40 对角占优, 160 对阶, 6 多步法, 222 多项式拟合, 96
二次插值, 59