简介面向数值计算课程与科研工程应用这份压缩包提供了雅可比迭代法和高斯-赛德尔迭代法的完整MATLAB实现帮助学习者快速解决大型线性方程组的近似求解问题。包内共5个文件全部为m格式源程序涵盖两种算法的函数定义与调用主程序便于直接运行、修改和对比迭代效果整个资源包仅3KB结构精简适合实验和作业场景。已有1967人学习这份资源既是课程设计的直接参考也能辅助理解两种迭代法在矩阵分裂、逐分量更新以及收敛速度上的差异。通过阅读和运行源码可掌握从系数矩阵分解、迭代公式编写到残差收敛判别的完整流程同时了解雅可比迭代适合并行计算、高斯-赛德尔迭代通常更快且更稳定的实际选型依据为后续处理大规模稀疏方程组打下扎实基础。1. 为什么线性方程组的迭代求解从雅可比和高斯-赛德尔开始求解稀疏线性方程组时迭代法比直接法更省内存这一点在真实工程里经常让新手意外几百阶矩阵直接消元还能撑住一旦到万阶以上一次直接法的内存和时间代价就会失控。雅可比迭代法与高斯-赛德尔迭代法正是迭代法里最基本、也最常被选作起点的两个方法。如果一个压缩包命名成这种里面通常跑的就是这两种方法的最小可运行实现。这篇笔记围绕同一个目标展开搞懂两种格式的差异在哪里用什么参数能让收敛快一点以及实际跑出问题时该去查哪里。2. 从矩阵分裂到谱半径雅可比与高斯-赛德尔为什么一个慢一个快2.1 雅可比与高斯-赛德尔的矩阵分裂与更新顺序两种方法都建立在同一个矩阵分裂上。把线性方程组 Ax b 的系数矩阵拆成 A D - L - U其中 D 是对角矩阵L 是严格下三角矩阵U 是严格上三角矩阵。这个拆分本身不改变方程的解改变的是怎么从当前近似解推进到下一步近似解。雅可比迭代法的格式是 x_new D^{-1}(b (LU)x_old)也就是下一轮迭代时所有分量都只用上一轮的值计算。这个特性让雅可比在并行环境里天然友好每个分量的计算互不依赖。高斯-赛德尔迭代法的格式是 x_new (D-L)^{-1}(b Ux_old)写成分量形式更直观求第 i 个分量时前 i-1 个分量已经更新过直接用新值后面的分量还用旧值。这种“边算边用”的串行更新让每一轮推进的信息量更大。两者的选择逻辑也因此不同追求可并行、实现简单、对初值鲁棒选雅可比追求串行环境下更快的收敛选高斯-赛德尔。一个容易误解的点是高斯-赛德尔在每个分量上比雅可比多做了更新但单轮的浮点运算代价没有质的差别差别主要体现在达到相同精度所需的轮数上。对比项雅可比迭代法高斯-赛德尔迭代法更新使用的值全部用上一轮的值已更新的分量用新值未更新的用旧值分量间耦合无天然可并行有串行依赖存储要求需要保存新旧两轮向量一个向量原地更新即可对称正定问题收敛轮数约为高斯-赛德尔的 2 倍更少2.2 收敛性判据谱半径、对角占优与收敛速度的直觉迭代法看起来是无限逼近实际是否收敛取决于迭代矩阵的谱半径。把两种方法统一写成 x_{k1} M x_k c 的形式收敛的充要条件是迭代矩阵 M 的谱半径 \rho(M) 1。雅可比对应的迭代矩阵是 D^{-1}(LU)高斯-赛德尔对应的是 (D-L)^{-1}U。工程里不会真的去手算谱半径但有两条硬性结论可以直接用如果 A 严格对角占优即每个对角元素的绝对值严格大于该行非对角元素绝对值之和那么雅可比和高斯-赛德尔都收敛如果 A 不可约且弱对角占优结论同样成立。注意“不可约”这个条件一维泊松离散矩阵恰好就是这种典型内层行对角占优达到临界但没有严格大于它依然能保证收敛。收敛速度的直觉来自谱半径和 1 的距离。若 \rho(M) 是 0.9每轮误差大约缩小到上一轮的 90%若 \rho(M) 是 0.999每轮后误差几乎不动。高斯-赛德尔看起来“更快”本质上是它的迭代矩阵谱半径在对称正定问题下通常约等于雅可比谱半径的平方0.999^2 比 0.999 离 1 更远所以收敛轮数接近减半。这个比例关系经常被误解成普适结论实际上只在一类问题里成立后文避坑章节会专门说。2.3 一个工程习惯先估算迭代矩阵的谱半径再选方法我一般会在动手写大矩阵之前先用一个小规模矩阵做谱半径估算。用幂法就能得到近似值不需要求特征值分解。给定系数矩阵 A先构造两个迭代矩阵再用随机向量反复乘观察范数增长或衰减的幅度。import numpy as np def build_jacobi_matrix(A): D np.diag(np.diag(A)) L -np.tril(A, -1) U -np.triu(A, 1) return np.linalg.inv(D) (L U) def spectral_radius_estimate(M, seed7, max_iter5000): rng np.random.default_rng(seed) v rng.normal(sizeM.shape[0]) # 幂法反复迭代每次归一化模长的变化趋势指向主特征值 for _ in range(max_iter): Mv M v norm np.linalg.norm(Mv, ordnp.inf) if norm 1e-15: return 0.0 v Mv / norm return float(np.linalg.norm(M v, ordnp.inf))这段代码里build_jacobi_matrix严格按 A D - L - U 的符号组装L和U都带负号。幂法每次把向量乘到迭代矩阵上并归一化若模长持续增长说明主特征值大于 1若衰减说明小于 1。注意幂法依赖于初始向量在主导特征向量方向上有分量所以我会换两个不同随机种子多跑一次结果一致才采信。谱半径估计的目的不是精确到小数点后几位而是判断这个矩阵适不适合做迭代求解以及大概需要多少轮。3. 把迭代写成能跑的代码简单但好用的雅可比与高斯-赛德尔实现3.1 用 numpy 实现两个迭代器收敛判定放在残差上而不是解上写迭代求解器的第一原则是收敛判据必须建立在残差上不能建立在相邻两步解的差值上。相邻两步解变化很小可能是因为矩阵条件数很差、解空间很平也可能因为真正收敛了这两者从“解差”里分不清。所以下面的实现里回归停机的判断用的是 ||b - Ax||。import numpy as np def jacobi(A, b, x0None, tol1e-8, maxit20000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n len(b) x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) diag np.diag(A) off A - np.diagflat(diag) # L U去掉对角后的剩余部分 for k in range(1, maxit 1): x_new (b - off x) / diag # 所有分量同时用旧值更新 r_norm np.linalg.norm(b - A x_new, ordnp.inf) if r_norm tol: return x_new, k, r_norm x x_new return x, maxit, np.linalg.norm(b - A x, ordnp.inf)diag是主对角元素组成的向量off是去掉对角后的矩阵off x相当于 (LU)x。(b - off x) / diag直接把向量除到对角元素上在 numpy 里是逐元素除法。每次更新后立刻用新向量算残差用无穷范数比较。注意这里如果对角元素有接近 0 的值会产生除零警告后文避坑部分会讲如何提前发现。高斯-赛德尔的实现比雅可比多一个内层循环因为每个分量更新依赖前面刚算出来的新值def gauss_seidel(A, b, x0None, tol1e-8, maxit20000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n len(b) x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) for k in range(1, maxit 1): x_old x.copy() # 必须保留旧值给未更新的分量使用 for i in range(n): s b[i] - A[i, :i] x[:i] - A[i, i1:] x_old[i1:] x[i] s / A[i, i] r_norm np.linalg.norm(b - A x, ordnp.inf) if r_norm tol: return x, k, r_norm return x, maxit, np.linalg.norm(b - A x, ordnp.inf)内层循环中A[i, :i] x[:i]用的是已更新的新分量A[i, i1:] x_old[i1:]用的是上一轮旧分量。x_old x.copy()这一步不能省因为 x 更新过程中会覆盖掉旧值后面的行还要拿旧值参与计算。几乎所有从雅可比改写成高斯-赛德尔翻车的人都栽在这一行上。这里还藏着两个实现细节。第一个是收敛判定放在每次完整扫描之后因为一轮内部分量更新并不代表整个方程已经近似满足第二个是maxit到期后返回当前结果但不抛异常我习惯把残差一并返回方便调用方判断是“收敛”还是“用满次数”。3.2 一维泊松测试算例同一套矩阵、右端项对比迭代轮数光看代码不跑算例很难建立对两种方法差异的直觉。一维泊松方程离散出来的三对角矩阵是最适合做对比的测试台结构简单、性质清晰、收敛情况可控。def poisson_1d(n): h2 (n 1) ** 2 A np.zeros((n, n)) for i in range(n): A[i, i] 2.0 * h2 if i 0: A[i, i - 1] -h2 if i n - 1: A[i, i 1] -h2 return A n 100 A poisson_1d(n) b np.ones(n) x_j, it_j, r_j jacobi(A, b, tol1e-6, maxit100000) x_g, it_g, r_g gauss_seidel(A, b, tol1e-6, maxit100000) print(Jacobi:, it_j, r_j) print(Gauss-Seidel:, it_g, r_g)在这个 n100 的算例里雅可比达到 tol1e-6 大约需要三万轮高斯-赛德尔大约只需要一万五千轮。这里的三对角矩阵对角元素是 2h²相邻元素是 -h²内部行对角占优只是临界成立不是严格占优它能保证收敛靠的是“不可约且弱对角占优”这条定理。把 n 增大到 400谱半径会更接近 1两种方法的轮数会按 n² 的规模往上翻这是线性收敛的正常表现不是程序挂死。跑完这个算例再看代码你会理解一个事实迭代轮数不是唯一标准。每轮的计算量同样重要雅可比每轮可以做向量化矩阵乘法高斯-赛德尔的内层循环在 Python 里是串行推导n 很大时单轮耗时可能明显更久。实际工程里谁更快永远是“轮数”和“每轮耗时”的乘积不能只看单方面。3.3 数值试验的 3 条读结果规则跑出数字之后不要急着认定结论先用三条规则过一遍。第一条用直接法兜底。对这个小规模算例直接调np.linalg.solve(A, b)拿到参考解和迭代结果比一下最大绝对误差。如果参考解与迭代解对不上多数情况不是迭代法不收敛而是手写代码里某个索引或符号写反了。第二条看残差下降曲线是否单调。打印每轮的||b - Ax_k||_∞如果稳定下降并进入线性区间说明迭代正常如果残差波动甚至上升基本可以断定矩阵性质不支持当前方法不要盲目加迭代次数。第三条检查初值和右端项的量级。很多“不收敛”其实是容差设得太苛刻比如矩阵和 b 的量级在 1e10却设了 1e-12 的绝对容差导致永远无法满足。先标准化再做迭代是减少这类误判的有效手段。4. 三个必调参数与收敛判定初值、容差、最大迭代次数4.1 迭代初值零向量不总是最差但也别只会给零从零向量出发是默认选择对很多对称正定问题来说这个选择并不会让收敛轮数显著变差。实际工程里更值得利用的是“上一个状态的解”。比如瞬态问题里上一时刻的解和当前时刻的解往往很接近把上一时刻的解作为初值迭代轮数可能缩减一个数量级这种收益比调任何容差都来得直接。另一个常见做法是用物理常识约束初值。温度场问题里初值取环境温度而不是全零结构位移问题里初值取上一载荷步结果都会让前几十轮迭代少走弯路。我一般会在函数入口允许传x0并默认置零这样调用方可以根据问题类型决定是否给一个更聪明的起点。如果初值给得太离谱比如物理上不可能出现的负温度或超大位移反而会让迭代前段出现大幅摆动。所以初值的关键不是“越接近就好”而是“不能显著违背物理约束”。4.2 残差判据与容差用 ||b-Ax_k||_∞ 替代相邻解差收敛判据的选择直接决定“这个解能不能用”。相邻两步解的差已经分析过不可靠更隐蔽的问题是绝对容差和相对容差。假设你的方程组是 x_1 x_2 1、x_1 - x_2 2残差降低到 1e-8 很容易但同样的容差放到边界元法生成的矩阵上可能跑一万轮也达不到。解决办法是使用相对残差判据||b - Ax_k|| / ||b||把量级差异消除掉。实际判断顺序是先算初始残差r0 ||b - Ax0||每轮算rk ||b - Ax_k||当rk / r0小于阈值时认定收敛。这样矩阵和右端项放大或缩小一百倍判定结果不受影响。容差的具体数值单精度场景通常取 1e-5 或 1e-6双精度场景取 1e-8 到 1e-10 已经足够。再往下压到 1e-12很多时候只是把迭代轮数白白拉长解的精度并不会跟着提高因为浮点舍入本身已经限制了可实现精度。4.3 最大迭代次数的设置逻辑线性收敛怎么预估轮数线性迭代的误差大致按几何级数衰减每轮误差约为上一轮的 \rho 倍。如果谱半径 \rho0.999每轮大约只缩小千分之一要让误差从 1 掉到 1e-6需要大约 1.4 万轮。这个估算公式是所需轮数约为ln(容差 / 初始误差) / ln(rho)用自然对数就算出来了。所以设置最大迭代次数之前先估算谱半径然后按这个公式给上限加一个安全系数。比如估算需要一万五千轮maxit就设成六万留足余量。如果残差曲线已经进入线性下降区间却因为maxit设小了被迫退出这是最可惜的浪费。还有一种情况是残差曲线看起来在下降但速度极慢比如从 1e-2 到 1e-4 花了一万轮到 1e-6 可能还要两万轮。这时候加maxit只能是备选方案更值得做的是换松弛法或预处理因为线性收敛就是这么慢硬等不是好策略。5. 迭代求解避坑5 个让我翻过车的典型场景5.1 现象迭代十几步就开始输出 NaN第一次跑雅可比迭代几下就满屏 NaN这是最经典的翻车现场。最常见的直接原因是系数矩阵对角线上存在零或接近零的元素迭代公式里x[i] s / A[i, i]这一步直接除出无穷大或非法值。更深层的原因往往是矩阵构造代码里索引错位比如有限差分里把对角元素写成了 0。解决办法分两步。第一步矩阵构造完后立刻打印np.diag(A)用眼睛确认对角线没有异常值第二步在迭代循环开头加一个保护判断检查对角元素绝对值是否小于某个下限比如1e-12。这看起来像多余的防御代码但能省掉大量查错时间。5.2 现象矩阵满足对角占优条件结果仍然收敛极慢有个同学拿一个三千阶的有限差分矩阵跑高斯-赛德尔残差每轮只下降一点点迭代两万轮还没到容差第一反应是代码写错了。检查发现矩阵确实满足弱对角占优且不可约理论上必然收敛问题出在谱半径和 1 的距离太近。这种情况下代码没有 bug纯属方法本身不适合当前矩阵。解决路径是换加速手段先算一下谱半径如果大于 0.999直接用 SOR 或更现代的 Krylov 子空间方法硬跑高斯-赛德尔属于把时间花在没有希望的路上。另一个值得做的改变是增大容差到物理问题能接受的尺度有时候 1e-6 是过度要求1e-3 已经满足工程需要。5.3 现象认为高斯-赛德尔一定更快却遇到反例“高斯-赛德尔一定比雅可比快”这个论断在不少书里被当作默认结论但实战里它会坑人。确实存在一些非对称矩阵雅可比收敛而高斯-赛德尔发散的反例理论上也早有证明两种方法的收敛性之间没有绝对的蕴含关系。遇到这种矩阵拿默认感性选择的高斯-赛德尔跑结果就是发散或震荡。解决方法是把“二选一”改成“先验后选”。对测试矩阵分别求两个迭代矩阵的谱半径哪个小于 1 且离 1 更远就用哪个。做完这一步再跑正式规模的问题就不会被“理所当然”的判断带沟里。幂法估算谱半径的代码在第 2.3 节已经给出直接复用即可。5.4 现象迭代提前退出得到的解却明显不对有人把收敛判据写成||x_{k1} - x_k|| tol结果解和真实值差了十万八千里程序却报“已收敛”。这正是前面强调过的坑相邻两步变化量小不等于残差小。矩阵条件数很差的时候解空间是狭长的迭代点在一个方向上几乎不动看起来“稳定”了但方程还没被满足。解决方式只有一个把所有停机条件替换成残差范数||b - A x_k||或相对残差。如果出于性能原因不想每轮都算一次矩阵乘向量至少也要每隔几轮算一次绝不能只看解的变化量。5.5 现象与直接法结果对不上半天找不到原因迭代收敛得正常残差也达标了把解和np.linalg.solve的结果一对发现若干分量完全不同。这种问题的根源往往不是迭代法本身而是手写代码里混用了新旧值。雅可比实现里错误地用到了本轮已更新的分量等于把雅可比写成了某种违规的高斯-赛德尔变形高斯-赛德尔实现里忘记保留旧值向量后面的分量的计算用了被覆盖的数据。排查手段是造一个 3 阶或 4 阶小矩阵把中间每一轮迭代的x打印出来逐行对照公式手算。小规模下能直接看到哪个分量被提前更新或错误保留比对着大矩阵挠头快得多。这也是我最常用的回归方式改完代码先跑小例子确认无误再放到大算例上。6. 从雅可比和高斯-赛德尔到 SOR一个值得掌握的提速技巧6.1 高斯-赛德尔的松弛改造高斯-赛德尔已经比雅可比快但还可以再推一把。SOR连续超松弛是在高斯-赛德尔基础上加一个松弛因子 \omega当 \omega1 时退回原方法\omega1 时相当于放大每轮修正量。实践里 \omega 在 1.2 到 1.8 之间常常能显著减少迭代轮数但 \omega 接近 2 时容易震荡发散需要试验着调。def sor(A, b, omega, x0None, tol1e-8, maxit20000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n len(b) x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) for k in range(1, maxit 1): for i in range(n): s b[i] - A[i, :i] x[:i] - A[i, i1:] x[i1:] x[i] (1 - omega) * x[i] omega * s / A[i, i] r_norm np.linalg.norm(b - A x, ordnp.inf) if r_norm tol: return x, k, r_norm return x, maxit, np.linalg.norm(b - A x, ordnp.inf)内层表达式里(1 - omega) * x[i]保留旧分量的一部分omega * s / A[i,i]是高斯-赛德尔的完整修正量两者按权重组合。\omega 扫描从一个简单方式开始先跑 \omega1 记录轮数再分别试 1.2、1.4、1.6、1.8看哪个轮数最少。一维泊松问题里 \omega 取到 1.8 附近时轮数可能比高斯-赛德尔少一个量级效果立竿见影。6.2 把直接法留在工具箱里作为验证后路接手任何迭代求解需求我的习惯都是先保留一个小规模的直接法求解路径用np.linalg.solve求参考解再和迭代结果比对确认一致后才把直接法从正式流程里拿掉。这样做的好处是迭代代码后续改动时能立即发现回归问题而不是等到跑大算例才发现结果飘了。我自己的经验是迭代法从来没有“配好参数就一劳永逸”这种好事矩阵一变、边界条件一变、网格拓扑一变收敛行为就可能重新洗牌。与其反复纠结单一方法不如把雅可比、高斯-赛德尔、SOR 三个函数放在同一个工具箱里需要时先跑小规模对比再决定用哪个。这套流程虽然朴素但能挡掉大部分低级的数值翻车希望帮到你。本文还有配套的精品资源点击获取
简介面向数值计算课程与科研工程应用这份压缩包提供了雅可比迭代法和高斯-赛德尔迭代法的完整MATLAB实现帮助学习者快速解决大型线性方程组的近似求解问题。包内共5个文件全部为m格式源程序涵盖两种算法的函数定义与调用主程序便于直接运行、修改和对比迭代效果整个资源包仅3KB结构精简适合实验和作业场景。已有1967人学习这份资源既是课程设计的直接参考也能辅助理解两种迭代法在矩阵分裂、逐分量更新以及收敛速度上的差异。通过阅读和运行源码可掌握从系数矩阵分解、迭代公式编写到残差收敛判别的完整流程同时了解雅可比迭代适合并行计算、高斯-赛德尔迭代通常更快且更稳定的实际选型依据为后续处理大规模稀疏方程组打下扎实基础。1. 为什么线性方程组的迭代求解从雅可比和高斯-赛德尔开始求解稀疏线性方程组时迭代法比直接法更省内存这一点在真实工程里经常让新手意外几百阶矩阵直接消元还能撑住一旦到万阶以上一次直接法的内存和时间代价就会失控。雅可比迭代法与高斯-赛德尔迭代法正是迭代法里最基本、也最常被选作起点的两个方法。如果一个压缩包命名成这种里面通常跑的就是这两种方法的最小可运行实现。这篇笔记围绕同一个目标展开搞懂两种格式的差异在哪里用什么参数能让收敛快一点以及实际跑出问题时该去查哪里。2. 从矩阵分裂到谱半径雅可比与高斯-赛德尔为什么一个慢一个快2.1 雅可比与高斯-赛德尔的矩阵分裂与更新顺序两种方法都建立在同一个矩阵分裂上。把线性方程组 Ax b 的系数矩阵拆成 A D - L - U其中 D 是对角矩阵L 是严格下三角矩阵U 是严格上三角矩阵。这个拆分本身不改变方程的解改变的是怎么从当前近似解推进到下一步近似解。雅可比迭代法的格式是 x_new D^{-1}(b (LU)x_old)也就是下一轮迭代时所有分量都只用上一轮的值计算。这个特性让雅可比在并行环境里天然友好每个分量的计算互不依赖。高斯-赛德尔迭代法的格式是 x_new (D-L)^{-1}(b Ux_old)写成分量形式更直观求第 i 个分量时前 i-1 个分量已经更新过直接用新值后面的分量还用旧值。这种“边算边用”的串行更新让每一轮推进的信息量更大。两者的选择逻辑也因此不同追求可并行、实现简单、对初值鲁棒选雅可比追求串行环境下更快的收敛选高斯-赛德尔。一个容易误解的点是高斯-赛德尔在每个分量上比雅可比多做了更新但单轮的浮点运算代价没有质的差别差别主要体现在达到相同精度所需的轮数上。对比项雅可比迭代法高斯-赛德尔迭代法更新使用的值全部用上一轮的值已更新的分量用新值未更新的用旧值分量间耦合无天然可并行有串行依赖存储要求需要保存新旧两轮向量一个向量原地更新即可对称正定问题收敛轮数约为高斯-赛德尔的 2 倍更少2.2 收敛性判据谱半径、对角占优与收敛速度的直觉迭代法看起来是无限逼近实际是否收敛取决于迭代矩阵的谱半径。把两种方法统一写成 x_{k1} M x_k c 的形式收敛的充要条件是迭代矩阵 M 的谱半径 \rho(M) 1。雅可比对应的迭代矩阵是 D^{-1}(LU)高斯-赛德尔对应的是 (D-L)^{-1}U。工程里不会真的去手算谱半径但有两条硬性结论可以直接用如果 A 严格对角占优即每个对角元素的绝对值严格大于该行非对角元素绝对值之和那么雅可比和高斯-赛德尔都收敛如果 A 不可约且弱对角占优结论同样成立。注意“不可约”这个条件一维泊松离散矩阵恰好就是这种典型内层行对角占优达到临界但没有严格大于它依然能保证收敛。收敛速度的直觉来自谱半径和 1 的距离。若 \rho(M) 是 0.9每轮误差大约缩小到上一轮的 90%若 \rho(M) 是 0.999每轮后误差几乎不动。高斯-赛德尔看起来“更快”本质上是它的迭代矩阵谱半径在对称正定问题下通常约等于雅可比谱半径的平方0.999^2 比 0.999 离 1 更远所以收敛轮数接近减半。这个比例关系经常被误解成普适结论实际上只在一类问题里成立后文避坑章节会专门说。2.3 一个工程习惯先估算迭代矩阵的谱半径再选方法我一般会在动手写大矩阵之前先用一个小规模矩阵做谱半径估算。用幂法就能得到近似值不需要求特征值分解。给定系数矩阵 A先构造两个迭代矩阵再用随机向量反复乘观察范数增长或衰减的幅度。import numpy as np def build_jacobi_matrix(A): D np.diag(np.diag(A)) L -np.tril(A, -1) U -np.triu(A, 1) return np.linalg.inv(D) (L U) def spectral_radius_estimate(M, seed7, max_iter5000): rng np.random.default_rng(seed) v rng.normal(sizeM.shape[0]) # 幂法反复迭代每次归一化模长的变化趋势指向主特征值 for _ in range(max_iter): Mv M v norm np.linalg.norm(Mv, ordnp.inf) if norm 1e-15: return 0.0 v Mv / norm return float(np.linalg.norm(M v, ordnp.inf))这段代码里build_jacobi_matrix严格按 A D - L - U 的符号组装L和U都带负号。幂法每次把向量乘到迭代矩阵上并归一化若模长持续增长说明主特征值大于 1若衰减说明小于 1。注意幂法依赖于初始向量在主导特征向量方向上有分量所以我会换两个不同随机种子多跑一次结果一致才采信。谱半径估计的目的不是精确到小数点后几位而是判断这个矩阵适不适合做迭代求解以及大概需要多少轮。3. 把迭代写成能跑的代码简单但好用的雅可比与高斯-赛德尔实现3.1 用 numpy 实现两个迭代器收敛判定放在残差上而不是解上写迭代求解器的第一原则是收敛判据必须建立在残差上不能建立在相邻两步解的差值上。相邻两步解变化很小可能是因为矩阵条件数很差、解空间很平也可能因为真正收敛了这两者从“解差”里分不清。所以下面的实现里回归停机的判断用的是 ||b - Ax||。import numpy as np def jacobi(A, b, x0None, tol1e-8, maxit20000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n len(b) x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) diag np.diag(A) off A - np.diagflat(diag) # L U去掉对角后的剩余部分 for k in range(1, maxit 1): x_new (b - off x) / diag # 所有分量同时用旧值更新 r_norm np.linalg.norm(b - A x_new, ordnp.inf) if r_norm tol: return x_new, k, r_norm x x_new return x, maxit, np.linalg.norm(b - A x, ordnp.inf)diag是主对角元素组成的向量off是去掉对角后的矩阵off x相当于 (LU)x。(b - off x) / diag直接把向量除到对角元素上在 numpy 里是逐元素除法。每次更新后立刻用新向量算残差用无穷范数比较。注意这里如果对角元素有接近 0 的值会产生除零警告后文避坑部分会讲如何提前发现。高斯-赛德尔的实现比雅可比多一个内层循环因为每个分量更新依赖前面刚算出来的新值def gauss_seidel(A, b, x0None, tol1e-8, maxit20000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n len(b) x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) for k in range(1, maxit 1): x_old x.copy() # 必须保留旧值给未更新的分量使用 for i in range(n): s b[i] - A[i, :i] x[:i] - A[i, i1:] x_old[i1:] x[i] s / A[i, i] r_norm np.linalg.norm(b - A x, ordnp.inf) if r_norm tol: return x, k, r_norm return x, maxit, np.linalg.norm(b - A x, ordnp.inf)内层循环中A[i, :i] x[:i]用的是已更新的新分量A[i, i1:] x_old[i1:]用的是上一轮旧分量。x_old x.copy()这一步不能省因为 x 更新过程中会覆盖掉旧值后面的行还要拿旧值参与计算。几乎所有从雅可比改写成高斯-赛德尔翻车的人都栽在这一行上。这里还藏着两个实现细节。第一个是收敛判定放在每次完整扫描之后因为一轮内部分量更新并不代表整个方程已经近似满足第二个是maxit到期后返回当前结果但不抛异常我习惯把残差一并返回方便调用方判断是“收敛”还是“用满次数”。3.2 一维泊松测试算例同一套矩阵、右端项对比迭代轮数光看代码不跑算例很难建立对两种方法差异的直觉。一维泊松方程离散出来的三对角矩阵是最适合做对比的测试台结构简单、性质清晰、收敛情况可控。def poisson_1d(n): h2 (n 1) ** 2 A np.zeros((n, n)) for i in range(n): A[i, i] 2.0 * h2 if i 0: A[i, i - 1] -h2 if i n - 1: A[i, i 1] -h2 return A n 100 A poisson_1d(n) b np.ones(n) x_j, it_j, r_j jacobi(A, b, tol1e-6, maxit100000) x_g, it_g, r_g gauss_seidel(A, b, tol1e-6, maxit100000) print(Jacobi:, it_j, r_j) print(Gauss-Seidel:, it_g, r_g)在这个 n100 的算例里雅可比达到 tol1e-6 大约需要三万轮高斯-赛德尔大约只需要一万五千轮。这里的三对角矩阵对角元素是 2h²相邻元素是 -h²内部行对角占优只是临界成立不是严格占优它能保证收敛靠的是“不可约且弱对角占优”这条定理。把 n 增大到 400谱半径会更接近 1两种方法的轮数会按 n² 的规模往上翻这是线性收敛的正常表现不是程序挂死。跑完这个算例再看代码你会理解一个事实迭代轮数不是唯一标准。每轮的计算量同样重要雅可比每轮可以做向量化矩阵乘法高斯-赛德尔的内层循环在 Python 里是串行推导n 很大时单轮耗时可能明显更久。实际工程里谁更快永远是“轮数”和“每轮耗时”的乘积不能只看单方面。3.3 数值试验的 3 条读结果规则跑出数字之后不要急着认定结论先用三条规则过一遍。第一条用直接法兜底。对这个小规模算例直接调np.linalg.solve(A, b)拿到参考解和迭代结果比一下最大绝对误差。如果参考解与迭代解对不上多数情况不是迭代法不收敛而是手写代码里某个索引或符号写反了。第二条看残差下降曲线是否单调。打印每轮的||b - Ax_k||_∞如果稳定下降并进入线性区间说明迭代正常如果残差波动甚至上升基本可以断定矩阵性质不支持当前方法不要盲目加迭代次数。第三条检查初值和右端项的量级。很多“不收敛”其实是容差设得太苛刻比如矩阵和 b 的量级在 1e10却设了 1e-12 的绝对容差导致永远无法满足。先标准化再做迭代是减少这类误判的有效手段。4. 三个必调参数与收敛判定初值、容差、最大迭代次数4.1 迭代初值零向量不总是最差但也别只会给零从零向量出发是默认选择对很多对称正定问题来说这个选择并不会让收敛轮数显著变差。实际工程里更值得利用的是“上一个状态的解”。比如瞬态问题里上一时刻的解和当前时刻的解往往很接近把上一时刻的解作为初值迭代轮数可能缩减一个数量级这种收益比调任何容差都来得直接。另一个常见做法是用物理常识约束初值。温度场问题里初值取环境温度而不是全零结构位移问题里初值取上一载荷步结果都会让前几十轮迭代少走弯路。我一般会在函数入口允许传x0并默认置零这样调用方可以根据问题类型决定是否给一个更聪明的起点。如果初值给得太离谱比如物理上不可能出现的负温度或超大位移反而会让迭代前段出现大幅摆动。所以初值的关键不是“越接近就好”而是“不能显著违背物理约束”。4.2 残差判据与容差用 ||b-Ax_k||_∞ 替代相邻解差收敛判据的选择直接决定“这个解能不能用”。相邻两步解的差已经分析过不可靠更隐蔽的问题是绝对容差和相对容差。假设你的方程组是 x_1 x_2 1、x_1 - x_2 2残差降低到 1e-8 很容易但同样的容差放到边界元法生成的矩阵上可能跑一万轮也达不到。解决办法是使用相对残差判据||b - Ax_k|| / ||b||把量级差异消除掉。实际判断顺序是先算初始残差r0 ||b - Ax0||每轮算rk ||b - Ax_k||当rk / r0小于阈值时认定收敛。这样矩阵和右端项放大或缩小一百倍判定结果不受影响。容差的具体数值单精度场景通常取 1e-5 或 1e-6双精度场景取 1e-8 到 1e-10 已经足够。再往下压到 1e-12很多时候只是把迭代轮数白白拉长解的精度并不会跟着提高因为浮点舍入本身已经限制了可实现精度。4.3 最大迭代次数的设置逻辑线性收敛怎么预估轮数线性迭代的误差大致按几何级数衰减每轮误差约为上一轮的 \rho 倍。如果谱半径 \rho0.999每轮大约只缩小千分之一要让误差从 1 掉到 1e-6需要大约 1.4 万轮。这个估算公式是所需轮数约为ln(容差 / 初始误差) / ln(rho)用自然对数就算出来了。所以设置最大迭代次数之前先估算谱半径然后按这个公式给上限加一个安全系数。比如估算需要一万五千轮maxit就设成六万留足余量。如果残差曲线已经进入线性下降区间却因为maxit设小了被迫退出这是最可惜的浪费。还有一种情况是残差曲线看起来在下降但速度极慢比如从 1e-2 到 1e-4 花了一万轮到 1e-6 可能还要两万轮。这时候加maxit只能是备选方案更值得做的是换松弛法或预处理因为线性收敛就是这么慢硬等不是好策略。5. 迭代求解避坑5 个让我翻过车的典型场景5.1 现象迭代十几步就开始输出 NaN第一次跑雅可比迭代几下就满屏 NaN这是最经典的翻车现场。最常见的直接原因是系数矩阵对角线上存在零或接近零的元素迭代公式里x[i] s / A[i, i]这一步直接除出无穷大或非法值。更深层的原因往往是矩阵构造代码里索引错位比如有限差分里把对角元素写成了 0。解决办法分两步。第一步矩阵构造完后立刻打印np.diag(A)用眼睛确认对角线没有异常值第二步在迭代循环开头加一个保护判断检查对角元素绝对值是否小于某个下限比如1e-12。这看起来像多余的防御代码但能省掉大量查错时间。5.2 现象矩阵满足对角占优条件结果仍然收敛极慢有个同学拿一个三千阶的有限差分矩阵跑高斯-赛德尔残差每轮只下降一点点迭代两万轮还没到容差第一反应是代码写错了。检查发现矩阵确实满足弱对角占优且不可约理论上必然收敛问题出在谱半径和 1 的距离太近。这种情况下代码没有 bug纯属方法本身不适合当前矩阵。解决路径是换加速手段先算一下谱半径如果大于 0.999直接用 SOR 或更现代的 Krylov 子空间方法硬跑高斯-赛德尔属于把时间花在没有希望的路上。另一个值得做的改变是增大容差到物理问题能接受的尺度有时候 1e-6 是过度要求1e-3 已经满足工程需要。5.3 现象认为高斯-赛德尔一定更快却遇到反例“高斯-赛德尔一定比雅可比快”这个论断在不少书里被当作默认结论但实战里它会坑人。确实存在一些非对称矩阵雅可比收敛而高斯-赛德尔发散的反例理论上也早有证明两种方法的收敛性之间没有绝对的蕴含关系。遇到这种矩阵拿默认感性选择的高斯-赛德尔跑结果就是发散或震荡。解决方法是把“二选一”改成“先验后选”。对测试矩阵分别求两个迭代矩阵的谱半径哪个小于 1 且离 1 更远就用哪个。做完这一步再跑正式规模的问题就不会被“理所当然”的判断带沟里。幂法估算谱半径的代码在第 2.3 节已经给出直接复用即可。5.4 现象迭代提前退出得到的解却明显不对有人把收敛判据写成||x_{k1} - x_k|| tol结果解和真实值差了十万八千里程序却报“已收敛”。这正是前面强调过的坑相邻两步变化量小不等于残差小。矩阵条件数很差的时候解空间是狭长的迭代点在一个方向上几乎不动看起来“稳定”了但方程还没被满足。解决方式只有一个把所有停机条件替换成残差范数||b - A x_k||或相对残差。如果出于性能原因不想每轮都算一次矩阵乘向量至少也要每隔几轮算一次绝不能只看解的变化量。5.5 现象与直接法结果对不上半天找不到原因迭代收敛得正常残差也达标了把解和np.linalg.solve的结果一对发现若干分量完全不同。这种问题的根源往往不是迭代法本身而是手写代码里混用了新旧值。雅可比实现里错误地用到了本轮已更新的分量等于把雅可比写成了某种违规的高斯-赛德尔变形高斯-赛德尔实现里忘记保留旧值向量后面的分量的计算用了被覆盖的数据。排查手段是造一个 3 阶或 4 阶小矩阵把中间每一轮迭代的x打印出来逐行对照公式手算。小规模下能直接看到哪个分量被提前更新或错误保留比对着大矩阵挠头快得多。这也是我最常用的回归方式改完代码先跑小例子确认无误再放到大算例上。6. 从雅可比和高斯-赛德尔到 SOR一个值得掌握的提速技巧6.1 高斯-赛德尔的松弛改造高斯-赛德尔已经比雅可比快但还可以再推一把。SOR连续超松弛是在高斯-赛德尔基础上加一个松弛因子 \omega当 \omega1 时退回原方法\omega1 时相当于放大每轮修正量。实践里 \omega 在 1.2 到 1.8 之间常常能显著减少迭代轮数但 \omega 接近 2 时容易震荡发散需要试验着调。def sor(A, b, omega, x0None, tol1e-8, maxit20000): A np.asarray(A, dtypefloat) b np.asarray(b, dtypefloat) n len(b) x np.zeros(n) if x0 is None else np.array(x0, dtypefloat) for k in range(1, maxit 1): for i in range(n): s b[i] - A[i, :i] x[:i] - A[i, i1:] x[i1:] x[i] (1 - omega) * x[i] omega * s / A[i, i] r_norm np.linalg.norm(b - A x, ordnp.inf) if r_norm tol: return x, k, r_norm return x, maxit, np.linalg.norm(b - A x, ordnp.inf)内层表达式里(1 - omega) * x[i]保留旧分量的一部分omega * s / A[i,i]是高斯-赛德尔的完整修正量两者按权重组合。\omega 扫描从一个简单方式开始先跑 \omega1 记录轮数再分别试 1.2、1.4、1.6、1.8看哪个轮数最少。一维泊松问题里 \omega 取到 1.8 附近时轮数可能比高斯-赛德尔少一个量级效果立竿见影。6.2 把直接法留在工具箱里作为验证后路接手任何迭代求解需求我的习惯都是先保留一个小规模的直接法求解路径用np.linalg.solve求参考解再和迭代结果比对确认一致后才把直接法从正式流程里拿掉。这样做的好处是迭代代码后续改动时能立即发现回归问题而不是等到跑大算例才发现结果飘了。我自己的经验是迭代法从来没有“配好参数就一劳永逸”这种好事矩阵一变、边界条件一变、网格拓扑一变收敛行为就可能重新洗牌。与其反复纠结单一方法不如把雅可比、高斯-赛德尔、SOR 三个函数放在同一个工具箱里需要时先跑小规模对比再决定用哪个。这套流程虽然朴素但能挡掉大部分低级的数值翻车希望帮到你。本文还有配套的精品资源点击获取
