SOR松弛因子微调，迭代次数从153暴涨到11828的背后

简介这是东北大学数值分析课程的一份完整实验报告面向数学系及软件等专业学生可用于学习线性方程组迭代法、撰写实验报告及参考MATLAB代码。资源包内仅含1个doc文档大小约186KB内容涵盖Jacobi迭代法、Gauss-Seidel迭代法与SOR松弛迭代法的原理、迭代公式、数据记录、结果分析及可直接运行的MATLAB实现。报告以典型三阶方程组为例在相同精度下对比三种方法的迭代次数与误差并针对SOR法给出ω从1.10到1.90的完整实验数据直观反映松弛因子对收敛速度的影响同时记录了不同初始向量下的迭代表现可帮助理解初值对迭代法收敛过程的作用。内容遵循东北大学理学院实验报告规范包含实验目的、原理、步骤、数据处理、结果分析等要素并预留指导教师签字栏可直接作为课程实验报告模板参考。目前已有2614人学习下载适合正在准备数值分析实验、复习迭代算法或需要MATLAB实现参考的同学。1. 一个ω从1.80调到1.90迭代次数为什么从153暴涨到11828SOR迭代法里松弛因子ω从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、UA [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 -(LU)Gauss-Seidel 取 M D - L、N USOR 取 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 - (LU)X于是 X D^{-1}(LU)X D^{-1}b。写成迭代式就是X_{k1} D^{-1}(LU)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)*(LU)*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)*(LU)*X_0 inv(D)*b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end注意D^{-1}(LU)和D^{-1}b其实跟迭代步无关理想实现应该在循环外只算一次但对3阶矩阵这种开销无所谓。真正重要的行为是每一步迭代X 的三个新分量都完全由 X_0 的旧值计算出来谁也不占谁的便宜。这就是“同时更新”天然可并行。坏处也在这里第一个分量 x1 的新值已经算出来了第二个方程里却还拿着旧值在算等于每轮都浪费了一次“信息复用”的机会。2.3 Gauss-Seidel新值刚算出来就被下一个方程拿去用Gauss-Seidel 的分量形式是算第 i 个新分量时小于 i 的角标全部用本轮新值大于 i 的角标才用旧值。写成方程组再整理就是X_{k1} (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 必须等 x1x3 必须等 x2没法直接并行。对3阶小矩阵这无所谓但对大规模并行求解这个差异会直接影响算法选型。2.4 SOR给Gauss-Seidel装上一个可调“油门”SOR 的矩阵形式是把 GS 的更新量再乘一个 ω然后加回旧值X_{k1} (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)*Dk*U)*X_0 k*inv(D-k*L)*b;当 ω1 时(1-ω)D ωU退化成 U公式就完整退回 Gauss-Seidelω1 叫欠松弛ω1 叫超松弛。所谓“松弛”字面意思是把这一步的更新量朝目标方向再推一把或收一把。推过头了不会立刻发散但会在目标附近来回震荡白白多跑很多轮——实验里 ω1.90 那 11828 次就是这么来的。三种方法写全了就是一张对照表方法迭代矩阵每步主要成本JacobiD^{-1}(LU)一次矩阵乘向量Gauss-Seidel(D-L)^{-1}U一次下三角回代SOR(D-ωL)^{-1}((1-ω)DωU)一次下三角回代三种迭代矩阵同构于“A M - N 分裂”这一个框架区别只在 M 拿走多少、N 剩下多少。这也是为什么同一份实验报告能一口气对比三种方法它们不是三个独立算法而是一个分裂框架下的三个参数化实例。3. 收敛速度的底层机制谱半径、渐近收敛因子与那11828次3.1 收敛的充要条件与“严格对角占优”这个安全垫把三种迭代都写成 X_{k1} B X_k f 之后误差服从 e_{k1} 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 74×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.10111.20141.30161.40191.50211.60321.70551.801531.9011828这组数据有两个特征值得琢磨前半个区间几乎是平缓的线性增长到 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]但把三组结果放在一起看会发现迭代次数几乎没怎么动初始向量JacobiGauss-SeidelSOR(ω1.25)SOR(ω1.50)SOR(ω1.75)[1;1;1]139142079[-1;-1;-1]149152187初值翻来覆去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 挂机跑实验回来发现循环还在转白白等了半小时——这种黑匣子式的试错本质上是把谱半径这个总开关误当成“调参游戏”。希望帮到你。本文还有配套的精品资源点击获取