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

告别龟速迭代用Python手把手实现一个简易多重网格求解器附完整代码在科学计算领域求解大型线性方程组是许多工程问题的核心挑战。传统迭代方法如Jacobi或Gauss-Seidel虽然实现简单但当面对高分辨率网格时收敛速度会急剧下降——这种现象被称为迭代停滞。想象一下在4K分辨率下处理流体模拟每个时间步都需要数万次迭代计算成本令人望而生畏。多重网格方法Multigrid Method正是为解决这一痛点而生。它巧妙利用了不同尺度网格对误差成分的过滤特性细网格擅长消除高频误差粗网格则能高效处理低频分量。通过在不同网格层次间传递信息多重网格可以实现接近O(N)的计算复杂度比传统方法的O(N²)有质的飞跃。本文将用Python从零构建一个可运行的两层V-cycle求解器带你穿透理论迷雾直达算法本质。1. 问题建模与基础工具我们选择二维泊松方程作为模型问题-∇²u f 在Ω内 u 0 在边界∂Ω上使用5点差分格式离散后得到线性方程组Auf。在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], formatcsr)1.2 Jacobi迭代器实现作为对比基准我们先实现经典Jacobi迭代def jacobi_iter(A, b, x0, max_iter100): 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[i1,j] fine[i,j-1] fine[i,j1] ) 0.0625*( fine[i-1,j-1] fine[i-1,j1] fine[i1,j-1] fine[i1,j1] ) 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 fine2.2 光滑迭代器优化采用带权Jacobi方法作为光滑器ω2/3时效果最佳def weighted_jacobi(A, b, x, omega2/3, iterations3): D_inv omega / A.diagonal() for _ in range(iterations): x D_inv * (b - A x) return x3. V-Cycle算法实现完整的V-Cycle流程如下预光滑在细网格上执行几次加权Jacobi迭代残差计算r b - Ax限制将残差限制到粗网格 r_coarse R(r)粗网格求解A_coarse e_coarse r_coarse插值将修正量插回细网格 e_fine P(e_coarse)修正解x e_fine后光滑再次执行光滑迭代def v_cycle(A, b, x, level, max_level3): 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), level1) # 插值并修正 error_fine prolongation(error_coarse.reshape(int(np.sqrt(len(error_coarse))),-1)).flatten() x error_fine # 后光滑 x weighted_jacobi(A, b, x) return x4. 性能对比与可视化我们构造一个已知解的问题进行测试N 64 # 网格尺寸 u_exact np.sin(np.pi*np.linspace(0,1,N2)[1:-1])[:,None] * np.sin(np.pi*np.linspace(0,1,N2)[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, labelJacobi) plt.semilogy(np.arange(10)*2, mg_res, o-, labelMultigrid (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, levels3): 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, tol1e-10) # 使用AMG作为求解器5.3 并行计算优化粗网格计算是天然的并行化机会from mpi4py import MPI def parallel_restriction(fine): comm MPI.COMM_WORLD rank comm.Get_rank() # 按区域划分网格并执行局部限制 # ... (具体实现取决于并行策略)在实际项目中多重网格常与以下技术配合使用预处理技术作为PCG方法的预条件子自适应网格根据解的特性动态调整网格密度混合精度粗网格计算使用低精度算术将上述代码扩展为3D版本时需要注意7点差分格式替代5点格式3D限制算子需要27点加权内存消耗随N³增长需要更积极的粗化策略一个实用的调试技巧在每层网格上输出残差范数观察误差成分如何在不同层级间转移。当高频误差在细网格上快速衰减低频误差在粗网格上显著降低时说明你的多重网格实现正在正确工作。