迭代法收敛速度之谜：谱半径如何让SOR慢到11828次？

简介：东北大学数值分析实验报告，围绕解线性方程组的迭代法展开，完整覆盖Jacobi迭代法、Gauss-Seidel迭代法与SOR迭代法三种经典解法。报告按学校实验规范组织，包含实验项目名称、目的、原理、内容、数据记录与处理、结果分析等环节，适合数值分析课程学生用于理解迭代法收敛条件、不同迭代格式的差异以及松弛因子对收敛速度的影响。附件为1个doc文档，大小186KB，内含可直接运行的MATLAB代码与实验数据表格，代码采用矩阵分解与误差阈值控制实现迭代过程，读者可对照运行并比较不同初始向量、不同ω取值下的迭代次数与误差变化。目前已有2614人浏览学习，适合需要完成实验报告或复习迭代法数值实现的在校学生与自学者。

1. 同一个方程组，Jacobi 14 次，SOR 11828 次：迭代法的收敛速度到底由谁决定？

一份数值分析实验报告里，解线性方程组的迭代法对比实验，藏着让人意外的一组数字：同一个系数矩阵 A=[4,-1,1;4,-8,1;-2,1,5]，同一个右端项 b=[7;-21;15]，都用初始向量 X0=[0;0;0]，误差阈值取 10^(-5)。Jacobi 迭代 14 次就收敛，Gauss-Seidel 迭代 9 次收敛，而 SOR 迭代法在松弛因子 ω=1.9 时足足迭代了 11828 次。这个结果和“SOR 是 Gauss-Seidel 的加速版”的直觉完全相反：一个精心调了参数的算法，在最接近收敛边界的区域，反而慢到几乎不收敛。

这件事值得搞懂，因为它揭示的不是算法写错了，而是迭代法收敛速度的底层逻辑：同样的方程组，不同迭代格式对应不同的迭代矩阵，谱半径大小直接决定每步误差缩减比例。这篇博文会从实验报告里的三种格式出发，拆开 D、L、U 的矩阵分裂，讲清楚为什么 ω 接近 2 会让迭代次数爆炸，并落到 MATLAB 代码里那些容易忽略的坑。

2. 三种迭代格式的底层逻辑：把 A 拆成 D、L、U，就是在选择哪些信息立刻用

2.1 统一视角：所有的迭代格式都来自同一个矩阵分裂

线性方程组 Ax=b，如果直接求逆，一次高斯消元即可，但大型稀疏矩阵不适合直接法。迭代法的基本思路是把 A 拆成A = M - N，其中 M 容易求逆，于是 Ax=b 等价于M x = N x + b，写成不动点形式：

x_{k+1} = M^{-1} N x_k + M^{-1} b

只要M^{-1} N的谱半径小于 1，不断代入就能逼近真解。Jacobi、Gauss-Seidel、SOR 的区别，全在 M 和 N 的选择上。

这里需要先解释实验代码里的符号约定。对一个矩阵 A，取出对角部分 D，严格下三角部分取负得到 L，严格上三角部分取负得到 U，于是A = D - L - U。注意，迭代法里的 L、U 与 LU 分解里的三角因子不是一回事，它们带着负号。实验代码里L=-tril(A,-1)、U=-triu(A,1)就是为了构造这个约定。如果不加负号，后面所有迭代矩阵的符号都会乱掉。

2.2 Jacobi：所有分量都用旧值，天然并行

Jacobi 选择M=D，N=L+U，得到迭代式：

x_{k+1} = D^{-1}(L+U)x_k + D^{-1}b

写成分量形式就是：

x_i^{(k+1)} = (b_i - Σ_{j≠i} A_ij x_j^{(k)}) / A_ii

每个分量更新时，等号右边只包含上一次迭代的 x_k 分量，所以所有分量可以同时计算。这种“同时更新”的格式在并行机上很友好，代价是它没有利用已经算出来的新值。实验里 Jacobi 用 14 次收敛，说明对这个问题它也不算慢，但也没快到哪去。

一个值得记住的反直觉点：Jacobi 并不是毫无用处的“笨方法”。对于很多对角占优矩阵，它的下降速度是稳定的，而且并行效率极高。实验报告里它比 Gauss-Seidel 慢，只是因为这个矩阵的下三角信息确实有价值，而 Jacobi 选择无视这些信息。

2.3 Gauss-Seidel：新值立即生效，单线程更占便宜

Gauss-Seidel 选择M=D-L，N=U，得到迭代式：

x_{k+1} = (D-L)^{-1}U x_k + (D-L)^{-1}b

