企业官网建设流程全解析

简介：这份资源是哈工大数值分析上机实验前四部分的Matlab程序代码合集，面向正在学习数值分析课程、需要动手实现经典数值方法的高校学生与自学者，可用于课程实验对照、算法复现与编程练习。压缩包共15个文件，以13个.m脚本为主体，另含1个.md说明文档与1个.asv自动备份文件，整体约8KB，体积轻量便于快速查阅。内容覆盖高斯列主元消去法、龙格现象的产生与克服、非线性方程求解、曲线拟合的最小二乘法等典型实验主题，每个实验配有多个脚本，便于对比不同实现思路与参数设置。已有154人学习下载，适合作为数值分析上机环节的参考代码，帮助读者理解算法流程、验证计算结果并积累Matlab编程经验，仅限学习参考使用。

1. 数值分析上机实验：一份能直接跑通的 Matlab 代码底稿

数值分析这门课最难受的地方，不是公式推不出来，而是上机实验里那些迭代过程——明明手算能收敛，写进 Matlab 却发散；明明课本上写着“取初值 x0 即可”，跑出来却是 NaN。这份哈工大数值分析上机实验前四部分的 Matlab 程序代码，解决的正是这个痛点：它把插值、数值积分、线性方程组迭代解法、非线性方程求根这几类经典实验的完整脚本都摆出来了，每个实验一个可独立运行的 .m 文件，参数、迭代次数、误差输出都写在明面上。适合两类人：一是正在上数值分析课、需要交上机报告的学生，二是想快速复习一遍数值方法实现细节的开发者。代码本身不复杂，但胜在结构清晰、注释到位，拿来当模板改比从零写省事得多。

2. 项目结构与技术选型：四个实验模块怎么组织

2.1 目录布局与文件命名

这份代码包按实验序号分目录，每个目录下是独立的 .m 脚本文件，命名基本遵循“实验内容_方法名.m”的规则。典型结构如下：

数值分析上机实验/ ├── 实验一_插值/ │ ├── lagrange_interp.m │ ├── newton_interp.m │ └── runge_demo.m ├── 实验二_数值积分/ │ ├── trapezoid_rule.m │ ├── simpson_rule.m │ └── romberg_integration.m ├── 实验三_线性方程组迭代/ │ ├── jacobi_iter.m │ ├── gauss_seidel.m │ └── sor_iter.m └── 实验四_非线性方程求根/ ├── bisection_method.m ├── newton_raphson.m └── secant_method.m

这种布局的好处是每个实验自包含，不依赖其他目录的函数，复制单个文件夹出去也能跑。缺点是公共工具函数（比如误差计算、迭代终止判断）在每个文件里重复实现了一遍。我一般会抽一个utils.m出来放公共逻辑，但作为教学代码，重复反而降低了阅读门槛——每个文件从头到尾能独立看懂，不用跳来跳去。

2.2 技术栈与运行环境

纯 Matlab 实现，没有用任何工具箱函数（比如interp1、integral、fsolve这些），所有算法都是手写的。这一点很关键：教学代码如果直接调内置函数，学生就看不到算法内部逻辑了。运行环境要求 Matlab R2016b 及以上即可，基本不挑版本。如果手头没有 Matlab，Octave 也能跑，但要注意两点：一是 Octave 的fprintf格式化输出和 Matlab 略有差异，二是部分脚本用了 Matlab 的inputParser做参数校验，Octave 对它的支持不完整，遇到报错就把那几行注释掉。

代码里没有用面向对象，全是函数式脚本，每个 .m 文件要么是一个函数定义，要么是一段可直接运行的脚本。这种写法对初学者友好，但也意味着变量作用域是全局的——跑完一个脚本再跑另一个，工作区里会残留上一轮的变量，偶尔会导致“明明改了参数结果没变”的乌龙。

2.3 四个实验模块的功能边界

