数学 - 微分方程数值解 - 第 3 章椭圆型方程的差分解法 - 3.2 五点差分格式 - Black_x

数学 - 微分方程数值解 - 第 3 章椭圆型方程的差分解法 - 3.2 五点差分格式

3.2 五点差分格式

3.2.1 五点差分格式的建立

(1) 建立差分格式

将区间 $[a,b]$ 做 $m$ 等分，记

h_{1} = \frac{b - a}{m}, x_{i} = a + i h_{1}, i = 0, \dots, m

$h_1 = \frac{b-a}{m}, \quad x_i = a + ih_1, \quad i=0,\cdots,m$

将区间 $[c,d]$ 做 $n$ 等分，记

h_{2} = \frac{d - c}{n}, y_{j} = c + j h_{2}, j = 0, \dots, n

$h_2 = \frac{d-c}{n}, \quad y_j = c + jh_2, \quad j=0,\cdots,n$

称 $h_1$ 为 $x$ 方向步长， $h_2$ 为 $y$ 方向步长。用两簇平行线 $x = x_i$ 与 $y=y_j$ 将区域 $\Omega$ 剖分为 $mn$ 个小矩形，称两簇直线的交点 $(x_i,y_j)$ 为网格结点。

所有结点构成集合

Ω_{h} = {(x_{i}, y_{j}) | 0 ⩽ i ⩽ m, 0 ⩽ j ⩽ n}

$\Omega_h = \{ (x_i, y_j) \, | \, 0 \leqslant i \leqslant m , \, 0 \leqslant j \leqslant n \}$

所有内结点构成集合

{\overset{˚}{Ω}}_{h} = {(x_{i}, y_{j}) | 1 ⩽ i ⩽ m - 1, 1 ⩽ j ⩽ n - 1}

$\mathring{\Omega}_h = \{ (x_i, y_j) \, | \, 1 \leqslant i \leqslant m-1, \, 1 \leqslant j \leqslant n-1\}$

所有边界结点构成集合

Γ_{h} = Ω_{h} ∖ {\overset{˚}{Ω}}_{h}

$\Gamma_h = \Omega_h \setminus \mathring{\Omega}_h$

为方便记

ω = {(i, j) | (x_{i}, y_{j}) \in {\overset{˚}{Ω}}_{h}}, γ = {(i, j) | (x_{i}, y_{j}) \in Γ_{h}},

$\omega = \{ (i, j) \, | \, (x_i, y_j) \in \mathring{\Omega}_h\}, \quad \gamma = \{ (i, j) \, | \, (x_i, y_j) \in \Gamma_h\},$

记

S_{h} = {v | v = {v_{i j} | \forall i, \forall j} 为 Ω_{h} 上的网格函数}

$S_h = \left\{ v \, | \, v = \{ v_{ij} \, | \, \forall i, \forall j \} \text{为} \Omega_h \text{上的网格函数}\right\}$

设 $v = \{ v_{ij} \, | \, 0 \leqslant i \leqslant m , \, 0 \leqslant j \leqslant n \} \in S_h$ ，引进如下记号：

D_{x} v_{i j} = \frac{1}{h_{1}} (v_{i + 1, j} - v_{i j}), D_{\bar{x}} v_{i j} = \frac{1}{h_{1}} (v_{i j} - v_{i - 1, j}),

$D_x v_{ij} = \frac{1}{h_1} (v_{i+1,j} - v_{ij}), \quad D_{\bar{x}} v_{ij} = \frac{1}{h_1} (v_{ij} - v_{i-1, j}),$

D_{y} v_{i j} = \frac{1}{h_{2}} (v_{i, j + 1} - v_{i j}), D_{\bar{y}} v_{i j} = \frac{1}{h_{2}} (v_{i j} - v_{i, j - 1}),

$D_y v_{ij} = \frac{1}{h_2} (v_{i,j+1} - v_{ij}), \quad D_{\bar{y}} v_{ij} = \frac{1}{h_2} (v_{ij} - v_{i, j-1}),$

δ_{x}^{2} v_{i j} = \frac{1}{h_{1}} (D_{x} v_{i j} - D_{\bar{x}} v_{i j}), δ_{y}^{2} v_{i j} = \frac{1}{h^{2}} (D_{y} v_{i j} - D_{\bar{y}} v_{i j})

$\delta_x^2 v_{ij} = \frac{1}{h_1} (D_x v_{ij} - D_{\bar{x}} v_{ij}), \quad \delta_y^2 v_{ij} = \frac{1}{h^2} (D_y v_{ij} - D_{\bar{y}} v_{ij})$

‖ v ‖_{\infty} = max_{\forall i, j} | v_{i j} |

$\| v \|_{\infty} = \max_{\forall i,j} |v_{ij}|$

在结点处考虑边值问题 $(3.1.1)$ ， $(3.1.2)$ 有

\begin{matrix} (3.2.1) & - [\frac{\partial^{2} u}{\partial x^{2}} (x_{i}, y_{j}) + \frac{\partial^{2} u}{\partial y^{2}} (x_{i}, y_{j})] = f (x_{i}, y_{j}), (i, j) \in ω \end{matrix}

$- \left[ \frac{\partial^2 u}{\partial x^2} (x_i, y_j) + \frac{\partial^2 u}{\partial y^2} (x_i, y_j) \right] = f(x_i, y_j), \quad (i,j) \in \omega \tag{3.2.1}$

\begin{matrix} (3.2.2) & u (x_{i}, y_{j}) = φ (x_{i}, y_{j}), (i, j) \in γ \end{matrix}

$u(x_i, y_j) = \varphi(x_i, y_j), \quad (i, j) \in \gamma \tag{3.2.2}$

