高级数值分析 第 3 章 常微分方程数值解法

本文由手写笔记《高级数值分析》扫描件的 LaTeX 转录稿改写而来，公式以红色保留原稿的红笔重点标记。

3.1 初值问题与离散化

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

$$y'=f(x,y),\qquad x_0\le x\le b,\qquad y(x_0)=y_0.$$

其积分形式为

$$y(x)=y_0+\int_{x_0}^{x}f\bigl(t,y(t)\bigr)\,\mathrm dt.$$

取 $N+1$ 个节点

$$a\le x_0<x_1<\cdots <x_N\le b, \qquad h_n=x_{n+1}-x_n,$$

分别记精确解为 $y(x_0),\ldots,y(x_N)$，数值解为 $y_0,\ldots,y_N$。数值格式一般写成

$$y_{n+1}=g\bigl(y_n,y_{n-1},\ldots,y_{n-k+1},x_n,h_n\bigr).$$

只含相邻一步未知量的格式称为单步法；含有前面 $k$ 个网格值的格式称为 $k$ 步法。若右端只含已知量，则为显式格式，否则为隐式格式。

3.1.1 局部截断误差与阶

将精确解代入数值格式，定义局部截断误差

$$T_{n+1}(h)=y(x_{n+1})- g\bigl(y(x_n),y(x_{n-1}),\ldots,y(x_{n-k+1}),x_n,h_n\bigr).$$

若

$$T_{n+1}(h)=O(h^{p+1}),$$

则称该方法具有 $p$ 阶局部精度（或 $p$ 阶）。

3.2 一致性、稳定性与收敛性

设一步格式写成

$$y_{n+1}=y_n+h\,\Phi(x_n,y_n,h),$$

其中 $\Phi$ 称为增量函数。若

$$\Phi(x,y,0)=f(x,y),$$

则称该格式与微分方程相容（或一致）。等价地，一致性可表述为

$$\lim_{h\to0}\max_{0\le n\le N} \left|\frac{T_{n+1}(h)}{h}\right|=0.$$

若存在与 $h$ 无关的常数 $K>0$ 和 $h_0>0$，使任意两组初值所产生的数值解满足

$$|y_n-\widetilde y_n| \le K\max_{0\le j\le n-1}|y_j-\widetilde y_j|, \qquad 0<h\le h_0,$$

则称该格式稳定。收敛性定义为

$$\lim_{h\to0}\max_{0\le n\le N} |y_n-y(x_n)|=0.$$

在通常的 Lipschitz 条件下，一致性与稳定性蕴含收敛性；若局部截断误差为 $O(h^{p+1})$，则全局误差为 $O(h^p)$。

3.3 常用单步方法

由积分形式对 $[x_n,x_{n+1}]$ 作不同求积近似，可得下列格式。

3.3.1 Euler 方法

向前 Euler 法（显式 Euler 法）：

$$y_{n+1}=y_n+h f(x_n,y_n).$$

向后 Euler 法（隐式 Euler 法）：

$$y_{n+1}=y_n+h f(x_{n+1},y_{n+1}).$$

梯形法：

$$y_{n+1}=y_n+\frac h2\bigl[f(x_n,y_n)+f(x_{n+1},y_{n+1})\bigr].$$

*改进 Euler 法（Heun 形式）*先作预测

$$\widetilde y_{n+1}=y_n+h f(x_n,y_n),$$

再校正为

$$y_{n+1}=y_n+\frac h2\bigl[f(x_n,y_n)+f(x_{n+1},\widetilde y_{n+1})\bigr].$$

3.4 绝对稳定性

考察测试方程

$$y'=\lambda y,\qquad \lambda\in\mathbb C,\quad \operatorname{Re}\lambda<0, \qquad y(0)=y_0.$$

将数值格式代入后得到差分方程。对一步法可写成

$$y_{n+1}=R(z)y_n,\qquad z=h\lambda,$$

其中 $R(z)$ 为稳定函数；对多步法则得到关于增长因子 $\eta$ 的特征方程。若所有特征根满足 $|\eta|<1$（单位圆上的根至多为单根），则称该格式在 $z=h\lambda$ 处绝对稳定。满足该条件的 $z$ 的集合称为绝对稳定域（边界由 $|R(z)|=1$ 给出）。

3.5 积分型方法

更一般地，可把一步格式写成局部积分公式