实验一插值模块包含拉格朗日插值和牛顿插值两种实现，附带一个龙格现象演示脚本，用等距节点插值 Runge 函数，画图对比不同节点数下的逼近效果。实验二数值积分覆盖梯形公式、辛普森公式和龙贝格积分，每个方法都输出了步长、积分值和误差估计。实验三线性方程组迭代解法实现了雅可比迭代、高斯-赛德尔迭代和逐次超松弛迭代（SOR），重点在于迭代矩阵的构造和收敛性判断。实验四非线性方程求根包含二分法、牛顿法和弦截法，每个方法都记录了迭代次数和最终误差。

四个模块之间没有调用关系，可以按任意顺序跑。但建议按实验序号来，因为后面的实验会用到前面实验里定义的辅助函数（比如实验四的牛顿法里用到了实验一里的差商计算逻辑，虽然是复制粘贴过去的）。

3. 核心算法拆解：从插值到求根的关键实现

3.1 拉格朗日插值与牛顿插值的代码对比

拉格朗日插值的核心是构造基函数然后加权求和。代码里是这样写的：

function y = lagrange_interp(x_nodes, y_nodes, x_query) % x_nodes: 插值节点横坐标（行向量） % y_nodes: 插值节点纵坐标（行向量） % x_query: 待求值的查询点（可以是向量） n = length(x_nodes); y = zeros(size(x_query)); for k = 1:n % 构造第 k 个拉格朗日基函数 L_k(x) Lk = ones(size(x_query)); for j = 1:n if j ~= k Lk = Lk .* (x_query - x_nodes(j)) / (x_nodes(k) - x_nodes(j)); end end y = y + y_nodes(k) * Lk; end end

这段代码的逻辑很直白：外层循环遍历每个节点，内层循环构造对应的基函数，最后累加。参数说明里x_query支持向量输入，是因为画图时需要一次性算出一整条曲线上的值。注意内层循环用了.*和./做逐元素运算，这是 Matlab 里处理向量化输入的标准写法，如果写成*和/会报矩阵维度错误。

牛顿插值的实现思路不同，它先构造差商表，再用秦九韶算法（Horner 法则）求值：

function y = newton_interp(x_nodes, y_nodes, x_query) n = length(x_nodes); % 构造差商表，第一列是函数值 divided_diff = zeros(n, n); divided_diff(:, 1) = y_nodes(:); for j = 2:n for i = j:n divided_diff(i, j) = (divided_diff(i, j-1) - divided_diff(i-1, j-1)) ... / (x_nodes(i) - x_nodes(i-j+1)); end end % 用 Horner 法则求值 y = divided_diff(n, n) * ones(size(x_query)); for k = n-1:-1:1 y = divided_diff(k, k) + (x_query - x_nodes(k)) .* y; end end

差商表的构造是这段代码的关键：divided_diff(i, j)表示从第i-j+1个节点到第i个节点的j-1阶差商。构造完成后，对角线元素就是牛顿插值多项式的系数。Horner 法则从最高次项开始往回卷，每一步做一次乘加，比直接展开多项式求值数值稳定性更好。实际跑的时候，如果节点数超过 20 个，拉格朗日插值的计算量会明显大于牛顿插值，因为前者是 O(n²) 的基函数构造，后者虽然差商表也是 O(n²)，但求值阶段是 O(n)。

3.2 龙贝格积分的递推实现

数值积分模块里最有意思的是龙贝格积分，它用递推的方式逐步提高精度。代码结构如下：

function [I, R] = romberg_integration(f, a, b, tol, max_iter) % f: 被积函数句柄 % a, b: 积分区间 % tol: 容差 % max_iter: 最大迭代次数 R = zeros(max_iter, max_iter); h = b - a; R(1, 1) = h / 2 * (f(a) + f(b)); % 梯形公式初值 for k = 2:max_iter h = h / 2; % 复化梯形递推 sum_f = 0; for i = 1:2^(k-2) sum_f = sum_f + f(a + (2*i - 1) * h); end R(k, 1) = R(k-1, 1) / 2 + h * sum_f; % Richardson 外推 for j = 2:k R(k, j) = R(k, j-1) + (R(k, j-1) - R(k-1, j-1)) / (4^(j-1) - 1); end % 收敛判断 if abs(R(k, k) - R(k-1, k-1)) < tol I = R(k, k); R = R(1:k, 1:k); return; end end I = R(max_iter, max_iter); warning('龙贝格积分未在最大迭代次数内收敛'); end