定义 $\Omega_h$ 上的网格函数

U = {U_{i j} | 0 ⩽ i ⩽ m, 0 ⩽ j ⩽ n}, U_{i j} = u (x_{i}, y_{j})

$U = \{U_{ij} \, | \, 0 \leqslant i \leqslant m , \, 0 \leqslant j \leqslant n \}, \quad U_{ij} = u(x_i, y_j)$

二阶导数用二阶中心差商来近似，由引理 $2.2.1$ 可知

\frac{\partial^{2} u}{\partial x^{2}} (x_{i}, y_{j}) = \frac{1}{h_{1}^{2}} [u (x_{i - 1}, y_{j}) - 2 u (x_{i}, y_{j}) + u (x_{i + 1}, y_{j})] - \frac{h_{1}^{2}}{12} \frac{\partial^{4} u}{\partial x^{4}} (ξ_{i j}, y_{j}) = δ_{x}^{2} U_{i j} - \frac{h_{1}^{2}}{12} \frac{\partial^{4} u}{\partial x^{4}} (ξ_{i j}, y_{j}), x_{i - 1} < ξ_{i j} < x_{i + 1}

$\frac{\partial^2 u}{\partial x^2} (x_i, y_j) = \frac{1}{h_1^2} \left[ u(x_{i-1}, y_j) - 2 u(x_{i}, y_j) + u(x_{i+1}, y_j)\right] - \frac{h_1^2}{12} \frac{\partial^4 u}{\partial x^4}(\xi_{ij}, y_j) = \delta_x^2 U_{ij} - \frac{h_1^2}{12} \frac{\partial^4 u}{\partial x^4}(\xi_{ij}, y_j), \quad x_{i-1} < \xi_{ij} < x_{i+1}$

\frac{\partial^{2} u}{\partial y^{2}} (x_{i}, y_{j}) = \frac{1}{h_{2}^{2}} [u (x_{i}, y_{j - 1}) - 2 u (x_{i}, y_{j}) + u (x_{i}, y_{j + 1})] - \frac{h_{2}^{2}}{12} \frac{\partial^{4} u}{\partial y^{4}} (x_{i}, η_{i j}) = δ_{y}^{2} U_{i j} - \frac{h_{2}^{2}}{12} \frac{\partial^{4} u}{\partial y^{4}} (x_{i}, η_{i j}), y_{j - 1} < η_{i j} < y_{j + 1}

$\frac{\partial^2 u}{\partial y^2} (x_i, y_j) = \frac{1}{h_2^2} \left[ u(x_{i}, y_{j-1}) - 2 u(x_{i}, y_j) + u(x_{i}, y_{j+1})\right] - \frac{h_2^2}{12} \frac{\partial^4 u}{\partial y^4}(x_i, \eta_{ij}) =\delta_y^2 U_{ij} - \frac{h_2^2}{12} \frac{\partial^4 u}{\partial y^4}(x_i, \eta_{ij}), \quad y_{j-1} < \eta_{ij} < y_{j+1}$

代入式 $(3.2.1)$ ，可得

\begin{matrix} (3.2.3) & - (δ_{x}^{2} U_{i j} + δ_{y}^{2} U_{i j}) = f (x_{i}, y_{j}) - \frac{h_{1}^{2}}{12} \frac{\partial^{4} u}{\partial x^{4}} (ξ_{i j}, y_{j}) - \frac{h_{2}^{2}}{12} \frac{\partial^{4} u}{\partial y^{4}} (x_{i}, η_{i j}), (i, j) \in ω \end{matrix}

$-(\delta_x^2 U_{ij} + \delta_y^2 U_{ij}) = f(x_i, y_j) - \frac{h_1^2}{12} \frac{\partial^4 u}{\partial x^4}(\xi_{ij}, y_j) - \frac{h_2^2}{12} \frac{\partial^4 u}{\partial y^4}(x_i, \eta_{ij}), \quad (i, j) \in \omega \tag{3.2.3}$

\begin{matrix} (3.2.4) & U_{i j} = φ (x_{i}, y_{j}), (i, j) \in γ \end{matrix}

$U_{ij} = \varphi(x_i, y_j), \quad (i,j) \in \gamma \tag{3.2.4}$

在上式中略去小量项

\begin{matrix} (3.2.5) & R_{i j} = - \frac{h_{1}^{2}}{12} \frac{\partial^{4} u}{\partial x^{4}} (ξ_{i j}, y_{j}) - \frac{h_{2}^{2}}{12} \frac{\partial^{4} u}{\partial y^{4}} (x_{i}, η_{i j}) \end{matrix}

$R_{ij} = - \frac{h_1^2}{12} \frac{\partial^4 u}{\partial x^4}(\xi_{ij}, y_j) - \frac{h_2^2}{12} \frac{\partial^4 u}{\partial y^4}(x_i, \eta_{ij}) \tag{3.2.5}$

用 $u_{ij}$ 代替 $U_{ij}$ 得到如下差分格式

\begin{matrix} (3.2.6) & - (δ_{x}^{2} u_{i j} + δ_{y}^{2} u_{i j}) = f (x_{i}, y_{j}), (i, j) \in ω \end{matrix}

$-(\delta_x^2 u_{ij} + \delta_y^2 u_{ij}) = f(x_i, y_j), \quad (i,j) \in \omega \tag{3.2.6}$

\begin{matrix} (3.2.7) & u_{i j} = φ (x_{i}, y_{j}), (i, j) \in γ \end{matrix}

$u_{ij} = \varphi(x_i, y_j), \quad (i,j) \in \gamma \tag{3.2.7}$

(2) 局部截断误差

