简介：面向数值计算学习者的MATLAB实现包，聚焦雅可比迭代法与高斯-赛德尔迭代法求解Ax=b型线性方程组。两种方法以迭代逼近替代直接求逆，尤其适合大型稀疏矩阵场景，便于理解迭代格式、收敛条件与算法差异。压缩包体积仅3KB，共含5个m文件，主程序和迭代函数分离，结构清晰，可直接运行或修改。已有1967人学习下载，适合数值分析课程实验、期末复习或工程应用入门。通过对比两种算法的更新规则，可直观体会高斯-赛德尔利用最新分量带来的加速效果；也可自行调整矩阵维度、初始向量与迭代精度，验证系数矩阵条件数对收敛速度的影响。这份资源将理论公式转化为可执行代码，能帮助初学者快速掌握迭代法编程实现与调试思路，是数值计算实践环节中一份轻量而完整的参考资料。
1. 迭代法求解线性方程组：一个 MATLAB 资源包能帮你省下多少调试时间
做数值计算分析时，凡是遇到大型稀疏线性方程组 Ax=b，直接法（高斯消去）往往不是首选——矩阵一大，消元产生的填充元会把稀疏结构毁掉，内存与耗时双双失控。迭代法才是这类场景的常规解法，而 Jacobi 雅可比迭代法和 Gauss-Seidel 高斯-赛德尔迭代法又是迭代法里最经典、最常被写进实验报告的一对。这套 MATLAB 资源包就是围绕这两个方法整理的：完整的迭代函数、收敛性分析脚本、若干典型算例。适合两类人：一类是正在上数值计算分析课程、需要交实验报告的学生；另一类是实际仿真里遇到大规模稀疏矩阵、想快速搭起迭代求解器的从业者。先说一个反直觉结论：Gauss-Seidel 并不保证永远比 Jacobi 快，存在某些矩阵让 Jacobi 收敛而 Gauss-Seidel 发散。这类边界知识，恰恰是调试时最容易翻车的地方。
2. Jacobi 迭代法的 MATLAB 实现：拆解 D、L、U 与收敛判定
2.1 迭代公式的由来：把 A 拆成三块
Jacobi 迭代法的出发点是把系数矩阵 A 拆成三部分：对角部分 D、严格下三角部分 L、严格上三角部分 U，即 A = D - L - U。这里的 L、U 取了负号，是为了把迭代格式写成矩阵乘法的标准形式——迭代式是 x^(k+1) = D^(-1)(b + (L+U)x^(k))。注意这里的 L、U 带负号，整理一下就得到最常用的分量形式：
x_i^(k+1) = (b_i - Σ_{j≠i} a_ij * x_j^(k)) / a_ii
理解这个公式的关键在于：Jacobi 迭代在更新第 i 个分量时，用的是上一轮 k 的全部旧值。这意味着每一轮内各个分量的更新互不依赖，天然适合并行计算。在 GPU 普及之后，Jacobi 因为并行友好又被重新重视起来，很多大规模稀疏矩阵的预处理器里都能看到它的影子。收敛性方面，迭代矩阵 M_J = D^(-1)(L+U) 的谱半径 ρ(M_J) < 1 是收敛的充要条件；而严格对角占优（每行对角元素绝对值大于其余元素绝对值之和）是工程上最常用的充分条件，实际矩阵能套上对角占优，基本就可以放心迭代了。
2.2 实现 jacobi.m：主函数设计与参数说明
常见的做法是把迭代封装成一个独立函数，输入系数矩阵 A、右端项 b、初始向量 x0、容差 tol、最大迭代次数 maxit，输出解 x、实际迭代步数 iter 和收敛标志 flag。下面是这份资源包里 jacobi.m 的核心实现：
function [x, iter, flag] = jacobi(A, b, x0, tol, maxit) % JACOBI Jacobi 迭代法求解 Ax=b % A: 系数矩阵（方阵）; b: 右端项列向量; x0: 初始解向量 % tol: 相邻迭代差容差; maxit: 最大迭代次数 % flag: 0 表示收敛, 1 表示达到 maxit 未收敛 D = diag(diag(A)); % 提取对角部分 L = -tril(A, -1); % 严格下三角取负 U = -triu(A, 1); % 严格上三角取负 M = D \ (L + U); % 迭代矩阵 M_J f = D \ b; % 常数项 x = x0; for k = 1:maxit x_new = M * x + f; % 一步雅可比迭代 err = norm(x_new - x, inf); % 无穷范数衡量变化量 x = x_new; if err < tol iter = k; flag = 0; return; end end iter = maxit; flag = 1; end这段代码的逻辑要点有三个。第一，迭代矩阵 M 用 D \ (L+U) 而不是 inv(D)*(L+U)，左除在 MATLAB 里的数值稳定性更好，且对对角矩阵来说左除就是逐元素除法，效率极高。第二，终止条件用的是相邻两轮解的无穷范数差，norm(x_new - x, inf) 表示所有分量里变化最大的那个绝对值，这个指标比二范数更严格、更容易暴露某个分量还在震荡的问题。第三，tol 和 maxit 是调试时最常动的两个参数：tol 设 1e-6 适合大部分教学算例，工程仿真里建议收紧到 1e-8；maxit 一般给 1000 起步，如果谱半径接近 1（比如 0.98 以上），1000 次可能还不够，得放大到 5000 甚至 10000。初始向量 x0 没把握就传 zeros(n,1)，迭代法对 x0 的要求不苛刻。
2.3 一个可复现算例：验证收敛与迭代轨迹
为了确认函数可用，资源包里附了一个 3×3 严格对角占优算例，我建议你拿到压缩包后先跑这个，再改自己的矩阵：
A = [4 -1 0; -1 4 -1; 0 -1 3]; b = [1; 2; 3]; x0 = zeros(3,1); tol = 1e-8; maxit = 1000; [x, iter, flag] = jacobi(A, b, x0, tol, maxit); fprintf('迭代次数: %d, 收敛标志: %d\n', iter, flag); disp('解 x ='); disp(x); % 验证残差 r = norm(b - A*x, 2); fprintf('残差 ||b - Ax||_2 = %.6e\n', r);这个矩阵每行对角元素（4、4、3）都严格大于该行其余元素绝对值之和（1、2、1），满足对角占优，理论上必收敛。跑完能看到 iter 大概在 20 次以内，残差在 1e-8 附近。如果残差下不去，先看是 flag 等于 1（达到最大步数没收敛）还是残差还在缓慢下降（谱半径离 1 太近）。典型的错误是以为加大 maxit 就能解决一切问题——如果谱半径大于 1，迭代矩阵会把误差分量不断放大，跑多少次都没用，这时候该换方法或重排方程组。
3. Gauss-Seidel 迭代法：一个语句的差距，收敛速度翻倍
3.1 即时更新思想：新值立刻投入使用
Gauss-Seidel 迭代法与 Jacobi 的唯一区别在于：计算第 i 个分量时，前 i-1 个分量已经用本轮的新值更新完毕，直接拿来参与求和；只有后面的分量仍用上一轮旧值。即迭代格式为 x^(k+1) = (D-L)^(-1)(b + Ux^(k))，分量形式是：
x_i^(k+1) = (b_i - Σ_{j<i} a_ij * x_j^(k+1) - Σ_{j>i} a_ij * x_j^(k)) / a_ii
这个“即时更新”让同一轮迭代的信息利用率更高。收敛时，Gauss-Seidel 的迭代步数通常是 Jacobi 的一半到三分之一，尤其是在对角占优很好的矩阵上，差距非常明显。存储方面，Gauss-Seidel 每轮只需一个解向量 x 就够了（新值直接覆盖旧值），而 Jacobi 在并行实现时需要 x_new 和 x 两份。表格里能看得更清楚：
| 对比项 | Jacobi | Gauss-Seidel |
|---|---|---|
| 收敛速度 | 较慢，步数约为 GS 的 1.2~3 倍 | 较快，多数场景优于 Jacobi |
| 存储需求 | 需要 x_new 与 x 两份向量 | 一份向量即可原地更新 |
| 并行性 | 各分量独立更新，天然并行 | 分量串行依赖，难并行 |
| 适用场景 | GPU 并行、教学演示 | 串行程序、常规工程计算 |
3.2 实现 gauss_seidel.m：与 Jacobi 的逐行差异
Gauss-Seidel 的矩阵形式实现只差了一行——迭代矩阵换成 (D-L) \ U，常数项换成 (D-L) \ b：
function [x, iter, flag] = gauss_seidel(A, b, x0, tol, maxit) % GAUSS_SEIDEL Gauss-Seidel 迭代法求解 Ax=b % 参数含义同 jacobi.m，额外要求 A 的对角元素均非零 D = diag(diag(A)); L = -tril(A, -1); U = -triu(A, 1); M = (D - L) \ U; % 迭代矩阵 M_GS f = (D - L) \ b; x = x0; for k = 1:maxit x_new = M * x + f; % 一步 GS 迭代 err = norm(x_new - x, inf); x = x_new; if err < tol iter = k; flag = 0; return; end end iter = maxit; flag = 1; end一个值得注意的细节：MATLAB 里 (D-L) \ U 用的是三角矩阵的左除算法，本质上做的是前代法，不会真的去求逆，所以即使 n 很大，这一步也不会出现显式逆矩阵的内存爆炸。如果你手头是稀疏矩阵，建议先把 A 存成 sparse 格式，再传给函数，三角左除会自动走稀疏路径，速度差上百倍都有可能。这份资源的代码里没有特意处理稀疏格式，但函数本身兼容，你在调用前加一行 A = sparse(A) 即可。
3.3 谁更快不一定：一组发散的边界案例
工程上普遍的经验是“Gauss-Seidel 总比 Jacobi 好”，我在实际项目里也确实 90% 的情况都验证了这个经验。但数值分析课上不会讲的反例是这样的——存在矩阵使 Jacobi 收敛而 Gauss-Seidel 发散。资源包里专门附了一组这样的矩阵构造脚本：
% 反例：Jacobi 收敛、Gauss-Seidel 发散 A = [1 2 -2; 1 1 1; 2 2 1]; b = [1; 1; 3]; x_jac = jacobi(A, b, zeros(3,1), 1e-6, 1000); x_gs = gauss_seidel(A, b, zeros(3,1), 1e-6, 1000); disp('Jacobi flag:'); disp(x_jac(3)); disp('GS flag:'); disp(x_gs(3)); % 分别计算两个迭代矩阵的谱半径 D = diag(diag(A)); L = -tril(A, -1); U = -triu(A, 1); rho_J = max(abs(eig(D \ (L + U)))); rho_GS = max(abs(eig((D - L) \ U))); fprintf('rho_J = %.4f, rho_GS = %.4f\n', rho_J, rho_GS);这个例子中 rho_J 小于 1 而 rho_GS 大于 1，直接导致两个方法一个收敛一个发散。这类反例的价值不在于让你记住这个具体矩阵，而在于提醒你：谱半径判断永远比“方法名听起来更高级”可靠。拿到新矩阵，先花 10 秒算两个迭代矩阵的谱半径，比盲目选方法省心得多。
4. 收敛性判定与参数选取：谱半径、对角占优与停机准则
4.1 谱半径计算与 eig 的精度陷阱
谱半径是判断迭代法收敛与否的最硬指标，定义是迭代矩阵特征值模的最大值。MATLAB 里直接算：
function rho = spectral_radius_of_jacobi(A) % 计算 Jacobi 迭代矩阵的谱半径 D = diag(diag(A)); L = -tril(A, -1); U = -triu(A, 1); M = D \ (L + U); rho = max(abs(eig(M))); end两个使用要点：第一，eig 会算出全部分特征值，对几百阶的小矩阵没问题；但上万阶的稀疏矩阵用 eig 会非常慢且吃内存，要改用 eigs(M, 1, 'largestabs') 只求模最大的那个特征值。第二，浮点误差会让特征值在 1 附近抖出假象，有时算出来 1.0000001，其实理论值恰好在 1 附近，所以判断线别设成 1，我一般用 0.999 作为“危险线”：大于 0.999 的，迭代会慢到让人失去耐心，建议直接换 SOR 或预处理。谱半径 0.5 和 0.95 的差距是数量级的——前者十几步收敛，后者可能几百上千步，这一步计算成本极低，却提前告诉你结果，别跳过。
4.2 对角占优判定与方程组重排
严格对角占优是收敛的充分非必要条件，工程里它最有用，因为判定起来一目了然。判定代码很简单：
function is_dd = is_strictly_diagonally_dominant(A) % 判断矩阵 A 是否严格对角占优 n = size(A, 1); is_dd = true; for i = 1:n sum_off = sum(abs(A(i, :))) - abs(A(i, i)); if abs(A(i, i)) <= sum_off is_dd = false; break; end end end如果你的矩阵不满足对角占优，先别急着换算法。常见做法是对方程组做行重排——目标是让每行绝对值最大的元素尽量落到对角位置。最简单的是贪心法：把每一行绝对值最大的元素所在列的索引提取出来，如果有冲突就按最大值降序处理。重排后的矩阵可能就满足对角占优了，Jacobi 和 GS 就能直接跑。注意重排等价于对原方程组左乘一个置换矩阵，解不变；但 Gauss-Seidel 的收敛性对行顺序敏感，重排后一定要重新算谱半径验证一遍，否则可能更糟。Jacobi 对行顺序不敏感，重排基本不影响收敛性。
4.3 停机准则：残差优先于相邻解差
默认的停机准则是相邻迭代步解向量的无穷范数差小于 tol，这个准则简单但有个隐蔽缺陷：收敛极慢时，相邻两步的变化可能已经小于 1e-8，但解离真实值还差着十万八千里。更可靠的做法是同时检查残差范数。资源包里的改进版循环长这样：
for k = 1:maxit x_new = M * x + f; x = x_new; res = norm(b - A * x, 2); % 残差的二范数 if res < tol_res % tol_res 通常比 tol 稍宽松 iter = k; flag = 0; return; end end参数选取上，我一般会给一套组合：tol_res = 1e-6 用于教学算例，1e-8 用于需要高精度结果的仿真；maxit 视矩阵规模而定，小矩阵给 1000，大规模稀疏矩阵给 n*100（n 是方程阶数）通常足够。注意 tol 取得太小会显著拉长迭代时间，谱半径接近 1 时尤其如此——已经知道 ρ=0.998 的情况下，硬要残差到 1e-10 是不现实的，需要考虑换方法。正确姿势是三步走：先求谱半径判断可行性，再定合理的残差目标，最后用残差下降曲线观察收敛趋势。
5. 常见问题与排查：五个让结果翻车的边界陷阱
5.1 迭代次数爆表但残差依然很大
现象：跑完 maxit 次迭代，flag 返回 1，打印出来的残差还在 1e-2 量级甚至更大。
原因：最常见的是迭代矩阵谱半径大于 1，误差分量被不断放大，怎么迭代都不会收敛；也可能是谱半径在 0.999 附近，收敛极慢，maxit 设的 1000 根本不够用。
解决：先用rho = max(abs(eig(D\(L+U))))算出谱半径。rho > 1 就别在 Jacobi 上浪费时间，直接换 Gauss-Seidel 或 SOR；rho 在 0.99 以上但小于 1，把 maxit 放大到 5000 并观察残差是否在对数尺度上线性下降——如果下降曲线是平的，说明收敛条件不成立，重排方程组或换方法才是正路。
5.2 Jacobi 收敛而 Gauss-Seidel 发散
现象：同一个矩阵，跑 Jacobi 十几步收敛，跑 Gauss-Seidel 却步数耗尽还发散，输出结果直接是 NaN。
原因：两个迭代矩阵的谱半径没有固定的序关系，存在 ρ(M_J) < 1 < ρ(M_GS) 的矩阵。单纯基于“GS 比 Jacobi 高级”做选择，就会踩进这个坑。
解决：换方法前各算一行谱半径，这比跑完整迭代快得多且不会翻车。如果 GS 发散而 Jacobi 收敛，用 Jacobi 就是最合理的选择，别想着“一定得用更高级的方法”。数值计算的第一原则是方法适配问题，而不是方法看起来多高深。
5.3 初始向量全部用零，却报“除零”错误
现象：A 存在零对角元，迭代公式里分母 a_ii 为 0，直接报错或得到 Inf。
原因：迭代法的分量更新公式要求所有对角元素非零，但现实工程矩阵某些对角元可能就是 0，比如网格剖分后某个节点恰好没有自联立项。这不是算法问题，是矩阵本身不适合直接迭代。
解决：先检查diag(A)是否全非零。有零就做行重排或列重排把非零元素移到对角；实在不行，把方程组两边同乘一个置换阵或做简单的行变换消去零对角。这个预处理步骤在工程上比调迭代参数更值得优先考虑。
5.4 大矩阵迭代一次要好几秒
现象：矩阵规模到几千阶以上，一次迭代循环耗时明显变长，总迭代时间完全不能接受。
原因：一是在用全稠密矩阵存储，迭代里的矩阵乘法 O(n^2) 的代价被放大了；二是循环体内部逐元素更新没用向量化，MATLAB 的 for 循环在这种场景性能很差。
解决：先把 A 转成 sparse 存储，再确保迭代函数里尽量用矩阵运算代替元素级循环。Jacobi 的一轮迭代可以写成x_new = M*x + f，M 本身也可以保持稀疏格式，这样单步代价从 O(n^2) 降到非零元个数级别。矩阵较大的情况下，别用 eig 求谱半径，改用 eigs 只求最大模特征值。
5.5 残差小但解明显不对
现象：残差norm(b - A*x, 2)已经压到 1e-8，但把解代回去肉眼可见不合理，比如解向量里有巨大的负值或者振荡值。
原因：矩阵病态，条件数很大。残差小只说明 Ax 逼近 b 做得好，不说明 x 逼近真实解；病态矩阵下，b 的微小扰动或中间过程的舍入误差会被条件数成倍放大。
解决：算一下cond(A)，超过 1e10 就要警惕。小幅扰动右端项 b（在千分之一量级），再迭代一次看解的变化幅度——变化巨大就确认病态。此时单纯调迭代法的参数没有意义，考虑用带预处理的 Krylov 子空间方法，或改用直接法配合高精度格式。迭代法不是万能的，病态矩阵是它的公认软肋。
6. 进阶：SOR 超松弛迭代与一个从实战中养成的习惯
6.1 从 Gauss-Seidel 到 SOR：松弛因子 ω 怎么选
SOR（逐次超松弛迭代）的思想很直接：把 Gauss-Seidel 一步的结果当作“候选方向”，用松弛因子 ω 做加权更新。公式是：
x^(k+1) = (1-ω) * x^(k) + ω * (GS一步的结果)
ω 在 0 到 2 之间，ω=1 时退化为普通 GS，ω>1 是超松弛，ω<1 是低松弛。超松弛的效果是让更新步子迈得更大，从而加速收敛——但步子太大会过冲震荡，步子太小又跑不动。工程上不存在解析的最优 ω，最常见做法是写个小脚本扫描：
omega_range = 0.5:0.1:1.9; best_omega = 1; best_iter = inf; for omega = omega_range [~, iter] = sor(A, b, x0, omega, 1e-8, 2000); if iter < best_iter best_iter = iter; best_omega = omega; end end fprintf('最优 omega = %.1f, 迭代 %d 步\n', best_omega, best_iter);sor.m 的实现并不复杂，核心就是把 GS 结果做个加权平均：
function [x, iter] = sor(A, b, x0, omega, tol, maxit) D = diag(diag(A)); L = -tril(A, -1); U = -triu(A, 1); x = x0; for k = 1:maxit x_gs = (D - L) \ (b - U*x); % 先算 GS 一步 x = x + omega * (x_gs - x); % 松弛加权 if norm(b - A*x, 2) < tol iter = k; return; end end iter = maxit; end扫描 0.5 到 1.9 的 ω，通常能看到迭代步数呈抛物线变化，最优 ω 在 1.3 到 1.8 之间不等。SOR 对 ω 敏感，但扫描成本低，值得做。
6.2 一个习惯：拿到新矩阵先走“三连”
有回做某个热传导模拟项目，40×40 网格离散出的矩阵不满足严格对角占优，我图省事直接把数据灌进 Jacobi，跑了 5000 步 flag 还是 1，浪费了一整晚查来查去，最后算了一眼谱半径——1.03，根本不收敛。从那以后，我每次拿到新方程组都强制走一遍三连流程：第一，检查对角占优；第二，求迭代矩阵谱半径；第三，跑 50 步提前看看残差下降曲线是不是线性下滑。整个过程不到 10 分钟，但省下来的调试时间是论小时算的。这份资源里已经包含了这三个步骤的脚本，建议你拿到后先按顺序跑一遍所有算例，再把自己的矩阵代进去，感受收敛曲线长什么样，以后遇到类似问题心里就有底了。希望帮到你。
本文还有配套的精品资源，点击获取