这里有两个关键点。第一是复化梯形的递推公式：R(k,1) = R(k-1,1)/2 + h*sum_f，它利用了上一次的梯形值和新增加的节点函数值，避免重复计算已有节点。第二是 Richardson 外推：R(k,j) = R(k,j-1) + (R(k,j-1)-R(k-1,j-1))/(4^(j-1)-1)，这个公式把梯形序列加速成辛普森序列、再加速成柯特斯序列，依此类推。参数tol控制的是对角线相邻两次的差值，不是绝对误差——这一点在写实验报告时容易被扣分，因为严格来说应该用后验误差估计。

实际跑的时候，如果被积函数在区间端点有奇性（比如1/sqrt(x)在 0 附近），梯形公式初值会算出 Inf，整个递推就崩了。这种情况需要先做变量替换把奇性消掉，代码里没有处理，属于使用时的注意事项。

3.3 雅可比迭代与 SOR 的收敛条件

线性方程组迭代解法模块里，雅可比迭代的实现最简洁：

function [x, iter] = jacobi_iter(A, b, x0, tol, max_iter) % A: 系数矩阵（要求对角占优） % b: 右端向量 % x0: 初始猜测 n = length(b); x = x0; D = diag(diag(A)); L_plus_U = A - D; for iter = 1:max_iter x_new = D \ (b - L_plus_U * x); if norm(x_new - x, inf) < tol x = x_new; return; end x = x_new; end warning('雅可比迭代未收敛'); end

这里用矩阵形式写迭代：x_new = D\(b - (L+U)x)，比逐元素循环快得多，也更贴近课本上的矩阵表达式。norm(x_new - x, inf)取的是无穷范数，即最大分量差，这是迭代法常用的终止准则。注意D \ (...)用的是左除，因为 D 是对角矩阵，Matlab 会自动识别并做逐元素除法，效率很高。

SOR 迭代在雅可比的基础上加了一个松弛因子omega：

function [x, iter] = sor_iter(A, b, x0, omega, tol, max_iter) n = length(b); x = x0; for iter = 1:max_iter x_old = x; for i = 1:n sigma = 0; for j = 1:n if j ~= i sigma = sigma + A(i, j) * x(j); end end x(i) = (1 - omega) * x_old(i) + omega * (b(i) - sigma) / A(i, i); end if norm(x - x_old, inf) < tol return; end end warning('SOR 迭代未收敛'); end

SOR 这里用的是逐元素更新，因为高斯-赛德尔类型的迭代必须用最新计算出的分量，没法写成纯矩阵形式。omega的取值很关键：omega=1退化成高斯-赛德尔，0<omega<1是低松弛，1<omega<2是超松弛。理论上的最优松弛因子需要知道雅可比迭代矩阵的谱半径，代码里没有自动计算，需要手动调参。我一般会先跑几个omega值看迭代次数，画个曲线找最小值。

4. 跑通与调试：环境配置和常见报错处理

4.1 运行前的准备工作

把压缩包解压到任意目录，然后在 Matlab 里把当前工作目录切到解压后的根目录，或者用addpath(genpath(pwd))把所有子目录加进搜索路径。每个实验目录下都有一个run_*.m脚本作为入口，直接运行即可。比如实验一：

% 在 Matlab 命令行中执行 cd('实验一_插值'); runge_demo; % 运行龙格现象演示

runge_demo.m会生成一张对比图，横轴是插值区间，纵轴是函数值和插值多项式值，不同节点数用不同颜色区分。如果图窗没弹出来，检查一下是不是在无头模式下运行（比如通过 SSH 连的服务器），那种情况下需要加-nodisplay参数或者改用saveas存成图片文件。

4.2 参数调整与结果验证

每个脚本开头都有一组可调参数，比如插值节点数、积分区间、迭代容差。以实验二为例：

% 实验二参数设置（在脚本开头） f = @(x) exp(x) .* cos(x); % 被积函数 a = 0; b = pi/2; % 积分区间 tol = 1e-8; % 容差 max_iter = 20; % 最大迭代次数