这里 M 是下三角矩阵，所以M^{-1}的过程本质上就是“从第一个方程到最后一个方程依次代入”。由于前一个分量刚被更新，后面分量马上就能用它，信息复用让同一轮迭代的信息传递更快。实验里它 9 次收敛，比 Jacobi 快。

但需要强调：Gauss-Seidel 并不总是优于 Jacobi。存在某些矩阵，Jacobi 收敛而 Gauss-Seidel 发散，也存在正好相反的矩阵。它不是 Jacobi 的“严格升级版”，而是一种改变了信息流向的格式。实验报告里这个矩阵下三角部分发挥了作用，所以我们看到 Gauss-Seidel 胜出，但这不能推广成“Gauss-Seidel 永远更快”。

2.4 SOR：给 Gauss-Seidel 加一个油门，但油门踩过头会翻车

SOR 的全称是逐次超松弛（Successive Over-Relaxation）。它在 Gauss-Seidel 的基础上引入权重 ω：先算一次 Gauss-Seidel 的增量，然后把这个增量乘上 ω，再加到当前解上。格式写成：

x_{k+1} = (D-ωL)^{-1}((1-ω)D+ωU)x_k + ω(D-ωL)^{-1}b

ω=1 时退化为 Gauss-Seidel；ω<1 是欠松弛，ω>1 是超松弛。从矩阵分裂的角度看，SOR 相当于在 M 里混合了 D 和 L 的比例，试图通过调整 M 的形状来压低迭代矩阵的谱半径。

实验报告里，当 ω 从 1.1 升到 1.9，迭代次数从 11 一路涨到 11828，说明这个“油门”不是越大越快。SOR 的收敛速度对 ω 极度敏感，真实世界里的最优 ω 往往需要靠估计或扫描得到，这也是它最“玄学”的地方。很多人以为 SOR 自带加速光环，实际上调不好会比 Jacobi 还慢。

2.5 一次手算迭代：看矩阵分裂如何决定新解

以 Jacobi 为例，初始 X0=[0;0;0]，第一次迭代等于直接把 b 的每个分量除以对角元：

X1 = [7/4; -21/(-8); 15/5] = [1.75; 2.625; 3]

精确解是[2;4;3]，所以第一次迭代后第一个分量误差 0.25，第二个分量误差 1.375，第三个分量已经精确。这个例子说明，迭代法的每一步都在把误差向零压缩，但压缩速度取决于迭代矩阵的结构。Gauss-Seidel 之所以能 9 步收敛，是因为它更新第二个分量时，会用上刚刚算出来的X1(1)；而 Jacobi 要等到下一轮才能用新值，信息滞后一轮。

2.6 三种格式对照表

这张表的价值在于：三种方法本质上是同一个矩阵分裂框架下的三个特例，差异只在 M 怎么选。理解了这一层，就不会再把“迭代次数”当作黑盒结果，而是能根据地推导出为什么某种格式更快。

3. 谱半径与松弛因子：为什么 ω 越接近 2，迭代次数越离谱

3.1 收敛的充要条件：迭代矩阵的谱半径必须小于 1

无论哪种迭代格式，都会得到x_{k+1}=B x_k + c。真解 x* 满足x*=B x*+c。两个式子相减，误差向量e_k=x_k-x*递推为e_{k+1}=B e_k，所以e_k=B^k e_0。收敛与否，取决于B^k是否趋于零矩阵。线性代数里有一个结论：B^k→0当且仅当 B 的所有特征值的模都小于 1，也就是谱半径ρ(B)<1。这就是迭代法最核心的判据。实验报告里没有算谱半径，但三种方法都能收敛，说明它们对应的迭代矩阵谱半径都小于 1。

一个更实用的推论是：谱半径越小，渐近收敛越快。如果ρ(B)=0.5，每步误差减半，10 步后误差缩到约千分之一；如果ρ(B)=0.99，每步只缩小 1%，同样缩到千分之一需要约 600 步。这能解释实验中 SOR(ω=1.9) 为什么崩溃：它的迭代矩阵谱半径一定非常接近 1。

3.2 把实验数据画成一条曲线：ω 从 1.1 到 1.9，谱半径在悄悄逼近 1

实验报告里记录了三组初始向量下的迭代次数。以 X0=[0;0;0] 为例：

这个趋势不是线性的。ω 从 1.1 到 1.5，迭代次数从 11 增加到 21，还算平缓；从 1.6 开始陡然变陡，到 1.9 几乎是垂直上升。如果把这个数据反过来推，SOR 迭代矩阵的谱半径随 ω 单调上升，在 ω=1.9 时可能已经达到 0.9999 级别。每一轮迭代误差只缩小万分之几，自然需要上万次。

