数值分析讲义（三）：常微分方程初值问题与刚性 Part I

… views … visitors … total

Read in English

建议先阅读 数值分析讲义（二）：数值积分。本篇是 Part I，只整理常微分方程（ordinary differential equation, ODE）初值问题中的基本数值方法、相容性、稳定性和收敛性；

---

3.1 引言

自然科学、工程技术和经济学中的许多应用都会导向常微分方程的初值问题。

初值问题

给定函数

\[f:[a,b]\times\mathbb{R}^n\to\mathbb{R}^n\]

以及初值 $y_0\in\mathbb{R}^n$。要求寻找函数

\[y:[a,b]\to\mathbb{R}^n,\]

使其导数 $y’$ 满足形如

\[y'(t)=f(t,y(t)),\qquad t\in[a,b]\]

的常微分方程，并且满足初始条件 $y(a)=y_0$。下面采用英语缩写 IVP（initial value problem）表示初值问题，简写为

\[\text{(IVP)}\qquad y'(t)=f(t,y(t)),\qquad t\in[a,b], \tag{3.1}\] \[y(a)=y_0. \tag{3.2}\]

在很多情形下，$t$ 表示时间，因此，把这类问题称为初值问题是自然的。

应用

运动方程（例如车辆动力学、行星运动）、反应动力学、电路仿真等。

下面的定理给出了 (IVP) 解的存在唯一性基础。

定理 3.1.1（存在唯一性定理）
设 $f:[a,b]\times\mathbb{R}^n\to\mathbb{R}^n$ 连续。此外，存在固定常数 $L>0$，使得

\[\|f(t,y)-f(t,z)\|\le L\|y-z\|\]

对所有 $t\in[a,b]$ 和 $y,z\in\mathbb{R}^n$ 成立（Lipschitz 条件）。则：

a) 根据 Picard-Lindelöf 定理，对每个 $y_0\in\mathbb{R}^n$，(IVP) 恰有一个解

\[y\in C^1([a,b];\mathbb{R}^n).\]

b) 若 $y,z$ 分别是初值 $y(a)=y_0$ 和 $z(a)=z_0$ 对应的解，则

\[\|y(t)-z(t)\| \le e^{L(t-a)}\|y_0-z_0\|, \qquad \forall t\in[a,b]. \tag{3.3}\]

证明可参见 Heuser [3] 或 Walter [7]。其中 b) 是 Gronwall 引理的一个推论。

注： b) 表明解连续依赖于初值 $y_0$。

3.1.1 数值方法的基本概念

为了数值求解 (IVP)，把区间 $[a,b]$ 划分为若干子区间：

\[t_j=a+jh,\qquad j=0,1,\ldots,N,\qquad h=\frac{b-a}{N}.\]

对 (IVP) 积分，并记

\[y_j=y(t_j),\qquad j=0,\ldots,N,\]

得到