称 $R_{ij}$ 为差分格式 $(3.2.5)$ 的局部截断误差，它表示为差分格式 $(3.2.5)$ 用精确解代替近似解后等式两边之差，即

R_{i j} = - \frac{1}{h_{1}^{2}} [u (x_{i - 1}, y_{j}) - 2 u (x_{i}, y_{j}) + u (x_{i + 1}, y_{j})] - \frac{1}{h_{2}^{2}} [u (x_{i}, y_{j - 1}) - 2 u (x_{i}, y_{j}) + u (x_{i}, y_{j + 1})] - f (x_{i}, y_{j})

$R_{ij} = -\frac{1}{h_1^2} [u(x_{i-1}, y_j) - 2 u(x_i, y_j) + u(x_{i+1}, y_j)] -\frac{1}{h_2^2} [u(x_{i}, y_{j-1}) - 2 u(x_i, y_j) + u(x_{i}, y_{j+1})] - f(x_i, y_j)$

记

\begin{matrix} (3.2.8) & M_{4} = max {max_{(x, y) \in \bar{Ω}} | \frac{\partial^{4} u}{\partial x^{4}} (x, y) |, max_{(x, y) \in \bar{Ω}} | \frac{\partial^{4} u}{\partial y^{4}} (x, y) |} \end{matrix}

$M_4 = \max \left\{ \max_{(x,y) \in \bar{\Omega}} \left| \frac{\partial^4 u}{\partial x^4} (x,y) \right|, \max_{(x,y) \in \bar{\Omega}} \left| \frac{\partial^4 u}{\partial y^4} (x,y) \right| \right\} \tag{3.2.8}$

则有

\begin{matrix} (3.2.9) & | R_{i j} | ⩽ \frac{(h_{1}^{2} + h_{2}^{2})}{12} M_{4}, 1 ⩽ i ⩽ m - 1, 1 ⩽ j ⩽ n - 1 \end{matrix}

$|R_{ij}| \leqslant \frac{(h_1^2 + h_2^2)}{12} M_4, \quad 1 \leqslant i \leqslant m-1, \quad 1 \leqslant j \leqslant n-1 \tag{3.2.9}$

局部截断误差可以反映差分格式对原方程的近似程度，由上式显然可知当 $M_4 < \infty$ 时，差分格式 $(3.2.5)$ 与微分方程 $(3.1.1)$ 是相容的。

3.2.2 差分格式的求解

差分格式 $(3.2.6)$ ， $(3.2.7)$ 是以 $(u_{ij})$ 可视为未知量的线性方程组，式 $(3.2.6)$ 可改写为

\begin{matrix} (3.2.10) & - \frac{1}{h_{2}^{2}} u_{i, j - 1} - \frac{1}{h_{1}^{2}} u_{i - 1, j} + 2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}}) u_{i j} - \frac{1}{h_{1}^{2}} u_{i + 1, j} - \frac{1}{h_{2}^{2}} u_{i, j + 1} = f (x_{i}, y_{j}), (i, j) \in ω \end{matrix}

$-\frac{1}{h_2^2} u_{i, j-1} - \frac{1}{h_1^2} u_{i-1, j} + 2 \left( \frac{1}{h_1^2} + \frac{1}{h_2^2} \right) u_{ij} - \frac{1}{h_1^2} u_{i+1, j} - \frac{1}{h_2^2} u_{i, j+1} = f(x_i, y_j), \quad (i,j) \in \omega \tag{3.2.10}$

我们考虑将原始的差分格式写成两种形式：线性方程组与矩阵方程。

(1) 线性方程组形式

记

\begin{aligned} u_{j} = [\begin{matrix} u_{1 j} \\ u_{2 j} \\ ⋮ \\ u_{m - 1, j} \end{matrix}], 0 ⩽ j ⩽ n \end{aligned}

$\begin{align*} \bm{u}_j = \begin{bmatrix} u_{1j} \\ u_{2j} \\ \vdots \\ u_{m-1,j} \end{bmatrix}, \quad 0 \leqslant j \leqslant n \end{align*}$

利用式 $(3.2.7)$ 可将式 $(3.2.10)$ 改写为

\begin{matrix} (3.2.11) & D u_{j - 1} + C u_{j} + D u_{j + 1} = f_{j}, 1 ⩽ j ⩽ n - 1 \end{matrix}

$D \bm{u}_{j-1} + C \bm{u}_j + D \bm{u}_{j+1} = \bm{f}_j, \quad 1 \leqslant j \leqslant n-1 \tag{3.2.11}$

其中，

\begin{aligned} C = [\begin{matrix} 2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}}) & - \frac{1}{h_{1}^{2}} \\ - \frac{1}{h_{1}^{2}} & 2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}}) & - \frac{1}{h_{1}^{2}} \\ ⋱ & ⋱ & ⋱ \\ - \frac{1}{h_{1}^{2}} & 2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}}) & - \frac{1}{h_{1}^{2}} \\ - \frac{1}{h_{1}^{2}} & 2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}}) \end{matrix}] \end{aligned}

$\begin{align*} C = \begin{bmatrix} 2 \left( \frac{1}{h_1^2} + \frac{1}{h_2^2} \right) & -\frac{1}{h_1^2} \\ -\frac{1}{h_1^2} & 2 \left( \frac{1}{h_1^2} + \frac{1}{h_2^2} \right) & -\frac{1}{h_1^2} \\ & \ddots & \ddots & \ddots \\ & & -\frac{1}{h_1^2} & 2 \left( \frac{1}{h_1^2} + \frac{1}{h_2^2} \right) & -\frac{1}{h_1^2} \\ & & & -\frac{1}{h_1^2} & 2 \left( \frac{1}{h_1^2} + \frac{1}{h_2^2} \right) \end{bmatrix} \end{align*}$

