简介：这份资源是哈工大数值分析课程上机实验前四部分的Matlab程序代码合集，面向正在学习数值分析、需要动手实现经典算法的本科生与自学者。内容覆盖高斯列主元消去法、龙格现象的产生与克服、非线性方程求解以及曲线拟合的最小二乘法等典型实验，每个实验均配有可运行的脚本，便于对照理论推导验证数值结果。压缩包共15个文件，以13个m脚本为主体，另含1个md说明文档和1个asv自动备份文件，整体约8KB，体积轻量、结构清晰，适合直接导入Matlab调试运行。目前已有150人学习下载。读者可借此熟悉矩阵分解、迭代法、插值与拟合等方法的编程实现，理解舍入误差与截断误差对精度的影响，并参考README梳理实验组织方式，为后续数值计算课程或科研编程打下基础。
1. 数值分析上机实验的 Matlab 代码到底该怎么用
很多人第一次拿到“数值分析上机实验 Matlab 程序代码”这类压缩包时，第一反应是解压、打开、运行、看结果，然后关掉。等到真正要交实验报告、要改参数、要解释误差来源的时候，才发现自己根本没看懂那段代码在干什么。这个标题对应的东西其实很具体：它是一套围绕数值分析前四个经典模块写成的 Matlab 脚本集合，通常覆盖插值、数值积分、线性方程组求解、非线性方程求根这几块。它解决的不是“从零学数值分析”的问题，而是“我已经知道算法名字，但不知道怎么在 Matlab 里把它跑出来、怎么调参数、怎么验证结果对不对”的问题。适合两类人：一类是正在做上机实验、需要一份能跑通的参考实现；另一类是想把课本公式翻译成可执行代码、但被矩阵索引和循环边界卡住的开发者。下面我按实际复现顺序，把这几块拆开讲清楚。
2. 插值与拟合：从拉格朗日到三次样条的代码落地
2.1 为什么插值实验最容易“看起来对、其实错”
插值实验的典型翻车点不是算法写错，而是节点选取和区间划分出了问题。拉格朗日插值在节点增多时会出现龙格现象，两端震荡剧烈，但如果你只打印中间几个点，结果看起来非常漂亮。我一般会先固定一组等距节点，再换切比雪夫节点做对比，这样一眼就能看出差异。Matlab 里实现拉格朗日插值不需要调工具箱，直接双重循环即可，但要注意基函数的分母在节点重合时会除零，实际写的时候要加一个极小量保护或者直接跳过相同节点。
% 拉格朗日插值：x为节点，y为函数值，xx为待求点 function yy = lagrange_interp(x, y, xx) n = length(x); yy = zeros(size(xx)); for k = 1:length(xx) L = ones(size(x)); for i = 1:n for j = 1:n if i ~= j L(i) = L(i) * (xx(k) - x(j)) / (x(i) - x(j)); end end end yy(k) = sum(y .* L); end end这段代码的逻辑很直白：对每个待求点，构造一组基函数 L_i，再和函数值做内积。参数上，x 和 y 长度必须一致，xx 可以是标量也可以是向量。实际跑的时候建议把 xx 设成 200 个等距点，这样画出来的曲线才平滑。如果发现曲线在端点附近剧烈抖动，那不是代码错了，是节点分布的问题，换切比雪夫节点就能缓解。
2.2 三次样条插值的边界条件怎么选
三次样条比拉格朗日稳定得多，但边界条件的选择会直接影响端点附近的表现。Matlab 自带 spline 函数默认使用非扭结边界条件，而很多教材要求自然边界条件，也就是二阶导在端点为零。如果你直接拿 spline 的结果去对教材答案，可能会发现端点附近对不上。我一般会自己写一个三弯矩法的实现，这样边界条件完全可控。
% 三弯矩法三次样条：自然边界条件 function [xx, yy] = cubic_spline_natural(x, y, n_query) n = length(x); h = diff(x); % 构造三对角方程组 mu = zeros(n-1,1); lambda = zeros(n-1,1); d = zeros(n-1,1); for i = 2:n-1 mu(i) = h(i-1) / (h(i-1) + h(i)); lambda(i) = h(i) / (h(i-1) + h(i)); d(i) = 6 / (h(i-1) + h(i)) * ((y(i+1)-y(i))/h(i) - (y(i)-y(i-1))/h(i-1)); end % 自然边界：M1 = Mn = 0 M = zeros(n,1); % 追赶法解三对角 a = [0; mu(2:n-1); 0]; b = [1; 2*ones(n-2,1); 1]; c = [0; lambda(2:n-1); 0]; rhs = [0; d(2:n-1); 0]; for i = 2:n m = a(i) / b(i-1); b(i) = b(i) - m * c(i-1); rhs(i) = rhs(i) - m * rhs(i-1); end M(n) = rhs(n) / b(n); for i = n-1:-1:1 M(i) = (rhs(i) - c(i) * M(i+1)) / b(i); end % 求值 xx = linspace(x(1), x(end), n_query); yy = zeros(size(xx)); for k = 1:length(xx) for i = 1:n-1 if xx(k) >= x(i) && xx(k) <= x(i+1) hi = h(i); t1 = (x(i+1) - xx(k))^3 / (6*hi) * M(i); t2 = (xx(k) - x(i))^3 / (6*hi) * M(i+1); t3 = (y(i) - M(i)*hi^2/6) * (x(i+1) - xx(k)) / hi; t4 = (y(i+1) - M(i+1)*hi^2/6) * (xx(k) - x(i)) / hi; yy(k) = t1 + t2 + t3 + t4; break; end end end end这里的关键参数是 mu 和 lambda，它们由相邻区间长度决定。自然边界条件把首尾弯矩设为零，如果你需要固定斜率边界，就要改 rhs 的首尾项。追赶法解三对角时，b 的第一个和最后一个元素必须是 1，对应边界方程。跑完之后建议把样条曲线和原函数画在一起，看最大误差出现在哪个区间，通常是在节点稀疏的地方。
2.3 拟合和插值的区别在实验报告里怎么体现
很多实验报告把插值和拟合混在一起写，其实两者的目标完全不同。插值要求曲线严格过点，拟合只要求整体误差最小。Matlab 里做最小二乘拟合用 polyfit 就够了，但要注意阶数不能太高，否则和拉格朗日一样会震荡。我一般会先画散点图，再试 1 到 3 阶，看残差平方和的变化趋势。如果阶数增加但残差下降不明显，就说明过拟合了。实验报告里最好把插值和拟合的结果放在同一张图上，用不同线型区分，这样对比一目了然。
3. 数值积分：复化梯形与辛普森的精度对比
3.1 复化梯形公式的实现和误差观察
数值积分实验的核心不是写出公式，而是理解步长和精度的关系。复化梯形公式最简单，但收敛阶只有二阶，也就是说步长减半，误差大约降到四分之一。我一般会写一个循环，把区间等分成 2 的幂次份，分别计算积分值，然后打印误差随步长的变化。
% 复化梯形积分 function I = trap_composite(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)); end参数 n 是区间份数，必须是正整数。实际跑的时候，我会把 n 取 10、20、40、80，然后观察误差是否按 1/n^2 下降。如果误差下降速度明显慢于理论值，那就要检查被积函数是不是在端点附近有奇异点，或者是不是振荡太剧烈导致采样不足。
3.2 辛普森公式的适用条件和常见误用
辛普森公式要求区间份数为偶数，因为它是每两个小区间做一个抛物线拟合。很多人直接拿奇数份数去套公式，结果最后一个区间被漏掉，误差反而比梯形还大。我一般会在函数入口加一个判断，如果 n 是奇数就自动加一。
% 复化辛普森积分 function I = simpson_composite(f, a, b, n) if mod(n, 2) ~= 0 n = n + 1; % 自动调整为偶数 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这段代码里，4 倍系数对应奇数索引点，2 倍系数对应偶数索引点。跑的时候建议和梯形公式用同样的 n，对比误差。通常辛普森在 n 较小时就能达到梯形 n 很大时的精度，这就是高阶收敛的优势。但要注意，如果被积函数只有一阶连续导数，辛普森的收敛阶会退化，这时候梯形反而更稳。
3.3 用龙贝格积分做精度验证
龙贝格积分是梯形公式的外推加速，适合用来验证前两种方法的结果。Matlab 里实现龙贝格表不需要太多代码，核心是不断把步长减半，然后用理查森外推。我一般会打印龙贝格表的前几行，看对角线元素是否快速收敛。如果对角线元素变化很小，说明已经达到机器精度，再算下去没有意义。
% 龙贝格积分表 function R = romberg_table(f, a, b, 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; for m = 2:k R(k,m) = R(k,m-1) + (R(k,m-1) - R(k-1,m-1)) / (4^(m-1) - 1); end end endmax_iter 一般取 6 到 8 就够了，再多会因为浮点舍入误差导致对角线不再下降。跑完之后看 R(k,k) 这一列，如果相邻两个值之差小于 1e-12，就可以认为收敛了。这个表在实验报告里很有说服力，能直接展示外推加速的效果。
4. 线性方程组求解：直接法与迭代法的选择边界
4.1 高斯消元为什么需要列主元
高斯消元是线性方程组实验的必做项，但很多人写出来的版本没有选主元，遇到主元接近零的矩阵就会翻车。我一般会在每步消元前找当前列绝对值最大的元素，交换行之后再消元。这样虽然增加了比较次数，但数值稳定性提升非常明显。
% 列主元高斯消元 function x = gauss_elim(A, b) n = length(b); Ab = [A, b]; for k = 1:n-1 [~, idx] = max(abs(Ab(k:n, k))); idx = idx + k - 1; if idx ~= k Ab([k, idx], :) = Ab([idx, k], :); end for i = k+1:n factor = Ab(i,k) / Ab(k,k); Ab(i, k:n+1) = Ab(i, k:n+1) - factor * Ab(k, k:n+1); end end x = zeros(n,1); x(n) = Ab(n,n+1) / Ab(n,n); for i = n-1:-1:1 x(i) = (Ab(i,n+1) - Ab(i,i+1:n) * x(i+1:n)) / Ab(i,i); end end这段代码里，max 找的是当前列从 k 到 n 的最大值，交换行之后再做消元。回代从最后一行开始，逐行往上算。实际跑的时候，建议构造一个希尔伯特矩阵做测试，这种矩阵条件数很大，不选主元很容易得到完全错误的结果。
4.2 雅可比迭代和高斯赛德尔迭代的收敛条件
迭代法适合大型稀疏矩阵，但收敛性依赖矩阵性质。雅可比迭代要求矩阵严格对角占优，高斯赛德尔的要求稍微宽松一些，但也不是所有矩阵都收敛。我一般会先算谱半径，如果谱半径大于等于 1，就直接换直接法，不要浪费时间。
% 高斯赛德尔迭代 function [x, iter] = gauss_seidel(A, b, x0, 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) = (b(i) - sigma) / A(i,i); end if norm(x - x_old, inf) < tol return; end end endtol 一般取 1e-6 到 1e-8，max_iter 取 1000 左右。注意高斯赛德尔在更新 x(i) 时用的是最新的 x 值，而雅可比用的是上一轮的值，这是两者最大的区别。如果迭代不收敛，先检查矩阵是不是对角占优，再看谱半径。
4.3 条件数对结果的影响怎么量化
条件数是衡量矩阵病态程度的指标，Matlab 里用 cond 函数就能算。条件数越大，右端项微小扰动引起的解的变化越大。我一般会在实验报告里加一个扰动实验：把 b 加一个 1e-6 的随机扰动，看解的变化幅度。如果变化超过 1e-3，就说明矩阵病态严重，直接法也要谨慎使用。
| 矩阵类型 | 条件数量级 | 建议方法 |
|---|---|---|
| 对角占优 | 1e0~1e2 | 迭代法或直接法 |
| 希尔伯特 | 1e4~1e10 | 列主元直接法 |
| 随机稀疏 | 1e2~1e5 | 迭代法加预处理 |
这张表是我自己总结的经验值，具体问题还要具体分析。条件数在 1e10 以上时，双精度浮点已经很难保证有效数字，这时候要考虑正则化或者换用更高精度。
5. 非线性方程求根：牛顿法为什么会在某些点上翻车
5.1 二分法和牛顿法的适用场景对比
二分法保证收敛但速度慢，牛顿法速度快但要求初值足够接近根，而且导数不能为零。我一般会先用二分法缩小范围，再用牛顿法加速。Matlab 里实现二分法很简单，关键是判断区间端点函数值是否异号。
% 二分法求根 function [root, iter] = bisection(f, a, b, tol, max_iter) if f(a) * f(b) > 0 error('区间端点函数值同号，无法保证有根'); end for iter = 1:max_iter c = (a + b) / 2; if f(c) == 0 || (b - a) / 2 < tol root = c; return; end if f(a) * f(c) < 0 b = c; else a = c; end end root = (a + b) / 2; endtol 一般取 1e-8，max_iter 取 100。二分法每次迭代区间减半，大约 30 次就能达到 1e-8 的精度。如果函数在区间内不连续，二分法会失效，这时候要先画图确认根的位置。
5.2 牛顿法的导数计算和收敛判据
牛顿法的迭代公式是 x_{k+1} = x_k - f(x_k)/f'(x_k)。导数可以用解析式，也可以用差商近似。我一般优先用解析导数，因为差商会引入额外误差。收敛判据用相邻两次迭代的差或者函数值的大小都可以，但要注意函数值很小不代表根很准，因为可能遇到平坦区域。
% 牛顿法求根 function [root, iter] = newton_root(f, df, x0, tol, max_iter) x = x0; for iter = 1:max_iter fx = f(x); dfx = df(x); if abs(dfx) < 1e-12 error('导数接近零，牛顿法失效'); end x_new = x - fx / dfx; if abs(x_new - x) < tol root = x_new; return; end x = x_new; end root = x; end这段代码里，导数保护是必须的，否则遇到导数为零的点会直接除零。实际跑的时候，建议把每次迭代的 x 和 f(x) 打印出来，观察收敛速度。如果发现函数值震荡不下降，可能是初值选得不好，或者根附近导数变化太剧烈。
5.3 弦截法作为牛顿法的替代方案
弦截法不需要解析导数，用两个初始点构造割线来近似导数。收敛阶大约是 1.618，比二分法快但比牛顿法慢。适合导数难以解析求出的情况。Matlab 实现时要注意两个初始点不能太近，否则割线斜率不稳定。
% 弦截法求根 function [root, iter] = secant_root(f, x0, x1, tol, max_iter) for iter = 1:max_iter f0 = f(x0); f1 = f(x1); if abs(f1 - f0) < 1e-12 error('割线斜率接近零'); end x_new = x1 - f1 * (x1 - x0) / (f1 - f0); if abs(x_new - x1) < tol root = x_new; return; end x0 = x1; x1 = x_new; end root = x1; endx0 和 x1 的选取会影响收敛速度，一般建议根据函数图像选在根的两侧。如果迭代过程中发现 x_new 跑出定义域，要加一个边界检查，否则函数求值会报错。
6. 避坑与排查：数值分析上机实验里最容易踩的五个坑
6.1 现象：插值曲线在端点剧烈震荡，原因：等距节点加高阶多项式
这是龙格现象的典型表现。很多人以为是代码写错了，反复检查循环和索引，其实问题出在节点分布上。解决方法是换用切比雪夫节点，或者改用分段低次插值。如果实验要求必须用等距节点，那就把多项式阶数控制在 5 以内，再高就不可控了。
6.2 现象：积分结果和理论值差很多，原因：被积函数在区间内有奇异点
复化梯形和辛普森都要求被积函数在区间上足够光滑。如果函数在端点附近趋于无穷，或者导数不连续，公式的收敛阶会退化。解决方法是把奇异点拆出来单独处理，或者换用自适应积分。Matlab 的 integral 函数自带自适应能力，可以用来做对比验证。
6.3 现象：迭代法不收敛，原因：矩阵不满足对角占优或谱半径大于 1
雅可比和高斯赛德尔都有收敛条件，不是所有矩阵都能用。遇到不收敛的情况，先算谱半径，如果大于等于 1，就换直接法。如果矩阵是对称正定的，可以考虑共轭梯度法，收敛性更有保证。不要硬调松弛因子，那只是治标不治本。
6.4 现象：牛顿法迭代几次后发散，原因：初值离根太远或导数接近零
牛顿法是局部收敛的，初值选得不好就会跑到别的根或者发散。解决方法是先用二分法或画图缩小范围，再让牛顿法接手。如果导数在根附近为零，牛顿法会失效，这时候换弦截法或者修正牛顿法。我一般会在代码里加一个迭代次数上限和发散检测，防止死循环。
6.5 现象：实验结果和教材答案对不上，原因：边界条件或归一化方式不同
数值分析教材里的例题往往有特定的边界条件和归一化约定。比如三次样条的自然边界和固定边界结果不同，最小二乘的基函数选择也会影响系数。解决方法是先确认教材用的什么条件，再调整代码。如果实在对不上，就把两种结果都列出来，分析差异来源，这本身就是实验报告的一部分。
7. 把代码变成自己的：参数扫描与误差阶验证的实操技巧
最后一章我想讲一个我反复用的技巧：参数扫描加误差阶验证。很多人跑完代码，得到一个结果就结束了，但真正能让你理解算法的方式是系统地改变参数，观察误差怎么变。具体做法是：对同一个问题，把步长或迭代次数按 2 的幂次递增，记录每次的误差，然后在双对数坐标下画误差随参数变化的曲线，用 polyfit 拟合斜率。如果斜率接近理论收敛阶，说明代码和理论一致；如果偏离，就要找原因。
以复化梯形为例，理论收敛阶是 2，也就是误差和 h^2 成正比。你可以这样写：
% 误差阶验证：复化梯形 f = @(x) exp(x) .* sin(x); a = 0; b = pi; exact = integral(f, a, b); ns = 2.^(2:8); errors = zeros(size(ns)); for i = 1:length(ns) I = trap_composite(f, a, b, ns(i)); errors(i) = abs(I - exact); end loglog(ns, errors, 'o-'); hold on; p = polyfit(log(ns), log(errors), 1); fprintf('拟合斜率: %.4f\n', p(1));这段代码里，ns 取 2 的幂次是为了让步长按 2 倍变化，这样在双对数图上点分布均匀。polyfit 拟合出来的斜率就是收敛阶的估计值。如果斜率接近 -2，说明梯形公式表现正常；如果斜率明显偏小，就要检查被积函数是不是不够光滑，或者积分区间是不是包含了奇异点。
同样的方法可以用在牛顿法上，记录每次迭代的误差，看收敛阶是不是接近 2。如果发现收敛阶只有 1 左右，可能是导数计算不准确，或者根附近函数性质不好。这种验证方式比单纯看结果有没有对更有说服力，因为它直接检验了算法的核心性质。
还有一个习惯：每次写完一个数值算法，我都会用一个小规模问题手工验算一遍。比如 2x2 的线性方程组，手算消元过程，再和代码结果对比。这样能快速发现索引错误或者符号错误。对于插值，我会用一次多项式去插两个点，看结果是不是直线。这些简单测试花不了几分钟，但能省下大量调试时间。
最后说一个我踩过的坑：不要迷信 Matlab 内置函数的结果。内置函数为了通用性，往往做了很多默认选择，比如 spline 的非扭结边界、integral 的自适应策略。这些默认值在大多数情况下没问题，但和教材要求不一致时就会让你困惑。我的做法是先用内置函数快速验证思路，再自己实现一遍，确保每个参数都可控。这样既保证了效率，又保证了理解深度。希望帮到你。
本文还有配套的精品资源，点击获取
