数学 - 微分方程数值解 - 第 4 章 抛物型方程的差分解法 - 4.5 Crank-Nicolson 格式

4.5 Crank-Nicolson 格式

本节对于定解问题 $(3.1.1) \sim (3.1.3)$ 建立一个具有 $O(\tau^2 + h^2)$ 精度的无条件稳定的差分格式。

注意，对各个符号取上标 $k+\frac{1}{2}$ 和取下标 $k+\frac{1}{2}$ 的意义可能各不相同，需要仔细甄别。

4.5.1 差分格式的建立

(1) 建立差分格式

我们记 $t_{k+\frac{1}{2}}=\frac{1}{2}(t_k + t_{k+1})$，在点 $(x_i,t_{k+\frac{1}{2}})$ 处考虑微分方程 $(4.1.1)$，有

$$\frac{\partial u}{\partial t}(x_i,t_{k+\frac{1}{2}}) - a \frac{\partial^2 u}{\partial x^2}(x_i, t_{k+\frac{1}{2}}) = f(x_i, t_{k+\frac{1}{2}}), \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1$$

应用公式

$$\frac{\partial^2 u}{\partial x^2}(x_i,t_{k+\frac{1}{2}}) = \frac{1}{2} \left[ \frac{\partial^2 u}{\partial x^2}(x_i,t_k) + \frac{\partial^2 u}{\partial x^2}(x_i,t_{k+1})\right]-\frac{\tau^2}{8} \frac{\partial^4 u}{\partial x^2 \partial t^2}(x_i,\zeta_{ik}), \quad \zeta_{ik} \in (t_k, t_{k+1})$$

可以得到

$$\frac{\partial u}{\partial t}(x_i, t_{k+\frac{1}{2}}) - \frac{1}{2} a \left[ \frac{\partial^2 u}{\partial x^2}(x_i,t_k) + \frac{\partial^2 u}{\partial x^2}(x_i,t_{k+1})\right] = f(x_i, t_{k+\frac{1}{2}}) - \frac{a \tau^2}{8} \frac{\partial^4 u}{\partial x^2 \partial t^2}(x_i,\zeta_{ik}), \quad \zeta_{ik} \in (t_k, t_{k+1})$$

利用式子

$$\frac{\partial^2 u}{\partial x^2}(x_i,t_k) = \delta_x^2 U_i^k - \frac{h^2}{12} \frac{\partial^4 u}{\partial x^4}(\xi_{ik},t_k), \quad x_{i-1} < \xi_{ik} < x_{i+1}$$

和式子

$$\frac{\partial u}{\partial t}(x_i, t_{k+\frac{1}{2}}) = \delta_t U_i^{k+\frac{1}{2}} - \frac{\tau^2}{24} \frac{\partial^3 u}{\partial t^3}(x_i, \eta_{ik}), \quad t_k < \eta_{ik} < t_{k+1}$$

得到

$$\delta_t U_i^{k+\frac{1}{2}} - a \delta_x^2 U_i^{k+\frac{1}{2}} = f(x_i, t_{k+\frac{1}{2}}) + R_{ik}^{(4)} \tag{4.5.1}$$

其中

$$\delta_t U_i^{k+\frac{1}{2}} = \frac{1}{\tau} (U_i^{k+1} - U_i^{k}), \quad \delta_x^2 U_i^{k+\frac{1}{2}} = \frac{1}{2}( \delta_x^2 U_i^{k+1} + \delta_x^2 U_i^{k})$$

$$R_{ik}^{(4)} = \tau^2 \left[ \frac{1}{24} \frac{\partial^3 u}{\partial t^3}(x_i, \eta_{ik}) - \frac{a}{8} \frac{\partial^4 u}{\partial x^2 \partial t^2}(x_i,\zeta_{ik})\right] - \frac{a h^2}{24} \left[ \frac{\partial^4 u}{\partial x^4}(\xi_{ik},t_k) + \frac{\partial^4 u}{\partial x^4}(\xi_{i,k+1},t_{k+1}) \right]$$

