1. 从直接法到迭代法为什么我们需要超松弛迭代在数值计算的世界里解线性方程组是绕不开的基石任务。无论是模拟物理场的有限元分析还是训练机器学习模型时的参数求解最终都常常归结为求解一个形如Ax b的方程组。对于小规模、稠密的矩阵我们习惯使用高斯消元法、LU分解这类直接法它们像一把精准的手术刀通过有限步的精确运算直接得到解。然而当矩阵A的规模膨胀到成千上万甚至更大并且是一个稀疏矩阵即绝大部分元素为零时直接法就显得力不从心了。原因有二一是计算复杂度高例如高斯消元法为 O(n³)内存消耗巨大二是直接法在消元过程中会破坏矩阵原有的稀疏结构产生大量非零的“填充元”导致内存爆炸。想象一下你要处理一个描述城市电网或社交网络的矩阵它可能99%的元素都是0用直接法去解相当于把一张极度稀疏的网强行填满这无疑是灾难性的。这时迭代法就登场了。它的核心思想不是去“硬算”精确解而是从一个初始猜测解x⁽⁰⁾出发通过一个迭代公式x⁽ᵏ⁺¹⁾ B x⁽ᵏ⁾ f构造一个序列 {x⁽ᵏ⁾}希望这个序列能收敛到方程组的真实解x*。雅可比迭代和高斯-赛德尔迭代是两种最基础的迭代法。简单来说雅可比迭代是“同时更新”即用第 k 步的所有分量去计算第 k1 步的所有分量而高斯-赛德尔迭代是“逐个更新”一旦计算出第 k1 步的某个新分量就立刻用它替换旧值参与后续分量的计算。高斯-赛德尔通常比雅可比收敛更快因为它利用了最新的信息。那么超松弛迭代法又是什么你可以把它理解为高斯-赛德尔迭代法的“加速增强版”。在实际计算中尤其是对于某些特定类型的矩阵如对角占优矩阵高斯-赛德尔迭代产生的序列收敛方向是正确的但步幅可能太“保守”导致收敛速度依然不够理想。超松弛迭代法引入了一个松弛因子 ω其核心操作是在得到高斯-赛德尔迭代的“预测值”后并不全盘接受而是将其与上一步的值做一个加权平均。这个 ω 就像汽车的油门和刹车的混合体当 ω1 时它就是标准的高斯-赛德尔迭代当 0ω1 时它是在“刹车”称为低松弛有时用于改善不收敛的迭代格式当 1ω2 时它是在“加油门”称为超松弛目的是放大修正量加速收敛。我们今天重点讨论的正是这个能显著提速的逐次超松弛迭代法。2. 逐次超松弛迭代法的数学原理与收敛性分析理解了迭代法的动机我们深入到SOR方法的数学内核。给定线性方程组Ax b其中A是一个 n×n 的非奇异矩阵。我们将A分解为三个部分A D - L - U。这里D是A的对角线元素构成的对角矩阵-L是A的严格下三角部分即对角线以下包含负号-U是A的严格上三角部分即对角线以上包含负号。这个分解是理解所有基本迭代法的钥匙。高斯-赛德尔迭代的矩阵形式可以写为(D - L) x⁽ᵏ⁺¹⁾ U x⁽ᵏ⁾ b。解出x⁽ᵏ⁺¹⁾的表达式为x⁽ᵏ⁺¹⁾ (D - L)⁻¹ (U x⁽ᵏ⁾ b)。这个公式清晰地展示了“逐个更新”的特性因为(D - L)是下三角矩阵其逆矩阵作用在向量上时计算x⁽ᵏ⁺¹⁾的第 i 个分量只会用到已经更新过的前 i-1 个新分量和未更新的后 n-i 个旧分量。现在我们引入松弛因子 ω。SOR方法的迭代公式可以从分量形式直观理解。对于第 i 个未知数我们先计算高斯-赛德尔迭代给出的新值预测值x_gs_i^(k1) (b_i - Σ_{j1}^{i-1} a_ij * x_j^(k1) - Σ_{ji1}^{n} a_ij * x_j^(k)) / a_ii然后我们不直接把x_i^(k1)设为x_gs_i^(k1)而是将其与上一步的值x_i^(k)做一个加权组合x_i^(k1) (1 - ω) * x_i^(k) ω * x_gs_i^(k1)将x_gs_i^(k1)的表达式代入就得到了SOR方法的分量形式迭代公式x_i^(k1) (1-ω) * x_i^(k) (ω / a_ii) * (b_i - Σ_{j1}^{i-1} a_ij * x_j^(k1) - Σ_{ji1}^{n} a_ij * x_j^(k))将其写成矩阵形式经过推导可得(D - ωL) x⁽ᵏ⁺¹⁾ [(1-ω)D ωU] x⁽ᵏ⁾ ωb因此SOR方法的迭代矩阵T_ω为T_ω (D - ωL)⁻¹ [(1-ω)D ωU]迭代格式为x⁽ᵏ⁺¹⁾ T_ω x⁽ᵏ⁾ (D - ωL)⁻¹ (ωb)。收敛性是迭代法的生命线。一个迭代法是否有效取决于迭代矩阵T的谱半径 ρ(T)即所有特征值模的最大值。收敛的充分必要条件是ρ(T) 1且谱半径越小收敛速度通常越快。对于SOR方法有以下关键结论必要条件要使SOR方法可能收敛松弛因子 ω 必须落在区间 (0, 2) 内。如果 ω 不在此区间方法必然发散。充分条件如果系数矩阵A是对称正定的那么对于任意 0ω2SOR方法都收敛。这是一个非常强大且实用的定理因为工程和科学计算中很多问题如泊松方程离散化产生的矩阵都是对称正定的。最优松弛因子 ω_opt对于一类特殊的矩阵如具有性质A的矩阵许多椭圆型偏微分方程离散后得到存在一个最优的松弛因子 ω_opt使得SOR方法的收敛速度达到最快谱半径最小。这个最优值通常大于1即超松弛。ω_opt 的理论值可以通过原矩阵的雅可比迭代矩阵的谱半径 ρ(J) 来计算ω_opt 2 / (1 sqrt(1 - ρ(J)²))。但在实际中ρ(J) 往往难以精确获得因此 ω_opt 的确定常常依赖于经验估计或试算。注意松弛因子 ω 的选择是SOR方法的灵魂也是最大的调参难点。ω1 退化为高斯-赛德尔ω 略大于1如1.1~1.5常能加速但 ω 过大接近2极易导致发散。对于非对称正定矩阵寻找 ω_opt 更像一门艺术。3. 算法实现从数学公式到可执行代码理论清晰之后我们着手实现。我们将分别用 MATLAB 和 Python 实现SOR迭代法并比较两者在代码风格和效率上的特点。我们将实现一个通用的函数sor_solver它接受系数矩阵A、右端向量b、松弛因子ω、初始猜测x0、最大迭代次数max_iter和容差tol作为输入返回近似解和迭代信息。3.1 MATLAB 实现与逐行解析MATLAB 以其简洁的矩阵语法著称非常适合原型算法开发。以下是完整的实现function [x, iter, rel_err_history] sor_solver(A, b, omega, x0, max_iter, tol) % 逐次超松弛迭代法求解 Ax b % 输入 % A: n x n 系数矩阵 % b: n x 1 右端向量 % omega: 松弛因子 (0 omega 2) % x0: n x 1 初始向量 % max_iter: 最大迭代次数 % tol: 相对误差容差 % 输出 % x: 近似解向量 % iter: 实际迭代次数 % rel_err_history: 每次迭代的相对误差历史记录可选用于分析收敛 n length(b); x x0(:); % 确保x是列向量 rel_err_history zeros(max_iter, 1); % 预分配内存记录误差 % 检查输入有效性 if omega 0 || omega 2 error(松弛因子 omega 必须在区间 (0, 2) 内。); end if size(A,1) ~ n || size(A,2) ~ n error(矩阵 A 的维度与向量 b 不匹配。); end % 主迭代循环 for iter 1:max_iter x_old x; % 保存上一步的解 for i 1:n % 计算 sigma1: 使用已更新的新分量 (i之前) sigma1 A(i, 1:i-1) * x(1:i-1); % 计算 sigma2: 使用尚未更新的旧分量 (i之后) sigma2 A(i, i1:n) * x_old(i1:n); % SOR迭代公式的核心计算 x(i) (1 - omega) * x_old(i) (omega / A(i, i)) * (b(i) - sigma1 - sigma2); end % 计算相对误差使用无穷范数 ||x_new - x_old||_inf / ||x_new||_inf rel_err norm(x - x_old, inf) / (norm(x, inf) eps); % 加eps防止除零 rel_err_history(iter) rel_err; % 收敛判断 if rel_err tol rel_err_history rel_err_history(1:iter); % 截断记录 fprintf(SOR迭代在 %d 次迭代后收敛相对误差为 %.4e。\n, iter, rel_err); return; end end % 如果达到最大迭代次数仍未收敛 warning(达到最大迭代次数 %d 仍未收敛当前相对误差为 %.4e。, max_iter, rel_err); end代码要点与避坑指南向量化与循环的取舍SOR方法的本质是顺序更新第i个分量依赖于前i-1个已更新的新分量这阻碍了完全的向量化。因此内层对i的循环是必须的。但在计算sigma1和sigma2时我们使用了矩阵切片乘法A(i, 1:i-1) * x(1:i-1)这依然是向量化操作比用内层循环累加快得多。误差范数的选择我们使用无穷范数norm(..., inf)计算相对误差。无穷范数是向量元素绝对值中的最大值它能敏锐地反映单个分量最大的变化对于判断收敛很直观。你也可以使用2-范数欧几里得范数但计算量稍大。防止除零在计算相对误差的分母norm(x, inf)时我们加了一个极小的数eps。这是为了防止在解非常接近零向量时分母为零导致计算错误。这是一个稳健性编程的小技巧。历史记录rel_err_history记录了每次迭代的误差这对于后续绘制收敛曲线、分析不同 ω 的影响至关重要。预分配数组zeros(max_iter, 1)是MATLAB性能优化的关键习惯能避免循环中动态调整数组大小带来的巨大开销。3.2 Python (NumPy) 实现与风格对比Python 结合 NumPy 库在科学计算领域与MATLAB分庭抗礼。其语法更通用生态更丰富。以下是等效的Python实现import numpy as np def sor_solver(A, b, omega, x0None, max_iter1000, tol1e-6): 逐次超松弛迭代法求解 Ax b。 参数 A : numpy.ndarray n x n 系数矩阵。 b : numpy.ndarray n x 1 右端向量。 omega : float 松弛因子应在 (0, 2) 区间内。 x0 : numpy.ndarray, 可选 初始猜测向量。默认为零向量。 max_iter : int, 可选 最大迭代次数。默认为1000。 tol : float, 可选 收敛容差基于无穷范数的相对误差。默认为1e-6。 返回 x : numpy.ndarray 近似解向量。 iter_count : int 实际迭代次数。 rel_err_history : list 每次迭代的相对误差记录。 n len(b) if x0 is None: x np.zeros_like(b, dtypenp.float64) else: x x0.astype(np.float64).copy() # 确保为浮点数并创建副本 x_old x.copy() rel_err_history [] # 输入验证 if omega 0 or omega 2: raise ValueError(f松弛因子 omega 必须在 (0, 2) 区间内当前为 {omega}。) if A.shape ! (n, n): raise ValueError(f矩阵 A 的形状 {A.shape} 与向量 b 的长度 {n} 不匹配。) # 主迭代循环 for iter_count in range(1, max_iter 1): x_old[:] x # 保存当前解到旧变量 for i in range(n): # 使用点积计算求和项更简洁 sigma np.dot(A[i, :i], x[:i]) np.dot(A[i, i1:], x_old[i1:]) # SOR迭代公式 x[i] (1 - omega) * x_old[i] (omega / A[i, i]) * (b[i] - sigma) # 计算相对误差无穷范数 rel_err np.linalg.norm(x - x_old, ordnp.inf) / (np.linalg.norm(x, ordnp.inf) np.finfo(float).eps) rel_err_history.append(rel_err) # 收敛判断 if rel_err tol: print(fSOR迭代在 {iter_count} 次迭代后收敛相对误差为 {rel_err:.4e}。) return x, iter_count, rel_err_history # 未收敛警告 print(f警告达到最大迭代次数 {max_iter} 仍未收敛最终相对误差为 {rel_err_history[-1]:.4e}。) return x, max_iter, rel_err_historyPython实现的关键差异与技巧数据类型管理Python是动态类型但NumPy数组有明确的dtype。我们使用x0.astype(np.float64).copy()来确保计算使用高精度的浮点数并避免对输入参数的意外修改。.copy()是防止视图共享导致错误的良好习惯。向量化计算与MATLAB类似我们使用np.dot进行向量内积来计算sigma。A[i, :i]和x[:i]的切片操作非常高效。注意Python的切片是左闭右开区间。误差计算中的微小量np.finfo(float).eps获取了双精度浮点数的机器精度作用同MATLAB中的eps用于稳定分母。历史记录存储Python中使用列表list来动态存储rel_err_history这比预分配NumPy数组更符合Python习惯因为列表的append操作在预期内是高效的。如果迭代次数非常多如10万预分配NumPy数组可能更有优势。输出与打印使用f-string进行格式化输出更现代、清晰。两种语言实现的对比小结语法MATLAB的矩阵操作语法更数学化如直接A(i, 1:i-1) * x(1:i-1)PythonNumPy更接近通用编程np.dot(A[i, :i], x[:i])。性能对于这种以逐分量循环为核心的算法两者性能在同等优化下相差不大。但对于可完全向量化的操作MATLAB的底层优化有时更激进。Python的优势在于其庞大的生态库如SciPy提供了更专业的迭代求解器scipy.sparse.linalg.spsolve。部署与生态MATLAB在控制系统、信号处理等传统工程领域有优势Python则在机器学习、数据科学和Web集成方面更胜一筹。4. 实战案例求解泊松方程模型问题为了验证我们的SOR求解器并直观感受 ω 的影响我们构造一个经典的测试问题二维泊松方程在正方形区域上的五点差分离散。这个问题产生的矩阵是对称正定的、稀疏的、对角占优的是迭代法的理想测试床。考虑方程-∇²u f在单位正方形 [0,1]×[0,1] 上边界条件为 u0。我们用均匀网格进行离散网格点数为 m×m内部网格点数为 n (m-1)×(m-1)。标准的五点差分格式会生成一个块三对角矩阵A其每行至多有5个非零元素。右端项b由源项 f 和边界条件决定。我们编写一个函数来生成这个问题的矩阵和右端项以 f1 为例即求解 -∇²u 1。import numpy as np from scipy.sparse import diags, csr_matrix from scipy.sparse.linalg import spsolve, norm as spnorm import time import matplotlib.pyplot as plt def generate_poisson_matrix(m): 生成二维泊松方程-∇²u f五点差分离散后的矩阵A稀疏格式和右端项bf1。 m: 每个方向的网格点数包括边界内部点数为 (m-2)*(m-2) 返回稀疏矩阵A_csr右端向量b n_inner m - 2 # 每个维度内部点数 n n_inner * n_inner # 总未知数个数 h 1.0 / (m - 1) # 网格间距 # 主对角线元素为4/h^2 main_diag np.full(n, 4.0 / (h * h)) # 上下次对角线元素为 -1/h^2对应x方向相邻点 off_diag_x np.full(n - 1, -1.0 / (h * h)) # 将每块最后一行对应网格最右边点的x方向连接置零 for i in range(n_inner): off_diag_x[(i1)*n_inner - 1] 0.0 # 远近次对角线元素为 -1/h^2对应y方向相邻点距离为n_inner off_diag_y np.full(n - n_inner, -1.0 / (h * h)) diagonals [main_diag, off_diag_x, off_diag_x, off_diag_y, off_diag_y] offsets [0, 1, -1, n_inner, -n_inner] A_sparse diags(diagonals, offsets, shape(n, n), formatcsr) b np.ones(n) # 源项 f 1 return A_sparse, b # 生成一个中等规模的问题 (31x31网格内部29x29841个未知数) m 31 A_sp, b generate_poisson_matrix(m) n A_sp.shape[0] print(f问题规模: {n} 个未知数矩阵非零元素: {A_sp.nnz}) # 使用SciPy的直接稀疏求解器得到“精确解”作为参考 print(使用SciPy直接稀疏求解器计算参考解...) x_exact spsolve(A_sp, b) # 这是一个高效的直接法求解器用于参考 # 测试我们的SOR求解器 omega_values [1.0, 1.2, 1.5, 1.7, 1.9] # 测试不同的松弛因子 max_iter 5000 tol 1e-8 x0 np.zeros(n) results {} for omega in omega_values: print(f\n--- 测试 ω {omega} ---) start_time time.time() x_sor, iter_count, err_history sor_solver(A_sp.toarray(), b, omega, x0, max_iter, tol) # 注意将稀疏矩阵转为稠密矩阵以适配我们的函数 elapsed_time time.time() - start_time # 计算与参考解的绝对误差无穷范数 abs_error np.linalg.norm(x_sor - x_exact, ordnp.inf) results[omega] { iterations: iter_count, time: elapsed_time, final_rel_err: err_history[-1], abs_error_to_exact: abs_error, err_history: err_history } print(f 迭代次数: {iter_count}, 耗时: {elapsed_time:.4f}秒, 与精确解绝对误差: {abs_error:.4e}) # 绘制收敛曲线 plt.figure(figsize(10, 6)) for omega, data in results.items(): plt.semilogy(data[err_history], labelfω{omega}, linewidth2) plt.xlabel(迭代次数, fontsize12) plt.ylabel(相对误差 (log scale), fontsize12) plt.title(SOR方法不同松弛因子下的收敛历史, fontsize14) plt.grid(True, whichboth, linestyle--, alpha0.7) plt.legend(fontsize11) plt.tight_layout() plt.show() # 寻找近似最优ω通过网格搜索 print(\n 网格搜索寻找近似最优ω ) search_omegas np.linspace(1.0, 1.95, 20) # 在1.0到1.95间搜索 best_omega 1.0 best_iter max_iter search_results [] for omega in search_omegas: _, iter_count, _ sor_solver(A_sp.toarray(), b, omega, x0, max_iter2000, tol1e-6) search_results.append((omega, iter_count)) if iter_count best_iter and iter_count 2000: # 确保是收敛的 best_iter iter_count best_omega omega print(fω{omega:.3f}, 迭代次数{iter_count}) print(f\n近似最优松弛因子 ω_opt ≈ {best_omega:.3f} 所需迭代次数 {best_iter})案例运行结果分析当你运行上述代码很可能会观察到以下现象ω1.0高斯-赛德尔收敛稳定但速度中等。ω≈1.2~1.7收敛速度显著加快迭代次数可能减少到高斯-赛德尔的一半甚至更少。曲线在图中下降得更陡峭。ω接近2如1.9收敛可能变得不稳定误差曲线可能出现震荡甚至无法在最大迭代次数内收敛。这印证了 ω 必须小于2的理论条件。最优 ω对于这个特定的泊松问题理论最优 ω_opt 可以通过其雅可比迭代矩阵的谱半径估算。我们的网格搜索通常会找到一个在1.7~1.8之间的值使得迭代次数最少。这个值大于1证实了超松弛的加速效果。实操心得对于生产环境中的大规模问题将稀疏矩阵转为稠密矩阵A_sp.toarray()是极其低效且耗内存的这里只是为了适配我们之前写的稠密矩阵求解器。一个真正实用的SOR求解器应该直接处理稀疏矩阵的存储格式如CSR在循环中只访问非零元素。这需要对代码中的求和部分sigma的计算进行重写利用稀疏矩阵的行索引和列索引指针。这是性能优化的关键一步也是从教学代码到工业级代码的飞跃。5. 性能调优、常见陷阱与进阶讨论掌握了基础实现后我们探讨如何让SOR方法更快、更稳并避开一些常见的坑。5.1 稀疏矩阵的高效实现如前所述对于大规模问题必须使用稀疏矩阵。修改我们的Python求解器以支持SciPy的CSR格式def sor_solver_sparse(A_csr, b, omega, x0None, max_iter1000, tol1e-6): 针对CSR稀疏矩阵优化的SOR求解器 n len(b) if x0 is None: x np.zeros_like(b, dtypenp.float64) else: x x0.astype(np.float64).copy() x_old x.copy() rel_err_history [] # 获取CSR格式的数据 data A_csr.data indices A_csr.indices indptr A_csr.indptr for iter_count in range(1, max_iter 1): x_old[:] x for i in range(n): # 遍历第i行的非零元素 sigma 0.0 row_start indptr[i] row_end indptr[i1] diag_val 0.0 for j_idx in range(row_start, row_end): j indices[j_idx] a_ij data[j_idx] if j i: diag_val a_ij elif j i: sigma a_ij * x[j] # 使用已更新的新值 else: # j i sigma a_ij * x_old[j] # 使用旧值 if diag_val 0: raise ValueError(f矩阵A的第{i}行对角线元素为零SOR迭代无法进行。) x[i] (1 - omega) * x_old[i] (omega / diag_val) * (b[i] - sigma) rel_err np.linalg.norm(x - x_old, ordnp.inf) / (np.linalg.norm(x, ordnp.inf) np.finfo(float).eps) rel_err_history.append(rel_err) if rel_err tol: print(f稀疏SOR迭代在 {iter_count} 次迭代后收敛。) return x, iter_count, rel_err_history print(f警告稀疏SOR达到最大迭代次数 {max_iter}。) return x, max_iter, rel_err_history这个版本只对非零元素进行操作内存占用和计算量大幅降低是处理大规模稀疏问题的正确方式。5.2 松弛因子 ω 的选取策略ω 的选择是SOR方法的精髓也是主要调参负担。理论指导对于具有性质A一致有序、2-循环相容的矩阵最优松弛因子公式为ω_opt 2 / (1 sqrt(1 - ρ(J)²))其中 ρ(J) 是雅可比迭代矩阵的谱半径。对于矩形区域上的泊松方程若网格点数为 M×N有近似公式ρ(J) ≈ cos(π/M) cos(π/N)。经验法则对于许多实际问题ω 在1.2 到 1.6之间往往能取得不错的效果。可以从1.1开始以0.1为步长增加观察收敛速度的变化。自适应松弛更高级的策略是使用自适应松弛因子在迭代过程中根据残差的变化动态调整 ω。但这需要更复杂的逻辑和监控。诊断工具绘制收敛曲线误差 vs. 迭代次数是最直观的诊断工具。如果曲线平滑下降说明 ω 合适如果剧烈震荡说明 ω 可能过大如果下降极其缓慢ω 可能偏小或矩阵性质不佳。5.3 预处理技术与SOR的结合单纯的SOR对于条件数很大的病态矩阵可能收敛很慢甚至不收敛。这时可以引入预处理技术。预处理的思想是找一个非奇异矩阵M使得新方程组M⁻¹Ax M⁻¹b的系数矩阵M⁻¹A的条件数更好或者特征值分布更集中从而加速迭代收敛。SOR本身可以看作是一种预处理其迭代矩阵依赖于(D - ωL)。更强大的预处理子如不完全LU分解 (ILU)常与Krylov子空间方法如GMRES, BiCGSTAB结合。虽然SOR通常不直接与这些外部预处理子结合但理解预处理的概念有助于你在更复杂的求解场景中选择合适的工具。5.4 常见陷阱与调试技巧不收敛检查 ω首先确认 ω 是否在 (0,2) 区间内。对于非对称正定矩阵即使 ω 在此区间也可能不收敛。检查矩阵性质确保矩阵是行对角占优或至少不可约对角占优。可以使用np.all(np.abs(np.diag(A)) np.sum(np.abs(A), axis1) - np.abs(np.diag(A)))来检查严格行对角占优。对于对称矩阵检查是否正定所有特征值大于0。检查代码仔细核对迭代公式特别是下三角部分用的是新值x[:i]还是旧值x_old[:i]。这是高斯-赛德尔与雅可比的核心区别也是常见错误点。收敛速度慢调整 ω进行网格搜索寻找更优的 ω。改进初始猜测如果有一个接近真实解的初始猜测x0可以大幅减少迭代次数。例如可以用上一次计算的结果作为下一次的初值或者在多尺度计算中用粗网格的解插值到细网格作为初值。考虑其他迭代法对于某些问题共轭梯度法对称正定或GMRES非对称可能比SOR快得多尤其是结合了好的预处理子之后。除零错误确保矩阵对角线元素a_ii没有零。在迭代前可以检查if np.any(np.diag(A) 0): ...。如果存在零对角元通常需要对矩阵进行行重排如使用对角线优势化或矩阵重排序技术或者选择不需要对角线元素的迭代法如雅可比法需要但高斯-赛德尔和SOR需要。性能瓶颈对于大规模问题务必使用稀疏矩阵格式并优化内层循环避免在Python层进行大量低效的循环。可以考虑使用Numba对核心循环进行即时编译或者用Cython/C编写扩展模块。对于超大规模问题SOR这类点迭代法可能依然不够快需要考虑更先进的Krylov子空间方法或代数多重网格法。