简介东北大学数值分析实验报告围绕线性方程组迭代解法完整呈现Jacobi迭代法、Gauss-Seidel迭代法与SOR迭代法的原理、步骤及MATLAB实现通过同一算例记录不同方法的迭代次数与误差直观对比收敛速度及松弛因子ω的影响面向数值分析课程学习者、需仿写实验报告的高校学生。资源为1个doc文档压缩包大小约186KB内含实验目的、实验原理、数据记录表、完整代码段与结果分析结构清晰便于直接参考或按需修改。已有2614人学习可作为数值分析实验课的有效参考资料。借助附带的MATLAB代码可重现三种迭代法从初始向量逐步逼近精确解的过程帮助理解迭代格式构造与收敛条件为后续更复杂的数值计算问题奠定基础。1. 数值分析实验报告不只有公式一份能直接跑通的迭代法对比资料同一个三阶方程组Jacobi 迭代 14 次才收敛Gauss-Seidel 只要 9 次而 SOR 方法把松弛因子 ω 调到 1.9 时迭代次数竟然涨到 11828 次。这个反直觉的结果来自一份某高校理学系数值分析课程的实验报告内容是把三种常用迭代法拆开对比附完整 MATLAB 源码、实验数据表格和结论推导属于拿到就能直接复现的课程文档资料。对正在写数值分析实验报告的学生来说它是现成的样板对想搞清楚迭代法收敛速度差异的从业者来说它把参数、数据、结论之间的对应关系摆得很清楚。适合两类人要交实验报告的学生和被 SOR 收敛性绕晕的开发者。2. 三种迭代法的收敛逻辑矩阵分解框架与松弛因子的真实边界三种方法共用同一个矩阵拆法。把系数矩阵 A 写成 A D − L − UD 是对角矩阵L、U 是严格下三角、严格上三角矩阵。注意这里 L、U 前面都带了负号这是为了把迭代公式写成「新解 迭代矩阵 × 旧解 常数」的标准形式。MATLAB 里三行拆完但符号和函数参数都有讲究。2.1 同一个 D-L-U 分解两种更新策略决定收敛节奏D diag(diag(A)); % 内层 diag 抽出对角线向量外层 diag 扩回对角矩阵 L -tril(A, -1); % 取严格下三角并取负-1 表示排除主对角线 U -triu(A, 1); % 取严格上三角并取负1 表示从主对角线上方开始diag 一个函数两种用法是 MATLAB 新手最容易懵的地方内层 diag(A) 把矩阵对角线抽成列向量外层 diag 再把向量扩成对角矩阵。tril(A,-1) 和 triu(A,1) 里的第二个参数控制从哪里开始取-1 和 1 都是为了把主对角线排除在外。负号不能省直接写 tril(A,-1) 不取负下面的迭代矩阵符号会整体错位收敛出的解直接带反号。拆完之后Jacobi 和 Gauss-Seidel 的差别就一句话Jacobi 把 D 单独留在左边其余全部当作旧解系数Gauss-Seidel 把 D−L 挪到左边系数里只留下 U。Jacobi X_{k1} D^{-1}(LU)X_k D^{-1}b Gauss-Seidel X_{k1} (D-L)^{-1} U X_k (D-L)^{-1}b看起来只差一个 L 挪到哪边实际差异在更新顺序。Jacobi 这一轮所有分量都用上一轮的值算天然适合并行Gauss-Seidel 每算出一个新分量立刻覆盖旧值后面的分量用本轮最新值串行结构让它在很多问题上收敛更快。这份报告里同一方程组、同一起点 [0;0;0]Jacobi 14 次、Gauss-Seidel 9 次就是这个差异的直接体现。还要交代一个前提实验用的 A 满足严格对角占优|4||−1||1|、|−8||4||1|、|5||−2||1| 全部成立所以两种方法对任意初始向量都收敛。实验里换了三组初始向量全部收敛和理论对得上。如果 A 不满足这个条件迭代矩阵谱半径可能大于 1那就不是 14 次还是 9 次的问题而是根本不收敛。2.2 松弛因子 ω 不是越大越快数据把反例摆得很直白SOR 是在 Gauss-Seidel 基础上加一个松弛因子 ω标准写法是 (D−ωL)X_{k1} ((1−ω)DωU)X_k ωb。ω1 时它退化成 Gauss-Seidelω1 叫超松弛ω1 叫欠松弛。很多教材把 SOR 讲成「加速版 Gauss-Seidel」但这份实验数据给出了非常直白的反例。ω1.101.201.301.401.501.601.701.801.90迭代次数1114161921325515311828初始向量 [0;0;0]、容差 1e-5 不变ω 从 1.1 涨到 1.9迭代次数从 11 一路涨到 11828后段几乎是垂直起飞。原因在于 SOR 的收敛范围是 0ω2但收敛快慢由 SOR 迭代矩阵的谱半径决定ω 越接近 2谱半径越贴向 1每轮误差只能衰减一点点自然要跑成千上万轮。1.9 那行在实验里基本等于死循环。更值得说的是最优 ω。对一致排序的矩阵有一套经典估算公式 ω_opt 2/(1√(1−ρ²))ρ 是 Jacobi 迭代矩阵 D^{-1}(LU) 的谱半径。这份报告的 A 不满足这套公式的全部前提但可以拿来粗算ρ 约在 0.48 附近ω_opt 大约 1.06实验里 ω1.10 恰好是整张表迭代次数最少的一列。理论估算和经验最优对得上这个细节值得写进报告的分析部分。结论就是ω 不是越大越快。对这类矩阵把 ω 调到实验覆盖的最小值 1.1 也才 11 次并没有压过 Gauss-Seidel 的 9 次。SOR 的「加速」红利需要矩阵结构支撑盲目超松弛只会收获一个天文数字的迭代次数。2.3 初始向量只改变圈数不改变收敛性报告做了三组初始向量的对照X[0;0;0]、X[1;1;1]、X[−1;−1;−1]。横着读这张表Gauss-Seidel 三行全是 9 次最稳Jacobi 稳定在 13 到 14 次SOR 随 ω 和初始向量明显摆动。初始向量JacobiGauss-SeidelSOR ω1.25SOR ω1.50SOR ω1.75[0;0;0]149152180[1;1;1]139142079[-1;-1;-1]149152187注意一个细节报告数据表里 [0;0;0] 那组 SOR 记录的是 15/21/80结果分析里 [1;1;1] 组写的是 14/20/79同一条件出现了两个版本这是抄录问题我在第 4 章单独讲。回到迭代本身迭代是否收敛由迭代矩阵谱半径决定初始向量只影响前几轮的误差分布所以改初值不改变「收敛或不收敛」的结论。差别体现在圈数上[1;1;1] 离真解 [2;4;3] 相对近第一轮误差初值小提前一两轮达标很正常。写结论时不要写成「初始向量决定收敛性」要说「在收敛前提下初始向量影响迭代次数」。3. MATLAB 代码逐行拆解从矩阵输入到迭代计数的完整流程这份实验报告把完整源码写在「实验内容及步骤」一节里MATLAB 就是这门课指定的开发语言。我按顺序拆标注清楚再抄避免把时间浪费在符号错误上。3.1 Jacobi 与 Gauss-Seidel对角/下三角/上三角矩阵怎么拆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); % 上三角矩阵取负 X_0 [0;0;0]; % 初始向量 X_1 inv(D)*(LU)*X_0 inv(D)*b; % 第一次迭代 dis abs(X_1 - X_0); % 相邻两次解的绝对差 p 1; for i 1:3 while dis(i) 10^(-5) % 检查第 i 个分量的变化量 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第一行按行输入矩阵分号分隔行、逗号分隔列b 用分号写成列向量这两行定义了整个实验唯一的算例真解是 [2;4;3]。拆矩阵的三行和 2.1 节完全对应注释里把取负标出来了抄代码时最容易漏的就是这个负号。X_0 是初值X_1 先算一步dis 存的是相邻两轮解向量的分量差。后面 for i1:3 嵌套 while dis(i)1e-5 的意思是逐个分量检查只要第 i 个分量的变化还没小于 1e-5就继续迭代。p 从 1 开始累加最后就是报告里的迭代次数。这里有个关键点while 检查的是「相邻两次迭代解的变化量」不是「与真解的误差」。真解 [2;4;3] 是已知的但代码故意不拿它做停机条件因为一般方程组根本没有真解可用。用相邻差代替真误差是迭代法的标准做法实验里算出的解和真解在小数点后四位完全一致靠的就是把相邻差压到 1e-5 以下。Gauss-Seidel 的代码结构一模一样只有迭代式两处不同X_0 [0;0;0]; X_1 inv(D-L)*U*X_0 inv(D-L)*b; p 1; for i 1:3 while dis(i) 10^(-5) X_1 inv(D-L)*U*X_0 inv(D-L)*b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end和 Jacobi 的唯一区别是 inv(D) 换成 inv(D−L)左边从 D 变成 D−L右边 U 前不再乘 D 的逆。矩阵拆法没变代码层面一次改动迭代次数从 14 降到 9这就是「用本轮刚算出的分量」换来的收敛加速。inv(D−L) 不用担心不可逆D−L 是下三角矩阵对角元全是 A 的主对角元非零可逆性有保证。3.2 SOR 迭代公式的 MATLAB 写法ω 怎么传参A [4,-1,1; 4,-8,1; -2,1,5]; b [7;-21;15]; 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; p 1; for i 1:3 while dis(i) 10^(-5) X_1 inv(D-k*L)*((1-k)*Dk*U)*X_0 k*inv(D-k*L)*b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end对照标准公式 (D−ωL)X_{k1} ((1−ω)DωU)X_k ωb代码把等式两边同时左乘 inv(D−kL)拆成两项inv(D-kL)((1-k)DkU)X_0 负责旧解部分kinv(D-kL)*b 负责被 ω 放大的常数项。变量名用 k 不用 omega是原作者的习惯跑实验时改这一行就行。改 ω 的方式是直接改 k。报告里 ω1.25、1.50、1.75 三组数据来自三次改 k 重跑1.10 到 1.90 的扫描则是每隔 0.1 改一次。手动改参数没问题但扫 9 个值就要改 9 次我在第 6 章会把它封装成函数到时候一行循环跑完整张表。3.3 误差判断与迭代次数统计while 逻辑里的小机关三段代码共用一个停机逻辑for i1:3 外层 while dis(i)1e-5 内层。这个写法的实际行为是先盯第 1 个分量把它迭代到变化小于 1e-5再盯第 2 个分量再盯第 3 个。p 是三个 while 累加的总圈数。隐患在于如果第 1 个分量早就稳了第 2 个还在抖代码会继续跑p 里有一部分圈数是白跑的。这个实验矩阵规模小、收敛快影响不大规模一上去迭代次数就会虚高。我一般会换成全局判据一步到位tol 1e-5; X_0 zeros(3,1); X_1 X_0; p 0; while true X_1 inv(D)*(LU)*X_0 inv(D)*b; % 换成对应方法的迭代式 p p 1; if max(abs(X_1 - X_0)) tol break; end X_0 X_1; endmax(abs(X_1−X_0)) 取整个解向量变化量的最大值任何一个分量没稳住都继续迭代全部稳定才停机。这样 p 的含义是严格的「全体分量都达标所需轮数」三种方法比较迭代次数才有可比性。注意这里的循环体和第 3.1 节不同它把 X_0 X_1 放在了判断之后避免最后一次迭代白记一次。4. 避坑指南四个把迭代结果带偏的常见问题这些坑有的是我第一遍跑代码就撞上的有的是对照数据表和结果分析时发现的。每一条按现象、原因、解决给全照着排查能省半天。4.1 收敛判据写成单分量判断迭代次数虚高现象同一套逻辑换初始向量后迭代次数忽高忽低和教材参考值对不上换成全局判据重跑数值发生了变化。原因for i1:3 加 while 的组合是按分量逐个判断前面已经收敛的分量不会让循环提前退出p 把多余的圈数全算进去了。这个实验里矩阵小误差不大但原理上就是不严谨。解决换成 max(abs(X_1−X_0))tol 的全局判据第 3.3 节的写法直接可用报告里所有迭代次数以全局判据重跑结果为准。4.2 循环里反复算 inv(D)矩阵阶数一高就卡死现象扫描 ω 的 9 个值要等很久把同一个脚本丢给 50 阶矩阵跑了十几分钟没出结果。原因inv(D)、inv(D−L) 是显式求逆代价 O(n³)还放在 while 里每圈都算。对角矩阵求逆本应是对角元取倒数D−L 的求逆也永远不该写成 inv。解决迭代前把矩阵分解结果算一次循环里只做矩阵乘法和加法MJ D \ (LU); % 左除等价于 inv(D)*(LU)但数值更稳 cb D \ b; while true X_1 MJ*X_0 cb; p p 1; if max(abs(X_1 - X_0)) tol break; end X_0 X_1; endMATLAB 官方也不推荐为解方程调 inv左除 \ 是更稳的做法并且只在循环外算一次。4.3 ω 越界或取到 2 附近结果直接爆炸现象ω1.90 时迭代次数是 11828把 ω 改成 2.0代码要么不收敛要么出现 NaN死循环成了常态。原因SOR 理论收敛要求 0ω2ω2 时迭代矩阵有特征值到达边界谱半径压不到 1 以下ω1.9 时谱半径已经很接近 1每轮误差衰减极慢11828 次就是这么来的。解决扫描 ω 前先算谱半径 rho max(abs(eig(M\N)))rho≥1 直接跳过扫描步长避开端点例如 1.05:0.05:1.95不要踩到 2.0。4.4 数据抄录不一致同一条件出现两个迭代次数现象ω1.25、初始向量 [1;1;1] 在报告数据表里是 15 次结果分析里是 14 次ω1.75 也有 80 和 79 两个版本。原因这类实验报告通常是先跑代码再填表格某一行重复跑一遍、或把不同初始向量的结果串行就出现了两个数。数据表和结论不是同一次运行的输出。解决以复现代码的输出为准写报告时把容差、ω 扫描步长、初始向量写进参数说明结论里加一句「以下数据基于 tol1e-5 的全局判据重跑」,让每个数字都可回溯。5. 实验数据怎么读从迭代次数表到站得住的结论数据和结论之间隔着一层解读。这份报告的数据表信息密度很高但很多人只抄数字不读数字写出的结论经不起追问。5.1 数据表里其实有三组对照方法、ω、初始向量把三组初始向量的结果并排放能看出三层规律。第一层方法对照Gauss-Seidel 三行全是 9 次最稳Jacobi 稳定在 13 到 14 次SOR 的值随 ω 和初始向量明显摆动。第二层ω 对照同一初始向量下ω 从 1.25 涨到 1.75迭代次数从十几涨到八十多单调变差。第三层初始向量对照[1;1;1] 离真解 [2;4;3] 最近整体少一两轮[−1;−1;−1] 离得最远SOR 在 ω1.75 时到了 87 次。这张表最值得注意的发现是即使 SOR 取到实验范围内的最优值也没能压过 Gauss-Seidel。也就是说这个算例里超松弛几乎不带来收益。这个结论比「SOR 收敛速度快」更接近真实数据写进报告会更有分量。5.2 把「收敛」说成「对该方程组收敛」结论的适用范围报告的结论「对于该方程组随着松弛因子 ω 取值增大迭代次数增加」带两个限定词该方程组、本次实验条件。它是现象描述没有上升到方法优劣。想把它说得更精确可以拆成三层第一层本算例中 Gauss-Seidel 的 9 次优于 Jacobi 的 14 次第二层本算例中 SOR 最优 ω 估算值约 1.06实测 ω1.10 时迭代最少11 次但仍然多于 Gauss-Seidel第三层结论只对这类严格对角占优、规模小、非对称的矩阵成立外推到别的方程组之前必须先算谱半径。写实验报告最常见的毛病是把话说满比如「SOR 收敛速度优于 Gauss-Seidel」,这份数据就不支持。把范围收缩到「在给定 A、b、容差和初始向量的条件下」结论才是硬的。评审看到这样的限定表述反而会更信任你的分析。5.3 报告六要素对应一条可复现链路实验目的对应你要回答的问题实验原理对应三种迭代公式和收敛条件内容及步骤对应可运行的 MATLAB 代码数据记录对应迭代次数与解向量表格结果与分析对应从表格推出的限定性结论。把这六要素串起来时把 tol、ω 扫描范围、初始向量统一写进参数说明评审和同学都能直接复现。这份文档资料的落地价值就在这里它不只是一段公式推导而是一条从矩阵输入到结论的完整链路每个环节都有数字和代码撑着。6. 把三份脚本改造成通用迭代求解器谱半径预判与参数扫描报告里的三次实验是手动改 k 重跑最后我把它封装成函数顺带输出谱半径。这样一来ω 扫描是一行循环跑之前就知道某个参数组合会不会收敛。function [x, iter, rho] iterSolve(A, b, method, omega, x0, tol, maxIter) % method: jacobi / gs / sor % omega: 仅 SOR 使用其他方法传 1 即可 n length(b); D diag(diag(A)); L -tril(A, -1); U -triu(A, 1); x x0(:); tol 1e-5; % 未传容差时用默认值 maxIter 5000; % 防止不收敛时死循环 switch lower(method) case jacobi M D; N L U; rhs b; case gs M D - L; N U; rhs b; case sor M D - omega*L; N (1-omega)*D omega*U; rhs omega*b; otherwise error(未知方法名); end rho max(abs(eig(M \ N))); % 迭代矩阵谱半径rho1 时多半不收敛 iter -1; for k 1:maxIter x_new M \ (N*x rhs); if max(abs(x_new - x)) tol x x_new; iter k; break; end x x_new; end end这个函数把三种方法的差异收敛到 switch 分支里Jacobi、Gauss-Seidel、SOR 的区别只剩 M、N、rhs 三个变量。rho 用 eig(M\N) 求迭代矩阵谱半径是预判收敛性的黑匣子开关。调用方式很简单A [4,-1,1; 4,-8,1; -2,1,5]; b [7;-21;15]; [x, it, rho] iterSolve(A, b, sor, 1.1, zeros(3,1), 1e-5, 5000); fprintf(x [%.4f; %.4f; %.4f], iter %d, rho %.4f\n, x, it, rho);复现报告里的 ω 扫描一行循环就够了for w 1.1:0.1:1.9 [x, it, rho] iterSolve(A, b, sor, w, zeros(3,1), 1e-5, 5000); fprintf(omega%.2f iter%d rho%.4f\n, w, it, rho); end跑出来的迭代次数和报告基本一致rho 随 ω 增大逐步逼近 1表格里那些几百上千的迭代次数在谱半径里早就写好了答案。从那以后我拿到任何新的系数矩阵都先跑一遍 eig(M\N) 看谱半径再决定用哪个方法、ω 从哪扫起不再盲目堆迭代次数。这份报告里那个 11828 次的教训比任何教材都直白松弛因子省不掉但可以提前算。希望帮到你。本文还有配套的精品资源点击获取