注意到初边值条件 $(4.1.2)$ 和 $(4.1.3)$，有

$$U_i^0 = \varphi(x_i), \quad 0 \leqslant i \leqslant m \tag{4.5.2}$$

$$U_0^k = \alpha(t_k), \quad U_m^k = \beta(t_k), \quad 1 \leqslant k \leqslant n \tag{4.5.3}$$

在 $(4.5.1) \sim (4.5.3)$ 中略去小量项 $R_{ik}^{(4)}$，并用 $u_{i}^k$ 代替 $U_i^k$，得到如下差分格式：

$$\delta_t u_i^{k+\frac{1}{2}} - a \delta_x^2 u_i^{k+\frac{1}{2}} = f(x_i, t_{k+\frac{1}{2}}) , \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1\tag{4.5.4}$$

$$u_i^0 = \varphi(x_i), \quad 0 \leqslant i \leqslant m \tag{4.5.5}$$

$$U_0^k = \alpha(t_k), \quad U_m^k = \beta(t_k), \quad 1 \leqslant k \leqslant n \tag{4.5.6}$$

称差分格式 $(4.5.4) \sim (4.5.6)$ 为 Crank-Nicolson 格式。

(2) 局部截断误差

称 $R_{ik}^{(4)}$ 为差分格式 $(4.5.4)$ 的局部截断误差，记

$$c_2 = \max \left\{ \frac{1}{24} \max_{0\leqslant x \leqslant 1\\0\leqslant t \leqslant T} \left| \frac{\partial^3 u}{\partial t^3}(x,t)\right|,\, \frac{a}{8} \max_{0\leqslant x \leqslant 1\\0\leqslant t \leqslant T} \left| \frac{\partial^4 u}{\partial x^2 \partial t^2}(x,t)\right|,\, \frac{a}{12} \max_{0\leqslant x \leqslant 1\\0\leqslant t \leqslant T} \left| \frac{\partial^4 u}{\partial x^4}(x,t)\right|,\, \right\} \tag{4.5.7}$$

则有

$$R_{ik}^{(4)} \leqslant c_2 (\tau^2 + h^2), \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1\tag{4.5.8}$$

4.5.2 差分格式的求解

记

$$u^k = (u_0^k, u_1^k, \cdots, u_{m-1}^k, u_m^k)$$

由式 $(4.5.5)$ 知第 $0$ 层上的值 $u^0$ 已知，设已经确定出了第 $k$ 层的值 $u^k$，则关于第 $k+1$ 层值的差分格式为

$$\delta_t u_i^{k+\frac{1}{2}} - a \delta_x^2 u_i^{k+\frac{1}{2}} = f(x_i, t_{k+\frac{1}{2}}), \quad 1 \leqslant i \leqslant m-1$$

$$u_0^{k+1} = \alpha(t_{k+1}), \quad u_m^{k+1} = \beta(t_{k+1})$$

可以写成如下矩阵形式

$$\begin{align*} \begin{bmatrix} 1 + r & -\frac{r}{2} \\ -\frac{r}{2} & 1 + r & -\frac{r}{2} \\ & \ddots & \ddots & \ddots \\ & & -\frac{r}{2} & 1 + r & -\frac{r}{2} \\ & & &-\frac{r}{2} & 1 + r \end{bmatrix} \begin{bmatrix} u_1^{k+1} \\ u_2^{k+1} \\ \vdots\\ u_{m-2}^{k+1} \\ u_{m-1}^{k+1} \end{bmatrix} = \begin{bmatrix} 1 - r & \frac{r}{2} \\ \frac{r}{2} & 1 - r & \frac{r}{2} \\ & \ddots & \ddots & \ddots \\ & & \frac{r}{2} & 1 - r & \frac{r}{2} \\ & & &\frac{r}{2} & 1 - r \end{bmatrix} \begin{bmatrix} u_1^{k} \\ u_2^{k} \\ \vdots\\ u_{m-2}^{k} \\ u_{m-1}^{k} \end{bmatrix} + \begin{bmatrix} \frac{r}{2}(u_0^k + u_0^{k+1}) + \tau f(x_1, t_{k+\frac{1}{2}}) \\ \tau f(x_2, t_{k+\frac{1}{2}}) \\ \vdots\\ \tau f(x_{m-2}, t_{k+\frac{1}{2}}) \\ \frac{r}{2}(u_m^k + u_m^{k+1}) + \tau f(x_{m-1}, t_{k+\frac{1}{2}}) \end{bmatrix} \end{align*} \tag{4.5.9}$$

