解线性方程组的三种迭代法原理与MATLAB实现全拆解

简介东北大学理学院数学系数值分析课程实验报告以解线性方程组的迭代法与迭代格式比较为主线系统介绍了Jacobi迭代法、Gauss-Seidel迭代法与SOR松弛迭代法。报告严格按学校规范撰写包含实验目的、原理、内容、数据记录和处理、结果分析等完整环节方便高校学生对照学习或完善课程实验报告。压缩包仅含1个doc文档大小186KB结构紧凑。文中给出了三种算法的MATLAB实现从矩阵分解、迭代公式到误差阈值控制均有对应代码并记录了不同初始向量与松弛因子ω下的迭代次数和收敛误差例如ω从1.10到1.90时迭代次数由11次增至11828次能直观展示松弛参数对收敛速度的影响。已有2614人学习下载可作为理解迭代法原理、查看算法实现细节及撰写实验报告的重要参考。1. 一份让你少熬夜的数值分析实验报告三种迭代法从原理到 MATLAB 代码全拆解如果你正在数值分析课上被“解线性方程组的迭代法”折腾大概率经历过这样的场景讲义上的公式一看就会自己写 MATLAB 一跑就死循环或者迭代次数忽大忽小根本不知道问题出在哪。这份实验报告正好是解药——它把 Jacobi 迭代法、Gauss-Seidel 迭代法和 SOR 松弛迭代法从原理、代码到实验数据完整整理了一遍连不同松弛因子下迭代次数从 11 次到 11828 次的对比都给你算好了。适合三类人正在做类似课程实验、需要复现数据写报告的同学准备面试或考研复试、想快速理解迭代法实现细节的人以及在实际工程里需要手写迭代求解器、想找个靠谱参考的开发者。这份资源不只是一堆代码更是一份能直接跑通、能照着改的完整实验样本。2. 实验设计解析线性方程组迭代法的收敛条件与矩阵分裂2.1 迭代法为什么可以解线性方程组直接法用高斯消元或 LU 分解一次性求出精确解迭代法则把一个Ax b的问题改写成固定点形式x Bx f然后从一个初始向量出发反复计算x_{k1} Bx_k f直到前后两步的差足够小。这里最关键的是收敛条件迭代矩阵B的谱半径必须小于 1也就是所有特征值的模都小于 1。如果谱半径大于等于 1迭代就会发散或者原地打转这也是很多同学实验报告里出现“死循环”的根本原因。这份实验报告选用的测试矩阵是三阶的A [4, -1, 1; 4, -8, 1; -2, 1, 5]; b [7; -21; 15];它不是严格对角占优矩阵第二行对角线是 -8但该行其他元素的绝对值之和是 5满足对角占优第一行对角线 4 比 2 大第三行 5 比 3 大实际上整体是严格对角占优的所以 Jacobi 和 Gauss-Seidel 都收敛。实验报告特意选了这样一个行列式不为零、迭代性质又不太极端的例子非常合适初学者观察不同方法的收敛速度差异。2.2 矩阵分解D、L、U 在 MATLAB 里怎么取迭代法把矩阵A拆成三部分对角部分D、严格下三角部分L和严格上三角部分U。注意这里L和U的方向取决于你自己的定义报告里的定义是D diag(diag(A)); % 对角矩阵 L -tril(A, -1); % 下三角矩阵取负 U -triu(A, 1); % 上三角矩阵取负为什么取负因为后续迭代公式里需要的是A D - L - U。tril(A, -1)返回 A 中主对角线以下的部分但不含主对角线triu(A, 1)返回主对角线以上的部分。如果你直接用tril(A, -1)得到的是一个“带正号的下三角”那么A D L U就会写错正负号。这个细节是很多实现翻车的起点后面避坑章节还会再提到。以这个测试矩阵为例D是对角线[4, -8, 5]L是[0,0,0; -4,0,0; 2,-1,0]的相反数U是[0,1,-1; 0,0,-1; 0,0,0]的相反数。写出迭代公式后就明白为什么 MATLAB 代码里到处是inv(D)、inv(D-L)。2.3 三种迭代格式的矩阵表达与收敛性Jacobi 迭代公式的核心是“用旧值计算新值”X_1 inv(D) * (L U) * X_0 inv(D) * b;把A D - L - U代入Ax b得到Dx (LU)x b所以迭代矩阵是inv(D)*(LU)。每次迭代都基于上一轮的整个向量X_0所以各个分量可以并行计算。Gauss-Seidel 迭代的关键是“用最新值计算下一个分量”X_1 inv(D - L) * U * X_0 inv(D - L) * b;这里把下三角部分并入左侧求逆等价于在计算第 i 个分量时已经使用本轮更新过的前 i-1 个分量。Gauss-Seidel 通常比 Jacobi 收敛快实验数据也印证了初始向量为零时Jacobi 需要 14 次Gauss-Seidel 只要 9 次。SOR 是 Gauss-Seidel 的加速版公式里多了一个松弛因子ω代码是k 1.75; % 松弛因子 omega X_1 inv(D - k*L) * ((1-k)*D k*U) * X_0 k * inv(D - k*L) * b;当ω1时 SOR 退化成 Gauss-Seidelω1叫超松弛通常用于加速收敛ω1叫低松弛有时用于改善稳定性。理论要求0 ω 2超出这个区间大概率发散。2.4 SOR 松弛因子的“玄学”取值不同迭代次数从 11 次变成 11828 次实验报告里最有意思的部分是 SOR 松弛因子扫描。初始向量为零向量时ω 值迭代次数1.10111.20141.30161.40191.50211.60321.70551.801531.9011828这个数据我第一次看到时也吓了一跳。1.9 和 2.0 只差 0.1迭代次数却从一百多次暴涨到上万次。原因在于 SOR 的最优松弛因子和系数矩阵的特征值分布强相关当ω接近上界 2 时迭代矩阵的谱半径会无限接近 1收敛速度变得极慢。这不是代码 bug而是数值方法本身的特性。这也是为什么实验报告的数据比单纯跑一次代码有价值——它帮你提前划出了“危险区间”。3. 核心代码拆解一个迭代循环从初始化到误差收敛的完整链路3.1 原始代码逐行解读矩阵、初值、迭代步报告的 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); % 上三角矩阵取负 X_0 [0;0;0]; % 初始向量 X_1 inv(D) * (L U) * 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) * (L U) * X_0 inv(D) * b; dis abs(X_1 - X_0); X_0 X_1; p p 1; end end这里p是迭代次数计数器。初始迭代执行了一次所以p从 1 开始。for i 1:3加while的组合很有意思它并不是一次性判断所有分量是否满足误差而是分别检查每个分量。我们稍后会说这种写法的隐患但先理解执行流程当第 1 个分量误差大于阈值时不断迭代直到它小于阈值然后跳出 while进入 i2再让第 2 个分量满足最后 i3。因为每轮迭代都会更新整个向量所以这个流程实际上会让所有分量最终都满足但可能比统一的范数判断多算一些迭代。3.2 误差判断的写法按分量判断和按范数判断的差别实验报告里的误差判断是while(dis(i) 10^(-5))只观察某个分量的绝对误差。常见做法是计算整个向量的无穷范数或 2 范数while max(abs(X_1 - X_0)) 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两种写法在大多数情况下结果一致但按分量判断存在一个潜在问题如果某个分量在迭代中波动比如先收敛又变坏外层 for 循环已经越过了它就不会再检查可能最终这个分量的误差其实略超阈值。对这个三阶测试矩阵而言不会出问题但换成更复杂的矩阵时就可能“表面收敛、实际没收敛”。所以我通常建议改用max(abs(...))或norm(...)作为退出条件。3.3 实验数据三种方法在不同初始向量下的迭代次数这份实验报告记录了三组初始向量的结果这部分数据可以直接用到你自己的报告里初始向量[0;0;0]Jacobi 14 次Gauss-Seidel 9 次SORω1.10 时 11 次ω1.50 时 21 次ω1.75 时 80 次初始向量[1;1;1]Jacobi 13 次Gauss-Seidel 9 次SORω1.25 时 14 次ω1.50 时 20 次ω1.75 时 79 次初始向量[-1;-1;-1]Jacobi 14 次Gauss-Seidel 9 次SORω1.25 时 15 次ω1.50 时 21 次ω1.75 时 87 次可以看出Gauss-Seidel 对初始向量不敏感始终 9 次收敛Jacobi 基本稳定在 13~14 次SOR 则对 ω 和初值都很敏感。这个数据告诉我们一个实际经验如果你的问题是线性的Gauss-Seidel 是性价比很高的选择SOR 的加速效果需要谨慎调参调不好反而更慢。3.4 一个隐藏问题inv 求逆的效率与数值稳定性报告代码里反复出现inv(D)、inv(D-L)、inv(D-k*L)。在小规模三阶矩阵下毫无问题但如果你把这个代码改造成大规模稀疏矩阵求解器inv会先把矩阵变成稠密矩阵内存和计算量瞬间爆炸。更稳的做法是改用左除运算符X_1 (D - L) \ (U * X_0) (D - L) \ b;\对稀疏矩阵会走 LU 分解或迭代求解数值稳定性也更好。这是从“能跑”到“能用在工程里”的关键一步二次开发章节会继续展开。4. 运行与复现在 MATLAB 里把三种解法跑通并画出收敛曲线4.1 环境准备MATLAB 或 Octave脚本文件怎么建最简单的方式是打开 MATLAB新建一个脚本命名为iterative_demo.m把三种方法的代码分三段贴进去然后点运行。如果没有 MATLAB 授权用 GNU Octave 也能跑通两者语法在这份报告里几乎完全兼容。需要注意的一点是文件名不要用中文也不要和内置函数重名比如别叫sum.m或plot.m。复制代码时最容易踩的坑是MATLAB 的矩阵换行用分号但如果你从网页复制可能把全角分号或空格粘进去。建议先检查矩阵行的分隔符是不是半角的;或换行。4.2 跑通 Jacobi 与 Gauss-Seidel注意看迭代次数 p 的累加逻辑运行 Jacobi 部分后命令窗口不会直接显示迭代次数需要自己加一行输出disp([Jacobi 迭代次数: , num2str(p)]);我在复现时习惯把p的初值改成 0然后让循环里的p p 1统一计数。因为原报告在 while 之前已经做了一次迭代p初始为 1这样最终得到的数字才和实验数据对得上。Gauss-Seidel 部分的思路完全相同只是迭代矩阵变成了inv(D-L)*U。跑完你应该看到p14和p9。4.3 批量测试不同松弛因子用 for 循环复现 SOR 实验数据实验报告里的 ω 扫描是手写一个个改的但复现时更推荐用循环。下面这段代码可以直接复现表里的数据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); omegas [1.10,1.20,1.30,1.40,1.50,1.60,1.70,1.80,1.90]; for k omegas 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; while(dis(1) 1e-5 || dis(2) 1e-5 || dis(3) 1e-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 fprintf(omega%.2f, 迭代次数%d\n, k, p); end注意这里我把 while 条件改成“任一分量不满足就继续”避免原始代码里for while结构可能带来的漏检问题。输出结果会和实验报告基本一致ω1.10 时 11 次ω1.90 时 11828 次。跑最后一条时要有耐心可能等几秒。4.4 修改初始向量和误差阈值观察收敛速度的变化把初始向量改成[1;1;1]只需要改X_0那一行。改成[-1;-1;-1]同理。误差阈值10^(-5)也可以改成10^(-8)你会看到迭代次数普遍增加。如果改得太小比如10^(-15)三阶矩阵用 double 浮点数就可能无法达到循环会进入无限震荡。这个现象可以当作实验报告里的“延伸思考”来写。如果想要更直观可以加一个误差记录数组把每次迭代的max(abs(X_1-X_0))存下来最后用semilogy画收敛曲线。这个技巧放到第六章。5. 避坑指南迭代法实验里最容易翻车的 5 个细节5.1 现象while 循环直接卡死跑几分钟都不停原因迭代矩阵谱半径大于等于 1或者误差阈值设成了不可达的数。最常见的是代码里L和U的正负号写反导致迭代矩阵根本不是预期那个。解决先打印D、L、U看符号再计算迭代矩阵的特征值比如max(abs(eig(inv(D)*(LU))))如果谱半径大于等于 1就不要再跑了。这是我做迭代法实验最深的血泪教训——不查谱半径就盲目跑循环等于让电脑陪你死等。5.2 现象Gauss-Seidel 结果和 Jacobi 几乎一样原因矩阵分解时没有取负导致D-L和Dtril(A,-1)混用。很多同学抄代码时把L -tril(A,-1)写成L tril(A,-1)那么D-L实际上变成了D 下三角根本不符合推导。解决统一用A D - L - U的约定打印A - (D - L - U)看是否为零矩阵。这个检查成本极低能避免一整晚的困惑。5.3 现象命令窗口报错“Matrix must be square”或“奇异矩阵”原因inv(D-k*L)在某个松弛因子下矩阵不可逆。理论上只要 A 非奇异D-kL通常也可逆但如果你把 ω 设成负数或大于 2这个矩阵可能变得病态。解决在循环里用 try-catch 捕获奇异矩阵错误或者改用左除\左除对近奇异矩阵会给出警告但至少能返回一个结果。另外检查一下是不是把矩阵的一行数据输错了比如[4,-1,1;4,-8,1;-2,1,5]中间的分号是否被替换成了逗号。5.4 现象松弛因子 1.9 时迭代次数跳到 11828 次这是 bug 吗不是 bug是 SOR 方法的固有现象。当 ω 接近上界 2 时迭代矩阵谱半径趋近 1收敛速度趋近于零。遇到这种情况不要慌张把它写进实验分析里反而是加分项。解决如果你想看到完整的收敛过程可以把误差阈值放大到10^(-3)迭代次数会显著减少也能观察到误差缓慢下降的曲线。如果只是验证代码正确性建议只测到 1.8。5.5 现象用 Octave 运行结果和 MATLAB 不一样原因Octave 的浮点数运算和 MATLAB 基本一致但某些内置函数如inv、eig底层实现有细微差异在接近收敛边界时可能表现出不同迭代次数。另外循环里的变量作用域和显示精度也可能不同。解决关键结果不要直接用比较而是用abs(p1 - p2) 1容忍。如果迭代次数差很多检查是不是用了不同的数据文件或没清空工作区。6. 二次开发把实验代码改造成通用的迭代法求解器6.1 封装成函数支持 method 参数和松弛因子原报告的代码是平铺的每换一个方法就要复制改一遍。更好的做法是封装成一个函数function [x, iter, err_history] iterative_solve(A, b, method, omega, tol, x0) % method: jacobi, gs, sor % omega: 仅 SOR 使用 % tol: 迭代终止阈值 % x0: 初始向量默认零向量 D diag(diag(A)); L -tril(A, -1); U -triu(A, 1); n size(A, 1); x x0; err_history []; iter 0; while true switch method case jacobi x_new D \ ((L U) * x) D \ b; case gs x_new (D - L) \ (U * x) (D - L) \ b; case sor x_new (D - omega*L) \ (((1-omega)*D omega*U) * x) ... omega * ((D - omega*L) \ b); otherwise error(未知方法); end err max(abs(x_new - x)); err_history [err_history; err]; x x_new; iter iter 1; if err tol || iter 10000 break; end end end这个函数把三种方法统一到了一套switch结构里最大的好处是测试不同方法时不需要反复改代码。我把原来的 while 条件改成了while true靠迭代上限和误差阈值双重退出避免死循环。6.2 增加误差记录与收敛曲线在函数里我把每次迭代的误差存进err_history这样调用后可以直接画收敛曲线[x, iter, his] iterative_solve(A, b, sor, 1.1, 1e-5, [0;0;0]); semilogy(1:length(his), his, -o); xlabel(迭代次数); ylabel(误差对数坐标);用semilogy画误差下降曲线能直观看到 SOR 在不同 ω 下收敛速度的差异。你会发现 ω1.9 的曲线下降得非常平缓这就是 11828 次的来源。6.3 验证方法用严格对角占优矩阵和谱半径判断改造完之后必须验证函数没有跑偏。我一般会做两步第一步构造一个严格对角占优矩阵比如三对角矩阵用直接法A\b的结果和迭代解对比第二步计算迭代矩阵的谱半径确认小于 1。另外再跑一遍原始实验数据确保 Jacobi 返回 14 次、Gauss-Seidel 返回 9 次。验证通过后再去调 ω就心中有数了。从那以后我每次做迭代法实验都会先写一个iterative_solve这样的通用函数再批量跑参数而不是复制粘贴三次代码。这个习惯帮我少走了很多弯路希望你也能直接用上。希望帮到你。本文还有配套的精品资源点击获取