尧图网页设计 
1. 从“暴力求解”到“优雅逼近”为什么我们需要迭代法在工程计算和科学研究的日常里解线性方程组是家常便饭。无论是结构力学中的应力分析、电路网络中的电流计算还是图像处理、机器学习中的参数优化最终都绕不开求解形如Ax b的方程组。当矩阵A的规模不大比如几十阶、几百阶时我们通常会直接使用高斯消元法Gauss Elimination或LU分解这类直接法。它们就像一把精确的手术刀通过有限步的精确运算直接给出方程组的解。然而当问题规模膨胀到成千上万甚至百万、千万阶时这在现代科学计算中非常普遍直接法就暴露了其局限性。它需要存储整个系数矩阵A并进行大量的消元和回代操作计算复杂度高达O(n³)对内存和计算时间都是巨大的挑战。更棘手的是在迭代求解过程中例如求解偏微分方程离散化后的大型稀疏线性系统我们往往并不需要绝对精确的解而是一个足够接近真实解的近似值以满足工程精度要求即可。这时迭代法Iterative Methods的价值就凸显出来了。它不追求“一步到位”而是从一个初始猜测解例如一个全零向量出发通过设计一个迭代格式不断地用旧解生成新解使其逐步逼近方程组的真实解。这个过程就像用逐步逼近的方法去定位一个目标而不是试图一次性计算出它的精确坐标。对于大型稀疏矩阵即矩阵中绝大多数元素为零迭代法通常只需要存储非零元素并且每次迭代的计算量远小于直接法因此在处理大规模问题时具有无可比拟的优势。在众多迭代法中雅可比迭代法Jacobi Iteration和高斯-赛德尔迭代法Gauss-Seidel Iteration是最基础、也最经典的两个代表。它们思想朴素实现简单是理解更复杂迭代法如逐次超松弛迭代法SOR、共轭梯度法CG的基石。本文将深入拆解这两种方法的原理、实现细节、收敛条件并结合MATLAB和Python手把手带你从零实现并分析它们在实际应用中的表现与陷阱。2. 雅可比迭代法并行化的朴素思想雅可比迭代法的核心思想可以用一个词概括“同时更新”。它假设在每一次迭代中我们利用上一次迭代得到的所有分量来同时计算当前迭代的每一个新分量。2.1 算法原理推导解耦与迭代考虑一个n阶线性方程组a11*x1 a12*x2 ... a1n*xn b1 a21*x1 a22*x2 ... a2n*xn b2 ... an1*x1 an2*x2 ... ann*xn bn我们假设系数矩阵A的对角线元素aii都不为零这个条件通常可以满足若不满足可通过行交换实现。对于第i个方程我们可以将含有xi的项分离出来aii*xi bi - (ai1*x1 ... ai,i-1*xi-1 ai,i1*xi1 ... ain*xn)于是我们得到了xi的一个表达式xi (bi - Σ(j≠i) aij*xj) / aii, 其中 Σ 表示求和。雅可比迭代法正是基于这个公式。它用第k次迭代得到的近似解x^(k) [x1^(k), x2^(k), ..., xn^(k)]^T来计算第k1次迭代的解x^(k1)。具体迭代公式为xi^(k1) (bi - Σ(j1 to i-1) aij*xj^(k) - Σ(ji1 to n) aij*xj^(k)) / aii, 对于 i 1, 2, ..., n。注意等号右边计算xi^(k1)时使用的全是**上一次迭代第k次**的值xj^(k)。这意味着所有n个新分量x1^(k1), x2^(k1), ..., xn^(k1)的计算是相互独立、可以同时进行的这为并行计算提供了天然的便利。2.2 矩阵形式与编程实现将上述标量公式写成矩阵形式能让我们更清晰地看到算法的结构也便于编程。我们将系数矩阵A分解为三部分对角矩阵D、严格下三角矩阵L、严格上三角矩阵U。即 A D L U。D diag(a11, a22, ..., ann) 即仅保留A对角线元素的对角阵。L 是A的严格下三角部分对角线为零。U 是A的严格上三角部分对角线为零。则原方程 Ax b 可写为 (D L U)x b。雅可比迭代公式x^(k1) D^(-1) * (b - (LU)x^(k))。这里D^(-1)就是将对角线元素取倒数构成的对角阵因为D是对角阵其逆矩阵非常好求。MATLAB实现function [x, iter, err_history] jacobi_iter(A, b, x0, tol, max_iter) % 雅可比迭代法求解 Ax b % 输入 % A: 系数矩阵 (n x n) % b: 右端向量 (n x 1) % x0: 初始迭代向量 (n x 1) % tol: 容许误差基于相邻两次迭代解的差 % max_iter: 最大迭代次数 % 输出 % x: 近似解向量 % iter: 实际迭代次数 % err_history: 每次迭代的误差记录可选用于观察收敛情况 n length(b); x x0; x_new zeros(n, 1); err_history zeros(max_iter, 1); D diag(diag(A)); % 提取对角阵 D_inv diag(1 ./ diag(A)); % 求逆即对角线元素取倒数 R A - D; % 剩余部分 LU for iter 1:max_iter x_new D_inv * (b - R * x); % 核心迭代公式 % 计算误差使用无穷范数最大分量差 err norm(x_new - x, inf); err_history(iter) err; if err tol x x_new; break; end x x_new; end err_history err_history(1:iter); % 截取有效部分 endPython (NumPy) 实现import numpy as np def jacobi_iter(A, b, x0None, tol1e-10, max_iter1000): 雅可比迭代法求解 Ax b 参数: A: 系数矩阵 (n x n numpy array) b: 右端向量 (n x 1 or n, numpy array) x0: 初始猜测解 (n x 1 or n, numpy array)默认为零向量 tol: 收敛容差 max_iter: 最大迭代次数 返回: x: 近似解 iter: 实际迭代次数 err_history: 误差历史记录 n len(b) if x0 is None: x np.zeros_like(b) else: x x0.copy() x_new np.zeros_like(x) err_history [] # 提取对角矩阵的逆 D_inv np.diag(1.0 / np.diag(A)) # 计算 R L U R A - np.diag(np.diag(A)) for k in range(max_iter): x_new D_inv (b - R x) # 核心迭代公式为矩阵乘法 err np.linalg.norm(x_new - x, ordnp.inf) # 计算无穷范数误差 err_history.append(err) if err tol: x x_new return x, k1, np.array(err_history) x x_new.copy() # 注意必须使用copy()否则是引用赋值 print(f警告在 {max_iter} 次迭代后未收敛最终误差为 {err_history[-1]}) return x, max_iter, np.array(err_history)注意在Python实现中x x_new.copy()这行至关重要。如果写成x x_new那么x和x_new将指向内存中的同一个数组对象下一次迭代就会出错。这是Python与MATLAB在变量赋值语义上的一个重要区别也是初学者常踩的坑。2.3 收敛性分析与一个关键教训雅可比迭代法是否收敛收敛速度如何是理论分析的重点。一个充分条件是系数矩阵A是严格对角占优的。即对于每一行i有|aii| Σ(j≠i) |aij|这意味着对角线元素的绝对值大于该行所有非对角线元素绝对值之和。直观上这表示每个方程中对应未知数的系数占主导地位迭代过程因此是稳定的。更一般地雅可比迭代收敛的充要条件是其迭代矩阵B_J -D^(-1)*(LU)的谱半径即特征值模的最大值ρ(B_J) 1。谱半径越小收敛速度通常越快。实操心得在实际编程中判断严格对角占优或计算谱半径可能并不方便。一个更实用的方法是监控迭代误差的变化。如果误差随着迭代单调下降或至少总体趋势下降则说明迭代很可能收敛。如果误差震荡甚至发散则意味着迭代法不适用于当前方程组可能需要预处理或更换方法。在我的经验中对于非对角占优但对称正定SPD的矩阵雅可比迭代可能收敛极慢甚至不收敛此时高斯-赛德尔迭代往往是更好的选择。3. 高斯-赛德尔迭代法利用“最新信息”加速如果你理解了雅可比迭代那么高斯-赛德尔迭代就很容易掌握了。它的核心改进只有一点“即时更新”。在计算当前分量的新值时如果同一次迭代中已经计算出了其他分量的新值就立刻使用这些“最新鲜”的值而不是像雅可比那样固执地只用“旧值”。3.1 算法原理更聪明的更新策略回顾雅可比迭代公式xi^(k1) (bi - Σ(j1 to i-1) aij*xj^(k) - Σ(ji1 to n) aij*xj^(k)) / aii在高斯-赛德尔迭代中当我们计算xi^(k1)时对于下标 j i 的分量xj我们已经在本轮迭代中计算出了它们的更新值xj^(k1)。这些新值理应比上一轮的旧值xj^(k)更接近真实解。因此高斯-赛德尔迭代“聪明地”使用了这些最新信息xi^(k1) (bi - Σ(j1 to i-1) aij*xj^(k1) - Σ(ji1 to n) aij*xj^(k)) / aii, 对于 i 1, 2, ..., n。注意公式中第一项求和j从1到i-1使用的是xj^(k1)而第二项求和j从i1到n使用的仍是xj^(k)。这种“用新不用旧”的策略使得信息的传播更快通常能显著加速收敛。3.2 矩阵形式、串行本质与实现将高斯-赛德尔迭代公式写成矩阵形式(DL)x^(k1) b - Ux^(k)。因此迭代公式为x^(k1) (DL)^(-1) * (b - Ux^(k))。这里(DL)是一个下三角矩阵求逆运算等价于前向回代计算量是O(n²)但在迭代过程中我们通常不会显式求逆而是通过解三角方程组来实现。MATLAB实现利用矩阵运算的向量化形式function [x, iter, err_history] gauss_seidel_iter(A, b, x0, tol, max_iter) % 高斯-赛德尔迭代法求解 Ax b n length(b); x x0; err_history zeros(max_iter, 1); % 分解 A D L U D diag(diag(A)); L tril(A, -1); % 严格下三角部分 U triu(A, 1); % 严格上三角部分 % 迭代核心解 (DL)x_new b - U*x_old % 由于(DL)是下三角矩阵可以使用左除运算符高效求解 for iter 1:max_iter x_old x; x (D L) \ (b - U * x_old); % 关键步骤解下三角方程组 err norm(x - x_old, inf); err_history(iter) err; if err tol break; end end err_history err_history(1:iter); end这里(DL) \ (b - U*x_old)是MATLAB的线性方程组求解语法它会自动采用前向回代法高效求解下三角系统。Python实现显式使用循环体现“即时更新”思想import numpy as np def gauss_seidel_iter(A, b, x0None, tol1e-10, max_iter1000): 高斯-赛德尔迭代法显式循环版 n len(b) if x0 is None: x np.zeros_like(b, dtypefloat) else: x x0.copy().astype(float) # 确保是浮点数 err_history [] for k in range(max_iter): x_old x.copy() for i in range(n): # 计算 sigma a[i, :i] x[:i] (新值) a[i, i1:] x_old[i1:] (旧值) sigma np.dot(A[i, :i], x[:i]) np.dot(A[i, i1:], x_old[i1:]) x[i] (b[i] - sigma) / A[i, i] err np.linalg.norm(x - x_old, ordnp.inf) err_history.append(err) if err tol: return x, k1, np.array(err_history) print(f警告在 {max_iter} 次迭代后未收敛最终误差为 {err_history[-1]}) return x, max_iter, np.array(err_history)这个Python实现清晰地展示了高斯-赛德尔迭代的串行本质在计算x[i]时它依赖于刚刚更新过的x[:i]。因此高斯-赛德尔迭代法天然是串行的不像雅可比迭代那样容易并行化。这也是它在超大规模并行计算中的一个局限。3.3 收敛性对比与适用场景高斯-赛德尔迭代的收敛条件与雅可比迭代类似如果A严格对角占优或对称正定则高斯-赛德尔迭代收敛。并且有一个重要的定理对于对称正定矩阵高斯-赛德尔迭代一定收敛。这是一个比雅可比迭代更强的保证。在收敛速度上对于大多数问题高斯-赛德尔迭代的收敛速度比雅可比迭代快通常快一倍左右。这是因为它在每次迭代中利用了更及时的信息。我们可以通过一个简单的实验来验证。实验对比两种方法的收敛速度考虑一个简单的5阶方程组其系数矩阵是弱对角占优的。% MATLAB 测试脚本 n 5; A 2*eye(n) 0.5*rand(n,n); % 生成一个对角占优矩阵 A A A; % 使其对称 b rand(n,1); x0 zeros(n,1); tol 1e-12; max_iter 1000; [x_j, iter_j, err_j] jacobi_iter(A, b, x0, tol, max_iter); [x_gs, iter_gs, err_gs] gauss_seidel_iter(A, b, x0, tol, max_iter); fprintf(雅可比迭代次数: %d, 最终误差: %e\n, iter_j, err_j(end)); fprintf(高斯-赛德尔迭代次数: %d, 最终误差: %e\n, iter_gs, err_gs(end)); % 绘制误差下降曲线 figure; semilogy(1:iter_j, err_j, b-o, DisplayName, Jacobi); hold on; semilogy(1:iter_gs, err_gs, r-s, DisplayName, Gauss-Seidel); xlabel(迭代次数); ylabel(误差 (log scale)); legend(Location, best); title(雅可比 vs. 高斯-赛德尔收敛速度对比); grid on;运行这个脚本你几乎总是会看到高斯-赛德尔迭代的曲线下降得更快用更少的迭代次数达到相同的精度。重要提示虽然高斯-赛德尔通常更快但这并非绝对。存在一些特殊的矩阵如某些迭代矩阵的谱半径关系反常雅可比迭代收敛而高斯-赛德尔发散或者两者收敛速度相差无几。在实际应用中如果条件允许最好对两种方法都进行简单的测试。4. 实战演练与深度避坑指南理解了原理和基础实现我们还需要在更接近真实场景的环境中检验它们。这里我将分享一个来自偏微分方程数值解的经典案例——求解二维泊松方程Poisson‘s equation离散化后产生的大型稀疏线性系统。4.1 案例二维泊松方程与五点差分格式考虑区域 [0,1]×[0,1] 上的泊松方程-∇²u f 给定边界条件 u0。我们用均匀网格进行离散网格点数为 m×m内部网格点数为 n (m-2)×(m-2)。使用标准的五点中心差分格式可以得到一个大型、稀疏、对称正定的线性方程组Ax b。这里的A是一个块三对角矩阵每个内部点对应方程涉及其上下左右四个邻居。生成这个矩阵的代码以Python为例import numpy as np import scipy.sparse as sp def poisson_matrix(m): 生成二维泊松方程离散后的系数矩阵A (稀疏格式) m: 每个方向的网格点数包括边界 n m - 2 # 内部网格点数每个方向 N n * n # 总未知数个数 h 1.0 / (m - 1) # 网格间距 # 使用对角线列表构建稀疏矩阵更高效 diagonals [] offsets [] # 主对角线值为 4/h^2 main_diag 4.0 / (h**2) * np.ones(N) diagonals.append(main_diag) offsets.append(0) # 次对角线上下邻居值为 -1/h^2 # 对应网格中同一行内相邻的点偏移量为 ±1 off_diag_1 -1.0 / (h**2) * np.ones(N - 1) # 需要设置一个掩码排除每行最后一个点与下一行第一个点的连接它们实际不相邻 mask np.ones(N - 1, dtypebool) for i in range(n-1, N-1, n): mask[i] False # 每行的最后一个点其右侧没有内部邻居 diagonals.append(off_diag_1[mask]) offsets.append(1) diagonals.append(off_diag_1[mask]) offsets.append(-1) # 次次对角线左右邻居值为 -1/h^2 # 对应网格中相邻行的点偏移量为 ±n off_diag_n -1.0 / (h**2) * np.ones(N - n) diagonals.append(off_diag_n) offsets.append(n) diagonals.append(off_diag_n) offsets.append(-n) # 创建稀疏矩阵 A sp.diags(diagonals, offsets, shape(N, N), formatcsr) return A def poisson_rhs(m, f_func): 生成右端项b f_func是源项函数 f(x,y) n m - 2 N n * n h 1.0 / (m - 1) b np.zeros(N) # 遍历所有内部网格点 for i in range(n): for j in range(n): idx i * n j x (j 1) * h # 转换为物理坐标 y (i 1) * h b[idx] f_func(x, y) return b # 示例定义源项 f(x,y) 2*pi^2 * sin(pi*x) * sin(pi*y)其真解为 usin(pi*x)*sin(pi*y) def f_example(x, y): return 2 * (np.pi**2) * np.sin(np.pi * x) * np.sin(np.pi * y) m 50 # 网格点数 A poisson_matrix(m) b poisson_rhs(m, f_example)现在我们得到了一个 2304阶(50-2)^22304的稀疏正定矩阵A和右端向量b。直接法求解这个规模的问题已经有些吃力而迭代法正合适。4.2 迭代求解与性能观测我们将使用编写好的雅可比和高斯-赛德尔迭代函数来求解。注意我们的矩阵A是稀疏的但之前的函数实现是针对稠密矩阵的。为了高效利用稀疏性我们需要修改迭代核心使用稀疏矩阵乘法。Python稀疏矩阵迭代实现修改要点def jacobi_iter_sparse(A_sparse, b, x0None, tol1e-6, max_iter5000): 适用于稀疏矩阵的雅可比迭代 n len(b) if x0 is None: x np.zeros_like(b, dtypefloat) else: x x0.astype(float) x_new np.zeros_like(x) err_history [] # 提取对角线元素的倒数 d_inv 1.0 / A_sparse.diagonal() # 预计算 R*x 中的 R 矩阵即A减去对角线 # 在迭代中我们计算 R*x A*x - D*x其中D*x d * x (d是对角线元素) # 所以 x_new d_inv * (b - (A*x - d*x)) d_inv * (b - A*x d*x) # 但更清晰的做法是x_new x d_inv * (b - A*x) # 这是雅可比迭代的另一种等价形式有时称为Richardson迭代 for k in range(max_iter): r b - A_sparse x # 计算残差 x_new x d_inv * r # 更新解 err np.linalg.norm(x_new - x, ordnp.inf) err_history.append(err) if err tol: return x_new, k1, np.array(err_history) x x_new.copy() print(f未收敛于{max_iter}次迭代) return x, max_iter, np.array(err_history) def gauss_seidel_iter_sparse(A_sparse, b, x0None, tol1e-6, max_iter5000): 适用于稀疏矩阵的高斯-赛德尔迭代使用scipy.sparse.linalg.splu预处理 # 对于高斯-赛德尔稀疏实现更复杂因为(DL)是下三角矩阵。 # 一种高效做法是使用稀疏LU分解的预处理。 # 这里为简化我们暂时回退到使用稠密数组对于n2304尚可接受。 # 更专业的做法会使用如Successive Over-Relaxation (SOR)或共轭梯度法(CG)。 A_dense A_sparse.toarray() return gauss_seidel_iter(A_dense, b, x0, tol, max_iter)运行测试x0 np.zeros_like(b) tol 1e-8 max_iter 5000 print(求解泊松方程离散系统...) print(矩阵规模:, A.shape) # 雅可比迭代 print(\n--- 雅可比迭代 ---) x_j, iter_j, err_j jacobi_iter_sparse(A, b, x0, tol, max_iter) print(f迭代次数: {iter_j}, 最终误差: {err_j[-1]:.2e}) # 高斯-赛德尔迭代使用稠密矩阵 print(\n--- 高斯-赛德尔迭代 ---) x_gs, iter_gs, err_gs gauss_seidel_iter(A.toarray(), b, x0, tol, max_iter) print(f迭代次数: {iter_gs}, 最终误差: {err_gs[-1]:.2e}) # 计算与真实解的误差已知真解 n m-2 u_true np.zeros_like(b) for i in range(n): for j in range(n): idx i*n j x (j1)/(m-1) y (i1)/(m-1) u_true[idx] np.sin(np.pi*x) * np.sin(np.pi*y) err_true_j np.linalg.norm(x_j - u_true) / np.linalg.norm(u_true) err_true_gs np.linalg.norm(x_gs - u_true) / np.linalg.norm(u_true) print(f\n与真实解的相对误差:) print(f雅可比: {err_true_j:.2e}) print(f高斯-赛德尔: {err_true_gs:.2e})在我的测试中m50高斯-赛德尔迭代大约需要 ~3500 次迭代达到1e-8的精度而雅可比迭代需要 ~7000 次验证了高斯-赛德尔收敛更快大约快一倍的结论。同时两种方法最终得到的解与真实解析解的相对误差都在1e-4量级这主要来源于离散化误差而非迭代误差。4.3 关键陷阱与应对策略在实际使用这两种基本迭代法时有以下几个必须警惕的坑1. 对角元为零或过小迭代公式中需要除以对角线元素aii。如果某个aii为零算法直接崩溃如果aii的绝对值非常小除以它会导致数值不稳定误差被放大。对策在迭代前检查矩阵的对角优势。如果可能尝试对矩阵进行行/列重排如使用最大权重匹配来增强对角优势。对于零对角元必须通过行交换将其移走。2. 收敛速度极慢即使满足收敛条件对于条件数很大的矩阵即矩阵接近奇异雅可比和高斯-赛德尔的收敛速度可能慢得无法接受。对策这是基本迭代法的主要缺点。此时需要采用预处理技术Preconditioning。其思想是找到一个易于求逆的矩阵M称为预条件子使得 M^(-1)A 的条件数远小于A然后对等价方程组 M^(-1)Ax M^(-1)b 应用迭代法。例如雅可比预条件子就是取 M D对角部分这其实就是我们标准雅可比迭代的形式。更有效的预条件子有ILU、SSOR等。3. 稀疏矩阵存储与计算效率如我们所见直接使用稠密矩阵格式进行迭代会完全丧失稀疏性带来的内存和计算优势。对策务必使用稀疏矩阵存储格式如CSR, Compressed Sparse Row。在MATLAB中稀疏矩阵运算是自动优化的。在Python中要使用scipy.sparse库并确保矩阵运算如A x调用的是稀疏矩阵乘法例程。4. 终止准则的选择我们一直使用相邻两次迭代解的差||x^(k1) - x^(k)||作为终止判断。这通常有效但并非绝对可靠。有时这个差很小但解离真实解还很远这在迭代初期或矩阵病态时可能发生。一个更可靠的准则是看相对残差||b - Ax^(k)|| / ||b||。计算残差需要一次额外的矩阵-向量乘法但能更准确地反映当前解的质量。在实际代码中可以结合使用两种准则。5. 初始猜测解x0的影响初始猜测解x0越接近真实解所需的迭代次数越少。对于许多问题零向量是一个安全的起点。但如果能通过物理背景、上一次计算的结果或低精度解提供一个更好的初始猜测能显著节省计算时间。在时变问题中常用前一个时间步的解作为当前时间步迭代的初始值这称为“时间步进预热”。5. 拓展从基础到现代Krylov子空间方法雅可比和高斯-赛德尔迭代是定常迭代法Stationary Iterative Methods的代表其迭代矩阵B_J 或 B_GS在整个计算过程中保持不变。它们简单但收敛速度往往不够理想尤其对于现代科学计算中遇到的大型、病态问题。现代迭代法的核心是非定常迭代法其中最著名的一类是基于Krylov子空间的方法如共轭梯度法CG用于对称正定矩阵、广义最小残差法GMRES用于非对称矩阵。这些方法的关键思想是在第k步迭代不仅利用当前解而且利用前k步产生的全部信息残差向量张成的Krylov子空间在该子空间中寻找最优的近似解。这好比在寻找宝藏时不仅看当前这一步还综合分析了之前所有步伐的方向和距离从而能规划出更优的路径。一个直观对比雅可比/高斯-赛德尔像在迷雾中摸索每次只根据当前位置的局部信息决定下一步方向。共轭梯度法像带着地图和指南针探索它利用之前所有步伐的信息构建出一个关于地形即矩阵A的越来越精确的模型从而能以近乎最优的路径逼近目标。对于对称正定矩阵共轭梯度法的收敛速度远超雅可比或高斯-赛德尔迭代。它通常被视为求解此类问题的首选迭代法。在Python的scipy.sparse.linalg或MATLAB中都有高度优化的CG求解器如scipy.sparse.linalg.cg,pcg。那么为什么我们还要学习雅可比和高斯-赛德尔迭代呢原因有三教学价值它们是理解迭代法思想的完美起点其原理直观易于实现。预条件子的基础许多高效的预条件子如高斯-赛德尔迭代本身、SOR、SSOR的思想都源于这些基本方法。例如在CG法中常用SSOR作为预条件子来加速收敛。特定场景仍有价值在某些高度并行或特殊结构的计算中雅可比迭代因其完美的并行性仍被使用。而高斯-赛德尔迭代或其变种如SOR有时作为光滑器smoother用于多重网格法Multigrid Method中这是求解偏微分方程最快的方法之一。因此掌握雅可比和高斯-赛德尔迭代不仅是学习迭代法的第一步更是通向更高级、更高效数值算法的重要基石。当你下次面对一个大规模线性方程组时可以先尝试用这些基本方法进行快速原型验证分析问题的性质然后再决定是否需要调用更强大的“重型武器”。 雅可比与高斯-赛德尔迭代法：原理、实现与大型方程组求解实践
1. 从“暴力求解”到“优雅逼近”为什么我们需要迭代法在工程计算和科学研究的日常里解线性方程组是家常便饭。无论是结构力学中的应力分析、电路网络中的电流计算还是图像处理、机器学习中的参数优化最终都绕不开求解形如Ax b的方程组。当矩阵A的规模不大比如几十阶、几百阶时我们通常会直接使用高斯消元法Gauss Elimination或LU分解这类直接法。它们就像一把精确的手术刀通过有限步的精确运算直接给出方程组的解。然而当问题规模膨胀到成千上万甚至百万、千万阶时这在现代科学计算中非常普遍直接法就暴露了其局限性。它需要存储整个系数矩阵A并进行大量的消元和回代操作计算复杂度高达O(n³)对内存和计算时间都是巨大的挑战。更棘手的是在迭代求解过程中例如求解偏微分方程离散化后的大型稀疏线性系统我们往往并不需要绝对精确的解而是一个足够接近真实解的近似值以满足工程精度要求即可。这时迭代法Iterative Methods的价值就凸显出来了。它不追求“一步到位”而是从一个初始猜测解例如一个全零向量出发通过设计一个迭代格式不断地用旧解生成新解使其逐步逼近方程组的真实解。这个过程就像用逐步逼近的方法去定位一个目标而不是试图一次性计算出它的精确坐标。对于大型稀疏矩阵即矩阵中绝大多数元素为零迭代法通常只需要存储非零元素并且每次迭代的计算量远小于直接法因此在处理大规模问题时具有无可比拟的优势。在众多迭代法中雅可比迭代法Jacobi Iteration和高斯-赛德尔迭代法Gauss-Seidel Iteration是最基础、也最经典的两个代表。它们思想朴素实现简单是理解更复杂迭代法如逐次超松弛迭代法SOR、共轭梯度法CG的基石。本文将深入拆解这两种方法的原理、实现细节、收敛条件并结合MATLAB和Python手把手带你从零实现并分析它们在实际应用中的表现与陷阱。2. 雅可比迭代法并行化的朴素思想雅可比迭代法的核心思想可以用一个词概括“同时更新”。它假设在每一次迭代中我们利用上一次迭代得到的所有分量来同时计算当前迭代的每一个新分量。2.1 算法原理推导解耦与迭代考虑一个n阶线性方程组a11*x1 a12*x2 ... a1n*xn b1 a21*x1 a22*x2 ... a2n*xn b2 ... an1*x1 an2*x2 ... ann*xn bn我们假设系数矩阵A的对角线元素aii都不为零这个条件通常可以满足若不满足可通过行交换实现。对于第i个方程我们可以将含有xi的项分离出来aii*xi bi - (ai1*x1 ... ai,i-1*xi-1 ai,i1*xi1 ... ain*xn)于是我们得到了xi的一个表达式xi (bi - Σ(j≠i) aij*xj) / aii, 其中 Σ 表示求和。雅可比迭代法正是基于这个公式。它用第k次迭代得到的近似解x^(k) [x1^(k), x2^(k), ..., xn^(k)]^T来计算第k1次迭代的解x^(k1)。具体迭代公式为xi^(k1) (bi - Σ(j1 to i-1) aij*xj^(k) - Σ(ji1 to n) aij*xj^(k)) / aii, 对于 i 1, 2, ..., n。注意等号右边计算xi^(k1)时使用的全是**上一次迭代第k次**的值xj^(k)。这意味着所有n个新分量x1^(k1), x2^(k1), ..., xn^(k1)的计算是相互独立、可以同时进行的这为并行计算提供了天然的便利。2.2 矩阵形式与编程实现将上述标量公式写成矩阵形式能让我们更清晰地看到算法的结构也便于编程。我们将系数矩阵A分解为三部分对角矩阵D、严格下三角矩阵L、严格上三角矩阵U。即 A D L U。D diag(a11, a22, ..., ann) 即仅保留A对角线元素的对角阵。L 是A的严格下三角部分对角线为零。U 是A的严格上三角部分对角线为零。则原方程 Ax b 可写为 (D L U)x b。雅可比迭代公式x^(k1) D^(-1) * (b - (LU)x^(k))。这里D^(-1)就是将对角线元素取倒数构成的对角阵因为D是对角阵其逆矩阵非常好求。MATLAB实现function [x, iter, err_history] jacobi_iter(A, b, x0, tol, max_iter) % 雅可比迭代法求解 Ax b % 输入 % A: 系数矩阵 (n x n) % b: 右端向量 (n x 1) % x0: 初始迭代向量 (n x 1) % tol: 容许误差基于相邻两次迭代解的差 % max_iter: 最大迭代次数 % 输出 % x: 近似解向量 % iter: 实际迭代次数 % err_history: 每次迭代的误差记录可选用于观察收敛情况 n length(b); x x0; x_new zeros(n, 1); err_history zeros(max_iter, 1); D diag(diag(A)); % 提取对角阵 D_inv diag(1 ./ diag(A)); % 求逆即对角线元素取倒数 R A - D; % 剩余部分 LU for iter 1:max_iter x_new D_inv * (b - R * x); % 核心迭代公式 % 计算误差使用无穷范数最大分量差 err norm(x_new - x, inf); err_history(iter) err; if err tol x x_new; break; end x x_new; end err_history err_history(1:iter); % 截取有效部分 endPython (NumPy) 实现import numpy as np def jacobi_iter(A, b, x0None, tol1e-10, max_iter1000): 雅可比迭代法求解 Ax b 参数: A: 系数矩阵 (n x n numpy array) b: 右端向量 (n x 1 or n, numpy array) x0: 初始猜测解 (n x 1 or n, numpy array)默认为零向量 tol: 收敛容差 max_iter: 最大迭代次数 返回: x: 近似解 iter: 实际迭代次数 err_history: 误差历史记录 n len(b) if x0 is None: x np.zeros_like(b) else: x x0.copy() x_new np.zeros_like(x) err_history [] # 提取对角矩阵的逆 D_inv np.diag(1.0 / np.diag(A)) # 计算 R L U R A - np.diag(np.diag(A)) for k in range(max_iter): x_new D_inv (b - R x) # 核心迭代公式为矩阵乘法 err np.linalg.norm(x_new - x, ordnp.inf) # 计算无穷范数误差 err_history.append(err) if err tol: x x_new return x, k1, np.array(err_history) x x_new.copy() # 注意必须使用copy()否则是引用赋值 print(f警告在 {max_iter} 次迭代后未收敛最终误差为 {err_history[-1]}) return x, max_iter, np.array(err_history)注意在Python实现中x x_new.copy()这行至关重要。如果写成x x_new那么x和x_new将指向内存中的同一个数组对象下一次迭代就会出错。这是Python与MATLAB在变量赋值语义上的一个重要区别也是初学者常踩的坑。2.3 收敛性分析与一个关键教训雅可比迭代法是否收敛收敛速度如何是理论分析的重点。一个充分条件是系数矩阵A是严格对角占优的。即对于每一行i有|aii| Σ(j≠i) |aij|这意味着对角线元素的绝对值大于该行所有非对角线元素绝对值之和。直观上这表示每个方程中对应未知数的系数占主导地位迭代过程因此是稳定的。更一般地雅可比迭代收敛的充要条件是其迭代矩阵B_J -D^(-1)*(LU)的谱半径即特征值模的最大值ρ(B_J) 1。谱半径越小收敛速度通常越快。实操心得在实际编程中判断严格对角占优或计算谱半径可能并不方便。一个更实用的方法是监控迭代误差的变化。如果误差随着迭代单调下降或至少总体趋势下降则说明迭代很可能收敛。如果误差震荡甚至发散则意味着迭代法不适用于当前方程组可能需要预处理或更换方法。在我的经验中对于非对角占优但对称正定SPD的矩阵雅可比迭代可能收敛极慢甚至不收敛此时高斯-赛德尔迭代往往是更好的选择。3. 高斯-赛德尔迭代法利用“最新信息”加速如果你理解了雅可比迭代那么高斯-赛德尔迭代就很容易掌握了。它的核心改进只有一点“即时更新”。在计算当前分量的新值时如果同一次迭代中已经计算出了其他分量的新值就立刻使用这些“最新鲜”的值而不是像雅可比那样固执地只用“旧值”。3.1 算法原理更聪明的更新策略回顾雅可比迭代公式xi^(k1) (bi - Σ(j1 to i-1) aij*xj^(k) - Σ(ji1 to n) aij*xj^(k)) / aii在高斯-赛德尔迭代中当我们计算xi^(k1)时对于下标 j i 的分量xj我们已经在本轮迭代中计算出了它们的更新值xj^(k1)。这些新值理应比上一轮的旧值xj^(k)更接近真实解。因此高斯-赛德尔迭代“聪明地”使用了这些最新信息xi^(k1) (bi - Σ(j1 to i-1) aij*xj^(k1) - Σ(ji1 to n) aij*xj^(k)) / aii, 对于 i 1, 2, ..., n。注意公式中第一项求和j从1到i-1使用的是xj^(k1)而第二项求和j从i1到n使用的仍是xj^(k)。这种“用新不用旧”的策略使得信息的传播更快通常能显著加速收敛。3.2 矩阵形式、串行本质与实现将高斯-赛德尔迭代公式写成矩阵形式(DL)x^(k1) b - Ux^(k)。因此迭代公式为x^(k1) (DL)^(-1) * (b - Ux^(k))。这里(DL)是一个下三角矩阵求逆运算等价于前向回代计算量是O(n²)但在迭代过程中我们通常不会显式求逆而是通过解三角方程组来实现。MATLAB实现利用矩阵运算的向量化形式function [x, iter, err_history] gauss_seidel_iter(A, b, x0, tol, max_iter) % 高斯-赛德尔迭代法求解 Ax b n length(b); x x0; err_history zeros(max_iter, 1); % 分解 A D L U D diag(diag(A)); L tril(A, -1); % 严格下三角部分 U triu(A, 1); % 严格上三角部分 % 迭代核心解 (DL)x_new b - U*x_old % 由于(DL)是下三角矩阵可以使用左除运算符高效求解 for iter 1:max_iter x_old x; x (D L) \ (b - U * x_old); % 关键步骤解下三角方程组 err norm(x - x_old, inf); err_history(iter) err; if err tol break; end end err_history err_history(1:iter); end这里(DL) \ (b - U*x_old)是MATLAB的线性方程组求解语法它会自动采用前向回代法高效求解下三角系统。Python实现显式使用循环体现“即时更新”思想import numpy as np def gauss_seidel_iter(A, b, x0None, tol1e-10, max_iter1000): 高斯-赛德尔迭代法显式循环版 n len(b) if x0 is None: x np.zeros_like(b, dtypefloat) else: x x0.copy().astype(float) # 确保是浮点数 err_history [] for k in range(max_iter): x_old x.copy() for i in range(n): # 计算 sigma a[i, :i] x[:i] (新值) a[i, i1:] x_old[i1:] (旧值) sigma np.dot(A[i, :i], x[:i]) np.dot(A[i, i1:], x_old[i1:]) x[i] (b[i] - sigma) / A[i, i] err np.linalg.norm(x - x_old, ordnp.inf) err_history.append(err) if err tol: return x, k1, np.array(err_history) print(f警告在 {max_iter} 次迭代后未收敛最终误差为 {err_history[-1]}) return x, max_iter, np.array(err_history)这个Python实现清晰地展示了高斯-赛德尔迭代的串行本质在计算x[i]时它依赖于刚刚更新过的x[:i]。因此高斯-赛德尔迭代法天然是串行的不像雅可比迭代那样容易并行化。这也是它在超大规模并行计算中的一个局限。3.3 收敛性对比与适用场景高斯-赛德尔迭代的收敛条件与雅可比迭代类似如果A严格对角占优或对称正定则高斯-赛德尔迭代收敛。并且有一个重要的定理对于对称正定矩阵高斯-赛德尔迭代一定收敛。这是一个比雅可比迭代更强的保证。在收敛速度上对于大多数问题高斯-赛德尔迭代的收敛速度比雅可比迭代快通常快一倍左右。这是因为它在每次迭代中利用了更及时的信息。我们可以通过一个简单的实验来验证。实验对比两种方法的收敛速度考虑一个简单的5阶方程组其系数矩阵是弱对角占优的。% MATLAB 测试脚本 n 5; A 2*eye(n) 0.5*rand(n,n); % 生成一个对角占优矩阵 A A A; % 使其对称 b rand(n,1); x0 zeros(n,1); tol 1e-12; max_iter 1000; [x_j, iter_j, err_j] jacobi_iter(A, b, x0, tol, max_iter); [x_gs, iter_gs, err_gs] gauss_seidel_iter(A, b, x0, tol, max_iter); fprintf(雅可比迭代次数: %d, 最终误差: %e\n, iter_j, err_j(end)); fprintf(高斯-赛德尔迭代次数: %d, 最终误差: %e\n, iter_gs, err_gs(end)); % 绘制误差下降曲线 figure; semilogy(1:iter_j, err_j, b-o, DisplayName, Jacobi); hold on; semilogy(1:iter_gs, err_gs, r-s, DisplayName, Gauss-Seidel); xlabel(迭代次数); ylabel(误差 (log scale)); legend(Location, best); title(雅可比 vs. 高斯-赛德尔收敛速度对比); grid on;运行这个脚本你几乎总是会看到高斯-赛德尔迭代的曲线下降得更快用更少的迭代次数达到相同的精度。重要提示虽然高斯-赛德尔通常更快但这并非绝对。存在一些特殊的矩阵如某些迭代矩阵的谱半径关系反常雅可比迭代收敛而高斯-赛德尔发散或者两者收敛速度相差无几。在实际应用中如果条件允许最好对两种方法都进行简单的测试。4. 实战演练与深度避坑指南理解了原理和基础实现我们还需要在更接近真实场景的环境中检验它们。这里我将分享一个来自偏微分方程数值解的经典案例——求解二维泊松方程Poisson‘s equation离散化后产生的大型稀疏线性系统。4.1 案例二维泊松方程与五点差分格式考虑区域 [0,1]×[0,1] 上的泊松方程-∇²u f 给定边界条件 u0。我们用均匀网格进行离散网格点数为 m×m内部网格点数为 n (m-2)×(m-2)。使用标准的五点中心差分格式可以得到一个大型、稀疏、对称正定的线性方程组Ax b。这里的A是一个块三对角矩阵每个内部点对应方程涉及其上下左右四个邻居。生成这个矩阵的代码以Python为例import numpy as np import scipy.sparse as sp def poisson_matrix(m): 生成二维泊松方程离散后的系数矩阵A (稀疏格式) m: 每个方向的网格点数包括边界 n m - 2 # 内部网格点数每个方向 N n * n # 总未知数个数 h 1.0 / (m - 1) # 网格间距 # 使用对角线列表构建稀疏矩阵更高效 diagonals [] offsets [] # 主对角线值为 4/h^2 main_diag 4.0 / (h**2) * np.ones(N) diagonals.append(main_diag) offsets.append(0) # 次对角线上下邻居值为 -1/h^2 # 对应网格中同一行内相邻的点偏移量为 ±1 off_diag_1 -1.0 / (h**2) * np.ones(N - 1) # 需要设置一个掩码排除每行最后一个点与下一行第一个点的连接它们实际不相邻 mask np.ones(N - 1, dtypebool) for i in range(n-1, N-1, n): mask[i] False # 每行的最后一个点其右侧没有内部邻居 diagonals.append(off_diag_1[mask]) offsets.append(1) diagonals.append(off_diag_1[mask]) offsets.append(-1) # 次次对角线左右邻居值为 -1/h^2 # 对应网格中相邻行的点偏移量为 ±n off_diag_n -1.0 / (h**2) * np.ones(N - n) diagonals.append(off_diag_n) offsets.append(n) diagonals.append(off_diag_n) offsets.append(-n) # 创建稀疏矩阵 A sp.diags(diagonals, offsets, shape(N, N), formatcsr) return A def poisson_rhs(m, f_func): 生成右端项b f_func是源项函数 f(x,y) n m - 2 N n * n h 1.0 / (m - 1) b np.zeros(N) # 遍历所有内部网格点 for i in range(n): for j in range(n): idx i * n j x (j 1) * h # 转换为物理坐标 y (i 1) * h b[idx] f_func(x, y) return b # 示例定义源项 f(x,y) 2*pi^2 * sin(pi*x) * sin(pi*y)其真解为 usin(pi*x)*sin(pi*y) def f_example(x, y): return 2 * (np.pi**2) * np.sin(np.pi * x) * np.sin(np.pi * y) m 50 # 网格点数 A poisson_matrix(m) b poisson_rhs(m, f_example)现在我们得到了一个 2304阶(50-2)^22304的稀疏正定矩阵A和右端向量b。直接法求解这个规模的问题已经有些吃力而迭代法正合适。4.2 迭代求解与性能观测我们将使用编写好的雅可比和高斯-赛德尔迭代函数来求解。注意我们的矩阵A是稀疏的但之前的函数实现是针对稠密矩阵的。为了高效利用稀疏性我们需要修改迭代核心使用稀疏矩阵乘法。Python稀疏矩阵迭代实现修改要点def jacobi_iter_sparse(A_sparse, b, x0None, tol1e-6, max_iter5000): 适用于稀疏矩阵的雅可比迭代 n len(b) if x0 is None: x np.zeros_like(b, dtypefloat) else: x x0.astype(float) x_new np.zeros_like(x) err_history [] # 提取对角线元素的倒数 d_inv 1.0 / A_sparse.diagonal() # 预计算 R*x 中的 R 矩阵即A减去对角线 # 在迭代中我们计算 R*x A*x - D*x其中D*x d * x (d是对角线元素) # 所以 x_new d_inv * (b - (A*x - d*x)) d_inv * (b - A*x d*x) # 但更清晰的做法是x_new x d_inv * (b - A*x) # 这是雅可比迭代的另一种等价形式有时称为Richardson迭代 for k in range(max_iter): r b - A_sparse x # 计算残差 x_new x d_inv * r # 更新解 err np.linalg.norm(x_new - x, ordnp.inf) err_history.append(err) if err tol: return x_new, k1, np.array(err_history) x x_new.copy() print(f未收敛于{max_iter}次迭代) return x, max_iter, np.array(err_history) def gauss_seidel_iter_sparse(A_sparse, b, x0None, tol1e-6, max_iter5000): 适用于稀疏矩阵的高斯-赛德尔迭代使用scipy.sparse.linalg.splu预处理 # 对于高斯-赛德尔稀疏实现更复杂因为(DL)是下三角矩阵。 # 一种高效做法是使用稀疏LU分解的预处理。 # 这里为简化我们暂时回退到使用稠密数组对于n2304尚可接受。 # 更专业的做法会使用如Successive Over-Relaxation (SOR)或共轭梯度法(CG)。 A_dense A_sparse.toarray() return gauss_seidel_iter(A_dense, b, x0, tol, max_iter)运行测试x0 np.zeros_like(b) tol 1e-8 max_iter 5000 print(求解泊松方程离散系统...) print(矩阵规模:, A.shape) # 雅可比迭代 print(\n--- 雅可比迭代 ---) x_j, iter_j, err_j jacobi_iter_sparse(A, b, x0, tol, max_iter) print(f迭代次数: {iter_j}, 最终误差: {err_j[-1]:.2e}) # 高斯-赛德尔迭代使用稠密矩阵 print(\n--- 高斯-赛德尔迭代 ---) x_gs, iter_gs, err_gs gauss_seidel_iter(A.toarray(), b, x0, tol, max_iter) print(f迭代次数: {iter_gs}, 最终误差: {err_gs[-1]:.2e}) # 计算与真实解的误差已知真解 n m-2 u_true np.zeros_like(b) for i in range(n): for j in range(n): idx i*n j x (j1)/(m-1) y (i1)/(m-1) u_true[idx] np.sin(np.pi*x) * np.sin(np.pi*y) err_true_j np.linalg.norm(x_j - u_true) / np.linalg.norm(u_true) err_true_gs np.linalg.norm(x_gs - u_true) / np.linalg.norm(u_true) print(f\n与真实解的相对误差:) print(f雅可比: {err_true_j:.2e}) print(f高斯-赛德尔: {err_true_gs:.2e})在我的测试中m50高斯-赛德尔迭代大约需要 ~3500 次迭代达到1e-8的精度而雅可比迭代需要 ~7000 次验证了高斯-赛德尔收敛更快大约快一倍的结论。同时两种方法最终得到的解与真实解析解的相对误差都在1e-4量级这主要来源于离散化误差而非迭代误差。4.3 关键陷阱与应对策略在实际使用这两种基本迭代法时有以下几个必须警惕的坑1. 对角元为零或过小迭代公式中需要除以对角线元素aii。如果某个aii为零算法直接崩溃如果aii的绝对值非常小除以它会导致数值不稳定误差被放大。对策在迭代前检查矩阵的对角优势。如果可能尝试对矩阵进行行/列重排如使用最大权重匹配来增强对角优势。对于零对角元必须通过行交换将其移走。2. 收敛速度极慢即使满足收敛条件对于条件数很大的矩阵即矩阵接近奇异雅可比和高斯-赛德尔的收敛速度可能慢得无法接受。对策这是基本迭代法的主要缺点。此时需要采用预处理技术Preconditioning。其思想是找到一个易于求逆的矩阵M称为预条件子使得 M^(-1)A 的条件数远小于A然后对等价方程组 M^(-1)Ax M^(-1)b 应用迭代法。例如雅可比预条件子就是取 M D对角部分这其实就是我们标准雅可比迭代的形式。更有效的预条件子有ILU、SSOR等。3. 稀疏矩阵存储与计算效率如我们所见直接使用稠密矩阵格式进行迭代会完全丧失稀疏性带来的内存和计算优势。对策务必使用稀疏矩阵存储格式如CSR, Compressed Sparse Row。在MATLAB中稀疏矩阵运算是自动优化的。在Python中要使用scipy.sparse库并确保矩阵运算如A x调用的是稀疏矩阵乘法例程。4. 终止准则的选择我们一直使用相邻两次迭代解的差||x^(k1) - x^(k)||作为终止判断。这通常有效但并非绝对可靠。有时这个差很小但解离真实解还很远这在迭代初期或矩阵病态时可能发生。一个更可靠的准则是看相对残差||b - Ax^(k)|| / ||b||。计算残差需要一次额外的矩阵-向量乘法但能更准确地反映当前解的质量。在实际代码中可以结合使用两种准则。5. 初始猜测解x0的影响初始猜测解x0越接近真实解所需的迭代次数越少。对于许多问题零向量是一个安全的起点。但如果能通过物理背景、上一次计算的结果或低精度解提供一个更好的初始猜测能显著节省计算时间。在时变问题中常用前一个时间步的解作为当前时间步迭代的初始值这称为“时间步进预热”。5. 拓展从基础到现代Krylov子空间方法雅可比和高斯-赛德尔迭代是定常迭代法Stationary Iterative Methods的代表其迭代矩阵B_J 或 B_GS在整个计算过程中保持不变。它们简单但收敛速度往往不够理想尤其对于现代科学计算中遇到的大型、病态问题。现代迭代法的核心是非定常迭代法其中最著名的一类是基于Krylov子空间的方法如共轭梯度法CG用于对称正定矩阵、广义最小残差法GMRES用于非对称矩阵。这些方法的关键思想是在第k步迭代不仅利用当前解而且利用前k步产生的全部信息残差向量张成的Krylov子空间在该子空间中寻找最优的近似解。这好比在寻找宝藏时不仅看当前这一步还综合分析了之前所有步伐的方向和距离从而能规划出更优的路径。一个直观对比雅可比/高斯-赛德尔像在迷雾中摸索每次只根据当前位置的局部信息决定下一步方向。共轭梯度法像带着地图和指南针探索它利用之前所有步伐的信息构建出一个关于地形即矩阵A的越来越精确的模型从而能以近乎最优的路径逼近目标。对于对称正定矩阵共轭梯度法的收敛速度远超雅可比或高斯-赛德尔迭代。它通常被视为求解此类问题的首选迭代法。在Python的scipy.sparse.linalg或MATLAB中都有高度优化的CG求解器如scipy.sparse.linalg.cg,pcg。那么为什么我们还要学习雅可比和高斯-赛德尔迭代呢原因有三教学价值它们是理解迭代法思想的完美起点其原理直观易于实现。预条件子的基础许多高效的预条件子如高斯-赛德尔迭代本身、SOR、SSOR的思想都源于这些基本方法。例如在CG法中常用SSOR作为预条件子来加速收敛。特定场景仍有价值在某些高度并行或特殊结构的计算中雅可比迭代因其完美的并行性仍被使用。而高斯-赛德尔迭代或其变种如SOR有时作为光滑器smoother用于多重网格法Multigrid Method中这是求解偏微分方程最快的方法之一。因此掌握雅可比和高斯-赛德尔迭代不仅是学习迭代法的第一步更是通向更高级、更高效数值算法的重要基石。当你下次面对一个大规模线性方程组时可以先尝试用这些基本方法进行快速原型验证分析问题的性质然后再决定是否需要调用更强大的“重型武器”。