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) 稳定性