简介：这是东北大学数值分析课程的一份完整实验报告，面向数学系及软件等专业学生，可用于学习线性方程组迭代法、撰写实验报告及参考MATLAB代码。资源包内仅含1个doc文档，大小约186KB，内容涵盖Jacobi迭代法、Gauss-Seidel迭代法与SOR松弛迭代法的原理、迭代公式、数据记录、结果分析及可直接运行的MATLAB实现。报告以典型三阶方程组为例，在相同精度下对比三种方法的迭代次数与误差，并针对SOR法给出ω从1.10到1.90的完整实验数据，直观反映松弛因子对收敛速度的影响；同时记录了不同初始向量下的迭代表现，可帮助理解初值对迭代法收敛过程的作用。内容遵循东北大学理学院实验报告规范，包含实验目的、原理、步骤、数据处理、结果分析等要素，并预留指导教师签字栏，可直接作为课程实验报告模板参考。目前已有2614人学习下载，适合正在准备数值分析实验、复习迭代算法或需要MATLAB实现参考的同学。

1. 一个ω从1.80调到1.90：迭代次数为什么从153暴涨到11828

SOR迭代法里，松弛因子ω从1.80调到1.90，表面看只是把“步子迈大”了一点点，可同一份代码、同一个10^(-5)的误差阈值下，迭代次数从153次暴涨到11828次。反过来，ω=1.10时这个数字只有11次。也就是说，一个看起来无足轻重的参数，能把同一个线性方程组的求解时间从“瞬间完成”拖到“肉眼可见卡住”。这正是数值分析里迭代法最容易被当成玄学的地方：Jacobi、Gauss-Seidel、SOR三种迭代法，拆开看不过是一行矩阵递推，合起来却牵涉矩阵分裂、谱半径、收敛性三大底层概念。本文要拆的这份某高校数值分析课程的实验报告，恰好用一套完整的MATLAB代码和三组实验数据，把这三件事串了起来，值得逐行读懂。

2. 一个矩阵分裂，三个递推式：Jacobi、Gauss-Seidel、SOR的数学同构

2.1 把A拆成D、L、U：负号约定不是玄学

实验报告的代码第一屏就是在拆矩阵。A是那个3阶方阵，b是右端项，然后一行一行取出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)都跳过了对角线，取出的是“严格”下三角和严格上三角；第二，取出来之后立刻乘了负号；第三，D、L、U 合起来满足 A = D - L - U，而不是直觉上的 A = D + L + U。

为什么要用减号？因为后面三种迭代法都可以写成“矩阵分裂”的形式：把 A 拆成 M - N，其中 M 要求可逆并且容易求逆。Jacobi 取 M = D、N = -(L+U)，Gauss-Seidel 取 M = D - L、N = U，SOR 取 M = (1/ω)(D - ωL)、N = (1/ω)((1-ω)D + ωU)。一旦按这个约定拆，三种迭代公式的右端项都干净得像教材原样，迭代矩阵的谱半径也直接落在我们关心的位置。如果反过来用 A = D + L + U 的拆法，代码里会到处冒出负号，稍不留神就把迭代矩阵写错。常见的实现思路就是像实验报告这样，tril/triu取完之后统一乘 -1，让公式和教材保持完全对齐。

2.2 Jacobi：三个分量同一轮齐步走

Jacobi 的推导很直接：由 A = D - L - U 得 DX = b - (L+U)X，于是 X = D^{-1}(L+U)X + D^{-1}b。写成迭代式就是：

X_{k+1} = D^{-1}(L+U)X_k + D^{-1}b

实验报告里的代码完全对应这个式子：

D = diag(diag(A)); % 对角矩阵 L = -tril(A,-1); % 严格下三角 U = -triu(A,1); % 严格上三角 X_0 = [0;0;0]; % 初始向量 X_1 = inv(D)*(L+U)*X_0 + inv(D)*b; % 第一步迭代 dis = abs(X_1 - X_0); % 相邻两步解的绝对差 p = 1; % 已迭代1次 for i = 1:3 while dis(i) >= 1e-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

注意D^{-1}(L+U)和D^{-1}b其实跟迭代步无关，理想实现应该在循环外只算一次，但对3阶矩阵这种开销无所谓。真正重要的行为是：每一步迭代，X 的三个新分量都完全由 X_0 的旧值计算出来，谁也不占谁的便宜。这就是“同时更新”，天然可并行。坏处也在这里：第一个分量 x1 的新值已经算出来了，第二个方程里却还拿着旧值在算，等于每轮都浪费了一次“信息复用”的机会。

2.3 Gauss-Seidel：新值刚算出来就被下一个方程拿去用

