数值分析讲义（三）：常微分方程初值问题与刚性 Part II
建议先阅读 数值分析讲义（三）：常微分方程初值问题与刚性 Part I。本篇继续整理常微分方程（ordinary differential equation, ODE）初值问题（initial value problem, IVP）中的刚性微分方程、稳定区域、A 稳定（A-stability）和 L 稳定（L-stability）。
3.2 刚性微分方程
在许多应用中，例如化学反应过程，以及偏微分方程的半离散化中，会出现刚性系统。虽然它们同样也是初值问题，但对于许多方法（不是所有方法），为了得到准确解，它们会迫使步长 $h$ 小到难以接受。
从一个由 $n$ 个常微分方程组成的初值问题出发。相应地，记 IVP$_n$ 为 $n$ 维 initial value problem：
\[\text{(IVP}_n\text{)}\qquad y'(t)=f(t,y(t)),\qquad t\in[a,b],\] \[y(a)=y_0,\]其中
\[f:[a,b]\times\mathbb{R}^n\to\mathbb{R}^n, \qquad y_0\in\mathbb{R}^n.\]文献中对“刚性系统”这一概念的定义并不完全统一。直观地说，刚性问题的解往往同时包含两个时间尺度：一部分变化较慢，另一部分会在很短时间内迅速衰减。
线性情形记为 LIVP$_n$（linear initial value problem）：
\[\text{(LIVP}_n\text{)}\qquad y'(t)=Ay(t)+c,\qquad t\in[a,b],\] \[y(a)=y_0,\]其中矩阵 $A\in\mathbb{R}^{n\times n}$，向量 $c\in\mathbb{R}^n$。
此外，设 $A\in\mathbb{R}^{n\times n}$ 可对角化，并有相应特征值 $\lambda_i$ 和特征向量 $v_i$。若 $y_P$ 是一个特解，则通解具有形式
\[y(t)=y_H(t)+y_P(t), \qquad y_H(t)=\sum_{i=1}^{n}C_i e^{\lambda_i t}v_i.\]若实部小于0
\[\operatorname{Re}(\lambda_i)<0,\qquad i=1,\ldots,n,\]则
\[\lim_{t\to\infty} y_H(t)\to 0,\]因此所有解都趋近于 $y_P$。其中，$y_H$ 中满足 $\operatorname{Re}(\lambda_i)\ll -1$ 的项衰减非常快，而满足 $\operatorname{Re}(\lambda_i)\not\ll -1$ 的项衰减明显较慢。也就是说，如果系统既有实部远小于 $0$ 的特征值，又有实部接近 $0$ 但仍为负的特征值，就会同时出现快、慢两个尺度；这类系统称为刚性系统，见定义 3.2.2。
刚性问题中，快衰减分量很快从解里消失，但它对应的特征值实部远小于 $0$，仍会限制显式方法可选的稳定步长。
例 3.2.1： 考虑问题
\[y'=Ay, \qquad y(0)=y_0:= \begin{pmatrix} C_1+C_2\\ C_1-C_2 \end{pmatrix},\]其中 $C_1,C_2\in\mathbb{R}$，并且
\[A= \begin{pmatrix} \frac{\lambda_1+\lambda_2}{2} & \frac{\lambda_1-\lambda_2}{2}\\ \frac{\lambda_1-\lambda_2}{2} & \frac{\lambda_1+\lambda_2}{2} \end{pmatrix}.\]$A$ 的特征值为 $\lambda_1,\lambda_2$，对应特征向量分别为
\[\begin{pmatrix}1\\1\end{pmatrix} \qquad\text{和}\qquad \begin{pmatrix}1\\-1\end{pmatrix}.\]例如当
\[\lambda_1=-1,\qquad \lambda_2=-1000\]时，解为
\[y(t) =C_1\begin{pmatrix}1\\1\end{pmatrix}e^{-t} +C_2\begin{pmatrix}1\\-1\end{pmatrix}e^{-1000t}.\]第二项在极短时间后几乎不再起作用。第一项起主导作用，并且在 $t\to\infty$ 时也趋于 0。对于一个合适的积分方法，人们期望它在不对步长作强限制的情况下给出满足
\[\lim_{j\to\infty}u_j=0\]的近似值 $u_j$。
然而，例如，如果使用显式 Euler 方法，那么由
\[u_0=y_0 =C_1\begin{pmatrix}1\\1\end{pmatrix} +C_2\begin{pmatrix}1\\-1\end{pmatrix}\]得到
\[u_1=(I+hA)u_0 =C_1(1+h\lambda_1)\begin{pmatrix}1\\1\end{pmatrix} +C_2(1+h\lambda_2)\begin{pmatrix}1\\-1\end{pmatrix},\]进而归纳得到
\[u_j =C_1(1+h\lambda_1)^j\begin{pmatrix}1\\1\end{pmatrix} +C_2(1+h\lambda_2)^j\begin{pmatrix}1\\-1\end{pmatrix}.\]若 $C_2\ne0$，则必须选择
\[|1+h\lambda_2|<1,\]也就是
\[-h\lambda_2=1000h<2,\]才能保证 $\lim_{j\to\infty}u_j=0$。合适的方法应尽可能对所有 $h>0$ 都确保这一点。
对多模态系统，步长必须同时让所有 $q_i=\lambda_i h$ 落在稳定区域内；因此最快衰减模态往往决定显式方法的最大步长。
因此，尽管解本身几乎不再变化，Euler 方法仍需要非常小的步长。这类微分方程称为刚性微分方程。形式化定义并不统一，下面的定义使用最为广泛。
定义 3.2.2： 如果 $A$ 的特征值实部非正，并且 $A$ 同时具有满足 $\operatorname{Re}(\lambda_i)\ll -1$ 的特征值，以及实部接近 $0$ 但仍为负的特征值，则称初值问题 (LIVP$_n$) 为刚性的。
下面转入刚性微分方程的数值处理。
为了推导刚性微分方程的一个简单模型方程，先考虑 $c=0$ 的 (LIVP$_n$)，即
\[y'=Ay, \qquad y(0)=y_0. \tag{3.9}\]设 $A$ 可对角化。那么存在 $M\in\mathbb{R}^{n\times n}$，使得
\[MAM^{-1}=\operatorname{diag}(\lambda_1,\ldots,\lambda_n),\]其中 $\lambda_1,\ldots,\lambda_n$ 是 $A$ 的特征值。令 $z=My$，则
\[z'=MAy=MAM^{-1}z =\operatorname{diag}(\lambda_1,\ldots,\lambda_n)z, \qquad z(0)=My_0.\]因此，$z=My$ 的分量 $z_i$ 满足
\[z_i'=\lambda_i z_i, \qquad z_i(0)=(My_0)_i. \tag{3.10}\]对刚性微分方程，还满足 $\operatorname{Re}(\lambda_i)\le0$，其中一些特征值的实部远小于 $0$，另一些特征值的实部接近 $0$ 但仍为负。
观察
如果一个数值方法对所有微分方程 (3.10) 表现良好，那么它通常也会适合原系统。
因此引出如下模型问题。
模型方程
\[y'=\lambda y, \qquad y(0)=1, \qquad \lambda\in\mathbb{C},\quad \operatorname{Re}(\lambda)<0. \tag{3.11}\]其解为
\[y(t)=e^{\lambda t},\]并且由于 $\operatorname{Re}(\lambda)<0$，有
\[\lim_{t\to\infty}y(t)=0. \tag{3.12}\]也就是说，解会根据 \(|\operatorname{Re}(\lambda)|\) 的大小以非常不同的速度衰减。为了使一个方法适合刚性微分方程，下面的要求被证明是有用的。
要求
对 (3.11) 数值得到的近似解，应尽可能好地反映精确解 $y(t)=e^{\lambda t}$ 的性质，特别是 (3.12)。
于是得到下面的定义。
定义 3.2.3（A 稳定，绝对稳定）
这里的 A 稳定来自 A-stability，也常译为绝对稳定。
若一个方法应用于模型问题 (3.11) 时，对每个步长 $h>0$ 都产生一个序列 ${u_j}_{j\in\mathbb{N}_0}$，满足
\[|u_{j+1}|\le |u_j|, \qquad \forall j\ge0,\]则称该方法绝对稳定，或 A 稳定。
对许多单步方法，应用到模型问题 (3.11) 时有关系
\[u_{j+1}=R(q)u_j, \qquad q=\lambda h,\]其中 $R:D\to\mathbb{C}$，$0\in D\subseteq\mathbb{C}$。
这里可以这样读：把一个单步法套到测试方程 $y’=\lambda y$ 上以后，一步更新通常不再需要保留完整的微分方程形式，而会化成“旧值乘以一个复数”的形式。这个复数就是 $R(q)$，称为放大因子；它只依赖于组合量 $q=\lambda h$，而不是分别依赖 $\lambda$ 和 $h$。因此，$\lambda$ 的衰减速度和步长 $h$ 会一起决定一次数值步到底是在缩小误差、保持幅值，还是放大误差。
这里的 $D$ 是 $R$ 有定义的复数范围。写成 $0\in D\subseteq\mathbb{C}$，只是说明 $R$ 至少在 $q=0$ 附近有意义，但不一定对所有复数 $q$ 都有定义。比如后面的隐式 Euler 有 $R(q)=1/(1-q)$，所以 $q=1$ 处没有定义；隐式梯形规则有 $R(q)=(1+q/2)/(1-q/2)$，所以 $q=2$ 处没有定义。稳定区域只会在 $R$ 有定义的地方讨论。
有了这个写法，稳定性判断就变得很直接：
\[u_j=R(q)^j u_0.\]如果 \(|R(q)|\le1\) ,一步不会放大当前数值；如果 $|R(q)|<1$，反复迭代后这个模态会衰减；如果 $|R(q)|>1$，即使精确解本来应该衰减，数值解也可能越算越大。稳定区域 $S$ 本质上就是所有满足 $|R(q)|\le1$ 的 $q=\lambda h$ 的集合。
例 3.2.4： 把显式 Euler 方法应用到模型问题 (3.11)，得到
\[u_{j+1} =u_j+h\lambda u_j =(1+\lambda h)u_j =(1+q)u_j,\]所以显式 Euler 方法的稳定函数为
\[R(q)=1+q.\]定义 3.2.5
$R$ 称为该单步方法的稳定函数。集合
$ S={q\in\mathbb{C}: |R(q)|\le1} $
称为该单步方法的稳定区域。
显然有
\[\text{A 稳定} \Longleftrightarrow |R(q)|\le1\quad \forall q\in\mathbb{C},\ \operatorname{Re}(q)<0\] \[\Longleftrightarrow S\supset\{q\in\mathbb{C}:\operatorname{Re}(q)<0\}.\]定义 3.2.6（L 稳定）
L-stability 是数值常微分方程文献中的约定名称。这里的 $L$ 通常不再像 A-stability 中的 A 那样展开成一个统一使用的英文短语；实际要记住的是它比 A 稳定更强，额外要求稳定函数在负实轴远端还要趋于 0。也就是说，对于 $q=\lambda h$ 中实部很大的负数，数值方法不仅不能放大，还应当把这种快速衰减模态强力压下去。
若一个方法 A 稳定，并且其稳定函数还满足
\[\lim_{q\to-\infty}R(q)=0,\]则称该方法 L 稳定。
A 稳定保证左半平面不放大；L 稳定还要求很快衰减的模态在数值上也被压到接近 0。
3.2.1 一些方法的稳定区域
显式 Euler 方法
把显式 Euler 方法应用到模型问题 (3.11)，得到
\[u_{j+1}=u_j+h\lambda u_j=(1+\lambda h)u_j,\]因此稳定函数为
\[R(q)=1+q.\]稳定区域为
\[S=\{q\in\mathbb{C}: |1+q|\le1\}.\]因此显式 Euler 方法不是 A 稳定的（例如取 $q=-1+2i$）。
注 3.2.7： 甚至可以证明，所有显式 Runge-Kutta 方法都不是 A 稳定的。
隐式 Euler 方法
隐式 Euler 方法对模型问题 (3.11) 给出
\[u_{j+1}=u_j+h\lambda u_{j+1},\]因此
\[u_{j+1}=\frac{1}{1-\lambda h}u_j.\]这给出稳定函数
\[R(q)=\frac{1}{1-q}, \qquad q\ne1,\]以及稳定区域
\[S=\{q\in\mathbb{C}: |1-q|\ge1\} \supset \{q\in\mathbb{C}: \operatorname{Re}(q)<0\}.\]因此隐式 Euler 方法是 A 稳定的，甚至是 L 稳定的。
隐式梯形规则
方法方程为
\[u_{j+1} =u_j+\frac{h}{2}\bigl(f(u_j)+f(u_{j+1})\bigr).\]对模型问题 (3.11)，得到
\[u_{j+1} =u_j+\frac{h}{2}\lambda(u_j+u_{j+1}),\]因此
\[u_{j+1} = \frac{1+\lambda h/2}{1-\lambda h/2}u_j.\]所以
\[R(q)=\frac{1+q/2}{1-q/2}, \qquad q\ne2,\]稳定区域为
\[S =\{q\in\mathbb{C}: |1+q/2|\le |1-q/2|\} =\{q\in\mathbb{C}: \operatorname{Re}(q)\le0\}.\]因此隐式梯形规则是 A 稳定的，但不是 L 稳定的，因为
\[\lim_{q\to-\infty}R(q)=-1.\]每个小图都画在 $q=\lambda h$ 的复平面上：横轴是 $\operatorname{Re}(q)$，纵轴是 $\operatorname{Im}(q)$。隐式 Euler 的稳定区域不是单纯的左半平面，而是圆盘 $|1-q|<1$ 的外部；隐式梯形规则的稳定区域才是左半平面 $\operatorname{Re}(q)\le0$。
隐式 Runge-Kutta 方法
隐式 Runge-Kutta 方法特别适合刚性微分方程。若 Butcher 表中的系数 $\alpha_{ij}$ 不构成严格下三角矩阵，就得到隐式 Runge-Kutta 方法。其方法方程为
\[k_i=k_i(t,u,h) := f\left( t+\gamma_i h,\, u+h\sum_{l=1}^{r}\alpha_{il}k_l \right), \qquad i=1,\ldots,r,\] \[\varphi(t,h;u) =\sum_{i=1}^{r}\beta_i k_i. \tag{3.13}\]注意这里求和到 $r$，而不是 $i-1$。
隐式 Runge-Kutta 方法是一个显式单步方法，只是各级 $k_i$ 由一个非线性方程组的解给出。实际上可以选择系数 $\alpha_{ij},\beta_i,\gamma_i$，使其成为一个 $L$ 稳定且阶数为 $p=2r$ 的方法。
显式 Runge-Kutta 可以逐级计算；隐式 Runge-Kutta 的 $k_i$ 通常互相依赖，因此每个时间步要先解出整组 stage。
Butcher 表
\[\begin{array}{c|ccccc} \gamma_1 & \alpha_{11} & \cdots & \cdots & \alpha_{1,r-1} & \alpha_{1,r}\\ \gamma_2 & \alpha_{21} & \cdots & \cdots & \alpha_{2,r-1} & \alpha_{2,r}\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ \gamma_r & \alpha_{r1} & \cdots & \cdots & \alpha_{r,r-1} & \alpha_{r,r}\\ \hline & \beta_1 & \beta_2 & \cdots & \beta_{r-1} & \beta_r \end{array}\]如果
\[\beta=(\beta_1,\ldots,\beta_r)^T\in\mathbb{R}^r\]且 $A=(\alpha_{ij})$ 是 $\alpha$ 系数组成的矩阵，则对模型方程 (3.11)，$r$ 级 Runge-Kutta 方法可按如下方式计算：
\[u_{j+1} = \left( 1+\lambda h\,\beta^T(I-\lambda h A)^{-1}\mathbf{1} \right)u_j\] \[=\left( 1+q\,\beta^T(I-qA)^{-1}\mathbf{1} \right)u_j,\]其中 $\mathbf{1}\in\mathbb{R}^r$ 是所有分量都为 1 的向量。由此可得
\[R(q) =1+q\,\beta^T(I-qA)^{-1}\mathbf{1} = \frac{\det(I-qA+q\mathbf{1}\beta^T)} {\det(I-qA)}.\]因此 $R(q)$ 是一个有理函数。
返回阅读 数值分析讲义（三）：常微分方程初值问题与刚性 Part I。
参考文献
- [1] P. Deuflhard and F. Bornemann. Numerische Mathematik II. de Gruyter, Berlin, 2002. 3.1.5.
- [2] P. Deuflhard and F. Hohmann. Numerische Mathematik I. de Gruyter, Berlin, 2008. 1.2.3.
- [3] H. Heuser. Gewöhnliche Differentialgleichungen. Teubner, Stuttgart, 1989. 3.1.
- [4] R. Plato. Numerische Mathematik kompakt. Vieweg Verlag, Braunschweig, 2000. 1.2.3, 6.3.2.
- [5] J. Stoer. Numerische Mathematik 1. Springer Verlag, Berlin, 1994. 1.2.3, 4.4.2.
- [6] W. Törnig and P. Spellucci. Numerische Mathematik für Ingenieure und Physiker 2. Springer Verlag, Berlin, 1990. 1.2.3.
- [7] W. Walter. Gewöhnliche Differentialgleichungen. Springer, Berlin, 1986. 3.1.
- [8] J. Werner. Numerische Mathematik 2. Vieweg Verlag, Braunschweig, 1992. 6.1.4.
英文缩写与术语说明
- ODE：ordinary differential equation，常微分方程。
- IVP：initial value problem，初值问题。
- LIVP：linear initial value problem，线性初值问题。
- RK4：fourth-order Runge-Kutta method，经典四阶 Runge-Kutta 方法。
- A-stability：A 稳定，也称绝对稳定。
- L-stability：L 稳定；除 A 稳定外，还要求稳定函数在负实轴远端趋于 0。
来源、版权与使用说明
本文主要参考 TU Darmstadt 信息学专业公开仓库中的数值分析基础课讲义： mathe3-script-2011-SoSe.pdf 原仓库包含 The Unlicense 授权说明。本文作为个人学习、翻译与知识整理用途发布，文中的中文表述、补充解释和图表重制不代表原作者或官方立场。 本文中的个人整理、中文表述、补充解释以及我重新制作的图表，可在注明作者与原文链接的前提下，用于非商业学习、交流和引用。由于本文部分内容基于 TU Darmstadt 公开讲义的翻译与整理，原始讲义及其中可能包含的材料仍应以其原作者、原仓库及相关授权说明为准。若需进行商业使用、系统转载、出版，或大规模改编，建议同时确认原始材料的授权状态。 如文中存在翻译、公式、术语或理解上的疏漏，或相关权利方认为内容使用不当，欢迎联系我指出，我会及时处理或删除。
AI feedback
AnonymousLoading AI feedback…