数学 - 微分方程数值解 - 第 4 章 抛物型方程的差分解法 - 4.2 向前 Euler 格式

4.2 向前 Euler 格式

为了用差分方法求解问题 \((4.1.1)\sim (4.1.3)\)，将如下求解区域剖分。

\[\Omega = \{(x,t) \, | \, 0\leqslant x \leqslant 1, 0\leqslant t \leqslant T\} \]

将区间 \([0,1]\) 做 \(m\) 等分，将区间 \([0,T]\) 做 \(n\) 等分，并记空间上的步长 \(h=1/m\)，时间上的步长 \(\tau = T/n\)，此外 \(x_i=ih\)，\(0\leqslant i \leqslant m\)，\(t_k=k\tau\)，\(0\leqslant k \leqslant n\)，用两簇平行直线

\[x = x_i, \quad 0 \leqslant i \leqslant m\\ t = t_k, \quad 0 \leqslant k \leqslant n \]

将 \(\Omega\) 分割成矩形的网格，记 \(\Omega_h=\{x_i \, | \, 0 \leqslant i \leqslant m\}\)，\(\Omega_{\tau}=\{t_k \, | \, 0 \leqslant k \leqslant n\}\)，\(\Omega_{h\tau}=\Omega_h \times \Omega_{\tau}\)，称 \((x_i,t_k)\) 为结点；称 \(x=0\)，\(x=1\) 以及 \(t=0\) 上的结点为边界结点，其他所有结点为内部结点；称在直线 \(t=t_k\) 上的所有结点 \(\{(x_i,t_k) \, | \, 0 \leqslant i \leqslant m\}\) 为第 \(k\) 层结点。

设 \(\{v_i^k \, | \, 0 \leqslant i \leqslant m, 0 \leqslant k \leqslant n\}\) 为 \(\Omega_{h \tau}\) 上的一个网格函数，引入以下记号：

4.2.1 差分格式的建立

定义 \(\{\Omega_{h\tau}\}\) 上的网格函数

\[U=\{U_i^k \, | \, 0 \leqslant i \leqslant m, 0\leqslant k \leqslant n\}, \quad U_i^k = u(x_i, t_k) \]