Gauss-Seidel 的分量形式是：算第 i 个新分量时，小于 i 的角标全部用本轮新值，大于 i 的角标才用旧值。写成方程组再整理，就是：

X_{k+1} = (D-L)^{-1}U X_k + (D-L)^{-1}b

实验报告里 GS 的循环体跟 Jacobi 几乎长得一模一样，只是把inv(D)换成了inv(D-L)。这里有个容易误读的点：矩阵形式的代码依然是“X_0 全部旧值算出 X_1”，看起来不像就地更新。但实际上(D-L)^{-1}右乘一个向量，等价于一次自下而上的三角回代，新算出的分量在同一个表达式的后续方程里立刻生效。换句话说，矩阵形式把分量式里“边算边用”的过程压缩进了一次回代里。

为什么 GS 通常更快？因为每一轮迭代里，已经更新的信息被后面的方程及时消化掉了，不像 Jacobi 要等到下一轮才用上。代价是这种更新天然串行：x2 必须等 x1，x3 必须等 x2，没法直接并行。对3阶小矩阵这无所谓，但对大规模并行求解，这个差异会直接影响算法选型。

2.4 SOR：给Gauss-Seidel装上一个可调“油门”

SOR 的矩阵形式是把 GS 的更新量再乘一个 ω，然后加回旧值：

X_{k+1} = (D-ωL)^{-1}((1-ω)D + ωU)X_k + ω(D-ωL)^{-1}b

实验代码里直接用 k 表示 ω：

k = 1.75; % 松弛因子 ω D = diag(diag(A)); % 对角矩阵 L = -tril(A,-1); % 严格下三角 U = -triu(A,1); % 严格上三角 X_0 = [0;0;0]; % 初始向量 X_1 = inv(D-k*L)*((1-k)*D+k*U)*X_0 + k*inv(D-k*L)*b;

当 ω=1 时，(1-ω)D + ωU退化成 U，公式就完整退回 Gauss-Seidel；ω<1 叫欠松弛，ω>1 叫超松弛。所谓“松弛”，字面意思是把这一步的更新量朝目标方向再推一把或收一把。推过头了不会立刻发散，但会在目标附近来回震荡，白白多跑很多轮——实验里 ω=1.90 那 11828 次就是这么来的。

三种方法写全了就是一张对照表：

三种迭代矩阵同构于“A = M - N 分裂”这一个框架，区别只在 M 拿走多少、N 剩下多少。这也是为什么同一份实验报告能一口气对比三种方法：它们不是三个独立算法，而是一个分裂框架下的三个参数化实例。

3. 收敛速度的底层机制：谱半径、渐近收敛因子与那11828次

3.1 收敛的充要条件，与“严格对角占优”这个安全垫

把三种迭代都写成 X_{k+1} = B X_k + f 之后，误差服从 e_{k+1} = B e_k，于是 e_k = B^k e_0。要让误差对任意初值都趋于0，核心条件只有一个：迭代矩阵 B 的谱半径 ρ(B) < 1。谱半径是 B 所有特征值模长的最大值，它直接决定误差每步被压缩多少。

谱半径不好算，所以实验课都会给一个更省事的充分条件：严格对角占优。也就是每一行，对角元绝对值的和必须大于该行所有非对角元绝对值之和。对报告里这个矩阵逐行检查：|4| > |-1| + |1| = 2；|-8| > |4| + |1| = 5；|5| > |-2| + |1| = 3。三行全满足，所以 Jacobi 和 Gauss-Seidel 的收敛性在理论上就有保证。这也解释了实验结论里“三种迭代法均收敛”不是运气，而是矩阵本身选得好：一个严格对角占优矩阵，天然就是演示迭代法收敛性最安全的教学道具。

SOR 的收敛性比前两者更挑矩阵。严格对角占优并不自动保证任意 ω∈(0,2) 的 SOR 都收敛，本实验里 ω 从 1.10 一路测到 1.90 确实都收敛，但不能把这个结论推广到所有矩阵。能做的最稳妥的验证，就是像实验报告这样逐点实测。

顺便验证一下解的精确值。报告表格里写的“精确解243”，其实是 [2;4;3] 的连写。带回原方程：4×2 - 1×4 + 1×3 = 7；4×2 - 8×4 + 1×3 = -21；-2×2 + 1×4 + 5×3 = 15。右端项完全吻合，说明这个例子选的是有唯一解且解很规则的3阶系统，便于学生在迭代次数上做对比。

3.2 为什么Gauss-Seidel比Jacobi快：信息复用