改完参数直接重新运行脚本即可。验证结果是否正确，最直接的办法是和解析解对比：exp(x)*cos(x)在[0, pi/2]上的积分有闭式解，算出来大约是 1.9052，如果龙贝格积分输出和这个值差太多，说明参数或实现有问题。另一个验证手段是看收敛阶：梯形公式是二阶收敛，步长减半误差应该降到约 1/4；辛普森公式是四阶，步长减半误差降到约 1/16。如果实测收敛阶对不上，大概率是代码里某个系数写错了。

4.3 常见报错与排查思路

跑这份代码最容易遇到的报错是“未定义函数或变量”。原因通常是当前工作目录不对，或者函数文件名和函数名不一致。Matlab 要求文件名必须和文件内第一个函数名相同，比如lagrange_interp.m里定义的函数必须叫lagrange_interp，改文件名会导致调用失败。

第二个高频报错是矩阵维度不匹配。插值函数里x_nodes和y_nodes必须是等长向量，如果从 Excel 或文本文件导入数据时行列搞反了，就会触发这个错误。排查方法是在报错行之前加一句size(x_nodes)和size(y_nodes)打印维度。

第三个坑是迭代法不收敛但代码不报错。雅可比迭代和高斯-赛德尔迭代都有收敛条件：系数矩阵必须严格对角占优，或者迭代矩阵谱半径小于 1。如果随便构造一个矩阵丢进去，可能迭代几百次还在震荡。代码里设了max_iter兜底，超了会打 warning，但不会中断。看到 warning 时先检查矩阵是否对角占优，不满足的话要么换方法（比如用直接法），要么做预处理。

5. 避坑与常见问题：跑这份代码时容易翻车的几个地方

5.1 现象：龙格现象演示图看起来“没问题”

原因：runge_demo.m默认只画了 5 个和 9 个节点的情况，这两个节点数下插值多项式在区间中部还能看，边缘的震荡不明显。很多人跑完这张图就以为等距插值没问题，实际上节点数加到 15 以上，边缘震荡会急剧放大。

解决：手动把节点数改成[5, 9, 15, 21]四组，重新跑一遍，把图窗放大看区间端点附近。会看到 21 个节点时插值曲线在两端几乎垂直上下跳，这就是龙格现象。理解了这个现象，才能明白为什么实验一后面要引入切比雪夫节点——代码包里没有实现切比雪夫节点插值，但可以自己加，把x_nodes = linspace(a, b, n)换成x_nodes = (a+b)/2 + (b-a)/2 * cos((2*(1:n)-1)*pi/(2*n))即可。

5.2 现象：龙贝格积分跑出来是 NaN

原因：被积函数在积分区间端点或内部有奇点，比如f = @(x) 1./sqrt(x)在x=0处是 Inf。梯形公式初值R(1,1) = h/2*(f(a)+f(b))直接算出 Inf，后续递推全是 NaN。

解决：先检查被积函数在[a,b]上的取值，用f(linspace(a,b,100))扫一遍看有没有 Inf 或 NaN。如果有奇点，做变量替换把奇点消掉，比如x = t^2可以把1/sqrt(x)变成2的常数被积函数。代码包里没有内置这个处理，需要自己改。

5.3 现象：SOR 迭代的迭代次数比高斯-赛德尔还多

原因：omega取值不当。SOR 的收敛速度对omega很敏感，omega偏离最优值哪怕 0.1，迭代次数可能翻倍。代码里默认omega=1.2，但这个值只对特定矩阵有效，换个矩阵就不灵了。

解决：写个小循环扫omega从 0.1 到 1.9，步长 0.1，记录每个omega下的迭代次数，画出来找最小值点。更严谨的做法是先算雅可比迭代矩阵的谱半径rho，最优松弛因子理论值是omega_opt = 2/(1+sqrt(1-rho^2))。代码包里没有算谱半径的函数，可以自己加一行rho = max(abs(eig(inv(D)*(L_plus_U))))。

5.4 现象：牛顿法求根迭代次数比预期多很多

