数学 - 微分方程数值解 - 第 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) 稳定性