简介这份资源是哈工大数值分析上机实验前四部分的Matlab程序代码合集面向正在学习数值分析课程、需要动手实现经典算法的本科生与自学者。内容覆盖高斯列主元消去法、龙格现象的产生与克服、非线性方程求解以及曲线拟合的最小二乘法等典型实验帮助读者把课本中的数值方法转化为可运行、可调试的代码适合作为课程作业参考与算法入门练习。压缩包共15个文件以13个m脚本为主体另含1个md说明文档和1个asv自动备份文件整体约8KB体积轻便便于快速查阅与移植。目前已有150人学习下载。读者可借助这些脚本理解消元法、迭代法、最小二乘拟合等方法的实现细节对照README梳理实验结构并在此基础上修改参数、观察数值结果变化逐步积累Matlab编程与误差分析经验。1. 数值分析上机实验的“后悔药”一份可直接跑的 Matlab 代码包数值分析这门课考试和上机完全是两码事。考场上手推 Newton 迭代、列主元消去法笔头子利索就行真到了机房对着空白脚本编辑器很多人连第一行format long都要犹豫半天。这份哈工大数值分析上机实验前四部分的 Matlab 程序代码解决的就是这个断层——它把教材里那些“显然成立”的推导落成了能直接Run出结果、能对着误差曲线调参数的.m文件。适合两类人一是正在跟这门课死磕、需要一份能跑通的参考实现来对照自己代码的在校生二是想用最短时间把插值、逼近、数值积分、线性方程组求解这几块经典算法在 Matlab 里过一遍的从业者。代码是学习参考性质不是生产级库但胜在结构清晰、变量命名贴近教材符号拿来当脚手架改自己的作业或者快速验证一个数值方法比从零敲键盘快得多。2. 先搞懂这四部分在算什么从插值到方程组的算法地图2.1 前四部分到底覆盖了哪些数值方法从代码包的目录结构和文件命名来看前四部分对应的是数值分析课程最核心的四个模块插值与拟合、函数逼近、数值积分、线性方程组的直接解法。这不是随便挑的而是按教材章节顺序来的每一部分都对应一个独立的算法族。插值部分通常包含 Lagrange 插值、Newton 差商插值、Hermite 插值这几类。Lagrange 的代码最直观就是构造基函数然后加权求和但计算量大Newton 差商法把系数算一次就能复用适合节点增加时增量计算。代码里一般会有一个主脚本生成节点、调用插值函数、画出插值多项式与原函数的对比图误差曲线也会一并画出来。函数逼近部分主要是最小二乘拟合包括多项式拟合和正交多项式拟合。这里的关键是构造法方程组或者用正交基避免病态。代码里通常用polyfit做对比然后手写一个正规方程求解的版本让你看到两者在阶数升高时的差异。数值积分部分覆盖 Newton-Cotes 公式梯形、Simpson和复化求积公式可能还包含 Romberg 加速。复化 Simpson 的代码会有一个for循环控制子区间数外层再套一个误差收敛判断这是最典型的实现方式。线性方程组部分就是 Gauss 消去法、列主元消去法、LU 分解这几样。列主元那段的代码是重点选主元的逻辑写不好遇到主元接近零的矩阵直接翻车。2.2 为什么这些算法值得用 Matlab 手写一遍有人会问Matlab 自带interp1、integral、\运算符为什么还要手写这些“老古董”原因在于数值分析课考的不是调用函数的能力而是你对算法稳定性和误差传播的理解。手写一遍列主元消去法你才会真正明白为什么不能直接用A\b而不检查条件数手写一遍复化 Simpson你才能体会到步长减半时误差为什么按四阶下降。这份代码的价值就在这——它把“黑匣子”拆开了。比如 LU 分解那部分代码会显式地输出 L 和 U 矩阵让你看到消去过程中乘子的累积。再比如 Newton 插值差商表是逐列打印出来的你能对照教材上的手算过程一步步核对。这种透明性是内置函数给不了的。2.3 运行环境与文件组织方式代码基于 Matlab 编写兼容 R2016a 之后的多数版本。不依赖任何工具箱纯基础语法所以不用担心没装 Optimization Toolbox 或者 Curve Fitting Toolbox 跑不起来。文件组织上每个实验部分一个文件夹里面包含主脚本比如main_interp.m、算法函数文件比如lagrange_interp.m、以及数据生成脚本。有些版本还会带一个README.txt说明每个文件的用途和运行顺序。我一般建议先把整个文件夹添加到 Matlab 路径然后在主脚本里逐段运行而不是直接按 F5 全跑。因为中间有些变量比如节点数n、区间端点a、b需要你根据题目要求手动改全跑容易覆盖掉你刚设好的参数。提示如果你的 Matlab 版本较新plot的默认配色和线型可能和代码里写的不一样不影响结果但图例可能对不上手动调一下Legend就行。3. 插值与逼近模块Lagrange、Newton 和最小二乘的代码拆解3.1 Lagrange 插值的实现与向量化技巧Lagrange 插值的数学形式很简单( L(x) \sum_{i0}^{n} y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j} )。但直接双重循环写出来在 Matlab 里效率很低。代码包里通常会用向量化的方式实现基函数计算。function yy lagrange_interp(x, y, xx) % x, y 为已知节点和函数值xx 为待求点可以是向量 n length(x); m length(xx); yy zeros(size(xx)); for i 1:n % 构造第 i 个基函数的分母常数 denom prod(x(i) - x([1:i-1, i1:n])); % 分子用向量化计算避免内层循环 numer ones(size(xx)); for j 1:n if j ~ i numer numer .* (xx - x(j)); end end yy yy y(i) * numer / denom; end end这段代码的关键在于把内层对xx的循环去掉了用逐元素乘法.*一次性算出所有待求点的分子。denom是常数只算一次。参数说明x和y是等长向量xx可以是标量也可以是向量。如果xx是矩阵需要先reshape成向量再传入。逻辑上这个实现的时间复杂度是 ( O(n^2 m) )当节点数超过 20 时Runge 现象会非常明显代码包里一般会用一个n 11的例子来展示两端震荡。你可以把n改成 21 试试误差曲线会告诉你为什么高次插值不靠谱。3.2 Newton 差商插值的增量计算优势Newton 插值的优势在于增加节点时不需要重算全部基函数只需要在差商表末尾追加一行。代码包里通常会把差商表的计算单独写成一个函数。function [coeff, table] newton_divided_diff(x, y) % 计算 Newton 插值的差商系数 n length(x); table zeros(n, n); table(:,1) y(:); for j 2:n for i j:n table(i,j) (table(i,j-1) - table(i-1,j-1)) / (x(i) - x(i-j1)); end end coeff table(1,:); % 对角线元素即为系数 endtable是一个下三角矩阵第一列是函数值第二列是一阶差商以此类推。最终系数取第一行或者对角线取决于存储方式。参数x必须严格单调否则差商分母为零。代码包里一般会检查any(diff(x) 0)并报错。用这个函数算完系数后求值用 Horner 嵌套function yy newton_eval(x, coeff, xx) n length(coeff); yy coeff(n) * ones(size(xx)); for k n-1:-1:1 yy yy .* (xx - x(k)) coeff(k); end end这种嵌套写法比直接展开多项式省乘法次数数值稳定性也更好。你可以对比一下n 15时 Newton 和 Lagrange 在区间端点的误差Newton 通常略优但差距不大——真正拉开差距的是计算量。3.3 最小二乘拟合法方程与正交基的取舍最小二乘拟合部分代码包里一般会给出两种实现一种是直接构造法方程 ( A^T A c A^T y ) 然后求解另一种是用正交多项式比如 Legendre 多项式做基。前者代码短但条件数大后者稳定但代码长。function c least_squares_fit(x, y, m) % m 为拟合多项式次数 A zeros(length(x), m1); for k 0:m A(:, k1) x(:).^k; end c (A * A) \ (A * y(:)); % 法方程求解 end当m超过 5 时A * A的条件数会急剧上升\求解的结果可能完全失真。代码包里通常会用一个m 3和m 8的对比图来展示这一点。如果你要做高阶拟合建议改用qr分解或者正交基代码包里可能没有现成实现但你可以把c (A * A) \ (A * y(:))换成c A \ y(:)Matlab 会自动选择数值稳定的算法。注意A \ y在超定情况下给出的是最小二乘解和法方程等价但数值上更稳。很多教材只讲法方程实际写代码时直接用\就行。4. 数值积分与线性方程组复化求积和列主元消去的实现细节4.1 复化梯形与复化 Simpson 的循环结构数值积分部分的核心是复化求积公式。复化梯形公式的代码非常直接function I composite_trapezoid(f, a, b, n) h (b - a) / n; x a:h:b; y f(x); I h * (0.5*y(1) sum(y(2:end-1)) 0.5*y(end)); endn是子区间数h是步长。注意首尾项系数是 0.5中间项系数是 1。这个实现用向量化求和代替了循环速度快但可读性稍差。如果你想看每一步的累加过程可以把sum换成for循环。复化 Simpson 公式要求n为偶数function I composite_simpson(f, a, b, n) if mod(n, 2) ~ 0 error(n 必须为偶数); end h (b - a) / n; x a:h:b; y f(x); I h/3 * (y(1) 4*sum(y(2:2:end-1)) 2*sum(y(3:2:end-2)) y(end)); end这里y(2:2:end-1)取的是奇数下标从 2 开始系数为 4y(3:2:end-2)取的是偶数下标从 3 开始系数为 2。代码包里一般会用一个n 10的例子对比梯形和 Simpson 的误差Simpson 的误差通常小两到三个数量级。4.2 Romberg 加速与外推表的生成Romberg 积分是在复化梯形的基础上做 Richardson 外推。代码包里通常用一个二维数组R存储外推表function [I, R] romberg_integration(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) 0.5 * R(k-1,1) 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); endtol是收敛容差max_iter是最大迭代次数。外推公式里的4^(j-1) - 1是 Romberg 算法的核心每做一次外推误差阶提高两阶。代码包里一般会打印出完整的R表让你看到对角线元素如何快速收敛。4.3 列主元 Gauss 消去法选主元的判断逻辑线性方程组部分列主元消去法是重点。直接 Gauss 消去法在遇到主元为零时会除零崩溃列主元通过行交换避免这个问题。function x gauss_elimination_pivot(A, b) n size(A, 1); Ab [A, b(:)]; % 增广矩阵 for k 1:n-1 % 选主元 [~, idx] max(abs(Ab(k:n, k))); idx idx k - 1; if Ab(idx, k) 0 error(矩阵奇异); end % 交换行 if idx ~ k Ab([k, idx], :) Ab([idx, k], :); end % 消去 for i k1:n factor Ab(i, k) / Ab(k, k); Ab(i, k:n1) Ab(i, k:n1) - factor * Ab(k, k:n1); end end % 回代 x zeros(n, 1); x(n) Ab(n, n1) / Ab(n, n); for i n-1:-1:1 x(i) (Ab(i, n1) - Ab(i, i1:n) * x(i1:n)) / Ab(i, i); end end选主元那三行是关键max(abs(Ab(k:n, k)))找当前列从第k行到最后一行的最大绝对值idx是相对下标加k-1转成绝对下标。交换行用 Matlab 的矩阵索引一行搞定。消去步骤里factor是乘子Ab(i, k:n1)一次性更新整行比逐元素循环快。参数说明A是系数矩阵b是右端向量。如果A接近奇异max返回的主元可能很小但不为零这时候结果虽然不报错但误差很大。代码包里一般会额外算一下条件数cond(A)并打印出来提醒你注意。4.4 LU 分解的紧凑格式与乘子存储LU 分解把A拆成下三角L和上三角U适合多次求解同系数矩阵不同右端项的场景。代码包里通常用紧凑格式把L和U存在同一个矩阵里function [L, U] lu_decomposition(A) n size(A, 1); LU A; for k 1:n-1 if LU(k,k) 0 error(主元为零需要列主元); end for i k1:n LU(i,k) LU(i,k) / LU(k,k); % 乘子存入下三角 LU(i,k1:n) LU(i,k1:n) - LU(i,k) * LU(k,k1:n); end end L eye(n) tril(LU, -1); U triu(LU); endtril(LU, -1)取严格下三角部分加上单位矩阵得到Ltriu(LU)取上三角包括对角线得到U。这个实现没有选主元所以只适用于顺序主子式全不为零的矩阵。代码包里一般会用一个对称正定矩阵做例子保证分解不会失败。提示如果你要处理一般矩阵把列主元的逻辑加进去同时记录置换矩阵P分解变成 ( PA LU )。代码包里可能没有这个版本但你可以基于上面的 Gauss 消去法自己改。5. 避坑与排查跑这份代码时最容易翻车的五个地方5.1 现象插值曲线在端点剧烈震荡误差大得离谱原因节点数取得太多用了高次多项式插值。等距节点下的高次插值会出现 Runge 现象这是理论上的必然不是代码写错了。解决把节点数降到 10 以内或者改用分段插值比如三次样条。如果题目要求必须用高次考虑改用 Chebyshev 节点代码包里可能没有现成实现但你可以把x linspace(a, b, n)换成x (ab)/2 (b-a)/2 * cos((2*(1:n)-1)*pi/(2*n))。5.2 现象复化 Simpson 的结果和梯形差不多没有明显精度提升原因n取了奇数或者n太小。Simpson 公式要求n为偶数如果代码里没做检查奇数n会导致索引错位结果完全错误。解决在调用前确认mod(n, 2) 0。如果n已经是偶数但精度不够把n翻倍再试。Simpson 的误差是 ( O(h^4) )n从 10 加到 20误差应该降到原来的 1/16 左右。5.3 现象Gauss 消去法报“矩阵奇异”但矩阵行列式明明不为零原因主元在消去过程中变成了零或者极小值直接除法导致数值爆炸。这是顺序 Gauss 消去法的固有缺陷不是矩阵本身奇异。解决改用列主元消去法。如果已经用了列主元还报错检查矩阵是不是接近奇异算一下cond(A)。条件数超过 ( 10^{12} ) 时双精度浮点数已经无法保证结果有效。5.4 现象最小二乘法拟合高阶多项式时系数完全不对原因法方程 ( A^T A ) 的条件数是 ( A ) 的条件数的平方高阶时病态严重。m 8时条件数可能已经到 ( 10^{10} ) 以上。解决不要直接解法方程用c A \ y让 Matlab 内部做 QR 分解。或者改用正交多项式基比如 Legendre 多项式代码包里可能没有但你可以用legendreP函数需要 Symbolic Toolbox或者自己递推生成。5.5 现象Romberg 积分不收敛迭代到最大次数还没达到容差原因被积函数在区间内不光滑或者有奇点。Romberg 外推假设函数足够光滑如果函数只有连续但不可导外推表不会快速收敛。解决先检查被积函数在区间内是否连续。如果有奇点考虑变量替换或者改用自适应积分。如果函数光滑但收敛慢把tol放宽到1e-6试试有时候是容差设得太苛刻。6. 进阶用法把这份代码改造成自己的实验报告生成器跑通代码只是第一步真正省时间的是把它改造成能自动生成实验报告的工具。我一般会做三件事参数化输入、自动绘图、结果导出。参数化输入的意思是把主脚本里硬编码的n、a、b、f全部提到脚本开头用一个结构体params管理。这样换一道题只需要改结构体不用翻遍整个文件。params.n 10; params.a 0; params.b 1; params.f (x) sin(x) x.^2; params.method simpson;自动绘图的部分把误差曲线和结果对比图封装成一个函数每次算完自动调用function plot_results(x, y_true, y_approx, title_str) figure; plot(x, y_true, b-, LineWidth, 1.5); hold on; plot(x, y_approx, r--, LineWidth, 1.5); legend(真实值, 近似值); title(title_str); xlabel(x); ylabel(y); grid on; % 误差图 figure; semilogy(x, abs(y_true - y_approx), k-); title([误差曲线 - , title_str]); xlabel(x); ylabel(绝对误差); grid on; end结果导出用writematrix或者save存成.mat文件方便后续写报告时直接引用数据。如果你要交 Word 版报告可以用actxserver调 Word 接口自动插入表格和图但这属于另一个话题了。改造点原始代码改造后收益参数管理硬编码在脚本里结构体统一管理换题不翻文件绘图手动调用 plot封装成函数自动出图省去重复代码结果保存无自动存 .mat 和 .csv写报告直接引用误差分析肉眼对比自动算最大误差和收敛阶数据更可信最后说一个我自己的习惯每次跑完一组参数我会把n、h、最大误差、收敛阶记在一个表格里跑三组不同n就能算出实际收敛阶。比如复化 Simpson 的理论收敛阶是 4如果实测只有 2那说明代码里有 bug 或者函数不光滑。这个验证步骤花不了两分钟但能帮你避开“看起来对但其实错”的坑。从那以后我每次交实验报告前都强制走一遍收敛阶验证再也没被助教打回来过。希望帮到你。本文还有配套的精品资源点击获取