同一初值 [0;0;0]、同一阈值 1e-5 下，Jacobi 用了 14 次，Gauss-Seidel 用了 9 次。差的这 5 次就是“信息复用”的收益。

用分量视角看：Jacobi 每一轮的三次更新都读旧值，GS 每算出一个新 x1，立刻被第二个方程拿去用；新 x1、x2 又立刻被第三个方程拿去用。相当于每一轮迭代里，GS 比 Jacobi 多消化了半轮的新信息，所以同样走一步，GS 的“有效信息量”更高。需要注意这并不绝对，教材上有 Jacobi 收敛而 Gauss-Seidel 不收敛的反例，但在严格对角占优的矩阵上，GS 的优势是典型表现。

3.3 ω的谱半径轨迹：误差按ρ^k衰减，尾部必然爆炸

实验报告里 [0;0;0] 初值下有一组完整的 ω 扫描数据：

这组数据有两个特征值得琢磨：前半个区间几乎是平缓的线性增长，到 1.80 之后开始失控，1.90 直接给出一根刺穿图表的天线。

原因是 SOR 迭代矩阵的谱半径随 ω 变化不是线性的。已知 Jacobi 迭代矩阵的谱半径可以算一下：对这个 A，ρ(B_J) ≈ 0.335。代入最优松弛因子公式 ω_opt = 2 / (1 + √(1 - ρ²))，得到 ω_opt ≈ 1.03。实验数据里 1.10 最快（11次），和理论最优位置完全对得上；一旦越过最优点，谱半径开始回升，ω→2 时谱半径趋近于1。

误差衰减是 e_k ≈ ρ^k e_0，想达到 1e-5，需要的步数粗略估计是 k ≈ ln(ε) / ln(ρ)。ρ=0.335 时这个值只有约 10 步出头，对应 Jacobi 的 14 次非常合理；而当 ρ 逼近 0.99 甚至 0.999 时，k 会从几百跳到上万，正是 153 次到 11828 次之间发生的跃变。这也是为什么迭代法实验里，“跑不完”往往不是死循环，而是谱半径已经贴到了1。

4. 避坑：改这份MATLAB代码之前，先看这四个细节

4.1 停机判断用的是相邻迭代差，不是残差

现象：代码里dis = abs(X_1 - X_0)，阈值 1e-5，看起来“误差”达标了，但如果你另外算一下残差，可能还没到 1e-5。

原因：dis衡量的是两次迭代之间解的变化量，不是方程本身的残差 ‖b - AX_k‖。收敛时相邻迭代差确实趋近于0，但反过来不一定成立。矩阵条件数大时，可能出现“解已经几乎不动，但残差还明显偏大”的情况，这时候按相邻差停机就会提前收工。

解决：做实验时至少加一行验证：norm(b - A*X_1, inf)，让残差也低于阈值再停止循环。更严谨的写法是把停机条件直接改成残差范数，而不是相邻解差。

4.2 for套while的怪异结构，等价于while all，不是提前停机

现象：初次读这段代码，很容易以为某个分量先达标循环就停了，甚至担心它会漏算其他分量。

原因：仔细看for i = 1:3里套的是三个串行的while dis(i) >= 1e-5。第一个 while 退出后，第二个、第三个 while 还会继续执行，而每一个 while 内部一旦进入，都对全部三个分量做一轮完整迭代。所以整体行为等价于“只要还有分量没达标就继续迭代”，也就是while any(dis >= 1e-5)的逐段展开。

解决：想要一眼可读的版本，直接写while max(dis) >= 1e-5，效果完全一样。实验报告里这个写法在3阶问题上没问题，但它的可读性确实容易被误解。

4.3 小矩阵里敢用inv，规模一上来就翻车

现象：3阶矩阵用inv(D)、inv(D-L)看起来很直观，算得也快；把同样的写法搬到一个几百阶的稀疏矩阵上，速度瞬间崩掉，还可能算出明显不对的解。

原因：inv是显式求逆。稀疏矩阵的逆通常不再稀疏，显式求逆会破坏稀疏性，存储和乘法成本都失控；此外，求逆本身的数值误差会被放大，直接乘到右端项上。

解决：循环外先把迭代矩阵算好一次，循环体内只做矩阵乘向量；需要解(D-L)\v时，用反斜杠运算符(D-L)\v做回代，而不是inv(D-L)*v。这在小矩阵上结果一致，在大矩阵上是能不能跑完的区别。

4.4 ω=1.90已经踩在悬崖边上，再往2走就是死循环

现象：1.80 是 153 次，1.90 是 11828 次。如果试 ω=1.99，这个实验基本跑不完；真取 ω=2，代码会永远不退出。