这里需要强调：这个规律仅对这个具体矩阵成立。SOR 并不总是 ω 越大越慢；对很多矩阵，存在一个大于 1 的最优 ω，谱半径先降后升。这个矩阵的最优 ω 非常接近 1，才造成单调变差。所以“SOR 一定能加速”是误解，加速的前提是选对 ω。

3.3 用 MATLAB 验证谱半径：三行代码看到本质

实验报告里的迭代次数是宏观结果，谱半径是微观原因。可以用几行 MATLAB 直接验证：

A=[4,-1,1;4,-8,1;-2,1,5]; D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B_j = D \ (L+U); % Jacobi 迭代矩阵 rho_j = max(abs(eig(B_j))); % Jacobi 谱半径 B_gs = (D-L) \ U; % Gauss-Seidel 迭代矩阵 rho_gs = max(abs(eig(B_gs))); % Gauss-Seidel 谱半径 omega = 1.9; B_sor = (D-omega*L) \ ((1-omega)*D + omega*U); rho_sor = max(abs(eig(B_sor))); % SOR(1.9) 谱半径

eig返回矩阵全部特征值，abs取模，max找最大值，得到的就是谱半径。运行后你会看到rho_j和rho_gs都小于 1，且rho_gs明显小于rho_j；而rho_sor在 ω=1.9 时非常接近 1。这正是实验数据背后那个“看不见的变量”。调试代码时，先算谱半径，比闷头迭代几千次要有用得多。

3.4 初始向量影响迭代次数，但没有想象中那么大

实验报告还做了三组初始向量：X0=[0;0;0]、[1;1;1]、[-1;-1;-1]。Jacobi 的迭代次数分别是 14、13、14；Gauss-Seidel 都是 9；SOR(ω=1.25) 分别是 15、14、15；ω=1.5 分别是 21、20、21。差别很稳定，基本在 1 次以内。

原因在哪里？误差e_0 = x_0 - x*是三个初始向量与精确解[2;4;3]的差。e_k = B^k e_0，如果e_0恰好落在某个特征方向上，这一方向的收敛速度就是相应特征值的模；但一般情况下e_0是多个特征方向的混合，最终收敛速度由谱半径决定，初始向量只影响前几步的误差分布。所以实验里换初始向量，迭代次数只会小幅抖动，不会改变数量级。这个观察对理解迭代法很有用：调初始向量不是加速收敛的有效手段，真正有效的武器是改迭代格式或调松弛因子。

4. MATLAB 代码拆解：用 diag/tril/triu 表达矩阵分裂，为什么别用 inv(D) 硬乘

4.1 从实验代码看 D、L、U 的构造

实验报告里的代码是这样开始的：

A=[4,-1,1;4,-8,1;-2,1,5]; b=[7;-21;15]; D=diag(diag(A)); % 对角矩阵 L=-tril(A,-1); % 下三角矩阵（严格下三角取负） U=-triu(A,1); % 上三角矩阵（严格上三角取负）

这里最关键的是负号。tril(A,-1)返回去掉主对角线后的严格下三角部分，triu(A,1)返回严格上三角部分。把它们取负，是因为迭代公式要写成A = D - L - U。如果直接令L=tril(A,-1)、U=triu(A,1)，那就变成A = D + L + U，后续迭代矩阵的符号会乱。实验代码里X_1=inv(D)*(L+U)*X_0+inv(D)*b，正因为L+U是非对角部分取负后的和，所以迭代矩阵是D^{-1}(L+U)，符号正好对。

这种写法在数学推导里非常清晰，但它把矩阵分裂直接表达成了“算完一个矩阵再乘向量”的过程。对小型教学矩阵没问题，但拿到大规模问题里，这种形式会诱导你犯两个错：一是显式求逆，二是把迭代矩阵完整算出来。数值分析课上强调的“矩阵分裂”，本质上是想让你用 M 的三角结构去解方程，而不是真的把M^{-1}这个矩阵造出来。

4.2 为什么说 inv(D) 是“形式解”，不是“实现解”

实验代码里把迭代直接写成inv(D)*(L+U)*X_0+inv(D)*b。这种写法在公式推导里非常清楚，但在数值计算里不推荐。一是inv(D)显式求逆的计算量是 O(n^3)，即使 D 是对角矩阵，求逆本身不贵，但写成inv(D)*矩阵会让代码失去“解线性方程组”的语义；二是对于三角矩阵 M，应该用左除M\b而不是先求逆再乘。MATLAB 的左除会选择合适的分解，数值更稳定，也不会额外构造逆矩阵。对于对角矩阵 D，D\b其实就是逐个分量除以对角线元素，和inv(D)*b结果一样，但语义更准确。