$$y_{n+1}=y_n+h\,\Phi(x_n,y_n,h),$$

其中 $\Phi$ 由对积分

$$\int_{x_n}^{x_{n+1}} f\bigl(x,y(x)\bigr)\,\mathrm dx$$

的数值求积得到。例如左端点、右端点和梯形求积分别给出向前 Euler、向后 Euler 和梯形法。若求积公式对光滑函数的误差为 $O(h^{p+1})$，则相应微分方程格式具有 $p$ 阶局部精度。

3.6 Runge–Kutta 方法

对积分形式

$$y(x_n+h)-y(x_n)=\int_{x_n}^{x_n+h}f\bigl(x,y(x)\bigr)\,\mathrm dx$$

作求积近似，可得 Runge–Kutta（简称 RK）方法。其一般的显式 $m$ 阶段形式为

$$\begin{aligned} y_{n+1}&=y_n+h\sum_{i=0}^{m-1}a_iK_i,\\ K_0&=f(x_n,y_n),\\ K_i&=f\left(x_n+c_i h, y_n+h\sum_{j=0}^{i-1}b_{ij}K_j\right), \qquad i=1,\ldots,m-1. \end{aligned}$$

一致性要求

$$\sum_{i=0}^{m-1}a_i=1, \qquad c_i=\sum_{j=0}^{i-1}b_{ij}.$$

例如二阶显式 RK（中点格式）为

$$\begin{aligned} y_{n+1}&=y_n+hK_1,\\ K_0&=f(x_n,y_n),\\ K_1&=f\left(x_n+\frac h2,y_n+\frac h2K_0\right). \end{aligned}$$

RK 方法的系数常用 Butcher 表表示：

$$\begin{array}{c|cccc} c_0&b_{00}\\ c_1&b_{10}&b_{11}\\ \vdots&\vdots&\vdots&\ddots\\ c_{m-1}&b_{m-1,0}&b_{m-1,1}&\cdots&b_{m-1,m-1}\\ \hline &a_0&a_1&\cdots&a_{m-1} \end{array}$$

（显式方法的对角线及其上方系数为零）。经典四阶 RK 方法的 Butcher 表为

$$\begin{array}{c|cccc} 0 &0\\ \frac12&\frac12\\ \frac12&0&\frac12\\ 1&0&0&1\\ \hline &\frac16&\frac13&\frac13&\frac16 \end{array}$$

相应计算公式为

$$\begin{aligned} y_{n+1}&=y_n+\frac h6\left(K_0+2K_1+2K_2+K_3\right),\\ K_0&=f(x_n,y_n),\\ K_1&=f\left(x_n+\frac h2,y_n+\frac h2K_0\right),\\ K_2&=f\left(x_n+\frac h2,y_n+\frac h2K_1\right),\\ K_3&=f(x_n+h,y_n+hK_2). \end{aligned}$$

经典四阶 RK 方法具有四阶精度，并且相容、稳定且收敛。

3.7 线性多步法

3.7.1 一般形式

线性 $k$ 步法的一般形式定义为

$$y_{n+k}=\sum_{i=0}^{k-1}\alpha_i y_{n+i} +h\sum_{i=0}^{k}\beta_i f_{n+i}, \qquad f_{n+i}=f(x_{n+i},y_{n+i}).$$

其中 $\alpha_i,\beta_i$ 为常数。若 $\beta_k=0$，则格式为显式；若 $\beta_k\ne0$，则为隐式。使用 $k$ 步法时，需要先由单步法求出 $k-1$ 个初始值。

将上式写成算子形式，定义第一特征多项式与第二特征多项式

$$\rho(\xi)=\xi^k-\sum_{i=0}^{k-1}\alpha_i\xi^i, \qquad \sigma(\xi)=\sum_{i=0}^{k}\beta_i\xi^i.$$

则线性多步法可记为

$$\rho(E)y_n=h\,\sigma(E)f_n,$$

其中 $E$ 为移位算子。

3.7.2 Adams 格式与 BDF 格式

Adams 格式只使用最近的函数值，常写为

$$y_{n+k}=y_{n+k-1}+h\sum_{i=0}^{k}\beta_i f_{n+i}.$$

当右端不含 $f_{n+k}$ 时得到 Adams–Bashforth 显式格式；含有 $f_{n+k}$ 时得到 Adams–Moulton 隐式格式。预测–校正方法（PCE）通常先用 Adams–Bashforth 预测，再用 Adams–Moulton 校正。Milne–Simpson 与 Hamming 方法也是常见的预测–校正格式，但 Milne 与 Simpson 格式的稳定性相对较差。

