简介东北大学数值分析实验报告围绕解线性方程组的迭代法展开完整覆盖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 统一视角所有的迭代格式都来自同一个矩阵分裂线性方程组 Axb如果直接求逆一次高斯消元即可但大型稀疏矩阵不适合直接法。迭代法的基本思路是把 A 拆成A M - N其中 M 容易求逆于是 Axb 等价于M x N x b写成不动点形式x_{k1} 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 选择MDNLU得到迭代式x_{k1} D^{-1}(LU)x_k D^{-1}b写成分量形式就是x_i^{(k1)} (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 选择MD-LNU得到迭代式x_{k1} (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_{k1} (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 三种格式对照表方法迭代矩阵 B更新策略X0[0;0;0] 迭代次数JacobiD^{-1}(LU)所有分量同时更新14Gauss-Seidel(D-L)^{-1}U逐个分量更新新值立即用9SOR(D-ωL)^{-1}((1-ω)DωU)在 Gauss-Seidel 上做松弛11~11828随 ω 变化这张表的价值在于三种方法本质上是同一个矩阵分裂框架下的三个特例差异只在 M 怎么选。理解了这一层就不会再把“迭代次数”当作黑盒结果而是能根据地推导出为什么某种格式更快。3. 谱半径与松弛因子为什么 ω 越接近 2迭代次数越离谱3.1 收敛的充要条件迭代矩阵的谱半径必须小于 1无论哪种迭代格式都会得到x_{k1}B x_k c。真解 x* 满足x*B x*c。两个式子相减误差向量e_kx_k-x*递推为e_{k1}B e_k所以e_kB^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.10111.20141.30161.40191.50211.60321.70551.801531.9011828这个趋势不是线性的。ω 从 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]; Ddiag(diag(A)); L-tril(A,-1); U-triu(A,1); B_j D \ (LU); % 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、14Gauss-Seidel 都是 9SOR(ω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]; Ddiag(diag(A)); % 对角矩阵 L-tril(A,-1); % 下三角矩阵严格下三角取负 U-triu(A,1); % 上三角矩阵严格上三角取负这里最关键的是负号。tril(A,-1)返回去掉主对角线后的严格下三角部分triu(A,1)返回严格上三角部分。把它们取负是因为迭代公式要写成A D - L - U。如果直接令Ltril(A,-1)、Utriu(A,1)那就变成A D L U后续迭代矩阵的符号会乱。实验代码里X_1inv(D)*(LU)*X_0inv(D)*b正因为LU是非对角部分取负后的和所以迭代矩阵是D^{-1}(LU)符号正好对。这种写法在数学推导里非常清晰但它把矩阵分裂直接表达成了“算完一个矩阵再乘向量”的过程。对小型教学矩阵没问题但拿到大规模问题里这种形式会诱导你犯两个错一是显式求逆二是把迭代矩阵完整算出来。数值分析课上强调的“矩阵分裂”本质上是想让你用 M 的三角结构去解方程而不是真的把M^{-1}这个矩阵造出来。4.2 为什么说 inv(D) 是“形式解”不是“实现解”实验代码里把迭代直接写成inv(D)*(LU)*X_0inv(D)*b。这种写法在公式推导里非常清楚但在数值计算里不推荐。一是inv(D)显式求逆的计算量是 O(n^3)即使 D 是对角矩阵求逆本身不贵但写成inv(D)*矩阵会让代码失去“解线性方程组”的语义二是对于三角矩阵 M应该用左除M\b而不是先求逆再乘。MATLAB 的左除会选择合适的分解数值更稳定也不会额外构造逆矩阵。对于对角矩阵 DD\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 \ ((LU)*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 \ ((LU)*x b)等价于D^{-1}((LU)xb)但避免了显式求逆。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 i1:3 到底在检查什么原实验报告里的判据写得比较绕disabs(X_1-X_0); p1; for i1:3 while(dis(i)10^(-5)) X_1inv(D)*(LU)*X_0inv(D)*b; disabs(X_1-X_0); X_0X_1; pp1; end end这个循环不是标准的“所有分量都达到误差阈值就停止”。它是先看第一个分量是否达标不达标就一直迭代到第一个分量达标然后看第二个分量如果第二个分量不达标再继续迭代直到第二个分量达标然后看第三个。问题是当迭代继续时第一个分量可能又超出阈值但此时外层 for 已经进入 i2不会再回头检查第一个分量。所以最终停下来的条件严格说只是“第三个分量达标时停止”而不是“所有分量都达标”。对这个三维方程组由于分量间耦合不严重结果碰巧还正确但换一个矩阵这个判据就会给出错误的迭代次数甚至死循环。正确的做法是用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\(LU))))检查。如果算出来大于 1就不要硬跑换直接法或换预处理。很多“程序死循环”的 bug其实是数学上就不收敛不是代码问题。5.2 翻车现场二判据写成分量逐个检查结果随矩阵规模变化现象沿用实验报告里for i1:3while(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 翻车现场四用相邻迭代步误差作为“精度”结果算出的解其实不满足 Axb现象代码里用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/(1sqrt(1-ρ(J)^2))其中 ρ(J) 是 Jacobi 迭代矩阵的谱半径。但实际中 ρ(J) 往往未知要估计它本身就需要迭代或特征值计算代价可能比你省下的迭代次数还高。所以工程里常见做法是先用 Gauss-Seidel 跑几十步估算谱半径再据此选一个 ω或者干脆做一次 ω 扫描。这就是设计取舍Jacobi 和 Gauss-Seidel 是零参数、鲁棒但可能慢SOR 有速度潜力但引入了一个需要调参的“油门”调不好就比不用还慢。6.2 如果让我重写这份实验报告我会在代码里同时输出每次迭代的残差、相邻步误差和迭代矩阵谱半径并画出 ω-迭代次数曲线。这样实验报告就不是“记录三个算法的迭代次数”而是“验证谱半径决定收敛速度”的证据链。调 SOR 时我会先用eig算谱半径再调 ω每改一次参数就先看 ρ 是否下降而不是直接看迭代次数。这个习惯是从一次踩坑里养成的有次我把 ω 从 1.7 调到 1.8程序跑了十分钟没停最后发现谱半径已经是 0.9996才知道 SOR 的“加速”是有条件的。从那以后我每次调 SOR 都强制先算一遍谱半径曲线再决定 ω。希望帮到你。本文还有配套的精品资源点击获取
