简介这份资源是东北大学理学院数学系数值分析课程的实验报告面向正在学习数值分析、需要完成迭代法实验或复习线性方程组求解的学生。报告围绕解线性方程组的迭代法与迭代格式比较展开完整记录了Jacobi、Gauss-Seidel与SOR三种方法的实验目的、原理、内容步骤、数据记录及结果分析并附有可直接运行的MATLAB代码涵盖矩阵分解、迭代控制与误差计算等关键环节。压缩包共1个doc文件约186KB按学校规范排版包含实验封皮、指导教师签字栏与统一版头便于直接参考或套用格式。资源已有2617人学习下载读者可从中获得三种迭代法的公式推导与实现细节、不同松弛因子ω对收敛速度影响的对比数据以及初始向量变化时的迭代次数记录适合作为课程实验报告撰写与算法验证的参考材料。1. 数值分析实验包三套迭代法 MATLAB 代码与完整报告一次拿到有个场景我印象很深期末前一周某实验室的学弟抱着一台装了 MATLAB 的笔记本过来说他能背出 Jacobi 的公式但一让他把A拆成D、L、U再写成循环就卡在inv(D)和误差判断上跑出来的迭代次数和老师给的参考答案对不上。这类问题不是数学没学懂而是缺一份能直接跑、能对着改的工程化样例。这份数值分析实验报告资源正好补这个缺口它把解线性方程组的三种迭代法——Jacobi、Gauss-Seidel、SOR——从实验目的、原理、数据记录到 MATLAB 实现完整串起来还附了不同松弛因子 ω 下的迭代次数对比。适合正在做数值分析课程实验的学生也适合想快速复习迭代法收敛行为、需要一份可复现 MATLAB 脚本的开发者。下面按“能干什么、怎么跑、怎么改”拆开讲。2. 实验包结构拆解报告模板、三种迭代法与数据记录怎么组织这份资源不是单纯的代码压缩包而是“报告文档 代码片段 实验数据”三件套。理解它的组织方式比直接复制代码更重要因为后面你要改参数、补实验、换方程组都得知道每一块对应哪里。2.1 报告文档的固定骨架从正文看报告遵循的是一套学院统一的实验报告规范结构大致是实验项目名称、实验目的和要求、实验原理、实验内容及步骤、实验数据记录和处理、实验结果与分析。这套骨架看起来是形式要求实际对写代码的人有约束作用——它逼你把“为什么这么迭代”和“迭代出来什么结果”分开写清楚。对应到技术内容上可以这样映射报告章节对应的技术内容对开发者的用处实验目的认识迭代法收敛含义、影响收敛速度的因素明确要观察的指标迭代次数、误差实验原理Jacobi、Gauss-Seidel、SOR 三种迭代格式代码里的矩阵拆分依据实验内容及步骤三段 MATLAB 代码可直接运行的核心资产实验数据记录解向量、迭代次数、误差表验证自己跑的结果对不对实验结果与分析不同 ω 下迭代次数变化调参和写结论的参考这张表建议先存下来。后面你如果要把报告改成自己的课程作业按这张表逐项替换内容比从空白文档硬写快得多。2.2 三种迭代法的矩阵拆分逻辑三份代码共用同一个拆分套路把系数矩阵A拆成对角阵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); % 上三角矩阵取主对角线以上并取负这里有个容易看漏的细节L和U前面都带了负号。原因是迭代格式通常写成A D - L - U的形式把移项后的负号提前吸收进L、U后面公式里就能直接写D*X (LU)*X b这类形式。如果你自己重写时忘了这个负号迭代会直接发散或者收敛到错误的值这是血泪经验里最常见的一类翻车。2.3 数据记录表透露的对比维度报告里的数据表记录了三个维度的信息解向量X1 X2 X3、迭代次数、误差阈值。精确解是(2, 4, 3)误差阈值统一取10^(-5)。这张表的价值在于它给了你一个“标准答案”你自己跑代码时如果 Jacobi 迭代次数不是 14 次左右、解不是逼近(2,4,3)说明代码或矩阵拆分出了问题。数据表还记录了 SOR 在不同 ω 下的表现比如 ω1.25 时 15 次、ω1.50 时 21 次、ω1.75 时 80 次。这组数据不是随便填的它对应了后面要讲的“松弛因子不是越大越好”这个结论。你复现时应该能跑出接近的数量级数量级对不上就要回头查循环终止条件。2.4 运行环境与文件准备这份资源的核心是 MATLAB 脚本运行门槛不高但有几个前置条件要确认。第一MATLAB 版本。代码只用了diag、tril、triu、inv、abs这些基础函数从 R2010a 到较新版本都能跑不依赖工具箱。如果你手头是 Octave也能跑但inv在病态矩阵上的数值行为可能略有差异实验规模下基本无感。第二文件组织。建议按实验建目录一个实验一个文件夹里面放脚本和数据numerical-lab/ ├── exp01_iterative/ │ ├── jacobi.m │ ├── gauss_seidel.m │ ├── sor.m │ └── result.md └── report/ └── exp01_report.docx这样组织的好处是后面做第二个实验时不会和第一个混在一起。报告文档单独放report目录代码和结果放实验目录改代码时不会误动报告。第三运行方式。每段代码可以单独存成.m文件在命令行运行也可以直接粘到 MATLAB 的编辑器里按 F5。建议存文件因为后面你要改 ω、改初始向量、改误差阈值存文件方便反复调。3. 三段核心代码逐行讲Jacobi、Gauss-Seidel、SOR 怎么跑起来这一章是整份资源最值钱的部分。三段代码结构相似但公式不同逐段拆开讲重点放在“这行在干什么、参数怎么改、改错会怎样”。3.1 Jacobi 迭代先算全量再更新Jacobi 的特点是每一轮迭代都用上一轮的完整X_0算新的X_1不中途更新分量。代码主体如下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)) 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、L、U的构造前面讲过注意负号。X_0是初始猜测实验里用了[0;0;0]后面分析部分还试了[1;1;1]和[-1;-1;-1]说明初始向量会影响迭代次数但不影响最终收敛解。核心迭代式inv(D)*(LU)*X_0 inv(D)*b就是 Jacobi 公式X D^{-1}(LU)X D^{-1}b的直译。误差判断用的是abs(X_1 - X_0)逐分量比较阈值10^(-5)。这里有个设计上的小问题for i 1:3外层循环配合while(dis(i) 10^(-5))实际效果是三个分量依次等到各自收敛但因为X_0在循环里持续更新后一个分量的while开始时前面的分量已经收敛整体能跑通只是迭代计数p的语义会变得有点模糊——它统计的是三个while累计的迭代步数不是严格意义上的“轮数”。复现时如果你发现p比预期大先检查这里。参数怎么改10^(-5)是精度改小到10^(-6)迭代次数会增加X_0换成别的初值迭代次数会变但解不变A、b换成你的方程组要保证D可逆也就是主对角线没有零元素。3.2 Gauss-Seidel 迭代用上新分量Gauss-Seidel 和 Jacobi 的差别在于算第i个分量时前面已经算出的新分量直接用上不等整轮结束。代码里体现为把inv(D)换成了inv(D-L)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-L)*U*X_0 inv(D-L)*b; % 注意是 D-L 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-L)*b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end关键改动只有一处inv(D)变成inv(D-L)。数学上这是因为 Gauss-Seidel 把D - L移到了等式左边迭代式变成X (D-L)^{-1}UX (D-L)^{-1}b。代码里L已经取过负所以写D-L正好对应D 原下三角这个符号关系如果搞反inv(D-L)可能变成奇异矩阵直接报错。实验数据里 Gauss-Seidel 迭代 9 次比 Jacobi 的 14 次少符合“用上新分量收敛更快”的预期。你复现时如果 Gauss-Seidel 反而比 Jacobi 慢大概率是L的负号处理错了或者while的终止条件写成了整轮判断。3.3 SOR 迭代松弛因子 ω 是唯一变量SOR 在 Gauss-Seidel 基础上引入松弛因子k代码里用k表示 ω迭代式变成X (D-kL)^{-1}((1-k)D kU)X k(D-kL)^{-1}b。代码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)*D k*U)*X_0 k*inv(D-k*L)*b; dis abs(X_1 - X_0); p 1; for i 1:3 while(dis(i) 10^(-5)) X_1 inv(D-k*L)*((1-k)*D k*U)*X_0 k*inv(D-k*L)*b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end endk是唯一需要调的参数。k1时 SOR 退化成 Gauss-Seidel0k1叫低松弛1k2叫超松弛。实验数据很能说明问题ω 从 1.10 一路加到 1.90迭代次数从 11 次涨到 11828 次说明这个方程组在 ω 偏大时收敛急剧变慢。这不是说 SOR 不好而是说松弛因子要选在合适区间盲目调大反而翻车。参数说明k的取值范围理论上要满足0k2才可能收敛超出这个范围基本发散inv(D-k*L)要求D-k*L可逆主对角线元素不能因为k的取值变成零。改k时建议从 1.0 开始以 0.05 或 0.1 为步长扫一遍记录迭代次数画个曲线就能看出最优区间。3.4 三段代码的公共坑inv与循环终止三段代码都用了inv求逆。实验规模是 3×3inv完全够用但如果你把方程组扩到几十阶inv的数值稳定性和速度都会成为问题。常见做法是把inv(D)*b换成D\b把inv(D)*(LU)*X_0换成D\((LU)*X_0)用左除代替显式求逆。这是后面二次开发时最值得先改的一处。另一个公共点是for i 1:3加while的嵌套。这个写法能跑但迭代计数p的语义不干净。更清晰的写法是把整个迭代包在一个while里用norm(X_1-X_0, inf) 10^(-5)做整体终止条件p就是真正的迭代轮数。实验报告里记录的 14 次、9 次这些数字用整体终止条件复现会更稳定。4. 避坑与排查迭代次数对不上、结果发散、ω 调崩了怎么办这一章按“现象 → 原因 → 解决”列几条实际会遇到的坑。这些不是理论上的可能性是复现这份资源时大概率会撞上的。4.1 现象Jacobi 跑出来迭代次数远大于 14原因通常有两个。一是L、U的负号漏了导致迭代矩阵变成D^{-1}(-L-U)而不是D^{-1}(LU)收敛速度完全变样。二是while的终止条件写成了dis(i) 10^(-5)但dis没有在循环内更新死循环或者提前退出。解决先打印D、L、U三个矩阵和手算结果对一遍。以实验矩阵为例D应该是diag(4,-8,5)L应该是[0,0,0; -4,0,0; 2,-1,0]U应该是[0,1,-1; 0,0,-1; 0,0,0]。对不上就先修拆分再谈迭代。4.2 现象Gauss-Seidel 报错“矩阵接近奇异”原因inv(D-L)里的D-L行列式接近零。在实验矩阵上不会发生但如果你换了方程组主对角线元素很小或者L的负号处理错D-L可能奇异。解决先算det(D-L)接近零就检查A的主对角线是否占优。迭代法收敛的充分条件之一是严格对角占优换方程组时先验证这个条件不满足就别硬跑换方法或者做预处理。4.3 现象SOR 在某个 ω 下迭代次数爆炸原因ω 超出了该方程组的有效收敛区间。实验数据里 ω1.90 时 11828 次就是典型表现。SOR 的收敛区间和矩阵谱半径有关不是所有方程组都在0ω2内都好用。解决不要一上来就试 1.75。从 ω1.0 开始步长 0.1 扫到 1.9记录每个点的迭代次数找出最小值对应的区间再细化。实验数据已经给了参考这个方程组在 ω1.10 附近 11 次比 Gauss-Seidel 的 9 次还多说明最优 ω 可能小于 1值得往 0.9 方向再扫。4.4 现象换了初始向量后迭代次数变了怀疑代码错了原因初始向量影响迭代起点迭代次数变化是正常现象不是 bug。实验数据里[0;0;0]时 Jacobi 14 次[1;1;1]时 13 次[-1;-1;-1]时 14 次波动不大但确实存在。解决判断代码对不对看最终解是否收敛到(2,4,3)而不是看迭代次数是否完全一致。迭代次数受初值、终止条件、浮点误差共同影响允许有小幅差异。4.5 现象复制代码到 Octave 后diag(diag(A))行为不同原因极少数旧版本 Octave 对diag嵌套调用的返回维度处理和 MATLAB 有差异导致D变成向量而不是矩阵。解决拆成两步写d diag(A); D diag(d);兼容性最好。这类兼容问题不影响算法本身但会让新手误以为公式错了。5. 二次开发与改造把实验代码变成可复用的迭代法工具箱这份资源作为课程实验是完整的但作为可复用代码还有距离。下面几条改造思路按投入产出比排序第一条最值得先做。5.1 把三段脚本抽成函数现在三段代码都是脚本矩阵和参数写死在文件里。改成函数后换方程组、换精度、换 ω 都只改调用参数function [x, iter] jacobi_solve(A, b, x0, tol, maxiter) % JACOBI_SOLVE Jacobi 迭代解线性方程组 % A: 系数矩阵 b: 右端向量 x0: 初始向量 % tol: 误差阈值 maxiter: 最大迭代次数 D diag(diag(A)); L -tril(A, -1); U -triu(A, 1); x x0; for iter 1:maxiter x_new D \ ((L U) * x b); % 用左除代替 inv if norm(x_new - x, inf) tol x x_new; return; end x x_new; end warning(达到最大迭代次数未收敛); end这个函数做了三处改进用D \ (...)代替inv(D)*(...)数值更稳用norm(..., inf)做整体终止判断iter语义清晰加了maxiter防止死循环。Gauss-Seidel 和 SOR 可以照这个模板改SOR 多一个omega参数内部把D换成D - omega*L把(LU)换成(1-omega)*D omega*U。5.2 加一个 ω 扫描脚本实验数据里 ω 从 1.10 到 1.90 的迭代次数是手工记录的改造后可以自动扫omegas 0.8:0.05:1.9; iters zeros(size(omegas)); for idx 1:length(omegas) [~, iters(idx)] sor_solve(A, b, [0;0;0], 1e-5, 10000, omegas(idx)); end plot(omegas, iters, -o); xlabel(\omega); ylabel(迭代次数);跑出来的曲线能直接看出最优 ω 区间。实验数据提示这个方程组在 ω 小于 1.1 时可能更好扫描范围往 0.8 扩一下结论会比原报告更完整。5.3 验证方法用已知解反推改完代码怎么确认没改坏构造一个已知解的方程组取x_true [1;2;3]令b A * x_true然后跑迭代看是否收敛到x_true。这个方法比对着实验数据猜迭代次数可靠得多因为解是确定的迭代次数只是过程指标。5.4 一个具体技巧把收敛判据和报告数据对齐原报告的迭代次数是在特定终止条件下得到的。你改造后如果发现次数对不上先别怀疑算法检查终止条件原代码是逐分量abs判断改造后是norm(..., inf)两者在接近收敛时行为接近但不完全等价。要对齐报告数据就临时把终止条件改回逐分量判断跑一遍对比确认差异来源后再决定用哪种。我自己的习惯是每次改完迭代类代码先拿一个 2×2 的小方程组手算两三步和代码输出逐步对对上了再上 3×3、再上更大规模。这个习惯帮我省过很多次“公式看着对、跑出来就是不对”的排查时间。这份实验包的价值也在这里它给了一个规模足够小、结果可手算验证的起点你在这个起点上改函数、加扫描、换方程组每一步都有参照。希望帮到你。本文还有配套的精品资源点击获取
原理与实现详解)