下面给出一个更“数值分析风”的 Jacobi 实现：

A = [4,-1,1; 4,-8,1; -2,1,5]; b = [7;-21;15]; D = diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); tol = 1e-5; x = zeros(3,1); % 初始向量 for k = 1:10000 x_new = D \ ((L+U)*x + b); % 逐分量为 (b_i - sum_{j≠i} A_ij x_j) / A_ii if norm(x_new - x, inf) < tol break; end x = x_new; end fprintf('迭代次数: %d, 解: %g %g %g\n', k, x_new);

D \ ((L+U)*x + b)等价于D^{-1}((L+U)x+b)，但避免了显式求逆。norm(x_new-x,inf)计算所有分量差的最大绝对值，比只检查单个分量可靠。这里的循环上限 10000 只是保险，正常几轮就能 break。

如果把这份代码改成 Gauss-Seidel，只需要把迭代式换成：

x_new = (D-L) \ (U*x + b);

因为(D-L)是下三角矩阵，左除会做前代法，按行依次代入，速度比显式求逆快，数值也更稳。SOR 则写成：

omega = 1.5; x_new = (D-omega*L) \ ((1-omega)*D + omega*U) * x + omega * ((D-omega*L) \ b);

每次迭代都要解一个下三角方程组，这正是 SOR 比 Jacobi 计算量大的原因之一：它用更多的单步计算换取了更少的迭代次数。实验报告里 9 次 vs 14 次，看着差距不大，但对更大规模的矩阵，这个差距会被谱半径放大到几个数量级。

4.3 原代码里的 for i=1:3 到底在检查什么

原实验报告里的判据写得比较绕：

dis=abs(X_1-X_0); p=1; for i=1:3 while(dis(i)>=10^(-5)) X_1=inv(D)*(L+U)*X_0+inv(D)*b; dis=abs(X_1-X_0); X_0=X_1; p=p+1; end end

这个循环不是标准的“所有分量都达到误差阈值就停止”。它是先看第一个分量是否达标，不达标就一直迭代到第一个分量达标；然后看第二个分量，如果第二个分量不达标，再继续迭代，直到第二个分量达标；然后看第三个。问题是，当迭代继续时，第一个分量可能又超出阈值，但此时外层 for 已经进入 i=2，不会再回头检查第一个分量。所以最终停下来的条件，严格说只是“第三个分量达标时停止”，而不是“所有分量都达标”。对这个三维方程组，由于分量间耦合不严重，结果碰巧还正确，但换一个矩阵，这个判据就会给出错误的迭代次数，甚至死循环。正确的做法是用norm(X_1-X_0,inf) < 1e-5作为整体判据。

由此可以看到，实验报告代码里“能跑出正确结果”并不等于“实现正确”。数值分析实验的常见坑，就是终止条件写得似是而非，结果数据看起来对，换个场景就崩。这也是下一章要展开的内容。

5. 避坑指南：迭代法实验里最常出现的四个翻车现场

5.1 翻车现场一：格式写对了，就是不收敛

现象：照着教材抄了 Jacobi 或 Gauss-Seidel 公式，初始向量也给了，结果迭代几千次都不收敛，或者数值溢出。

原因：迭代矩阵谱半径大于等于 1。很多直接法里“矩阵可逆”就够了，但迭代法要求更苛刻，通常需要矩阵具备某种性质，比如严格对角占优或对称正定。如果原始矩阵不满足这些性质，即使用正确的公式也可能发散。实验里的矩阵 A=[4,-1,1;4,-8,1;-2,1,5] 是严格对角占优的：第一行 |4|>|-1|+|1|，第二行 |-8|>|4|+|1|，第三行 |5|>|-2|+|1|。这个条件保证了 Jacobi 收敛，也是实验能顺利跑完的前提。

解决：在写迭代法之前，先算一下迭代矩阵的谱半径。比如用rho = max(abs(eig(D\(L+U))))检查。如果算出来大于 1，就不要硬跑，换直接法或换预处理。很多“程序死循环”的 bug，其实是数学上就不收敛，不是代码问题。

5.2 翻车现场二：判据写成分量逐个检查，结果随矩阵规模变化

现象：沿用实验报告里for i=1:3+while(dis(i)>=tol)的写法，在三维问题上结果正常，一旦把代码改成 n 维，迭代次数突然变得很奇怪，或者有些维度还没收敛就停了。

