数学 - 微分方程数值解 - 第 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$ 的情况，此时必定存在