\begin{aligned} D = [\begin{matrix} - \frac{1}{h_{2}^{2}} \\ - \frac{1}{h_{2}^{2}} \\ ⋱ \\ - \frac{1}{h_{2}^{2}} \\ - \frac{1}{h_{2}^{2}} \end{matrix}], f_{j} = [\begin{matrix} f (x_{1}, y_{j}) + \frac{1}{h_{1}^{2}} φ (x_{0}, y_{j}) \\ f (x_{2}, y_{j}) \\ ⋮ \\ f (x_{m - 2}, y_{j}) \\ f (x_{m - 1}, y_{j}) + \frac{1}{h_{1}^{2}} φ (x_{m}, y_{j}) \end{matrix}] \end{aligned}

$\begin{align*} D = \begin{bmatrix} -\frac{1}{h_2^2} \\ & -\frac{1}{h_2^2} \\ & & \ddots \\ & & & -\frac{1}{h_2^2} \\ & & & &-\frac{1}{h_2^2} \end{bmatrix}, \quad \bm{f}_j = \begin{bmatrix} f(x_1, y_j) + \frac{1}{h_1^2} \varphi(x_0, y_j) \\ f(x_2, y_j) \\ \vdots \\ f(x_{m-2}, y_j) \\ f(x_{m-1}, y_j) + \frac{1}{h_1^2} \varphi(x_m, y_j) \end{bmatrix} \end{align*}$

将式 $(3.2.11)$ 进一步改写成

\begin{matrix} (3.2.12) & \begin{aligned} [\begin{matrix} C & D \\ D & C & D \\ ⋱ & ⋱ & ⋱ \\ D & C & D \\ D & C \end{matrix}] [\begin{matrix} u_{1} \\ u_{2} \\ ⋮ \\ u_{n - 2} \\ u_{n - 1} \end{matrix}] = [\begin{matrix} f_{1} - D u_{0} \\ f_{2} \\ ⋮ \\ f_{n - 2} \\ f_{n - 1} - D u_{n} \end{matrix}] \end{aligned} \end{matrix}

$\begin{align*} \begin{bmatrix} C & D \\ D & C & D \\ & \ddots & \ddots & \ddots \\ & & D & C & D \\ & & & D & C \end{bmatrix} \begin{bmatrix} \bm{u}_1 \\ \bm{u}_2 \\ \vdots \\ \bm{u}_{n-2} \\ \bm{u}_{n-1} \end{bmatrix} = \begin{bmatrix} \bm{f}_1 - D \bm{u}_0 \\ \bm{f}_2 \\ \vdots \\ \bm{f}_{n-2} \\ \bm{f}_{n-1} - D \bm{u}_n \end{bmatrix} \end{align*} \tag{3.2.12}$

上述线性方程组的系数矩阵是一个三对角块矩阵，每一行至多有 $5$ 个非零元素，通常称这种绝大多数元素为 $0$ 的矩阵为稀疏矩阵。常用迭代法求解以大型稀疏矩阵为系数矩阵的线性方程组。

(2) 矩阵方程形式

可以将式 $(3.2.6)$ 写成

\begin{aligned} \frac{1}{h_{1}^{2}} [\begin{matrix} - 1 & 2 & - 1 \end{matrix}] [\begin{matrix} u_{i - 1, j} \\ u_{i j} \\ u_{i + 1, j} \end{matrix}] + \frac{1}{h_{2}^{2}} [\begin{matrix} u_{i, j - 1} & u_{i j} & u_{i, j + 1} \end{matrix}] [\begin{matrix} - 1 \\ 2 \\ - 1 \end{matrix}] \end{aligned} = f (x_{i}, y_{j})

$\begin{align*} \frac{1}{h_1^2} \begin{bmatrix} -1 & 2 & -1 \end{bmatrix} \begin{bmatrix} u_{i-1,j} \\ u_{ij} \\ u_{i+1, j} \end{bmatrix} + \frac{1}{h_2^2} \begin{bmatrix} u_{i,j-1} & u_{ij} & u_{i, j+1} \end{bmatrix} \begin{bmatrix} -1 \\ 2 \\ -1 \end{bmatrix} \end{align*} = f(x_i, y_j)$

进一步写成

\begin{matrix} (3.2.13) & \frac{1}{h_{1}^{2}} A_{m - 1} X + \frac{1}{h_{2}^{2}} X A_{n - 1} = f + \frac{1}{h_{1}^{2}} u_{u b} + \frac{1}{h_{2}^{2}} u_{l r} \end{matrix}

$\frac{1}{h_1^2} A_{m-1} X + \frac{1}{h_2^2} X A_{n-1} = f + \frac{1}{h_1^2} u_{ub} + \frac{1}{h_2^2} u_{lr} \tag{3.2.13}$

其中，

\begin{aligned} A_{k} = {[\begin{matrix} 2 & - 1 \\ - 1 & 2 & - 1 \\ ⋱ & ⋱ & ⋱ \\ - 1 & 2 & - 1 \\ - 1 & 2 \end{matrix}]}_{k \times k} X = [\begin{matrix} u_{11} & u_{12} & \dots & u_{1, n - 1} \\ u_{21} & u_{22} & \dots & u_{2, n - 1} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ u_{m - 1, 1} & u_{m - 1, 2} & \dots & u_{m - 1, n - 1} \end{matrix}] \end{aligned}