后向差分格式（backward differentiation formula，BDF）可写为

$$y_{n+k}=\sum_{i=0}^{k-1}\alpha_i y_{n+i}+h\beta_k f_{n+k}.$$

常用的 Gear–BDF 系数如下：

$$\begin{array}{c|ccccc} \text{格式}&\alpha_0&\alpha_1&\alpha_2&\alpha_3&\beta_k\\ \hline \mathrm{BDF1}&1&-1&&&1\\ \mathrm{BDF2}&\frac43&-\frac13&&&\frac23\\ \mathrm{BDF3}&\frac{18}{11}&-\frac9{11}&\frac2{11}&&\frac6{11}\\ \mathrm{BDF4}&\frac{48}{25}&-\frac{36}{25}&\frac{16}{25}&-\frac3{25}&\frac{12}{25} \end{array}$$

3.7.3 多步法的性能与 Dahlquist 等价定理

对线性多步法，定义

$$\rho(\xi)=\xi^k-\sum_{i=0}^{k-1}\alpha_i\xi^i, \qquad \sigma(\xi)=\sum_{i=0}^{k}\beta_i\xi^i.$$

由定义，格式一致（相容）的充要条件为

$$\rho(1)=0, \qquad \rho'(1)=\sigma(1).$$

零稳定性由根条件刻画：若 $\rho(\xi)=0$ 的全部根满足

$$|\xi|\le1,$$

且单位圆上的根均为单根，则称该多步法零稳定；否则称其不零稳定。根条件保证舍入误差不会随步数无限放大。

定理 3.1（Dahlquist 等价定理）

在线性多步法的一般假设下，方法收敛当且仅当它一致且零稳定。特别地，一致性与零稳定性共同保证

$$\max_{0\le n\le N}|y_n-y(x_n)|\longrightarrow0 \qquad(h\to0).$$

若方法满足根条件，则其对初值扰动和局部误差稳定；若再满足一致性，则由 Dahlquist 等价定理可知方法收敛。通常局部截断误差为 $O(h^{p+1})$ 时，全局误差为 $O(h^p)$。

3.8 绝对稳定性与刚性问题

对测试方程

$$y'=\lambda y,$$

将线性多步格式代入，并记 $\bar h=h\lambda$，得到

$$y_{n+k}=\sum_{i=0}^{k-1}\alpha_i y_{n+i} +\bar h\sum_{i=0}^{k}\beta_i y_{n+i}.$$

若设 $y_n=\eta^n$，则增长因子满足特征方程

$$\rho(\eta)-\bar h\,\sigma(\eta)=0, \qquad \rho(\eta)=\eta^k-\sum_{i=0}^{k-1}\alpha_i\eta^i, \quad \sigma(\eta)=\sum_{i=0}^{k}\beta_i\eta^i.$$

当给定 $\bar h\in\mathbb C$ 时，若该方程的全部根满足

$$|\eta_i(\bar h)|\le 1,$$

且单位圆上的根为单根，则称格式在 $\bar h$ 处绝对稳定；所有这类 $\bar h$ 构成绝对稳定域。

以显式 Euler 法为例，格式化为

$$y_{n+1}=(1+\bar h)y_n,$$

故其绝对稳定条件为

$$|1+\bar h|\le1.$$

这说明当 $\lambda$ 为负实数时，步长必须满足 $-2\le h\lambda\le0$；复数情形则对应复平面上以 $-1$ 为圆心、半径为 $1$ 的圆盘。

3.9 常微分方程组

一阶常微分方程组可写成

$$\begin{cases} y_1'=f_1(x,y_1,\ldots,y_m),\\ y_2'=f_2(x,y_1,\ldots,y_m),\\ \ \vdots\\ y_m'=f_m(x,y_1,\ldots,y_m), \end{cases} \qquad\Longleftrightarrow\qquad \begin{cases} \boldsymbol y'=\boldsymbol f(x,\boldsymbol y),\\ \boldsymbol y(x_0)=\boldsymbol y_0. \end{cases}$$

特别地，$m$ 阶方程