由式 $(4.5.9)$ 可知在每一个时间层式 $(4.5.4) \sim (4.5.6)$ 在每一个时间层均为一个三对角线性方程组，具体算法可采用追赶法求解。

4.5.3 差分格式解的先验估计式

定理 4.5.1

考虑差分方程组

$$\delta_t u_i^{k+\frac{1}{2}} - a \delta_x^2 u_i^{k+\frac{1}{2}} = f_i^{k+\frac{1}{2}}, \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \tag{4.5.10}$$

$$u_i^0 = \varphi_i, \quad 1 \leqslant i \leqslant m-1 \tag{4.5.11}$$

$$u_0^k = 0, \quad u_m^k = 0, \quad 0 \leqslant k \leqslant n \tag{4.5.12}$$

设 $\{v_i^k \, | \, 0 \leqslant i \leqslant m, \, 0 \leqslant k \leqslant n\}$ 为上述差分方程组的解，则对任意步长比 $r$，有

$$\|u^k\|^2 \leqslant \|u^0\|^2 + \frac{1}{12a} \tau \sum_{l=0}^{k-1} \|f^{l+\frac{1}{2}}\|^2, \quad 0 \leqslant k \leqslant n \tag{4.5.13}$$

$$|u^k|_1^2 \leqslant |u^0|_1^2 + \frac{1}{2a} \tau \sum_{l=0}^{k-1} \|f^{l+\frac{1}{2}}\|^2, \quad 0 \leqslant k \leqslant n \tag{4.5.14}$$

$$\|u^k\|_{\infty} \leqslant \frac{1}{2} \left( |u^0|_1^2 + \frac{1}{2a} \tau \sum_{l=0}^{k-1} \|f^{l+\frac{1}{2}}\|^2 \right)^{1/2}, \quad 0 \leqslant k \leqslant n \tag{4.5.15}$$

证明：

证毕。

4.5.4 差分格式解的存在性与唯一性

定理 4.5.2

差分格式 $(4.5.4) \sim (4.5.6)$ 的解存在且唯一。

证明：差分格式 $(4.5.4) \sim (4.5.6)$ 对应的矩阵形式为式 $(4.5.9)$，易见其系数矩阵是严格对角占优的，因而存在唯一解。

证毕。

4.5.5 差分格式解的收敛性与稳定性

(1) 收敛性

给出差分格式的收敛性的相关定理。

定理 4.5.3

设 $\{u(x,t) \, | \, 0\leqslant x \leqslant 1,\, 0 \leqslant t \leqslant T\}$ 为定解问题 $(3.1.1) \sim (3.1.3)$ 的解，$\{u_i^k \, | \, 0 \leqslant i \leqslant m, \, 0 \leqslant k \leqslant n\}$ 为差分格式 $(4.5.4) \sim (4.5.6)$ 的解，则对任意步长比 $r$，有

$$\max_{1 \leqslant i \leqslant m-1} |u(x_i,t_k) - u_i^k| \leqslant \frac{c_2}{2}\sqrt{\frac{T}{2a}}(\tau^2 + h^2), \quad 1 \leqslant k \leqslant n$$

其中，$c_2$ 由式 $(4.5.7)$ 定义。

证明：

证毕。

(2) 稳定性