$\begin{align*} A_k = \begin{bmatrix} 2 & -1 \\ -1 & 2 & -1 \\ & \ddots & \ddots & \ddots \\ & & -1 & 2 & -1 \\ & & & -1 & 2 \end{bmatrix}_{k \times k} \quad X = \begin{bmatrix} u_{11} & u_{12} & \cdots & u_{1,n-1} \\ u_{21} & u_{22} & \cdots & u_{2,n-1} \\ \vdots & \vdots & \ddots & \vdots \\ u_{m-1,1} & u_{m-1,2} & \cdots & u_{m-1,n-1} \end{bmatrix} \end{align*}$

\begin{aligned} u_{l r} = [\begin{matrix} u_{10} & 0 & \dots & 0 & u_{1, n} \\ u_{20} & 0 & \dots & 0 & u_{2, n} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ u_{m - 1, 0} & 0 & \dots & 0 & u_{m - 1, n} \end{matrix}], u_{u b} = [\begin{matrix} u_{01} & u_{02} & \dots & u_{0, n - 1} \\ 0 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 0 \\ u_{m 1} & u_{m 2} & \dots & u_{m, n - 1} \end{matrix}] \end{aligned}

$\begin{align*} u_{lr} = \begin{bmatrix} u_{10} & 0 & \cdots & 0 & u_{1,n} \\ u_{20} & 0 & \cdots & 0 & u_{2,n} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ u_{m-1,0} & 0 & \cdots & 0 & u_{m-1,n} \end{bmatrix}, \quad u_{ub} = \begin{bmatrix} u_{01} & u_{02} & \cdots & u_{0,n-1} \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \\ u_{m1} & u_{m2} & \cdots & u_{m,n-1} \end{bmatrix} \end{align*}$

\begin{aligned} f = [\begin{matrix} f_{11} & f_{12} & \dots & f_{1, n - 1} \\ f_{21} & f_{22} & \dots & f_{2, n - 1} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ f_{m - 1, 1} & f_{m - 1, 2} & \dots & f_{m - 1, n - 1} \end{matrix}], f_{i j} = f (x_{i}, y_{j}) \end{aligned}

$\begin{align*} f = \begin{bmatrix} f_{11} & f_{12} & \cdots & f_{1,n-1} \\ f_{21} & f_{22} & \cdots & f_{2,n-1} \\ \vdots & \vdots & \ddots & \vdots \\ f_{m-1,1} & f_{m-1,2} & \cdots & f_{m-1,n-1} \end{bmatrix}, \quad f_{ij} = f(x_i, y_j) \end{align*}$

式 $(3.2.13)$ 可写成矩阵方程形式

\begin{matrix} (3.2.14) & A X + X B = C \end{matrix}

$AX + XB = C \tag{3.2.14}$

A = \frac{1}{h_{1}^{2}} A_{m - 1}, B = \frac{1}{h_{2}^{2}} A_{n - 1}, C = f + \frac{1}{h_{1}^{2}} u_{u b} + \frac{1}{h_{2}^{2}} u_{l r}

$A = \frac{1}{h_1^2} A_{m-1}, \quad B = \frac{1}{h_2^2} A_{n-1}, \quad C = f + \frac{1}{h_1^2} u_{ub} + \frac{1}{h_2^2} u_{lr}$

(3) 性质

记式 $(3.2.12)$ 的系数矩阵为 $G$ ，即

\begin{aligned} G = [\begin{matrix} C & D \\ D & C & D \\ ⋱ & ⋱ & ⋱ \\ D & C & D \\ D & C \end{matrix}] \end{aligned}

$\begin{align*} G = \begin{bmatrix} C & D \\ D & C & D \\ & \ddots & \ddots & \ddots \\ & & D & C & D \\ & & & D & C \end{bmatrix} \end{align*}$

我们可以发现，线性方程组与矩阵方程是等价的，它们之间具体的关系可用下述定理描述。

定理 3.2.1

对式 $(3.2.11)$ 有

$AX + XB = C \iff \left[ \text{Kron} (I, A) + \text{Kron} (B^T, I) \right] \text{vec}(X) = \text{vec}(C)$
其中第一个单位矩阵的阶数等于矩阵 $X$ 的列数；第二个单位矩阵的阶数等于矩阵 $X$ 的行数。 $\text{vec}$ 表示将矩阵用列拉直的方式变为列向量。 $\text{Kron}$ 为两矩阵的 $\text{kronecker}$ 积。

表示为线性方程组时，可以发现系数矩阵不仅是一个稀疏矩阵，还具有良好的代数性质。

定理 3.2.2

式 $(3.2.9)$ 的系数矩阵 $G$ 为对称正定矩阵。

证明：系数矩阵 $G$ 显然是一个对称矩阵。只需再证明该矩阵正定。

证毕。

(4) 求解的数值算法

由 $(3.2.9)$ 的线性方程组形式可以看出，该线性方程组的系数矩阵是一个三对角块矩阵，每一行至多有 $5$ 个非零元素，通常称这种绝大多数元素为零的矩阵为稀疏矩阵。常用迭代法求解以大型稀疏矩阵为系数矩阵的线性方程组。

$\text{Jacobi}$ 迭代法 对 $k=0,1,2,\cdots$ ，计算

u_{i j}^{(k + 1)} = [f (x_{i}, y_{j}) + \frac{1}{h_{2}^{2}} u_{i, j - 1}^{(k)} + \frac{1}{h_{1}^{2}} u_{i - 1, j}^{(k)} + \frac{1}{h_{1}^{2}} u_{i + 1, j}^{(k)} + \frac{1}{h_{2}^{2}} u_{i, j + 1}^{(k)}] / [2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}})], i = 1, \dots, m - 1, j = 1, \dots, n - 1

