块迭代法解离散Poisson方程的五点差分格式线性系统
我将为您讲解如何使用块迭代法（特别是块Jacobi和块Gauss-Seidel）求解由二维Poisson方程离散化产生的特殊线性系统。这个系统在科学计算中非常常见，具有块三对角结构。
一、问题背景与描述
考虑二维Poisson方程：
\[-\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) = f(x, y), \quad (x, y) \in \Omega = (0,1) \times (0,1) \]
在边界 \(\partial \Omega\) 上给定Dirichlet边界条件。
在矩形区域上使用均匀网格进行离散：
- \(x\)方向：网格点 \(x_i = i \cdot h\)，\(i = 0,1,\dots,m+1\)，网格间距 \(h = \frac{1}{m+1}\)
- \(y\)方向：网格点 \(y_j = j \cdot h\)，\(j = 0,1,\dots,n+1\)（为简化常取 \(n = m\)）
使用五点中心差分格式，在内部点 \((x_i, y_j)\) 处，离散方程为：
\[\frac{4u_{i,j} - u_{i-1,j} - u_{i+1,j} - u_{i,j-1} - u_{i,j+1}}{h^2} = f_{i,j} \]
其中 \(u_{i,j} \approx u(x_i, y_j)\)，\(f_{i,j} = f(x_i, y_j)\)。
线性系统形式：
将所有内部点按字典序排列（先固定 \(j\)，遍历所有 \(i\)），得到一个大型稀疏线性系统：
\[A\mathbf{u} = \mathbf{b} \]
其中：
- \(\mathbf{u}\) 是未知向量，包含所有内部点的 \(u_{i,j}\)
- \(A\) 是 \(N \times N\) 矩阵（\(N = m \times n\)）
- \(A\) 具有块三对角结构：主对角块为 \(m \times m\) 的三对角矩阵 \(T\)，次对角块为单位矩阵 \(I\) 的倍数
具体地，将网格按 \(y\) 方向（行方向）分块，则 \(A\) 可写为：
\[A = \begin{bmatrix} T & -I & & & \\ -I & T & -I & & \\ & \ddots & \ddots & \ddots & \\ & & -I & T & -I \\ & & & -I & T \end{bmatrix} \]
其中 \(T\) 是 \(m \times m\) 的三对角矩阵：
\[T = \begin{bmatrix} 4 & -1 & & & \\ -1 & 4 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 4 & -1 \\ & & & -1 & 4 \end{bmatrix} \]
\(I\) 是 \(m \times m\) 的单位矩阵。系数 \(h^{-2}\) 已被吸收到右端项 \(\mathbf{b}\) 中。
二、块迭代法的基本思想
传统点迭代法（如Jacobi、Gauss-Seidel）每次更新一个未知数 \(u_{i,j}\)。块迭代法则将未知数分组，每次更新一组未知数（一个子向量）。对于Poisson方程，自然的分组方式是：将同一行 \(y = y_j\) 上的所有 \(m\) 个未知数 \(u_{1,j}, u_{2,j}, \dots, u_{m,j}\) 作为一个块。
设 \(\mathbf{u}_j = [u_{1,j}, u_{2,j}, \dots, u_{m,j}]^\top\)，则线性系统可重写为：
\[\begin{aligned} T\mathbf{u}_1 - I\mathbf{u}_2 &= \mathbf{b}_1 \\ -I\mathbf{u}_{j-1} + T\mathbf{u}_j - I\mathbf{u}_{j+1} &= \mathbf{b}_j, \quad j = 2,\dots,n-1 \\ -I\mathbf{u}_{n-1} + T\mathbf{u}_n &= \mathbf{b}_n \end{aligned} \]
其中 \(\mathbf{b}_j\) 是第 \(j\) 行对应的右端项（包含 \(f_{i,j}\) 和边界条件贡献）。
三、块Jacobi迭代法
算法步骤：
-
初始化：选取初始猜测 \(\mathbf{u}^{(0)} = [ (\mathbf{u}_1^{(0)})^\top, \dots, (\mathbf{u}_n^{(0)})^\top ]^\top\)，常取零向量或某种插值。
-
迭代格式：
对于第 \(k+1\) 次迭代，解 \(n\) 个独立的 \(m \times m\) 三对角线性系统：
\[ T \mathbf{u}_j^{(k+1)} = \mathbf{b}_j + \mathbf{u}_{j-1}^{(k)} + \mathbf{u}_{j+1}^{(k)}, \quad j = 1,\dots,n \]
其中约定 \(\mathbf{u}_0^{(k)} = \mathbf{u}_{n+1}^{(k)} = \mathbf{0}\)（对应于边界条件）。
-
矩阵形式解释：
将 \(A\) 分解为 \(A = D_B + L_B + U_B\)，其中：- \(D_B = \text{blockdiag}(T, T, \dots, T)\) 是块对角矩阵
- \(L_B\) 是块下三角部分（仅包含 \(-I\) 块）
- \(U_B\) 是块上三角部分
块Jacobi迭代的矩阵形式为：
\[ D_B \mathbf{u}^{(k+1)} = \mathbf{b} - (L_B + U_B) \mathbf{u}^{(k)} \]
由于 \(D_B\) 是块对角矩阵，这相当于对每个块 \(j\) 独立求解 \(T \mathbf{u}_j^{(k+1)} = \text{右端项}\)。
-
求解子问题：
每个子问题 \(T \mathbf{u}_j = \mathbf{r}_j\) 是一个三对角系统，可用高效且稳定的追赶法（Thomas算法）在 \(O(m)\) 时间内求解，总成本为 \(O(mn) = O(N)\) 每迭代步。 -
收敛性：
- 块Jacobi迭代的收敛速度通常比点Jacobi快，因为每次更新利用了同一行内未知数之间的强耦合（通过 \(T\)）。
- 收敛条件依赖于 \(T\) 的特征值。由于 \(T\) 对称正定，且块Jacobi迭代矩阵的谱半径小于点Jacobi，故收敛更快。
- 但与点Jacobi类似，收敛速度仍较慢，对于大规模问题需要较多迭代步。
四、块Gauss-Seidel迭代法
算法步骤：
- 迭代格式：
按照 \(j = 1,2,\dots,n\) 的顺序更新块，每次使用最新可用的相邻块值：
\[ T \mathbf{u}_j^{(k+1)} = \mathbf{b}_j + \mathbf{u}_{j-1}^{(k+1)} + \mathbf{u}_{j+1}^{(k)}, \quad j = 1,\dots,n \]
注意：
- 当 \(j=1\) 时，\(\mathbf{u}_{0}^{(k+1)} = \mathbf{0}\)
- 当 \(j=n\) 时，\(\mathbf{u}_{n+1}^{(k)} = \mathbf{0}\)
- 矩阵形式解释：
块Gauss-Seidel对应于分裂 \(A = (D_B + L_B) + U_B\)，迭代格式为：
\[ (D_B + L_B) \mathbf{u}^{(k+1)} = \mathbf{b} - U_B \mathbf{u}^{(k)} \]
由于 \((D_B + L_B)\) 是块下三角矩阵，可顺序求解每个块系统。
-
求解过程：
从 \(j=1\) 开始，解 \(T \mathbf{u}_1^{(k+1)} = \mathbf{b}_1 + \mathbf{u}_{2}^{(k)}\)（因为 \(\mathbf{u}_{0}^{(k+1)} = \mathbf{0}\)）。
然后对 \(j=2\)，解 \(T \mathbf{u}_2^{(k+1)} = \mathbf{b}_2 + \mathbf{u}_{1}^{(k+1)} + \mathbf{u}_{3}^{(k)}\)，其中 \(\mathbf{u}_{1}^{(k+1)}\) 已更新。
如此继续，直到 \(j=n\)。 -
优点：
- 通常比块Jacobi收敛更快，因为利用了最新信息。
- 每迭代步的计算成本与块Jacobi相同（仍为 \(O(N)\)），但所需迭代步数更少。
五、收敛性与实用技巧
-
收敛条件：
- 对于Poisson方程，系数矩阵 \(A\) 对称正定、不可约对角占优，保证块Jacobi和块Gauss-Seidel均收敛。
- 块Gauss-Seidel的渐进收敛速度约是块Jacobi的两倍（对于模型问题）。
-
加速技术——块逐次超松弛（块SOR）：
引入松弛参数 \(\omega\)，将块Gauss-Seidel更新与外推结合：
\[ T \tilde{\mathbf{u}}_j = \mathbf{b}_j + \mathbf{u}_{j-1}^{(k+1)} + \mathbf{u}_{j+1}^{(k)} \]
\[ \mathbf{u}_j^{(k+1)} = (1-\omega) \mathbf{u}_j^{(k)} + \omega \tilde{\mathbf{u}}_j \]
最优参数 \(\omega_{\text{opt}}\) 可基于矩阵谱半径估计，显著加快收敛。
-
并行性：
- 块Jacobi具有天然并行性：所有块可同时求解。
- 块Gauss-Seidel是顺序的，但可通过红黑排序或多色排序实现并行化（将块分组，组内并行，组间顺序）。
-
预条件子应用：
块迭代法本身可作为Krylov子空间方法（如CG、GMRES）的预条件子。例如，块Jacobi预条件子：\(M = D_B\)，求解 \(M^{-1} A \mathbf{u} = M^{-1} \mathbf{b}\)，其中 \(M^{-1}\) 的应用即解多个独立的三对角系统。
六、算法总结
块迭代法解离散Poisson方程的核心在于：
- 利用问题的自然块结构（按网格行分组），将大规模系统分解为多个较小、易于求解的三对角系统。
- 块Jacobi：简单、并行，但收敛较慢。
- 块Gauss-Seidel：顺序执行，收敛更快，可通过技巧并行化。
- 高效求解子问题：每个块系统使用 \(O(m)\) 的追赶法求解，总成本线性于网格点数。
- 适用性：不仅限于Poisson方程，也适用于其他分离变量产生的块三对角系统。
这种方法在科学与工程计算中广泛应用，特别是在需要中等精度解、且希望避免复杂并行算法的情况下，提供了简单有效的求解途径。