告别龟速迭代：用Python手把手实现一个简易多重网格求解器（附完整代码）

告别龟速迭代：用Python手把手实现一个简易多重网格求解器（附完整代码）

在科学计算领域，求解大型线性方程组是许多工程问题的核心挑战。传统迭代方法如Jacobi或Gauss-Seidel虽然实现简单，但当面对高分辨率网格时，收敛速度会急剧下降——这种现象被称为"迭代停滞"。想象一下在4K分辨率下处理流体模拟，每个时间步都需要数万次迭代，计算成本令人望而生畏。

多重网格方法（Multigrid Method）正是为解决这一痛点而生。它巧妙利用了不同尺度网格对误差成分的"过滤"特性：细网格擅长消除高频误差，粗网格则能高效处理低频分量。通过在不同网格层次间传递信息，多重网格可以实现接近O(N)的计算复杂度，比传统方法的O(N²)有质的飞跃。本文将用Python从零构建一个可运行的两层V-cycle求解器，带你穿透理论迷雾，直达算法本质。

1. 问题建模与基础工具

我们选择二维泊松方程作为模型问题：

-∇²u = f 在Ω内 u = 0 在边界∂Ω上

使用5点差分格式离散后，得到线性方程组Au=f。在256×256网格上，传统迭代方法可能需要10⁵次迭代，而多重网格有望在数十次循环内收敛。

1.1 初始化计算环境

import numpy as np import matplotlib.pyplot as plt from scipy.sparse import diags def poisson_matrix(N): """生成N×N网格的泊松方程系数矩阵""" diag = 4*np.ones(N*N) off_diag = -1*np.ones(N*N-1) off_diag[N-1::N] = 0 # 每N个元素跳过一个 return diags([off_diag, diag, off_diag], [-1,0,1], format='csr')

1.2 Jacobi迭代器实现

作为对比基准，我们先实现经典Jacobi迭代：

def jacobi_iter(A, b, x0, max_iter=100): x = x0.copy() D_inv = 1 / A.diagonal() residuals = [] for _ in range(max_iter): r = b - A @ x x += D_inv * r residuals.append(np.linalg.norm(r)) return x, residuals

关键观察：Jacobi在前几十次迭代中迅速降低残差，但随后进入缓慢收敛阶段——这正是低频误差主导的信号。

2. 多重网格核心组件

2.1 网格转换算子

限制算子（Restriction）将细网格残差传递到粗网格：