$u_{ij}^{(k+1)} = \left[ f(x_i, y_j) + \frac{1}{h_2^2} u_{i, j-1}^{(k)} + \frac{1}{h_1^2} u_{i-1, j}^{(k)} + \frac{1}{h_1^2} u_{i+1, j}^{(k)} + \frac{1}{h_2^2} u_{i, j+1}^{(k)} \right] \Big/ \left[ 2\left( \frac{1}{h_1^2} + \frac{1}{h_2^2}\right) \right], \quad i=1,\cdots,m-1, \quad j=1,\cdots,n-1$

$\text{Gauss-Seidel}$ 迭代法 对 $k=0,1,2,\cdots$ ，计算

u_{i j}^{(k + 1)} = [f (x_{i}, y_{j}) + \frac{1}{h_{2}^{2}} u_{i, j - 1}^{(k + 1)} + \frac{1}{h_{1}^{2}} u_{i - 1, j}^{(k + 1)} + \frac{1}{h_{1}^{2}} u_{i + 1, j}^{(k)} + \frac{1}{h_{2}^{2}} u_{i, j + 1}^{(k)}] / [2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}})], i = 1, \dots, m - 1, j = 1, \dots, n - 1

$u_{ij}^{(k+1)} = \left[ f(x_i, y_j) + \frac{1}{h_2^2} u_{i, j-1}^{(k+1)} + \frac{1}{h_1^2} u_{i-1, j}^{(k+1)} + \frac{1}{h_1^2} u_{i+1, j}^{(k)} + \frac{1}{h_2^2} u_{i, j+1}^{(k)} \right] \Big/ \left[ 2\left( \frac{1}{h_1^2} + \frac{1}{h_2^2}\right) \right], \quad i=1,\cdots,m-1, \quad j=1,\cdots,n-1$

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

我们可以应用极值原理来给出差分方程组的解的先验估计式。

设 $u = \{u_{ij} \, | \, 0 \leqslant i \leqslant m, 0 \leqslant j \leqslant n\}$ 为 $\Omega_h$ 上的网格函数，简记

(L_{h} u)_{i j} = - (δ_{x}^{2} u_{i j} + δ_{y}^{2} u_{i j}), (i, j) \in ω

$(L_h u)_{ij} = -(\delta_x^2 u_{ij} + \delta_y^2 u_{ij}), \quad (i,j) \in \omega$

定理 3.2.3

设 $u = \{u_{ij} \, | \, 0 \leqslant i \leqslant m, 0 \leqslant j \leqslant n\}$ 为 $\Omega_h$ 上的网格函数，如果我们有如下式子成立

$(L_h u)_{ij} \leqslant 0, \quad (i,j) \in \omega$
则我们有

$\max_{(i,j) \in w} u_{ij} \leqslant \max_{(i,j) \in \gamma} u_{ij}$

证明：用反证法。设

max_{(i, j) \in w} u_{i j} > max_{(i, j) \in γ} u_{i j}

$\max_{(i,j) \in w} u_{ij} > \max_{(i,j) \in \gamma} u_{ij}$

令 $\max_{(i,j) \in w} u_{ij} = M$ ，则一定存在 $(i_0, j_0) \in \omega$ 使得 $u_{i_0, j_0} = M$ ，且 $u_{i_0 - 1, \,j_0}$ 、 $u_{i_0 + 1, \,j_0}$ 、 $u_{i_0, \,j_0-1}$ 、 $u_{i_0, \,j_0+1}$ 中至少有一个值小于 $M$ 。

因此我们可得

\begin{aligned} (L_{h} u)_{i_{0}, j_{0}} & = 2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}}) u_{i_{0}, j_{0}} - \frac{1}{h_{1}^{2}} (u_{i_{0} - 1, j_{0}} + u_{i_{0} + 1, j_{0}}) - \frac{1}{h_{2}^{2}} (u_{i_{0}, j_{0} - 1} + u_{i_{0}, j_{0} + 1}) \\ > 2 (\frac{1}{h_{1}^{2}} + \frac{1}{h_{2}^{2}}) M - \frac{1}{h_{1}^{2}} (M + M) - \frac{1}{h_{2}^{2}} (M + M) = 0 \end{aligned}

$\begin{align*} (L_h u)_{i_0, \, j_0} & = 2 \left( \frac{1}{h_1^2} + \frac{1}{h_2^2} \right) u_{i_0, \, j_0} - \frac{1}{h_1^2} \left( u_{i_0 - 1, \,j_0} + u_{i_0 + 1, \,j_0} \right) - \frac{1}{h_2^2} \left( u_{i_0, \,j_0-1} + u_{i_0, \,j_0+1} \right) \\ & > 2 \left( \frac{1}{h_1^2} + \frac{1}{h_2^2} \right) M - \frac{1}{h_1^2} (M + M) - \frac{1}{h_2^2} (M + M) = 0 \end{align*}$

与条件矛盾。

证毕。

下面利用上述的极值原理给出差分格式解的先验估计式。

定理 3.2.4

给定差分方程组

$-(\delta_x^2 u_{ij} + \delta_y^2 u_{ij}) = f_{ij}, \quad (i,j) \in \omega$
$u_{ij} = \varphi_{ij}, \quad (i,j) \in \gamma$
设 $\{u_{ij}\}$ 为上述差分方程组的解，则有

$\max_{(i,j) \in \omega} | u_{ij} | \leqslant \max_{(i,j) \in \gamma} | \varphi_{ij} | + \frac{1}{4} \left[ \left( \frac{b-a}{2}\right)^2 + \left( \frac{d-c}{2} \right)^2 \right] \max_{(i,j)\in \omega} |f_{ij}|$
上式中 $\left( \frac{b-a}{2}\right)^2 + \left( \frac{d-c}{2} \right)^2$ 为 $\Omega$ 的外接圆半径的平方。