原因：SOR 的迭代矩阵在 ω→2 时谱半径趋近于1，误差每步几乎不衰减，1e-5 可能要百万步才能摸到；ω=2 时谱半径达到1，误差不再被压缩，循环条件永远满足。

解决：选 ω 之前先估算最优点，比如用 Jacobi 谱半径代入公式 ω_opt = 2/(1+√(1-ρ²))，把实验范围缩在最优值附近。靠近2的区域不是“更快的捷径”，而是“收敛性正在消失的悬崖”。

5. 数据背后的设计：三组初值和一串ω，到底在验证什么

5.1 三组初值的迭代次数几乎一样：线性收敛的初值不敏感性

实验报告特意测了三组初始向量：[0;0;0]、[1;1;1]、[-1;-1;-1]，但把三组结果放在一起看，会发现迭代次数几乎没怎么动：

初值翻来覆去，GS 始终是 9 次，Jacobi 在 13 和 14 之间浮动，SOR 的差别也在个位数。这不是巧合，而是线性迭代的固有性质：误差 e_k = B^k e_0，初值只影响前面的常数系数，真正决定每步压缩多少的是谱半径 ρ。初值差几步，本质是前几轮各自的路径差异，迭代次数的大头完全由 ρ 决定。

实验故意设计三组初值，恰恰是想让观察者自己看出“初值对迭代次数影响很小”，从而把注意力从“调初值”移到“调迭代格式和松弛因子”上来。理解了这一点，以后遇到“换个初值迭代次数暴涨”的情况，第一反应就不该是调初值，而是检查迭代矩阵的谱半径是不是已经逼近1。

5.2 ω序列从1.10排到1.90：用一次实验“画”出谱半径曲线

1.10 到 1.90 每隔 0.1 取一个点，九个点连起来就是 SOR 谱半径右半支的采样。前段 11、14、16、19、21 几乎是线性爬坡，说明谱半径还在温和回升；到 1.70 的 55、1.80 的 153，斜率突然变陡；1.90 的 11828 则直接掉进“谱半径无限逼近1”的区域。

这种排点方式本身就是实验设计的核心：先放一个接近理论最优的 1.10，再逐步向危险区 2 靠拢，用迭代次数当“探针”去反推谱半径的形状。所以数据表不是记录一堆无意义的数字，而是在做一次绕开特征值计算的数值实验。

5.3 这份文档资料的真正用法：可复跑、可改造的迭代法基准代码

这份材料既是规范的实验报告模板，也是一份能直接跑的 MATLAB 迭代法基准代码。拿到手之后，最值得做的不是照着抄一遍，而是改三个地方观察行为变化：

把阈值从 1e-5 改成 1e-8，三种方法的迭代次数会同步上升，但 J、GS、SOR 之间的相对差距基本不变，可以验证“阈值只改变截止位置，不改变收敛速度的排序”。把矩阵 A 的某个对角元改小，破坏严格对角占优条件，比如把 A(1,1) 从 4 改成 1，立刻能看到 Jacobi 和 GS 不再保证收敛，实验会翻车或死循环。固定 ω=1.90 完整跑一次，用 11828 次循环直观感受“谱半径接近1”意味着什么。

这三组改造跑完，收敛性理论就不再是书上的定理，而是自己手里的实验现象。

6. 设计取舍：迭代法不是万能药，谱半径才是总开关

6.1 直接法与迭代法的分界线

回到最根本的问题：这个3阶方程组用A\b一行就能解出来，为什么还要研究迭代法？答案是规模。直接法的复杂度是 O(n^3)，几百阶就明显吃力；迭代法每步只是矩阵乘向量或三角回代，配合稀疏存储，单步成本是 O(nnz)，对大规模稀疏系统优势巨大。但迭代法的总成本不是固定值，它取决于谱半径——碰到 ρ=0.999 这种病态情况，一万步叠加起来反而比直接法还贵。

所以工程上的取舍很清楚：稠密小矩阵，直接法省心；稀疏大矩阵，迭代法才有意义；而迭代法能不能用，先看谱半径，不看迭代格式的“高级感”。

6.2 我做迭代实验的固定动作

从那以后，我拿到一个线性方程组，第一件事永远是先算迭代矩阵的谱半径，再决定要不要用迭代法，然后才轮得到选 GS 还是 SOR、ω 取多少。跑完循环也一定会加一行残差验证，而不是看到相邻误差达标就收工。有一次我图省事，没检查谱半径就把 ω 设成 1.9 挂机跑实验，回来发现循环还在转，白白等了半小时——这种黑匣子式的试错，本质上是把谱半径这个总开关误当成“调参游戏”。希望帮到你。

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