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

4.4 Richardson 格式

向前 Euler 格式和向后 Euler 格式的收敛阶均为 \(O(\tau + h^2)\)，记关于时间步长 \(\tau\) 是一阶的，关于空间步长 \(h\) 是二阶的。为了使得收敛阶提高到 \(O(\tau^2 + h^2)\)，一个很自然的想法是用关于时间 \(t\) 的一阶中心差商去近似关于时间的导数。

4.4.1 差分格式的建立

在结点处考虑定解问题 \((4.1.1)\)，有

\[\frac{\partial u}{\partial t}(x_i, t_k) - a \frac{\partial^2 u}{\partial x^2}(x_i, t_k) = f(x_i, t_k), \quad 1 \leqslant i \leqslant m-1, \quad 1 \leqslant k \leqslant n-1 \tag{4.4.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) = D_{\hat{t}} U_i^k - \frac{\tau^2}{6} \frac{\partial^3 u}{\partial t^3}(x_i,\tilde{\eta}_{ik}), \quad t_{k-1} < \tilde{\eta}_{ik} < t_{k+1} \]

代入到式 \((4.4.1)\) 得到

\[D_{\hat{t}} U_i^k - a \delta_x^2 U_i^k = f(x_i, t_k) + \frac{\tau^2}{6} \frac{\partial^3 u}{\partial t^3}(x_i,\tilde{\eta}_{ik}) - \frac{a h^2}{12} \frac{\partial^4 u}{\partial x^4} (\xi_{ik}, t_k), \quad 1 \leqslant i \leqslant m-1, \quad 1 \leqslant k \leqslant n-1 \tag{4.4.2} \]

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

\[U_i^0 = u(x_i, 0) = \varphi(x_i), \quad 0 \leqslant i \leqslant m \tag{4.4.3} \]

\[U_0^k = \alpha(t_k), \quad U_m^k = \beta(t_k), \quad 1 \leqslant k \leqslant n \tag{4.4.4} \]

在 \((4.3.2) \sim (4.3.4)\) 中略去小量项

\[R_{ik}^{(3)} = \frac{\tau^2}{6} \frac{\partial^3 u}{\partial t^3}(x_i, \tilde{\eta}_{ik}) - \frac{a h^2}{12} \frac{\partial^4 u}{\partial x^4} (\xi_{ik}, t_k) \]

并用 \(u_i^k\) 代替 \(U_i^k\)，得到如下差分格式

\[D_{\hat{t}} u_i^k - a \delta_x^2 u_i^k = f(x_i, t_k), \quad 1 \leqslant i \leqslant m-1, \quad 1 \leqslant k \leqslant n-1 \tag{4.4.5} \]

\[u_i^0 = \varphi(x_i), \quad 0 \leqslant i \leqslant m \tag{4.4.6} \]

\[u_0^k = \alpha(t_k), \quad u_m^k = \beta(t_k), \quad 1 \leqslant k \leqslant n \tag{4.4.7} \]

称差分格式 \((4.4.5) \sim (4.4.7)\) 为 Richardson 格式。称 \(R_{ik}^{(3)}\) 为差分格式 \((4.4.5)\) 的局部截断误差。显然 Richardson 格式是一个显式格式。

4.4.2 差分格式的求解

差分格式 \((4.4.5) \sim (4.4.7)\) 可写成

\[u_i^{k+1} = 2r (u_{i-1}^k - 2 u_i^k + u_{i+1}^k) + u_i^{k-1} + 2 \tau f(x_i, t_k), \quad 1 \leqslant i \leqslant m-1, \quad 1 \leqslant k \leqslant n-1 \]

应用 Richardson 格式计算第 \(k+1\) 层上的值时，需要用到第 \(k\) 层和第 \(k-1\) 层上的值，它是一个三层格式。实际计算时，当第 \(0\) 层和第 \(1\) 层上的值已知时，可依此求出第 \(2\) 层、第 \(3\) 层、\(\cdots\)，直到第 \(n\) 层的值，现第 \(0\) 层的值已由 \((4.4.6)\) 给出，第 \(1\) 层的值如何给出呢？

由方程 \((4.1.1)\) 以及初值条件 \((4.1.2)\)，有