$$y^{(m)}=f\bigl(x,y,y',\ldots,y^{(m-1)}\bigr)$$

可令

$$y_1=y,\quad y_2=y',\quad\ldots,\quad y_m=y^{(m-1)},$$

化为

$$\begin{cases} y_1'=y_2,\\ y_2'=y_3,\\ \ \vdots\\ y_{m-1}'=y_m,\\ y_m'=f(x,y_1,\ldots,y_m), \end{cases}$$

初值为

$$y_1(x_0)=y_0,\quad y_2(x_0)=y_0',\quad\ldots,\quad y_m(x_0)=y_0^{(m-1)}.$$

3.10 刚性常微分方程

对方程组 $\boldsymbol y'=\boldsymbol f(x,\boldsymbol y)$，记 Jacobi 矩阵

$$J=\frac{\partial\boldsymbol f}{\partial\boldsymbol y}.$$

若 $J$ 的全部特征值满足 $\operatorname{Re}\lambda_i<0$，而且不同特征值的衰减尺度相差很大，则称该问题为刚性 ODE。通常用刚性指标

$$S=\frac{\max_i|\operatorname{Re}\lambda_i|} {\min_i|\operatorname{Re}\lambda_i|}$$

衡量刚性程度。刚性问题的稳定性往往由最大的负特征值限制步长，而解的长期行为又由最小的衰减率决定，因此显式方法可能要求远小于精度所需的步长；此时宜采用隐式 Euler、BDF 等具有较大稳定域的方法。

3.11 积分方程形式与配点法

初值问题的积分形式为

$$y(x)=y_0+\int_{x_0}^{x}f\bigl(t,y(t)\bigr)\,\mathrm dt.$$

取基函数 $\phi_0,\ldots,\phi_m$，用

$$y_m(x)=\sum_{j=0}^{m}a_j\phi_j(x)$$

近似 $y(x)$，定义残量

$$R(x;\boldsymbol a)=y_m(x)-y_0- \int_{x_0}^{x}f\bigl(t,y_m(t)\bigr)\,\mathrm dt.$$

在配置点 $x_k$ 上令 $R(x_k;\boldsymbol a)=0$，得到关于系数 $\boldsymbol a$ 的非线性方程组。展开写成

$$R(x_k;\boldsymbol a)= \sum_{j=0}^{m}a_j\phi_j(x_k)-y_0- \int_{x_0}^{x_k}f\!\left(t,\sum_{j=0}^{m}a_j\phi_j(t)\right)\mathrm dt=0, \qquad k=0,\ldots,m.$$

该方法称为配点法；积分可用数值求积计算。

3.12 Galerkin 方法

Galerkin 方法取试验函数 $p_0,\ldots,p_m$，要求残量对每个试验函数正交：

$$\int_{x_0}^{b}p_i(x)R(x;\boldsymbol a)\,\mathrm dx=0, \qquad i=0,\ldots,m.$$

代入残量得到

$$\sum_{j=0}^{m}a_j\int_{x_0}^{b}p_i(x)\phi_j(x)\,\mathrm dx -y_0\int_{x_0}^{b}p_i(x)\,\mathrm dx -\int_{x_0}^{b}p_i(x)\int_{x_0}^{x} f\!\left(t,\sum_{j=0}^{m}a_j\phi_j(t)\right)\mathrm dt\,\mathrm dx=0.$$

记

$$G_{ij}=\int_{x_0}^{b}p_i(x)\phi_j(x)\,\mathrm dx, \qquad s_i=\int_{x_0}^{b}p_i(x)\,\mathrm dx,$$

则在 $f$ 线性或经线性化后可写成矩阵方程

$$G\boldsymbol a=\boldsymbol s;$$

一般非线性情形则形成关于 $\boldsymbol a$ 的非线性方程组，可用 Newton 法等求解。

3.12.0.0.1 残量最小化

还可直接令

$$\min_{\boldsymbol a}\lVert R\rVert_2^2 =\min_{\boldsymbol a}\int_{x_0}^{b}|R(x;\boldsymbol a)|^2\,\mathrm dx, \qquad \frac{\partial}{\partial a_i}\lVert R\rVert_2^2=0,$$

这给出与 Galerkin 正交条件等价的正规方程。对于向量方程组，以上过程逐分量进行，或直接以向量残量

$$\boldsymbol R(x)=\boldsymbol y_m'(x)-\boldsymbol f\bigl(x,\boldsymbol y_m(x)\bigr)$$

构造相应的配点、Galerkin 或最小二乘离散系统。