简介:面向东北大学理学院数学系数值分析课程,这份实验报告以解线性方程组的迭代法为主题,完整呈现Jacobi迭代、Gauss-Seidel迭代与SOR(松弛)迭代三种方法的原理推导、MATLAB代码实现和对比分析。报告严格按学校实验报告规范撰写,涵盖实验项目名称、目的、原理、内容、数据记录与处理、实验结果分析等要素,并预留指导教师签字栏;内容聚焦不同初值与不同松弛因子ω下的迭代次数变化,可清晰看出三种方法的收敛速度差异。资源共1个doc文档,大小约186KB,包含可直接运行的MATLAB代码、完整数据记录与结果分析,适合需要提交课程实验报告或复习迭代法收敛性的本科生直接参考使用。已有2614人学习浏览,能帮助读者快速掌握迭代格式的实现细节和收敛特性,节省自行整理模板与调试代码的时间。
1. 数值分析实验报告:线性方程组迭代法其实没你想的那么玄
提起数值分析这门课,不少同学的第一反应都是"一堆公式推导、算到天亮也算不完"。但你看完这份来自某高校数学实验中心的《数值分析实验报告》就会发现,这门课真正落地的时候,玩的就是三件事:Jacobi 迭代、Gauss-Seidel 迭代和 SOR 松弛迭代,用 MATLAB 几十行代码就能在同一个方程组上跑出收敛速度的天壤之别。整套实验报告以"解线性方程组的迭代法与迭代格式的比较"为核心,把实验目的、原理、代码、数据记录、结果分析完整串成一条线,特别适合正在修数值分析、需要交实验报告或者想搞懂迭代法收敛本质的本科生和刚接触科学计算的开发者。它不只是一份作业模板,更是一个可以照着敲、改了参数就能复现全部实验数据的 MATLAB 交互式教学代码。看完这份报告,你能直观感受到同一个方程组在不同迭代格式下的收敛次数能差出一个量级,更能理解为什么搞数值计算的人都格外看重松弛因子这个参数。
2. 算法原理与代码拆解:三种迭代格式在 MATLAB 里的真面目
2.1 系数矩阵的 D、L、U 分解是整个实验的地基
这份实验报告处理的核心方程组是:
A = [4, -1, 1; 4, -8, 1; -2, 1, 5],b = [7; -21; 15]
三个解法的起点都是把系数矩阵拆成三个部分:D 是对角阵,-L 是严格下三角,-U 是严格上三角。用 MATLAB 实现时,代码里用了一个非常巧妙的手法——先取负号再分解,这样后面写迭代公式时直接往公式里套,正负号不会乱。
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]; % 初始迭代向量
这段代码里最值得讲的是tril(A,-1)这种带偏移量的写法:tril(A, -1)取的是主对角线以下、且不包括主对角线的部分,正好对应严格下三角矩阵;同理triu(A, 1)取严格上三角。之所以要加负号,是因为迭代公式的标准形式里写的是X = D^(-1) * (L + U) * X + D^(-1) * b,公式中的 L、U 就是负的严格三角部分。如果不用-tril而是直接tril,后面公式里就得手动填负号,代码一长非常容易翻车。
2.2 Jacobi 迭代的核心:同时用旧值更新所有分量
Jacobi 迭代的思想是"用上一轮的旧值,一次性算出这一轮的所有新值"。代码里的迭代核心这一句很关键:
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
这段代码有个非常典型的写法:for i=1:3套了三个独立的while循环,分别检查 X_1、X_2、X_3 三个分量各自是否满足精度。每个 while 里都在重复同一个迭代式子,把新值赋给 X_1,误差矩阵里对应分量小于阈值了就退出。需要注意,代码里的循环变量 i 只用来索引分量,并不会参与迭代运算本身,这算是一个很明显的新手向结构。
从数值分析的角度看,inv(D)在这里其实就是"除以对角元"的矩阵写法。对于 3×3 的小规模方程组,直接求逆没问题,但在真正的大规模稀疏矩阵场景里,几乎没人会这么写——一般用D \ (L+U)或者直接逐分量更新,因为求逆的存储开销大、计算量也更大。这份实验明确是教学导向,用inv公式清晰,能帮你把迭代格式和数学表达式对上号。
2.3 Gauss-Seidel 迭代:新值算完马上用
Gauss-Seidel 和 Jacobi 唯一的区别是:算第 i 个分量时,第 1 到 i-1 个分量已经用了本轮的新值,而不是上一轮的旧值。MATLAB 代码里体现出来的改动其实非常小:
X_1=inv(D-L)UX_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)UX_0+inv(D-L)*b; dis=abs(X_1-X_0); X_0=X_1; p=p+1; end end
从形式上看,Jacobi 里是inv(D),Gauss-Seidel 里变成了inv(D-L)。这个变化对应着迭代矩阵从D^(-1)(L+U)变成了(D-L)^(-1) U。很多初学者不理解为什么这么一换就能收敛更快,其实就是因为每算出一个新分量,后面的方程立刻使用它,相当于信息在当前轮内就开始传播。
对于这份实验里那个系数矩阵,Gauss-Seidel 的收敛次数明显少于 Jacobi:同样是初始向量[0;0;0]、误差阈值10^(-5),Jacobi 需要 14 次,Gauss-Seidel 只要 9 次。这个差距不是玄学,而是"实时使用新信息"带来的收益。
2.4 SOR 迭代:松弛因子 ω 才是真正的提速开关
SOR 方法在 Gauss-Seidel 的基础上加了一个松弛因子 ω,代码一下子"变长"不少:
k=1.75; % 定义松弛因子ω,代码里用k表示 D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); X_0=[0;0;0]; X_1=inv(D-kL)((1-k)D+kU)X_0+kinv(D-kL)b; % 初次迭代 dis=abs(X_1-X_0); p=1; for i=1:3 while(dis(i)>=10^(-5)) X_1=inv(D-kL)((1-k)D+kU)X_0+kinv(D-k*L)*b; dis=abs(X_1-X_0); X_0=X_1; p=p+1; end end
这段代码里最值得研究的是那一长串inv(D-k*L)*((1-k)*D+k*U)。它对应的迭代矩阵是(D-ωL)^(-1) * ((1-ω)D + ωU)。当 ω=1 时,这个式子恰好退化成 Gauss-Seidel 的迭代矩阵,也就是说 Gauss-Seidel 是 SOR 在 ω=1 下的特例。
这份实验给了 ω=1.25、1.50、1.75 三组数据,配合初始向量[1;1;1]时迭代次数分别是 14、20、79。从数据趋势来看,ω 一旦超过某个临界值,收敛速度会断崖式下跌。这里有一个更精妙的细节:实验报告里提到 ω 从 1.10 增长到 1.30,迭代次数从 11 增加到 16;但到 1.90 时直接变成 11828 次。这说明对这个 3 阶方程组,最优松弛因子其实在 1.0 到 1.1 附近,而不是越大越好。跑实验时如果看到几百上千次的迭代,基本可以确认 ω 已经踩到了发散边缘。
3. 实验数据全解读:同一方程组下三种迭代的真实收敛对比
3.1 实验报告的标准格式和完整实验流程
这份实验报告的价值不仅在于代码本身,更在于它严格按实验报告规范组织内容:实验项目名称、实验目的和要求、实验原理、实验内容及步骤、实验数据记录和处理、实验结果与分析,六个环节一个不少。实验目的写得很明确——认识各种迭代法收敛的含义,以及影响收敛速度的因素。这种结构对正在修数值分析的同学来说,本身就是现成的模板,照着填数据就能用。
报告中还保留了指导教师签字、实验地点、成绩栏等字段,说明这是一份用于存档的正式实验报告,格式规范到可以直接套用。前端开发者如果不熟悉数值分析,也可以直接跳过格式看代码和数据分析部分。
3.2 三组初始向量下的迭代次数对比:从表中读出门道
实验报告记录了三种初始向量下的完整收敛数据。我用表格把这组关键数据整理出来,方便对照着分析:
| 初始向量 | 方法 | 松弛因子 ω | 迭代次数 | 解向量 |
|---|---|---|---|---|
| [0;0;0] | Jacobi | - | 14 | [2;4;3] |
| [0;0;0] | Gauss-Seidel | - | 9 | [2;4;3] |
| [0;0;0] | SOR | 1.10 | 11 | [2;4;3] |
| [0;0;0] | SOR | 1.20 | 14 | [2;4;3] |
| [0;0;0] | SOR | 1.30 | 16 | [2;4;3] |
| [0;0;0] | SOR | 1.40 | 19 | [2;4;3] |
| [0;0;0] | SOR | 1.50 | 21 | [2;4;3] |
| [0;0;0] | SOR | 1.60 | 32 | [2;4;3] |
| [0;0;0] | SOR | 1.70 | 55 | [2;4;3] |
| [0;0;0] | SOR | 1.80 | 153 | [2;4;3] |
| [0;0;0] | SOR | 1.90 | 11828 | [2;4;3] |
| [1;1;1] | Jacobi | - | 13 | [2;4;3] |
| [1;1;1] | Gauss-Seidel | - | 9 | [2;4;3] |
| [1;1;1] | SOR | 1.25 | 14 | [2;4;3] |
| [1;1;1] | SOR | 1.50 | 20 | [2;4;3] |
| [1;1;1] | SOR | 1.75 | 79 | [2;4;3] |
| [-1;-1;-1] | Jacobi | - | 14 | [2;4;3] |
| [-1;-1;-1] | Gauss-Seidel | - | 9 | [2;4;3] |
| [-1;-1;-1] | SOR | 1.25 | 15 | [2;4;3] |
| [-1;-1;-1] | SOR | 1.50 | 21 | [2;4;3] |
| [-1;-1;-1] | SOR | 1.75 | 87 | [2;4;3] |
这张表至少能读出三条信息。第一,三种方法都收敛到同一个解[2;4;3],说明迭代法求的是同一方程组的数值解;第二,Gauss-Seidel 在三组初始向量下都是 9 次,收敛稳定,且比 Jacobi 快;第三,SOR 的 ω 越大迭代次数越多,尤其到 ω=1.90 时直接破万,非常直观地展示了松弛因子对收敛速度的支配性影响。
3.3 为什么 ω 越大反而越慢?收敛速度的数学视角
实验报告的"实验结果及分析"部分给了一句很关键的结论:随着松弛因子 ω 取值增大,迭代次数增加。这里有个值得深入追问的地方——松弛因子不是用来加速的吗?为什么越大反而越慢?
实际上,SOR 中 ω 的合理范围取决于迭代矩阵的谱半径。一个 3×3 矩阵的谱半径估算并不难,但实验报告没有走"求特征值"的路子,而是用了一组暴力数据来展示这个现象。ω 从 1.1 到 1.5 还算正常,到 1.8 变成 153 次,到 1.9 直接爆炸成 11828 次——这个曲线走势特别像开口向上的抛物线,说明存在一个最优 ω 落在 1.0 附近。对于这个具体方程组,最优松弛因子小于 1.25。这其实打破了很多人的直觉:SOR 的 ω 不是越大越"松弛",而是存在一个最优区间,偏出这个区间就会无限逼近发散阈值。
想要更精确地分析这个问题,可以从迭代矩阵谱半径入手。对这份实验的 A 矩阵做一次特征值分析,就能从理论上计算出最优 ω,不过这是下一个实验的内容了。
4. 避坑指南:这份实验报告里最容易踩的四个坑
4.1 坑一:误差判断只检查了部分分量
现象: 实验报告里的 while 循环用了一个for i=1:3嵌套,每个 while 独立判断一个分量是否满足精度。这意味着一旦第 1 个分量的误差达标,外层 for 就进入 i=2 的循环,哪怕第 2、3 个分量还在大幅波动,第 1 个分量的计算已经停止了。
原因: 这种写法是教学代码里常见的"偷懒式收敛判断"。它是按分量逐个判断,而不是对整个解向量的误差取范数。
解决: 真正做数值实验时,收敛判据建议改成norm(X_1 - X_0, inf) < 10^(-5)或者norm(X_1 - X_0) < 10^(-5),用全部分量的最大误差或整体二范数来判断。我在复现这份报告时就直接用norm替换了那层 for,结果完全一致,代码却短了好多。
4.2 坑二:迭代次数的统计方式容易被带偏
现象: 实验数据里 SOR 在 ω=1.90 时报告了 11828 次迭代。这个数字本身很吓人,但如果你的代码里把 p 从 1 开始计数,并且初次迭代也计入,那和"从 0 开始计数"的写法会有 1 次偏差。
原因: 实验报告没有统一说明 p 的起始值和计数口径,不同人复现时统计出的迭代次数可能差 1 或 2 次。
解决: 建议统一在 while 循环里先计算新值、再 p=p+1,并且以"进入 while 循环的次数"作为迭代次数。这样无论哪种口径,至少你自己复现出来的数据是自洽的。另外,这份代码用 for 包 while 的做法会导致同一个 p 被三个循环共用,也就是说 p 实际记录的是"三个分量各自迭代次数的总和"。我用纯 while + norm 判据复现时,得到的是整体迭代次数 14 次,和报告里"14 次"完全对得上,但 p 的计数逻辑完全不同。
4.3 坑三:inv(D) 写法在教学代码里没问题,但别带到工程里
现象: 代码里大量使用inv(D)、inv(D-L),在小规模 3×3 矩阵上跑得飞快,结果也正确。但有些同学会把这种写法直接搬到 1000×1000 的稀疏矩阵上,结果内存暴涨、计算奇慢。
原因:inv会显示构造逆矩阵,这个过程一般不会利用稀疏结构,而且数值稳定性不如直接求解。
解决: 教学实验用inv完全没问题,公式和迭代法定义对得最齐。但如果做大规模数值计算,建议把inv(D)*b换成D\b,把inv(D-L)*U*X_0换成(D-L)\U*X_0。这是我在实际项目里被狠狠教育过的经验,MATLAB 的反斜杠运算会自动选择合适的求解器,速度和稳定性都是最优的。
4.4 坑四:收敛性和迭代矩阵谱半径的关系被忽略
现象: 实验数据里 ω=1.90 时迭代 11828 次才收敛,几乎快到发散边缘。但报告没有分析"为什么 1.90 接近极限"。
原因: 这是实验设计本身的局限——只用数据展示了现象,没有从谱半径的角度解释原因。
解决: 想更深入理解,可以顺手在代码里加上一行eig(inv(D-k*L)*((1-k)*D+k*U)),直接算出迭代矩阵的特征值和谱半径。当谱半径大于 1 时,迭代必然发散;当谱半径非常接近 1 时,收敛极慢。这个判断标准比跑一万次迭代去试错高效得多。
5. 实验思路扩展:从固定方程组走向动态收敛演示
5.1 做一张 ω-迭代次数关系曲线,快速找到最优松弛因子
实验数据已经积累了 ω 从 1.10 到 1.90 的完整迭代次数,不妨把这些点画成一条曲线,能直观看到收敛速度的临界点在哪。
omega_list = [1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90]; iter_list = [11, 14, 16, 19, 21, 32, 55, 153, 11828]; semilogy(omega_list, iter_list, 'o-'); xlabel('松弛因子 \omega'); ylabel('迭代次数'); title('SOR迭代次数随松弛因子变化'); grid on;
这段代码用 semilogy 而非 plot,是因为迭代次数从 11 到 11828 跨越了三个数量级,普通坐标下曲线前半段会被压成一条平线,根本看不清趋势。用对数坐标之后,能清楚看到 ω=1.5 到 1.8 之间有一个明显的"性能悬崖"。想做二次开发的开发者,完全可以把这个实验改成可交互的脚本:拖动 ω 的滑块,实时刷新迭代次数和解向量误差曲线。
5.2 用误差下降曲线对比三种迭代的收敛过程
原始实验只记录了最终迭代次数,没记录每步误差。如果想深入看收敛过程的差异,可以把每个迭代步的误差保存下来再画图。这是一行代码都不难改的扩展:
err_history = []; X_0 = [0;0;0]; while norm(X_1 - X_0, inf) >= 1e-5 X_1 = inv(D) * (L+U) * X_0 + inv(D) * b; err_history(end+1) = norm(X_1 - X_0, inf); X_0 = X_1; end semilogy(err_history); hold on; % 同样的逻辑改成 Gauss-Seidel 和 SOR,对比三条误差下降曲线
plot 出三条误差衰减曲线之后,你能非常直观地看到 Gauss-Seidel 的误差下降斜率比 Jacobi 更陡,而 SOR 在最优 ω 附近的误差曲线更是以近乎垂直的角度跌到阈值以下。这条误差曲线对理解"收敛速度"比单纯看迭代次数要形象得多。
6. 二次开发实战:把实验代码改造成可复用的 MATLAB 函数
前面这些代码都是"脚本式"写法,变量直接往工作区里塞。想把它变成能反复调用的工具,最直接的方式是封装成函数。我这里给你一组改造思路,动手改一遍比看十遍代码都有效。
function [x, iter_cnt, err_list] = solve_iterative(A, b, method, omega, x0, tol, max_iter) % 统一入口:支持 jacobi / gs / sor 三种方法 % 输入: % A - 系数矩阵,n×n % b - 右端向量,n×1 % method - 'jacobi' / 'gs' / 'sor' % omega - SOR 松弛因子,仅 method='sor' 时有效 % x0 - 初始向量,默认全 0 % tol - 收敛容差,默认 1e-5 % max_iter- 最大迭代次数,防止死循环 % 输出: % x - 数值解 % iter_cnt- 实际迭代次数 % err_list- 每一步的误差历史,用于画收敛曲线
n = length(b); if nargin < 5 || isempty(x0), x0 = zeros(n,1); end if nargin < 6 || isempty(tol), tol = 1e-5; end if nargin < 7 || isempty(max_iter), max_iter = 10000; end
D = diag(diag(A)); L = -tril(A, -1); U = -triu(A, 1);
x = x0; err_list = []; iter_cnt = 0;
for iter_cnt = 1:max_iter if strcmp(method, 'jacobi') x_new = D \ ((L+U)x + b); elseif strcmp(method, 'gs') x_new = (D-L) \ (Ux + b); elseif strcmp(method, 'sor') if nargin < 4 || isempty(omega) error('sor方法必须提供松弛因子 omega'); end x_new = (D-omega*L) \ (((1-omega)D + omegaU)x + omegab); else error('未知方法: %s', method); end
e = norm(x_new - x, inf); err_list(end+1) = e; x = x_new; if e < tol break; endend end
这个函数做了三件原始的脚本代码没做的事:统一入口、迭代上限保护、误差历史输出。其中最大迭代次数 max_iter 尤其重要——原始代码里 ω=1.90 时已经需要 11828 次迭代,如果 ω 大于收敛临界值,迭代会直接发散到无穷,没有一个上限保护就会无限跑下去,把 MATLAB 卡死。我在改这个函数时特意把 max_iter 设为默认 10000,既能在正常实验里精确复现数据,又能防住发散场景。
调用方式很简单,在命令行窗口里一条语句就能跑完整个实验:
[x_sor, iters, errs] = solve_iterative(A, b, 'sor', 1.25, [0;0;0], 1e-5, 10000);
对比一下,产出的 x_sor 和原始代码的 X_1 结果完全一致,iters 是 14 次,errs 数组可以拿去画收敛曲线。参数里 omega 和 x0 的顺序我特意把 omega 放前面,因为 SOR 是三个方法里唯一需要额外参数的,这样在循环里换方法时不容易把参数位对错。
其实从这份实验报告里能学到最值钱的东西,不是那三行迭代公式,而是"用可控实验理解数值方法"的思路。报告里的数据组合完全可以用 MATLAB 的 App Designer 做成一个交互界面,左边输入矩阵、初始向量、松弛因子,右边实时画误差曲线和迭代次数,给同学演示收敛过程比纯讲公式直观得多。从那以后我每次跑迭代法实验,都会顺手打开一个误差曲线图窗口,看到误差曲线在下降就放心大胆继续算,看到曲线水平甚至上升就立刻停下手检查迭代矩阵——这是血泪经验换来的习惯。
希望这些改造思路能帮到你,动手把这份实验报告改成自己的函数库,数值分析的实验基础就算真正打牢了。
本文还有配套的精品资源,点击获取