原因：初始猜测离根太远，或者函数在根附近导数接近零。牛顿法的收敛是局部的，初值选得不好可能发散或者收敛到别的根。代码里默认初值写死在脚本开头，换个函数忘了改初值就会出问题。

解决：先用二分法粗定位根的大致区间，再把区间中点作为牛顿法的初值。代码包里实验四的bisection_method.m和newton_raphson.m是分开的，可以手动串起来用：先跑二分法拿到根所在区间，再把区间端点代入牛顿法。另外检查一下导数表达式有没有写错，数值求导和解析求导差太多也会导致迭代异常。

5.5 现象：换台电脑跑，图窗里的中文变成乱码

原因：Matlab 的默认字体在不同操作系统下不一样，Windows 下中文正常，Linux 或 macOS 下可能显示成方框。代码里的title、xlabel、ylabel用了中文字符串。

解决：在脚本开头加一句set(0, 'DefaultAxesFontName', 'SimHei')（Windows）或set(0, 'DefaultAxesFontName', 'PingFang SC')（macOS）。更省事的办法是把图里的中文标签改成英文，反正实验报告里图注可以用中文写，图本身用英文标签不影响理解。

6. 二次开发：把教学代码改成能用的工具函数

这份代码作为教学材料是合格的，但直接拿来当工程工具用还有距离。最值得做的改造是把每个算法封装成带参数校验和错误处理的函数，统一输入输出接口。比如把插值函数改成支持“给定节点和查询点，返回插值值和误差估计”的形式：

function [y, err_est] = interp_unified(method, x_nodes, y_nodes, x_query) % method: 'lagrange' 或 'newton' % 返回插值结果和基于差商的误差估计 switch lower(method) case 'lagrange' y = lagrange_interp(x_nodes, y_nodes, x_query); err_est = []; % 拉格朗日插值不易给出误差估计 case 'newton' [y, dd] = newton_interp_with_diff(x_nodes, y_nodes, x_query); % 用 n+1 阶差商乘以节点多项式估计误差 n = length(x_nodes); omega = ones(size(x_query)); for k = 1:n omega = omega .* (x_query - x_nodes(k)); end err_est = abs(dd(n+1)) * abs(omega); % 需要额外一个节点算 n+1 阶差商 otherwise error('不支持的插值方法: %s', method); end end

这个改造的关键点在于误差估计：牛顿插值天然能给出余项表达式R(x) = f[x0,...,xn,x] * omega(x)，其中f[x0,...,xn,x]是 n+1 阶差商。实际实现时需要多给一个节点来算这个差商，或者用相邻两个 n 阶差商的差来近似。代码包里原来的newton_interp没有返回差商表，改造时要把它作为第二个输出参数暴露出来。

另一个值得做的改造是给迭代法加自动收敛性诊断。现在的代码只会在超迭代次数时打 warning，但不会告诉你为什么不收敛。可以加一段诊断逻辑：先算迭代矩阵的谱半径，如果大于等于 1 就直接报错并提示“该矩阵不满足收敛条件，建议改用直接法或做预处理”；如果小于 1 但接近 1，提示“收敛速度较慢，建议调整松弛因子或换方法”。谱半径的计算对大规模矩阵有开销，可以设一个阈值，矩阵维度超过 500 就跳过诊断直接跑。

验证改造是否成功，我一般用三个测试用例：一个对角占优矩阵（应该快速收敛）、一个对称正定但非对角占优矩阵（高斯-赛德尔可能收敛但雅可比发散）、一个病态矩阵（应该触发诊断提示）。三个用例都跑通，说明封装没问题。

还有一个容易忽略的点：原代码里所有函数都没有做输入维度检查。比如jacobi_iter里如果b是列向量而x0是行向量，x_new - x会触发维度广播，算出来的范数可能是错的但不报错。改造时在函数开头加一句assert(iscolumn(b) && iscolumn(x0), 'b 和 x0 必须是列向量')，能省掉很多调试时间。

从那以后我每次把教学代码往工程代码改，都强制先写三个测试用例再动主逻辑——教学代码的“能跑”和工程代码的“跑对”之间，差的就是这几行断言和诊断。希望帮到你。

本文还有配套的精品资源，点击获取