证明：我们记

C = max_{(i, j) \in ω} | f_{i j} |

$C = \max_{(i,j) \in \omega} |f_{ij}|$

证毕。

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

定理 3.2.5

差分格式 $(3.2.6) \sim (3.2.7)$ 的解存在且唯一。

证明：由于差分格式 $(3.2.6)$ 与 $(3.2.7)$ 是线性的，考虑其齐次方程组

\begin{matrix} (3.2.15) & - (δ_{x}^{2} u_{i j} + δ_{y}^{2} u_{i j}) = 0, (i, j) \in ω \end{matrix}

$-(\delta_x^2 u_{ij} + \delta_y^2 u_{ij}) = 0, \quad (i,j) \in \omega \tag{3.2.15}$

\begin{matrix} (3.2.16) & u_{i j} = 0, (i, j) \in γ \end{matrix}

$u_{ij} = 0, \quad (i,j) \in \gamma \tag{3.2.16}$

设 $\|u\|_{\infty} = M > 0$ ，则由 $(3.2.16)$ 知，存在 $(i,j) \in \omega$ 使得 $|u_{i_0, \, j_0}| = M$ ，且 $|u_{i_0-1, \, j_0}|$ 、 $|u_{i_0+1, \, j_0}|$ 、 $|u_{i_0, \, j_0-1}|$ 、 $|u_{i_0, \, j_0+1}|$ 中至少有一个小于 $M$ ，考虑式 $(3.2.15)$ 中 $(i,j) = (i_0,j_0)$ 的等式，有

(\frac{2}{h_{1}^{2}} + \frac{2}{h_{2}^{2}}) u_{i_{0}, j_{0}} = \frac{1}{h_{1}^{2}} (u_{i_{0} - 1, j_{0}} + u_{i_{0} + 1, j_{0}}) + \frac{1}{h_{2}^{2}} (u_{i_{0}, j_{0} - 1} + u_{i_{0}, j_{0} + 1})

$\left( \frac{2}{h_1^2} + \frac{2}{h_2^2} \right) u_{i_0, \, j_0} = \frac{1}{h_1^2} (u_{i_0 - 1, \, j_0} + u_{i_0 + 1, \, j_0}) + \frac{1}{h_2^2} \left( u_{i_0, \, j_0-1} + u_{i_0, \, j_0+1} \right)$

将上式两边取绝对值可得

(\frac{2}{h_{1}^{2}} + \frac{2}{h_{2}^{2}}) M ⩽ \frac{1}{h_{1}^{2}} (| u_{i_{0} - 1, j_{0}} | + | u_{i_{0} + 1, j_{0}} |) + \frac{1}{h_{2}^{2}} (| u_{i_{0}, j_{0} - 1} | + | u_{i_{0}, j_{0} + 1} |) < (\frac{2}{h_{1}^{2}} + \frac{2}{h_{2}^{2}}) M

$\left( \frac{2}{h_1^2} + \frac{2}{h_2^2} \right) M \leqslant \frac{1}{h_1^2} ( |u_{i_0 - 1, \, j_0}| + |u_{i_0 + 1, \, j_0}| ) + \frac{1}{h_2^2} \left( |u_{i_0, \, j_0-1}| + |u_{i_0, \, j_0+1}| \right) < \left( \frac{2}{h_1^2} + \frac{2}{h_2^2} \right) M$

矛盾。故 $M=0$ ，因而差分格式的解是存在且唯一的。

证毕。

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

(1) 收敛性

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

定理 3.2.6

设 $\{u(x,y) \, | \, a \leqslant x \leqslant b, \, c \leqslant y \leqslant d\}$ 是定解问题 $(3.1.1)$ 和 $(3.1.2)$ 的解， $\{u_{ij} \, | \, 0 \leqslant i \leqslant m, \, 0 \leqslant j \leqslant n\}$ 为差分格式 $(3.2.5)$ 和 $(3.2.6)$ 的解，则有

$\max_{(i,j) \in \omega} \left| u(x_i, y_j) - u_{ij} \right| \leqslant \frac{M_4}{48} \left[ \left( \frac{b-a}{2} \right)^2 + \left( \frac{d-c}{2} \right)^2 \right] (h_1^2 + h_2^2)$
其中 $M_4$ 由式 $(3.2.8)$ 定义。

证明：我们记

e_{i j} = u (x_{i}, y_{j}) - u_{i j}, (i, j) \in ω \cup γ

$e_{ij} = u(x_i, y_j) - u_{ij}, \quad (i,j) \in \omega \cup \gamma$

将式 $(3.2.3)$ 、 $(3.2.4)$ 分别与式 $(3.2.6)$ 、 $(3.2.7)$ 相减，得误差方程为

\begin{matrix} (3.2.17) & - (δ_{x}^{2} e_{i j} + δ_{y}^{2} e_{i j}) = R_{i j}, (i, j) \in ω \end{matrix}

$-(\delta_x^2 e_{ij} + \delta_y^2 e_{ij}) = R_{ij}, \quad (i,j) \in \omega \tag{3.2.17}$

\begin{matrix} (3.2.18) & e_{i j} = 0, (i, j) \in γ \end{matrix}

$e_{ij} = 0, \quad (i,j) \in \gamma \tag{3.2.18}$

其中 $R_{ij}$ 由式 $(3.2.5)$ 定义，由式 $(3.2.9)$ 可知，有

max_{(i, j) \in ω} | R_{i j} | ⩽ \frac{M_{4}}{12} (h_{1}^{2} + h_{2}^{2})

$\max_{(i,j) \in \omega} |R_{ij}| \leqslant \frac{M_4}{12} (h_1^2 + h_2^2)$