def restriction(fine): """全加权限制算子，细网格→粗网格""" coarse = np.zeros((fine.shape[0]//2, fine.shape[1]//2)) for i in range(1, fine.shape[0]-1): for j in range(1, fine.shape[1]-1): if i%2 == 1 and j%2 == 1: coarse[i//2,j//2] = 0.25*fine[i,j] + 0.125*( fine[i-1,j] + fine[i+1,j] + fine[i,j-1] + fine[i,j+1] ) + 0.0625*( fine[i-1,j-1] + fine[i-1,j+1] + fine[i+1,j-1] + fine[i+1,j+1] ) return coarse

插值算子（Prolongation）反向传递修正量：

def prolongation(coarse): """双线性插值，粗网格→细网格""" fine = np.zeros((coarse.shape[0]*2, coarse.shape[1]*2)) fine[1::2,1::2] = coarse[1:,1:] # 中心点 # 水平插值 fine[1::2,::2] = 0.5*(coarse[1:,:-1] + coarse[1:,1:]) # 垂直插值 fine[::2,1::2] = 0.5*(coarse[:-1,1:] + coarse[1:,1:]) # 对角插值 fine[::2,::2] = 0.25*(coarse[:-1,:-1] + coarse[:-1,1:] + coarse[1:,:-1] + coarse[1:,1:]) return fine

2.2 光滑迭代器优化

采用带权Jacobi方法作为光滑器（ω=2/3时效果最佳）：

def weighted_jacobi(A, b, x, omega=2/3, iterations=3): D_inv = omega / A.diagonal() for _ in range(iterations): x += D_inv * (b - A @ x) return x

3. V-Cycle算法实现

完整的V-Cycle流程如下：

1. 预光滑：在细网格上执行几次加权Jacobi迭代
2. 残差计算：r = b - Ax
3. 限制：将残差限制到粗网格 r_coarse = R(r)
4. 粗网格求解：A_coarse e_coarse = r_coarse
5. 插值：将修正量插回细网格 e_fine = P(e_coarse)
6. 修正解：x += e_fine
7. 后光滑：再次执行光滑迭代

def v_cycle(A, b, x, level, max_level=3): if level == max_level: # 最粗网格直接求解 return np.linalg.solve(A.toarray(), b) # 预光滑 x = weighted_jacobi(A, b, x) # 计算残差并限制 residual = b - A @ x coarse_residual = restriction(residual.reshape(int(np.sqrt(len(b))),-1)).flatten() # 粗网格修正 A_coarse = poisson_matrix(int(np.sqrt(len(coarse_residual)))) error_coarse = v_cycle(A_coarse, coarse_residual, np.zeros_like(coarse_residual), level+1) # 插值并修正 error_fine = prolongation(error_coarse.reshape(int(np.sqrt(len(error_coarse))),-1)).flatten() x += error_fine # 后光滑 x = weighted_jacobi(A, b, x) return x

4. 性能对比与可视化

我们构造一个已知解的问题进行测试：

N = 64 # 网格尺寸 u_exact = np.sin(np.pi*np.linspace(0,1,N+2)[1:-1])[:,None] * np.sin(np.pi*np.linspace(0,1,N+2)[1:-1]) f = 2*np.pi**2 * u_exact # -∇²u = f A = poisson_matrix(N) b = f.flatten() x0 = np.random.rand(N*N) # Jacobi迭代测试 _, jacobi_res = jacobi_iter(A, b, x0, 100) # 多重网格测试 mg_res = [] x = x0.copy() for cycle in range(10): x = v_cycle(A, b, x, 0) mg_res.append(np.linalg.norm(b - A @ x))

绘制收敛曲线：

plt.semilogy(jacobi_res, label='Jacobi') plt.semilogy(np.arange(10)*2, mg_res, 'o-', label='Multigrid (per cycle)') plt.xlabel('Iterations/cycles') plt.ylabel('Residual norm') plt.legend() plt.show()

典型结果会显示：

• Jacobi迭代在前50步快速收敛，之后几乎停滞
• 每个V-Cycle相当于约20次Jacobi迭代的效果
• 多重网格在6-8个循环后即可达到1e-6精度

5. 进阶优化方向

5.1 完全多重网格（FMG）

初始猜测采用从粗网格插值的解：

def full_multigrid(A, b, levels=3): if levels == 0: return np.linalg.solve(A.toarray(), b) # 在粗网格上求解 A_coarse = poisson_matrix(int(np.sqrt(len(b))//2)) b_coarse = restriction(b.reshape(int(np.sqrt(len(b))),-1)).flatten() x_coarse = full_multigrid(A_coarse, b_coarse, levels-1) # 插值并执行V-Cycle x_init = prolongation(x_coarse.reshape(int(np.sqrt(len(x_coarse))),-1)).flatten() return v_cycle(A, b, x_init, 0)

5.2 代数多重网格（AMG）

对于非结构化网格问题，可以基于矩阵模式自动构造粗网格：

from pyamg import ruge_stuben_solver def amg_solve(A, b): ml = ruge_stuben_solver(A) # 构建AMG层级 return ml.solve(b, tol=1e-10) # 使用AMG作为求解器

5.3 并行计算优化

粗网格计算是天然的并行化机会：

from mpi4py import MPI def parallel_restriction(fine): comm = MPI.COMM_WORLD rank = comm.Get_rank() # 按区域划分网格并执行局部限制 # ... (具体实现取决于并行策略)

在实际项目中，多重网格常与以下技术配合使用：

• 预处理技术：作为PCG方法的预条件子
• 自适应网格：根据解的特性动态调整网格密度
• 混合精度：粗网格计算使用低精度算术

将上述代码扩展为3D版本时，需要注意：

1. 7点差分格式替代5点格式
2. 3D限制算子需要27点加权
3. 内存消耗随N³增长，需要更积极的粗化策略

一个实用的调试技巧：在每层网格上输出残差范数，观察误差成分如何在不同层级间转移。当高频误差在细网格上快速衰减，低频误差在粗网格上显著降低时，说明你的多重网格实现正在正确工作。