\[y_{j+1} =y_j+\int_{t_j}^{t_{j+1}}y'(t)\,dt =y_j+\int_{t_j}^{t_{j+1}}f(t,y(t))\,dt. \tag{3.4}\]

右端积分无法精确计算，因为 $y(t)$ 未知。因此用插值型求积近似该积分，由此得到计算近似值

\[u_j\approx y(t_j),\qquad j=1,\ldots,N,\qquad u_0=y_0\]

的数值算法。误差

\[e_j=y(t_j)-u_j\]

称为离散化误差。

3.1.2 一些重要方法

显式 Euler 方法

若用矩形规则近似 (3.4) 中的积分，并采用左端点作为支撑点，即

\[\int_{t_j}^{t_{j+1}}f(t,y(t))\,dt \approx h f(t_j,y_j),\]

则得到显式 Euler 方法：

\[u_0:=y_0,\] \[u_{j+1}:=u_j+h f(t_j,u_j), \qquad j=0,\ldots,N-1. \tag{3.5}\]

也可以把差商

\[\frac{y(t_{j+1})-y(t_j)}{h}\]

看作 $y’(t)$ 的近似，因此它大约等于 $f(t_j,y_j)$。

方法示意图：显式 Euler 方法只使用左端点斜率

显式 Euler 的特点是一步可以直接算出，但只用左端点的一条切线，步长较大时误差会迅速累积。

隐式 Euler 方法

如果用右端点 $t_{j+1}$ 作为支撑点的矩形规则近似积分，则得到隐式 Euler 方法：

\[u_0:=y_0,\] \[u_{j+1}:=u_j+h f(t_{j+1},u_{j+1}), \qquad j=0,\ldots,N-1. \tag{3.6}\]

这里要注意，每走一步，都必须先从这个方程中解出 $u_{j+1}$。

方法示意图：隐式 Euler 方法使用右端点斜率

隐式 Euler 的斜率取在右端点，稳定性更好；代价是每一步通常要解一个非线性方程。

若用梯形规则近似 (3.4) 中的积分，则得到

\[u_{j+1} =u_j+\frac{h}{2} \left(f(t_j,u_j)+f(t_{j+1},u_{j+1})\right).\]

右端依赖于 $u_{j+1}$，因此该方法是隐式的。

方法示意图：隐式梯形规则取两端斜率的平均

梯形规则可以看成把左、右端点的斜率平均，因此比单纯 Euler 步更对称，但仍是隐式方法。

Heun 方法，第一个二阶 Runge-Kutta 方法

在右端用显式 Euler 步

\[u_{j+1}=u_j+h f(t_j,u_j)\]

替代 $u_{j+1}$，就得到 Heun 方法，也就是第一个二阶 Runge-Kutta 方法（Heun, 1900）：

\[u_0=y_0,\] \[u_{j+1} =u_j+\frac{h}{2} \left( f(t_j,u_j) +f(t_{j+1},u_j+h f(t_j,u_j)) \right), \qquad j=0,\ldots,N-1.\]

该方法也可以写成

\[u_{j+1}=u_j+\frac{h}{2}(k_1+k_2),\]

其中

\[k_1=f(t_j,u_j), \qquad k_2=f(t_{j+1},u_j+hk_1).\]

方法示意图：Heun 方法先预测，再用端点斜率修正

Heun 方法仍是显式的：红线从精确起点出发，所以画成起点切线；蓝线是在红色预测点处采样的向量场方向，不表示黑色精确曲线的切线。

改进 Euler 方法（二阶 Runge-Kutta 方法）

若用矩形规则近似积分，并用 Euler 步

\[u_j+\frac{h}{2}f(t_j,u_j)\]

近似 $u_{j+1/2}$，则得到改进 Euler 方法/二阶 Runge-Kutta 二阶方法（Runge, 1895）：

\[u_0=y_0,\] \[u_{j+1} =u_j+h f\left(t_j+\frac{h}{2},\, u_j+\frac{h}{2}f(t_j,u_j)\right), \qquad j=0,\ldots,N-1.\]

该方法也可以写成

\[u_{j+1}=u_j+hk_2,\]

其中

\[k_1=f(t_j,u_j), \qquad k_2=f\left(t_j+\frac{h}{2},u_j+\frac{h}{2}k_1\right).\]

方法示意图：改进 Euler 方法使用中点斜率

改进 Euler 方法不取右端预测斜率，而是在预测中点采样向量场方向，再用这个方向走完整一步；绿色线不应理解为黑色曲线某一点的切线。

经典四阶 Runge-Kutta 方法（RK4，fourth-order Runge-Kutta method）

最后，如果应用 Simpson 规则，并用 Taylor 展开适当地替代 $u_{j+1/2}$ 和 $u_{j+1}$，经过一些计算可得到非常精确且常用的经典四阶 Runge-Kutta 方法（下文简称 RK4）：

\[u_0=y_0,\] \[u_{j+1} =u_j+\frac{h}{6}(k_1+2k_2+2k_3+k_4), \qquad j=0,\ldots,N-1,\]

其中

\[k_1=f(t_j,u_j),\] \[k_2=f\left(t_j+\frac{h}{2},u_j+\frac{h}{2}k_1\right),\] \[k_3=f\left(t_j+\frac{h}{2},u_j+\frac{h}{2}k_2\right),\] \[k_4=f(t_{j+1},u_j+hk_3).\]

方法示意图：RK4 组合四个采样斜率

RK4 的精度来自多次向量场采样：$k_1$ 从起点切线出发，$k_2,k_3,k_4$ 是预测点处的采样方向，最终以 $k_1+2k_2+2k_3+k_4$ 组合。

例： 考虑一个串联电路：线圈电感为 $L$，电阻为 $R$，电容器电容为 $C$。电流可由下面的二阶线性微分方程描述：

\[LC\,I''(t)+RC\,I'(t)+I(t)=0.\]

为了使用数值方法，必须把该系统转化为一阶系统。为此令

\[y_1(t)=I(t), \qquad y_2(t)=I'(t).\]

得到系统

\[y'(t) = \begin{pmatrix} y_1'(t)\\ y_2'(t) \end{pmatrix} =f(t,y) = \begin{pmatrix} y_2(t)\\ -\frac{R}{L}y_2(t)-\frac{1}{LC}y_1(t) \end{pmatrix},\]

其中使用了

\[y_2'(t)=I''(t) =-\frac{R}{L}I'(t)-\frac{1}{LC}I(t).\]

现在取

\[L=C=1,\qquad R=\frac14,\]

于是

\[d=\frac14,\qquad k=1,\qquad \omega=\sqrt{1-\frac{1}{64}}.\]

初值取为

\[y_1(0)=I(0)=1,\qquad y_2(0)=I'(0)=1.\]

图 3.1： 振荡电路微分方程

\[I''(t)+\frac14 I'(t)+I(t)=0\]

在初值 $I(0)=I’(0)=1$ 下的解与近似；左图取 $n=50$，右图取 $n=100$。图中比较了显式 Euler、Heun、RK4 和精确解。

图 3.1：振荡电路微分方程的解与近似（左：$n=50$；右：$n=100$）

这张图按原讲义图 3.1 的题注组织为左右两栏：左栏 $n=50$，右栏 $n=100$。步数增加后，Heun 和 RK4 明显贴近精确解；显式 Euler 对该振荡问题的误差仍更明显。

解析解可以写成

\[I(t) =e^{-\frac18 t} \left(A\cos(\omega t)+B\sin(\omega t)\right),\]

其中

\[\omega =\sqrt{1-\left(\frac18\right)^2} =\frac{3\sqrt{7}}{8}.\]

由初值 $I(0)=1$ 可得 $A=1$（令 $t=0$）。由初值 $I’(0)=1$，通过求导并令 $t=0$ 得到

\[B=\frac{1+\frac18}{\omega} =\frac{3\sqrt{7}}{7}.\]

完整解为

\[I(t) =e^{-\frac18 t} \left( \cos(\omega t)+\frac{3\sqrt{7}}{7}\sin(\omega t) \right).\]

图 3.1 表明，用 Heun 方法和 RK4 得到的近似相当好。相反，显式 Euler 方法误差明显更大。

3.1.3 收敛性和相容性

现在考察上述方法的实际可用性和精度。这些方法可以写成一般形式

\[u_0=y_0,\] \[u_{j+1} =u_j+h\varphi(t_j,h;u_j,u_{j+1}), \qquad j=0,\ldots,N-1. \tag{3.7}\]

定义 3.1.2
函数 $\varphi(t,h;u,v)$ 在 (3.7) 中称为方法函数。如果 $\varphi$ 不依赖于 $v$，则称该方法为显式方法，否则称为隐式方法。

德语原文中的记号是 Die Funktion $\varphi(t,h;u,v)$。

这里 $u$ 表示当前步左端点的状态，在 (3.7) 中对应 $u_j$；$v$ 表示当前步右端点的状态，在 (3.7) 中对应 $u_{j+1}$。显式方法不真正使用 $v$，隐式方法则会让 $v$ 出现在方法函数中，因此需要通过方程一起求出 $u_{j+1}$。

量

\[\tau(t,h) = \frac{1}{h} \left( y(t+h)-y(t) -h\varphi(t,h;y(t),y(t+h)) \right),\]

其中 $h>0$、$t\in[a,b-h]$，称为方法 (3.7) 对 (IVP) 在位置 $t$ 处的局部截断误差或相容性误差。

换言之，它等于把精确解代入方法后得到的缺陷再除以 $h$。

定义 3.1.3
若存在常数 $C>0$ 和 $\bar h>0$，使得

\[\|\tau(t,h)\|\le Ch^p\]

对所有 $0<h\le\bar h$ 和所有 $t\in[a,b-h]$ 成立，则称方法 (3.7) 对 (IVP) 具有 $p$ 阶相容性。

性质示意图：相容性看的是一步局部缺陷

相容性不是先看很多步之后的总误差，而是问：如果这一步从精确解 $y(t)$ 出发，数值公式在 $t+h$ 处偏离精确解多少。

若存在常数 $K>0$，使得

\[\|\varphi(t,h;u,v)-\varphi(t,h;\tilde u,\tilde v)\| \le K\bigl(\|u-\tilde u\|+\|v-\tilde v\|\bigr)\]

对所有 $t\in[a,b]$ 和 $u,v,\tilde u,\tilde v\in\mathbb{R}^n$ 成立，则称方法 (3.7) 稳定。

性质示意图：稳定性看的是扰动是否被温和传播

稳定性在这里不是说单步误差为零，而是说方法函数 $\varphi(t,h;u,v)$ 对 $u$、$v$ 的小扰动不会产生不受控制的斜率变化；这正是后面把局部相容性误差转化为全局收敛性估计时需要的条件。

若存在常数 $M>0$、$H>0$，使得

\[\|e_j\| =\|y(t_j)-u_j\| \le Mh^p,\]

对 $j=0,\ldots,N$ 以及所有

\[h=\frac{b-a}{N}\le H\]

成立，则称方法 (3.7) 具有 $p$ 阶收敛性。

性质示意图：收敛性看的是网格整体误差

收敛性关心的是所有网格点上的全局离散化误差 $e_j=y(t_j)-u_j$。如果方法收敛，步长变小后整条数值轨道会贴近精确轨道。

例 3.1.4（显式 Euler 方法）： Euler 方法具有 1 阶相容性。

证明：设

\[f\in C^1([a,b]\times\mathbb{R}^n;\mathbb{R}^n)\]

且 $y$ 是 $y’=f(t,y)$ 的解。那么

\[y'\in C^1([a,b];\mathbb{R}^n),\]

所以

\[y\in C^2([a,b];\mathbb{R}^n).\]

Taylor 展开按分量给出，对某些 $\xi_i\in[0,1]$ 有

\[y(t+h) =y(t)+y'(t)h +\frac12\bigl(y_i''(t+\xi_i h)\bigr)_{1\le i\le n}h^2\] \[=y(t)+f(t,y(t))h +\frac12\bigl(y_i''(t+\xi_i h)\bigr)_{1\le i\le n}h^2.\]

因此

\[\|\tau(t,h)\| = \left\| \frac1h\bigl(y(t+h)-y(t)-h f(t,y(t))\bigr) \right\|\] \[=\frac12 \left\| \bigl(y_i''(t+\xi_i h)\bigr)_{1\le i\le n} \right\|h \le Ch.\]

这里 $C$ 是 $y’’$ 在 $[a,b]$ 上的一个上界常数。因此显式 Euler 方法具有 1 阶相容性。

3.1.4 一个收敛性定理

现在考虑显式单步方法的一个基本收敛性定理。

定理 3.1.5
设

\[y\in C^1([a,b];\mathbb{R}^n)\]

是 (IVP) 的解。方法 (3.7) 具有 $p$ 阶相容性并且稳定。则该方法具有 $p$ 阶收敛性。更精确地说，存在 $H>0$，使得全局离散化误差满足

\[\|e_j\| =\|y(t_j)-u_j\| \le \frac{e^{4K|t_j-a|}-1}{4K}\,2C h^p,\]

对 $j=0,\ldots,N$ 和所有

\[h=\frac{b-a}{N}\le H\]

成立。

证明（供感兴趣的读者参考）： 令

\[y_j=y(t_j),\qquad e_j=y_j-u_j,\qquad j=0,\ldots,N.\]

根据方法 (3.7) 和局部离散化误差的定义，对 $j=0,\ldots,N-1$ 有

\[u_{j+1}=u_j+h\varphi(t_j,h;u_j,u_{j+1}),\] \[y_{j+1} =y_j+h\varphi(t_j,h;y_j,y_{j+1}) +h\tau(t_j,h).\]

第二式减去第一式得到

\[e_{j+1} =e_j +h\left( \varphi(t_j,h;y_j,y_{j+1}) -\varphi(t_j,h;u_j,u_{j+1}) \right) +h\tau(t_j,h).\]

令

\[0<h=\frac{b-a}{N}\le\bar h.\]

由于 $t_j\in[a,b-h]$，相容性条件给出

\[\|\tau(t_j,h)\|\le Ch^p.\]

结合方法的稳定性和三角不等式，得到

\[\|e_{j+1}\| \le (1+hK)\|e_j\|+hK\|e_{j+1}\|+hCh^p.\]

取 $0<H\le\bar h$ 足够小，使得 $HK\le\frac12$。那么对所有

\[0<h=\frac{b-a}{N}\le H\]

都有

\[\|e_{j+1}\| \le \frac{1+hK}{1-hK}\|e_j\|+2Ch^{p+1} \le (1+4hK)\|e_j\|+2Ch^{p+1}.\]

下面的引理再结合 $e_0=0$，给出

\[\|e_{j+1}\| \le \frac{e^{4K|t_{j+1}-a|}-1}{4K}\,2C h^p.\]

定理得证。

为了完成证明，还需要下面的离散 Gronwall 引理，用来估计误差累积。

引理 3.1.6： 对数 $L>0$、$a_j\ge0$、$h_j>0$ 和 $b\ge0$，若

\[a_{j+1}\le(1+h_jL)a_j+h_jb, \qquad j=0,1,\ldots,n-1,\]

则

\[a_j \le \frac{e^{Lt_j}-1}{L}b +e^{Lt_j}a_0, \qquad t_j:=\sum_{i=0}^{j-1}h_i.\]

证明（供感兴趣的读者参考）： 当 $j=0$ 时命题显然成立。归纳步骤 $j\to j+1$ 由下式给出：

\[\begin{aligned} a_{j+1} &\le (1+h_jL) \left( \frac{e^{Lt_j}-1}{L}b +e^{Lt_j}a_0 \right) +h_jb\\ &\le \left( \frac{e^{L(t_j+h_j)}-1-h_jL}{L} +h_j \right)b +e^{L(t_j+h_j)}a_0\\ &= \frac{e^{Lt_{j+1}}-1}{L}b +e^{Lt_{j+1}}a_0. \end{aligned}\]

3.1.5 显式 Runge-Kutta 方法

可以通过推广 RK4 方法中的思想来构造高相容性阶的方法。

$r$ 级显式 Runge-Kutta 方法

这里选择方法函数

\[k_i(t,u,h)=k_i := f\left( t+\gamma_i h,\, u+h\sum_{j=1}^{i-1}\alpha_{ij}k_j \right), \qquad i=1,\ldots,r,\] \[\varphi(t,h;u) =\sum_{i=1}^{r}\beta_i k_i. \tag{3.8}\]

其中 $k_i=k_i(t,u,h)$ 称为第 $i$ 级。为了紧凑描述显式 Runge-Kutta 方法，通常把系数写成一个表，称为：

Butcher 表

\[\begin{array}{c|ccccc} \gamma_1 & 0 \\ \gamma_2 & \alpha_{21} & 0 \\ \gamma_3 & \alpha_{31} & \alpha_{32} & 0 \\ \vdots & \vdots & \vdots & \ddots & \ddots \\ \gamma_r & \alpha_{r1} & \cdots & \cdots & \alpha_{r,r-1} & 0 \\ \hline & \beta_1 & \beta_2 & \cdots & \beta_{r-1} & \beta_r \end{array}\]

Butcher 表例子

显式 Euler 方法：

\[\begin{array}{c|c} 0 & 0\\ \hline & 1 \end{array}\]

改进 Euler 方法：

\[\begin{array}{c|cc} 0 & 0 & 0\\ \frac12 & \frac12 & 0\\ \hline & 0 & 1 \end{array}\]

Heun 方法：

\[\begin{array}{c|cc} 0 & 0 & 0\\ 1 & 1 & 0\\ \hline & \frac12 & \frac12 \end{array}\]

用这种方法可以构造任意相容性阶 $p$ 的方法。为此必须把级数 $r$ 选得足够大。对局部截断误差作 Taylor 展开后，会得到关于系数的方程。

由 Taylor 展开可证明下面的定理。

定理 3.1.7
考虑带有方法函数 (3.8) 的 Runge-Kutta 方法 (3.7)，并满足

\[\gamma_i=\sum_{j=1}^{r}\alpha_{ij}, \qquad i=1,\ldots,r.\]

它对每个右端项

\[f\in C^p([a,b]\times\mathbb{R})\]

具有如下相容性阶，当且仅当系数满足相应方程：

$p=1$：系数满足

\[\sum_{i=1}^{r}\beta_i=1.\]

$p=2$：系数还满足

\[\sum_{i=1}^{r}\beta_i\gamma_i=\frac12.\]

$p=3$：系数还满足以下两个方程

\[\sum_{i=1}^{r}\beta_i\gamma_i^2=\frac13,\] \[\sum_{i,j=1}^{r}\beta_i\alpha_{ij}\gamma_j=\frac16.\]

$p=4$：系数还满足以下四个方程

\[\sum_{i=1}^{r}\beta_i\gamma_i^3=\frac14,\] \[\sum_{i,j=1}^{r}\beta_i\gamma_i\alpha_{ij}\gamma_j=\frac18,\] \[\sum_{i,j=1}^{r}\beta_i\alpha_{ij}\gamma_j^2=\frac{1}{12},\] \[\sum_{i,j,k=1}^{r}\beta_i\alpha_{ij}\alpha_{jk}\gamma_k=\frac{1}{24}.\]

证明略。详细推导可参见 Deuflhard 和 Bornemann [1]。

---

参考文献

• [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，初值问题。
• RK4：fourth-order Runge-Kutta method，经典四阶 Runge-Kutta 方法。

来源、版权与使用说明

本文主要参考 TU Darmstadt 信息学专业公开仓库中的数值分析基础课讲义： mathe3-script-2011-SoSe.pdf 原仓库包含 The Unlicense 授权说明。本文作为个人学习、翻译与知识整理用途发布，文中的中文表述、补充解释和图表重制不代表原作者或官方立场。 本文中的个人整理、中文表述、补充解释以及我重新制作的图表，可在注明作者与原文链接的前提下，用于非商业学习、交流和引用。由于本文部分内容基于 TU Darmstadt 公开讲义的翻译与整理，原始讲义及其中可能包含的材料仍应以其原作者、原仓库及相关授权说明为准。若需进行商业使用、系统转载、出版，或大规模改编，建议同时确认原始材料的授权状态。 如文中存在翻译、公式、术语或理解上的疏漏，或相关权利方认为内容使用不当，欢迎联系我指出，我会及时处理或删除。