原因：如上一章分析，这个判据实际上只保证最后一个被检查的分量达标，而不是整体误差。n 变大后，前面的分量可能因为没有后续检查而“放飞自我”，循环却已结束。

解决：统一用norm(x_new - x_old, inf) < tol作为终止条件，或者用while norm(x_new-x_old,inf) >= tol。inf范数就是最大分量差，简单直观。如果坚持用分量判据，也要写成all(abs(x_new-x_old) < tol)，而不是嵌套 for 和 while。

5.3 翻车现场三：SOR 的 ω 调大后迭代次数暴涨，以为是程序错误

现象：SOR 把 ω 从 1.5 调到 1.8，迭代次数从 21 跳到 153；从 1.8 到 1.9，直接变成 11828。很多人第一反应是inv(D-k*L)哪里写错了。

原因：迭代次数暴涨不是 bug，而是谱半径随 ω 增大趋近于 1。SOR 的迭代矩阵在 ω 接近 2 时，谱半径会趋于 1，导致收敛速度趋近于零。这符合 SOR 理论：松弛因子的有效范围是 (0,2)，但谱半径的分布随矩阵而定。对这个矩阵，最优 ω 就在 1 附近，ω 越大越糟糕。

解决：不要盲目加大 ω。实验前先画一条“ω-迭代次数”曲线，或者直接计算 SOR 迭代矩阵的特征值，观察谱半径随 ω 的变化。找到曲线最低点再定 ω。另外，SOR 的迭代矩阵格式里用了 k，注意它和 ω 是同一个参数，别在循环里混用。

5.4 翻车现场四：用相邻迭代步误差作为“精度”，结果算出的解其实不满足 Ax=b

现象：代码里用norm(X_1-X_0, inf)<1e-5作为终止条件，最后输出 X=[2.00004;4.0000;3.0000]，看起来误差很小，但把 X 代回原方程，残差norm(A*X-b)其实是 1e-4 甚至更大。

原因：迭代法常用的停机准则是“相邻两步之差足够小”，但两步之差小不代表当前解与真解误差小，更不代表残差小。特别当迭代矩阵谱半径接近 1 时，两步之差会被压缩，但距离真解还可能很远。实验报告里记录的是“迭代解法”在停止时与精确解的差，这个差异是 10^(-5) 量级，但如果你只看相邻步差，可能会误判精度。

解决：如果实验报告要求验证精度，就在收敛后加一行res = norm(A*x_new - b, inf)，把残差也打印出来。严谨的话，把停机准则写成“残差小于允许值”或“同时满足残差和步长阈值”。在做数值分析实验时，这一步能帮你发现很多“看似收敛实则没收敛”的情况。实验报告里那个 2.00004 其实也提示了这一点：解的实际误差并没有达到 10^(-6)，只是相邻步误差达到了阈值。

6. 设计取舍：用通用性换速度，还是用参数估计换稳健性？

6.1 为什么最优 ω 很少被直接算出来

SOR 的理论很漂亮：对对称正定矩阵，最优松弛因子是ω_opt = 2/(1+sqrt(1-ρ(J)^2))，其中 ρ(J) 是 Jacobi 迭代矩阵的谱半径。但实际中 ρ(J) 往往未知，要估计它本身就需要迭代或特征值计算，代价可能比你省下的迭代次数还高。所以工程里常见做法是先用 Gauss-Seidel 跑几十步，估算谱半径，再据此选一个 ω，或者干脆做一次 ω 扫描。这就是设计取舍：Jacobi 和 Gauss-Seidel 是零参数、鲁棒但可能慢；SOR 有速度潜力，但引入了一个需要调参的“油门”，调不好就比不用还慢。

6.2 如果让我重写这份实验报告

我会在代码里同时输出每次迭代的残差、相邻步误差和迭代矩阵谱半径，并画出 ω-迭代次数曲线。这样实验报告就不是“记录三个算法的迭代次数”，而是“验证谱半径决定收敛速度”的证据链。调 SOR 时，我会先用eig算谱半径，再调 ω；每改一次参数，就先看 ρ 是否下降，而不是直接看迭代次数。

这个习惯是从一次踩坑里养成的：有次我把 ω 从 1.7 调到 1.8，程序跑了十分钟没停，最后发现谱半径已经是 0.9996，才知道 SOR 的“加速”是有条件的。从那以后，我每次调 SOR 都强制先算一遍谱半径曲线，再决定 ω。希望帮到你。

本文还有配套的精品资源，点击获取