应用定理 $3.2.4$ 可得

max_{(i, j) \in ω} | e_{i j} | ⩽ \frac{1}{4} [{(\frac{b - a}{2})}^{2} + {(\frac{d - c}{2})}^{2}] max_{(i, j) \in ω} | R_{i j} | ⩽ \frac{M_{4}}{48} [{(\frac{b - a}{2})}^{2} + {(\frac{d - c}{2})}^{2}] (h_{1}^{2} + h_{2}^{2})

$\max_{(i,j) \in \omega} |e_{ij}| \leqslant \frac{1}{4} \left[ \left( \frac{b-a}{2}\right)^2 + \left( \frac{d-c}{2} \right)^2 \right] \max_{(i,j) \in \omega} |R_{ij}| \leqslant \frac{M_4}{48} \left[ \left( \frac{b-a}{2} \right)^2 + \left( \frac{d-c}{2} \right)^2 \right] (h_1^2 + h_2^2)$

证毕。

(2) 稳定性

假设在应用差分格式 $(3.2.6)$ 和 $(3.2.7)$ 时，计算右端函数 $f(x_i, y_j)$ 有误差 $f_{ij}$ ，计算边界值有误差 $\varphi_{ij}$ ，考虑如下差分格式

\begin{matrix} (3.2.19) & - (δ_{x}^{2} v_{i j} + δ_{y}^{2} v_{i j}) = f (x_{i}, y_{j}) + f_{i j}, (i, j) \in ω \end{matrix}

$-(\delta_x^2 v_{ij} + \delta_y^2 v_{ij}) = f(x_i, y_j) + f_{ij}, \quad (i,j) \in \omega \tag{3.2.19}$

\begin{matrix} (3.2.20) & v_{i j} = φ (x_{i}, y_{j}) + φ_{i j}, (i, j) \in γ \end{matrix}

$v_{ij} = \varphi(x_i, y_j) + \varphi_{ij}, \quad (i,j) \in \gamma \tag{3.2.20}$

设 $\{v_{ij}\}$ 为上述差分格式的解，我们记

ε_{i j} = v_{i j} - u_{i j}, (i, j) \in ω \cup γ

$\varepsilon_{ij} = v_{ij} - u_{ij}, \quad (i,j) \in \omega \cup \gamma$

将式 $(3.2.19)$ 、 $(3.2.20)$ 与式 $(3.2.6)$ 、 $(3.2.7)$ 相减，得到

\begin{matrix} (3.2.21) & - (δ_{x}^{2} ε_{i j} + δ_{y}^{2} ε_{i j}) = f_{i j}, (i, j) \in ω \end{matrix}

$-(\delta_x^2 \varepsilon_{ij} + \delta_y^2 \varepsilon_{ij}) = f_{ij}, \quad (i,j) \in \omega \tag{3.2.21}$

\begin{matrix} (3.2.22) & ε_{i j} = φ_{i j}, (i, j) \in γ \end{matrix}

$\varepsilon_{ij} = \varphi_{ij}, \quad (i,j) \in \gamma \tag{3.2.22}$

应用定理 $3.2.4$ 可得

max_{(i, j) \in ω} | ε_{i j} | ⩽ max_{(i, j) \in γ} | φ_{i j} | + \frac{1}{4} [{(\frac{b - a}{2})}^{2} + {(\frac{d - c}{2})}^{2}] max_{(i, j) \in ω} | f_{i j} |

$\max_{(i,j) \in \omega} |\varepsilon_{ij}| \leqslant \max_{(i,j) \in \gamma} |\varphi_{ij}| + \frac{1}{4} \left[ \left( \frac{b-a}{2}\right)^2 + \left( \frac{d-c}{2} \right)^2 \right] \max_{(i,j) \in \omega} |f_{ij}|$

由上式可知，当 $|\varphi_{ij}|$ 和 $|f_{ij}|$ 是小量时， $\max_{(i,j) \in \omega} |\varepsilon_{ij}|$ 也为小量。此时我们称差分格式 $(3.2.6)$ 、 $(3.2.7)$ 关于边界值和右端函数是稳定的。

同时我们称得到的式 $(3.2.21)$ 、 $(3.2.22)$ 为摄动误差方程组，可见它的形式与差分格式 $(3.2.6)$ 、 $(3.2.7)$ 是一样的，综合以上，我们可以得到如下定理。

定理 3.2.7

差分格式 $(3.2.6)$ 、 $(3.2.7)$ 在下述意义下关于边界值和右端函数是稳定的，设 $\{u_{ij}\}$ 为

$-(\delta_x^2 u_{ij} + \delta_y^2 u_{ij}) = f_{ij}, \quad (i,j) \in \omega$
$u_{ij} = \varphi_{ij}, \quad (i,j) \in \gamma$
的解，则有

$\max_{(i,j) \in \omega} |u_{ij}| \leqslant \max_{(i,j) \in \gamma} |\varphi_{ij}| + \frac{1}{4} \left[ \left( \frac{b-a}{2}\right)^2 + \left( \frac{d-c}{2} \right)^2 \right] \max_{(i,j) \in \omega} |f_{ij}|$

证明：见上段关于摄动误差方程组的叙述。

证毕。

3.2.6 Richardson 外推法（待定，不知道是否该写这节）

？？

posted on 2022-03-28 10:09 Black_x 阅读(1899) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型，支持深度思考和联网搜索！
· 基于 Docker 搭建 FRP 内网穿透开源项目（很简单哒）
· 25岁的心里话
· ollama系列01：轻松3步本地部署deepseek，普通电脑可用
· 按钮权限的设计及实现

Black_x