在结点处考虑微分方程 \((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 0 \leqslant k \leqslant n-1 \tag{4.2.1} \]

将

\[\begin{align*} \frac{\partial^2 u}{\partial x^2}(x_i,t_k) & = \frac{1}{h^2} \left[ u(x_{i-1},t_k) - 2 u(x_i, t_k) + u(x_{i+1},t_k)\right] - \frac{h^2}{12} \frac{\partial^4 u}{\partial x^4}(\xi_{ik},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} \end{align*} \]

和

\[\begin{align*} \frac{\partial u}{\partial t}(x_i,t_k) & = \frac{1}{\tau} \left[ u(x_i,t_{k+1}) - u(x_i, t_k) \right] - \frac{\tau}{2} \frac{\partial^2 u}{\partial t^2}(x_{i},\eta_{ik}) \\ & = D_t U_i^k - \frac{\tau}{2} \frac{\partial^2 u}{\partial t^2}(x_{i},\eta_{ik}), \quad t_{k} < \eta_{ik} < t_{k+1} \end{align*} \]

代入 \((4.2.1)\) 得到

\[D_t U_i^k - a \delta_x^2 U_i^k = f(x_i, t_k) + \frac{\tau}{2} \frac{\partial^2 u}{\partial t^2}(x_{i},\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 0 \leqslant k \leqslant n-1 \tag{4.2.2} \]

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

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

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

在 \((4.2.2) \sim (4.2.4)\) 中略去小量项

\[R_{ik}^{(1)} = \frac{\tau}{2} \frac{\partial^2 u}{\partial t^2}(x_{i},\eta_{ik}) - \frac{a h^2}{12} \frac{\partial^4 u}{\partial x^4}(\xi_{ik},t_k) \tag{4.2.5} \]

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

\[D_t u_i^k - a \delta_x^2 u_i^k = f(x_i, t_k), \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \tag{4.2.6} \]

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

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

称差分格式 \((4.2.6)\sim (4.2.8)\) 为向前 \(\text{Euler}\) 格式，称 \(R_{ik}^{(1)}\) 为差分格式 \((4.2.6)\) 的局部截断误差。记

\[c_1 = \max \left\{ \frac{1}{2} \max_{0\leqslant x \leqslant 1\\0\leqslant t \leqslant T} \left| \frac{\partial^2 u}{\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.2.9} \]

则有

\[|R_{ik}^{(1)}| \leqslant c_1 (\tau + h^2), \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \tag{4.2.10} \]

4.2.2 差分格式的求解

记 \(r\) 称为步长比，表示为

\[r = \frac{a \tau}{h^2} \]

利用上式，可将差分格式 \((4.2.6)\) 写成

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

由于差分格式 \((4.2.6) \sim (4.2.8)\) 是显式的，有时也称差分格式 \((4.2.6) \sim (4.2.8)\) 为古典显式格式，可写成矩阵形式

\[\begin{align*} \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-2r & r \\ r & 1-2r & r \\ & \ddots & \ddots & \ddots \\ & & r & 1-2r & r \\ & & & r & 1-2r \\ \end{bmatrix} \begin{bmatrix} u_1^{k} \\ u_2^{k} \\ \vdots \\ u_{m-2}^{k} \\ u_{m-1}^{k} \end{bmatrix} + \begin{bmatrix} \tau f(x_1, t_k) + r u_0^k \\ \tau f(x_2, t_k) \\ \vdots \\ \tau f(x_{m-2}, t_k) \\ \tau f(x_{m-1}, t_k) + r u_m^k \end{bmatrix} \end{align*} \]

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

定理 4.2.1

考虑如下差分方程

\[D_t u_i^k - a \delta_x^2 u_i^k = g_i^k, \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \tag{4.2.11} \]

\[u_i^0 = \varphi_i, \quad 0 \leqslant i \leqslant m \tag{4.2.12} \]

\[u_0^k = 0, \quad u_m^k = 0, \quad 1 \leqslant k \leqslant n \tag{4.2.13} \]

设 \(u = \{u_i^k \, | \, 0\leqslant i \leqslant m, \, 0\leqslant k \leqslant n \}\) 为上述差分方程的解，则当步长比 \(r \leqslant 1/2\) 时，有

\[\|u^k\|_{\infty} \leqslant \|\varphi\|_{\infty} + \tau \sum_{l=0}^{k-1} \|g^l\|_{\infty}, \quad 0 \leqslant k \leqslant n \]

其中，\(\|g^l\|_{\infty} = \max_{1\leqslant i \leqslant m-1} |g_i^l|\)。

证明：将式 \((4.2.11)\) 写成

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

则有

\[|u_i^{k+1}| \leqslant (1-2r) \|u^k\|_{\infty} + r (\|u^k\|_{\infty} + \|u^k\|_{\infty}) + \tau \|g^k\|_{\infty} = \|u^k\|_{\infty} + \tau \|g^k\|_{\infty}, \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \]

于是

\[\|u^{k+1}\|_{\infty} \leqslant \|u^k\|_{\infty} + \tau \|g^k\|_{\infty}, \quad 0 \leqslant k \leqslant n-1 \]

递推可得

\[\|u^k\|_{\infty} \leqslant \|\varphi\|_{\infty} + \tau \sum_{l=0}^{k-1} \|g^l\|_{\infty}, \quad 0 \leqslant k \leqslant n \]

证毕。

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

有差分格式求解的矩阵表达形式可以看出，对时间上第 \(k+1\) 层的值可由第 \(k\) 层的值显式表示出来。即若已知第 \(k\) 层的值 \(\{u_i^k \, | \, 0 \leqslant i \leqslant m\}\)，则由上式就可直接得到第 \(k+1\) 层上的值 \(\{u_i^{k+1} \, | \, 0 \leqslant i \leqslant m\}\)，因而对于任意步长比 \(r\)，均是唯一可解的。即既有存在性，也有唯一性。

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

(1) 收敛性

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

定理 4.2.2

设 \(\{u(x,t) \, | \,0 \leqslant x\leqslant 1,\, 0 \leqslant t \leqslant T\}\) 为定解问题的 \((4.1.1) \sim (4.1.3)\) 的解，\(\{u_i^k \, | \, 0 \leqslant i\leqslant m,\, 0 \leqslant k \leqslant n\}\) 为差分格式 \((4.2.6) \sim (4.2.8)\) 的解，则当 \(r \leqslant 1/2\) 时，有

\[\max_{0 \leqslant i \leqslant m} |u(x_i, t_k) - u_i^k| \leqslant c_1 T(\tau + h^2), \quad 0 \leqslant k \leqslant n \]

其中，\(c_1\) 由 \((4.2.9)\) 定义。

证明：我们记误差

\[e_i^k = u(x_i, t_k) - u_i^k, \quad 0 \leqslant i\leqslant m, \quad 0 \leqslant k \leqslant n \]

将 \((4.2.2) \sim (4.2.4)\) 分别与 \((4.2.6) \sim (4.2.8)\) 相减，得到误差方程

\[D_t e_i^k - a \delta_x^2 e_i^k = R_{ik}^{(1)}, \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \]

\[e_i^0 = 0, \quad 0 \leqslant i \leqslant m \]

\[u_0^k = 0, \quad u_m^k = 0, \quad 1 \leqslant k \leqslant n \]

应用定理 \(4.2.1\)，并注意到 \((4.2.10)\)，当 \(r \leqslant 1/2\) 时，有

\[\|e_k\|_{\infty} \leqslant \tau \sum_{l=0}^{k-1} \max_{1 \leqslant i \leqslant m-1} |R_{il}^{(1)}| \leqslant \tau c_1 k (\tau + h^2) \leqslant c_1 T (\tau + h^2), \quad 0 \leqslant k \leqslant n \]

证毕。

(2) 稳定性

如果在应用差分格式 \((4.2.6) \sim (4.2.8)\) 时，计算右端函数 \(f(x_i, t_k)\) 有误差 \(g_i^k\)，计算初值 \(\varphi(x_i)\) 有误差 \(\psi_i\)，则实际得到的是如下差分方程的解

\[D_t v_i^k - a \delta_x^2 v_i^k = f(x_i, t_k) + g_i^k, \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \tag{4.2.14} \]

\[v_i^0 = \varphi(x_i) + \psi_i, \quad 0 \leqslant i \leqslant m \tag{4.2.15} \]

\[v_0^k = \alpha(t_k), \quad v_m^k = \beta(t_k), \quad 1 \leqslant k \leqslant n \tag{4.2.16} \]

令

\[\varepsilon_i^k = v_i^k - u_i^k, \quad 0 \leqslant i \leqslant m, \quad 0 \leqslant k \leqslant n \]

将 \((4.2.6) \sim (4.2.8)\) 与 \((4.2.14) \sim (4.2.16)\) 相减，可得摄动方程组

\[D_t \varepsilon_i^k - a \delta_x^2 \varepsilon_i^k = g_i^k, \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \tag{4.2.17} \]

\[\varepsilon_i^0 = \psi_i, \quad 0 \leqslant i \leqslant m \tag{4.2.18} \]

\[\varepsilon_0^k = 0, \quad \varepsilon_m^k = 0, \quad 1 \leqslant k \leqslant n \tag{4.2.19} \]

应用定理 \(4.2.1\)，当 \(r \leqslant 1/2\) 时，有

\[\|\varepsilon^k\|_{\infty} \leqslant \|\psi\|_{\infty} + \tau \sum_{l=0}^{k-1} \|g^l\|_{\infty}, \quad 1 \leqslant k \leqslant n \]

上式说明当 \(\|\psi\|_{\infty}\) 和 \(\tau \sum_{l=0}^{n-1} \|g^l\|_{\infty}\) 很小时，误差 \(\max_{1 \leqslant k \leqslant n} \|\varepsilon^k\|_{\infty}\) 也很小。

上述结果可陈述如下：

定理 4.2.3

当 \(r \leqslant 1/2\) 时，差分格式 \((4.2.6) \sim (4.2.8)\) 关于初值和右端项在下述意义下是稳定的：考虑差分方程组

\[D_t u_i^k - a \delta_x^2 u_i^k = f_i^k, \quad 1 \leqslant i \leqslant m-1, \quad 0 \leqslant k \leqslant n-1 \]

\[u_i^0 = \psi_i, \quad 0 \leqslant i \leqslant m \]

\[u_0^k = 0, \quad u_m^k = 0, \quad 1 \leqslant k \leqslant n \]

设 \(\{u_i^k \, | \, 0 \leqslant i \leqslant m, \, 0 \leqslant k \leqslant n\}\) 为上述差分方程组的解，则有

\[\|u^k\|_{\infty} \leqslant \|u^0\|_{\infty} + \tau \sum_{l=0}^{k-1} \|f^l\|_{\infty}, \quad 1 \leqslant k \leqslant n \]

下面考虑 \(r > 1/2\) 的情况，此时必定存在