数值分析 笔记 8

常微分方程

形如：

$$g(t, y, y',\cdots, y^{(k)}) = 0$$

的方程称为常微分方程，其中 $y = y(t)$ 是单变量函数，ODE 的阶数定义为未知函数的最高阶导数次数

如果最高阶导数可以分离出来，则得到显式常微分方程：

$$y^{(k)}(t) = f(t, y, y', \cdots, y^{(k - 1)})$$

对于显式 ODE，我们可以令 $u_i(t) = y^{(i-1)}(t)$，得到 $k$ 个新变量，并构造得到一个一阶 ODE 方程组：

$$\begin{cases} u_1'(t) = u_2(t) \\ u_2'(t) = u_3(t) \\ \dots \\ u_k'(t) = f(t, u_1, u_2, \cdots, u_k) \end{cases}$$

这样我们只用考虑一阶 ODE 的求解即可

初值问题

一般来说，如果不给定初始条件，ODE 会有无穷多个解，为了得到确定的解，会给出在 $t = t_0$ 时候的函数值 $y_0$，这种问题称为初值问题

而给定初始值有扰动的情况下，解得到的函数会有偏差，我们利用这个偏差来定义初值问题的稳定性：

• 如果 $t\to\infty$ 的时候 $\Delta y \in (c_1, c_2)$，即自变量充分大时偏差有界，则称该问题稳定
• 如果 $\lim\limits_{t\to\infty}\Delta y = \infty$，即无穷远处偏差发散为无穷，则称该问题不稳定
• 如果 $\lim\limits_{t\to\infty}\Delta y = 0$，即无穷远处偏差趋近于 $0$，则称该问题渐进稳定

对于非线性 ODE 来说，其稳定性分析比线性情况更复杂$一个重要的方法是局部线性化，即在特定点附近用线性 ODE 近似原方程：

非线性 ODE 的局部稳定性分析

考虑一般形式的非线性 ODE：

$$y' = f(t, y)$$

在任意点 $(t_c, y_c)$ 处进行泰勒展开：

$$f(t, y) = f(t_c, y_c) + \dfrac{\partial f}{\partial t}\bigg|_{(t_c,y_c)}(t - t_c) + \boldsymbol{J} \cdot (y - y_c) + \cdots$$

其中雅可比矩阵 $\boldsymbol{J} = \dfrac{\partial f}{\partial y}\bigg|_{(t_c,y_c)}$ 是函数对状态变量的偏导数矩阵，具体分析步骤如下：

欧拉法

对于一阶初值问题，大部分情况下其没有解析解，因此一种常见的手段是步进地计算函数值，即分别计算：

$$\begin{align*} t_0,\quad t_1,\quad t_2 \,\,\cdots\,\, t_n\,\,\cdots \end{align*}$$

处的函数值，其中 $h_n = t_{n + 1} - t_n$ 被称为步长，计算方法为：

$$y_{n + 1} = G(y_{n + 1}, y_n, y_{n - 1}, \dots, y_{n - k})$$

其中 $G$ 是近似方程，分类包括：

• 若 $k = 0$ 则称为单步法，反之为多步法
• 若 $G$ 与 $y_{n + 1}$ 无关，则称为显格式方法，反之为隐格式方法

当求解 $y' = f(t, y)$ 的时候，一种显式单步法为：

$$y_{n + 1} = y_n + h_nf(t_n, y_n)$$

这被称之为欧拉法

推导是容易的，例如使用向前差分公式：

$$f(t_n, y_n) = y'(t_n) \approx \dfrac{y_{n + 1} - y_n}{h_n}$$

或者使用数值积分的方式：

$$y_{n + 1} - y_{n} = \int_{t_n}^{t_{n + 1}}y'(s)ds = \int_{t_n}^{t_{n + 1}}f(s, y(s))ds \approx h_nf(t_n, y_n)$$

稳定性

在计算函数近似值的过程中，我们需要考虑误差在递推计算过程中的传播过程，即欧拉法的稳定性，我们定义：

对于欧拉法来说，其 $\delta_n$ 在下一个节点上产生的误差为：

$$\begin{align*} \delta_{n + 1} &= \delta_n + h_n[f(t_n, y_n + \delta_n) - f(t_n, y_n)] \\ &= \delta_n + h_n\Big[\delta_n\dfrac{\partial f}{\partial y}(t_n, y_n) + O(\delta^{2}_n)\Big] \\ &= \delta_n \Big[1 + h_n\dfrac{\partial f}{\partial y}(t_n, y_n)\Big] + O(\delta^{2}_n) \end{align*}$$

因此当 $|1 + h_n\dfrac{\partial f}{\partial y}(t_n, y_n)| \leq 1$ 的时候，这种方法是稳定的

局部截断误差

在计算

$$y_{n + 1} = G(y_{n + 1}, y_n, y_{n - 1}, \dots, y_{n - k})$$

的时候，假设 $0 \leq i \leq k$ 的时候，$y_{n - i} = y(t_{n - i})$，即这些数据是准确的，则定义局部截断误差为：

$$l_{n + 1} = y(t_{n + 1}) - y_{n + 1}$$

对于模型问题来说，解为 $y = y_0 e^{\lambda(t - t_0)}$，局部截断误差为：

$$\begin{align*} l_{n + 1} &= y(t_{n})e^{h_n\lambda} - [y(t_n) + h_n\lambda y(t_n)] \\ &= y(t_n)(e^{h_n\lambda} - 1 - h_n\lambda) = O(h^{2}_n) \end{align*}$$

并且可以证明，任意欧拉法的局部截断误差都是 $O(h^2)$ 的

如果某种解法的局部截断误差为 $O(h^{p + 1})$，则我们称这个方法具有 $p$ 阶准确度

相容性

对于一般的显式单步法：

$$y_{n + 1} = y_n + h_n\varphi(t_n, y_n, h_n)$$

其满足相容性当且仅当其准确度阶数不小于 1，等价于：

$$\varphi(t, y, 0) = f(t, y)$$

向后欧拉法与梯形法

在欧拉法的推导中，一种是使用向前差分公式进行的，如果将这个公式换为向后差分公式或中心差分公式，则得到向后欧拉法或梯形法，对应的公式分别为：

$$\begin{align*} y_{n + 1} &= y_n + h_nf(t_{n + 1}, y_{n + 1}) \\ y_{n + 1} &= y_n + \dfrac{1}{2}h_n[f(t_{n}, y_{n}) + f(t_{n + 1}, y_{n + 1})] \end{align*}$$

这两种方法都是隐式单步法

可以证明的是，这两种方法的效果都远高于向前欧拉法，以稳定的模型问题即 $\lambda \leq 0$ 为例：

• 向后欧拉法：

  $$y_{n + 1} = y_{n} + h_n\lambda y_{n + 1} \Rightarrow y_{n + 1} = \dfrac{1}{1 - h_n\lambda}y_n$$

  则误差的传播：

  $$|\delta_{n + 1}| = \dfrac{1}{|1 - h_n\lambda|}|\delta_n|$$

  由于 $\lambda \leq 0, h_n > 0$，因此一定有 $|\delta_{n + 1}| \leq |\delta_n|$
• 梯形法：

  $$y_{n + 1} = y_{n} + \dfrac{1}{2}h_n(\lambda y_n + \lambda y_{n + 1}) \Rightarrow y_{n + 1} = \dfrac{2 + h_n\lambda}{2 - h_n\lambda}y_n$$

  则误差的传播：

  $$|\delta_{n + 1}| = |\dfrac{2 + h_n\lambda}{2 - h_n\lambda}||\delta_n|$$

  由于 $\lambda \leq 0, h_n > 0$，因此也一定有 $|\delta_{n + 1}| \leq |\delta_n|$

这说明，无论步长如何，这两种方法都是稳定的，这被称为 无条件稳定，并且可以证明这对于任意的方程都是成立的

准确度方面：

• 向后欧拉法具有 $1$ 阶准确度
• 梯形法具有 $2$ 阶准确度

Runge-Kutta 方法

该方法源自欧拉法推导过程中的数值积分方法，即：

$$y(t_{n + 1}) = y(t_{n}) + \int_{t_{n}}^{t_{n + 1}}f(s, y(s))ds$$

将这个积分展开成一般的机械求积公式并适当变形，得到：

$$\begin{align} y_{n + 1} = y_{n} + h_n\sum\limits_{i = 1}^{r}c_if(t_n + \lambda_i h_n, y(t_n + \lambda_i h_n)) \end{align}$$

其中 $\lambda_{i} \in [0, 1]$ 是积分节点的相对距离因子，$r$ 为区间划分的数量

很显然，如果我们令 $r = 1, c_1 = 1, \lambda_1 = 0$，则得到了欧拉法

而对于更一般的情况，最大的问题在与我们只知道这个函数在 $t_n$ 处的值，而无法得知其在 $t_n + \lambda_ih_n$ 处的值，于是我们希望逐步使用已知的值来近似得到积分节点处的值，这便得到了一般的 Runge-Kutta 公式，具体来说：

下面介绍几种特殊的 R-K 公式

2 级 R-K 方法

可以得到，其一般的公式为：

$$\begin{cases} y_{n + 1} = y_n + h_n(c_1k_1 + c_2k_2) \\ k_1 = f(t_n, y_n) \\ k_2 = f\big(t_n + \lambda_2 h_n, y_n + \lambda_2h_nk_1\big) \\ \end{cases}$$

则一共需要确定 $c_1, c_2, \lambda_2$ 三个参数的值

下面以模型问题为例，研究应该如何确定这三个参数，代入公式：

$$\begin{align*} y_{n+1} &= y_n + h_nc_1\lambda y_n + h_nc_2\lambda(y_n +\lambda_2 h_n\lambda y_n) \\ &= y_n(1 + (c_1 + c_2)h_n\lambda + c_2\lambda_2(h_n\lambda)^{2}) \end{align*}$$

而对 $y(t_{n+1})$ 作 Taylor 展开可以得知，这种方法最多有 2 阶准确度，即：

$$l_{n + 1} = O(h^3)$$

此时要求：

$$\begin{cases} c_1 + c_2 = 1 \\[1em] c_2\lambda_2 = \dfrac{1}{2} \end{cases}$$

这个方程有无穷多组解，比较经典的解是：

1. Heun 方法（改进的欧拉法）：$c_1 = c_2 = 1/2, \lambda_2 = 1$，对应公式为：

$$\begin{cases} y_{n + 1} = y_n + \dfrac{h_n}{2}(k_1 + k_2) \\[1em] k_1 = f(t_n, y_n) \\[1em] k_2 = f\big(t_n + h_n, y_n + h_nk_1\big) \end{cases}$$

2. 中点公式：$c_1 = 0, c_2 = 1, \lambda_2 = 1/2$，对应公式为：

$$\begin{cases} y_{n + 1} = y_n + h_nk_2 \\[1em] k_1 = f(t_n, y_n) \\[1em] k_2 = f\big(t_n + \dfrac{h_n}{2}, y_n + \dfrac{h_n}{2}k_1\big) \end{cases}$$

这种方法其实就是用中点矩形法计算机械求积公式

2 级公式的稳定条件是：

$$h_n \leq \frac{-2}{\lambda}$$

3 级 R-K 方法

和 2 级的方法类似，3 级最多有 3 阶准确度，也能得到无数种，两种经典的解为：

1. 3 阶 Kutta 公式：

$$\begin{cases} y_{n + 1} = y_n + \dfrac{h_n}{6}(k_1 + 4k_2 + k_3) \\[1em] k_1 = f(t_n, y_n) \\[1em] k_2 = f(t_n + \dfrac{h_n}{2}, y_n + \dfrac{h_n}{2}k_1) \\[1em] k_3 = f(t_n + h_n, y_n - h_nk_1 + 2h_nk_2) \end{cases}$$

2. 3 阶 Ralston 公式：

$$\begin{cases} y_{n + 1} = y_n + \dfrac{h_n}{9}(2k_1 + 3k_2 + 4k_3) \\[1em] k_1 = f(t_n, y_n) \\[1em] k_2 = f(t_n + \dfrac{h_n}{2}, y_n + \dfrac{h_n}{2}k_1) \\[1em] k_3 = f(t_n + \dfrac{3h_n}{4}, y_n + \dfrac{3h_n}{4}k_2) \end{cases}$$

3 级公式的稳定条件是：

$$h_n \leq \frac{-2.51}{\lambda}$$

4 级 R-K 方法

和 2 级的方法类似，4 级最多有 4 阶准确度，也能得到无数种，一种经典的解为：

$$\begin{cases} y_{n + 1} = y_n + \dfrac{h_n}{6}(k_1 + 2k_2 + 2k_3 + k_4) \\[1em] k_1 = f(t_n, y_n) \\[1em] k_2 = f(t_n + \dfrac{h_n}{2}, y_n + \dfrac{h_n}{2}k_1) \\[1em] k_3 = f(t_n + \dfrac{h_n}{2}, y_n + \dfrac{h_n}{2}k_2) \\[1em] k_4 = f(t_n + h_n, y_n + h_nk_3) \end{cases}$$

对于这种方法，如果我们令 $f(t, y) = f(t)$，也即这个常微分方程是一个求积分的问题，则这种方法就等价于 Simpson 公式

4 级公式的稳定条件是：

$$h_n \leq \frac{-2.78}{\lambda}$$

多步法

多步法是指在公式：

$$y_{n + 1} = G(y_{n + 1}, y_n, y_{n - 1}, \dots, y_{n - k})$$

中，$k \geq 1$ 的方法，一般来说，求解初值问题的方法是线性多步法，固定步长为 $h$，则 $m$ 步对应的公式可以写成：

$$y_{n + 1} = \sum\limits_{i = 1}^{m}\alpha_i y_{n + 1 - i} + h\sum\limits_{i = 0}^{m}\beta_i f(t_{n + 1 - i}, y_{n + 1 - i})$$

通过 Taylor 展开等方法，可以证明，在达到 $p$ 阶准确度时，需要让 Taylor 展开的前 $p + 1$ 个系数为零，其中前三个对应的方程为：

$$\begin{cases} \sum\limits_{i = 1}^{m} \alpha_{i} = 1 \\[1em] \sum\limits_{i = 1}^{m}i\alpha_{i} = \sum\limits_{i = 0}^{m} \beta_i \\[1em] \dfrac{1}{2}\sum\limits_{i = 1}^{m}i^{2}\alpha_{i} = \sum\limits_{i = 1}^{m} i\beta_i \\[1em] \end{cases}$$

其中前两个方程保证了相容性

由于这种方法比较繁琐，一种更简单的保证有 $p$ 阶准确度的方式是：
若 $y(t) = t^{i}\quad (i = 0, 1, \dots, p)$ 均满足 $l_{n + 1} = 0$，则满足条件的参数列对应的多步法一定至少有 $p$ 阶准确度

Adams 公式

一种常用的多步法公式：

$$y_{n + 1} = y_n + h\sum\limits_{i = 0}^{m}\beta_i f(t_{n + 1 - i}, y_{n + 1 - i})$$

则利用 Taylor 展开可以推出，当 $\beta_0 = 0$ 的时候，其参数的系数矩阵为 Vandermonde 矩阵的转置，因此一定是非奇异的，有唯一解

可以证明：当 $\beta_0 = 0$ 时，其一定至少有 $m$ 阶准确度，反之一定具有至少 $m + 1$ 阶准确度

由于求解 Vandermonde 矩阵为系数的方程组求解比较困难，一种求系数的方式是，将数值积分公式中的被积函数 $f(s, y(s))$ 用 Langrange 插值，之后直接积分，并且根据系数对应关系得到 $\beta_i$

利用这种方法，可以得到几种常用的 Adams 公式，例如：

1. 显式 4 阶 Adams-Bashforth 公式：

   $$y_{n + 1} = y_{n} + \dfrac{h}{24}(55f_n - 59f_{n - 1} + 37f_{n - 2} - 9f_{n - 3})$$

2. 隐式 4 阶 Adams-Bashforth 公式：

   $$y_{n + 1} = y_{n} + \dfrac{h}{24}(9f_{n + 1} + 19f_{n} - 5f_{n - 1} + f_{n - 2})$$

Adams 参数表为：

几种显式公式的系数

几种隐式公式的系数

其中稳定阈值为 $q$ 代表该方法在 $(q, 0)$ 上稳定，而误差常数则是指局部截断误差中的首项系数