简介这份资源是哈工大数值分析上机实验前四部分的Matlab程序代码合集面向正在学习数值分析课程、需要动手实现经典算法的本科生与自学者。内容覆盖高斯列主元消去法、龙格现象的产生与克服、非线性方程求解以及曲线拟合的最小二乘法等典型实验帮助读者把课本中的数值方法落到可运行的代码上。压缩包共15个文件以13个m脚本为主体另含1个md说明文档和1个asv自动备份文件整体约8KB体积轻巧便于快速查阅与移植。目前已有150人学习下载。读者可借助这些脚本理解算法流程、对照实验结果并在此基础上改写参数、验证误差与精度为课程上机与后续数值计算实践提供直接参考。1. 数值分析上机实验从公式到能跑通的 Matlab 代码很多人学数值分析考试能算特征值、能背 Newton 迭代公式一到上机就卡住——不是不会数学是不知道公式怎么变成矩阵、循环和收敛判据。这个标题指向的是一套围绕数值分析上机实验前四部分的 Matlab 程序代码覆盖插值、数值积分、线性方程组求解、非线性方程求根这几类经典任务。它解决的不是数学懂不懂而是代码怎么写、参数怎么设、结果怎么验证。适合正在做数值分析课程上机、需要一份可对照参考实现的人也适合已经工作、想快速把数值方法落进 Matlab 脚本的工程师。下面按我实际写这类代码的顺序把每个实验拆到能直接复现。2. 插值与拟合拉格朗日、牛顿插值的 Matlab 落地插值是数值分析上机第一部分最常见的题目。核心诉求是给一组离散点构造一个多项式穿过它们再在任意点求值。很多人第一版代码能跑但点一多就数值不稳定或者求值点写死导致复用困难。2.1 拉格朗日插值的向量化写法拉格朗日插值的公式直白每个基函数是一个连乘最后加权求和。新手容易写成双重循环点一多就慢。我一般用向量化把基函数矩阵一次算出来。function yq lagrange_interp(x, y, xq) % x, y: 已知节点与函数值长度相同 % xq: 待求值点可以是向量 n length(x); m length(xq); L ones(n, m); % 基函数矩阵第 i 行对应第 i 个基函数 for i 1:n for j 1:n if j ~ i L(i,:) L(i,:) .* (xq - x(j)) / (x(i) - x(j)); end end end yq y * L; % 加权求和y 是行向量 end逻辑说明外层i遍历基函数内层j遍历除自己以外的节点逐点乘上(xq - x(j))/(x(i) - x(j))。L是n×m矩阵最后y * L一次得到所有求值点的结果。参数上x和y必须等长xq可以是标量或向量。注意x里不能有重复节点否则分母为零。2.2 牛顿插值的差商表与 Horner 式求值牛顿插值比拉格朗日更适合增量加点因为差商表可以复用。差商表是上三角结构第一列是函数值之后每列是低一阶差商。function [c, yq] newton_interp(x, y, xq) % 返回差商系数 c 和求值结果 yq n length(x); D zeros(n, n); D(:,1) y(:); % 第一列是函数值 for j 2:n for i j:n D(i,j) (D(i,j-1) - D(i-1,j-1)) / (x(i) - x(i-j1)); end end c diag(D); % 对角线就是牛顿系数 % Horner 式求值 yq c(n) * ones(size(xq)); for k n-1:-1:1 yq yq .* (xq - x(k)) c(k); end end逻辑说明差商表D的第j列由前一列相邻两项之差除以节点间距得到。对角线元素D(i,i)就是牛顿插值系数。求值用嵌套乘法从最高次系数开始每次乘(xq - x(k))再加下一个系数。参数上x必须严格递增或递减否则差商公式里的分母符号会乱。xq支持向量ones(size(xq))保证维度一致。2.3 龙格现象与节点选择的取舍等距节点做高次插值会出现龙格现象两端剧烈振荡。这不是代码 bug是数学性质。常见做法是改用切比雪夫节点或者分段低次插值。如果上机题目没指定节点我一般先用等距跑一遍看振荡再换切比雪夫节点对比。n 10; k 0:n; x_cheb cos((2*k 1) * pi / (2*(n1))); % 切比雪夫节点落在 [-1,1] x_cheb sort(x_cheb); % 排序后用于差商 y_cheb 1 ./ (1 25 * x_cheb.^2); % 经典龙格函数 xq linspace(-1, 1, 500); [~, yq] newton_interp(x_cheb, y_cheb, xq);参数说明n是多项式次数节点数n1。切比雪夫节点在区间两端更密能压制振荡。注意x_cheb必须排序后再传给newton_interp因为差商公式依赖节点顺序。这一步不做结果会完全错属于血泪经验。3. 数值积分复化梯形、复化 Simpson 与精度对比数值积分上机的核心是给定被积函数和区间用离散采样逼近积分值并控制误差。题目通常要求比较不同方法的收敛阶所以代码要能输出误差随步长变化的表。3.1 复化梯形公式的实现与误差阶验证复化梯形把区间分成n等份每份用梯形面积近似。公式简单但收敛只有二阶。function I trap_comp(f, a, b, n) % 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逻辑说明x是等距节点y是对应函数值。梯形公式里首尾各取一半权重中间全权重。参数上n越大精度越高但误差按O(h^2)下降。验证误差阶时取n 10, 20, 40, 80算误差再取对数比。f (x) exp(x) .* sin(x); a 0; b pi; exact (exp(pi) - 1) / 2; % 解析积分值 ns [10 20 40 80 160]; for n ns err abs(trap_comp(f, a, b, n) - exact); fprintf(n%4d err%.3e\n, n, err); end运行后误差大约每翻倍n降为四分之一说明二阶收敛。如果误差不按这个规律降先检查f是否支持向量输入exp(x).*sin(x)里的点乘不能漏。3.2 复化 Simpson 公式与偶数区间约束Simpson 公式用二次多项式近似收敛四阶但要求子区间数为偶数。function I simpson_comp(f, a, b, n) % 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逻辑说明奇数下标节点权重为 4偶数下标不含首尾权重为 2。参数上n为偶数时节点数n1为奇数正好满足 Simpson 的配对要求。如果题目给的n是奇数要么报错要么自动加一。我一般直接报错避免静默改参数导致结果对不上。3.3 两种方法的误差对比表与选型建议把两种方法放同一张表里看收敛速度选型就清楚了。子区间数 n复化梯形误差复化 Simpson 误差10约 1e-2约 1e-520约 2.5e-3约 6e-740约 6e-4约 4e-880约 1.6e-4约 2e-9表格里的数量级是典型值具体取决于被积函数。选型上如果函数光滑且要求高精度优先 Simpson如果函数只有一阶连续导数或者采样点已经固定梯形更稳。注意 Simpson 对高频振荡函数不一定占优因为四阶收敛的前提是四阶导数有界。4. 线性方程组直接法与迭代法的 Matlab 实现线性方程组求解是数值分析上机的重头戏通常要求实现 Gauss 消元、LU 分解、Jacobi 迭代、Gauss-Seidel 迭代并比较直接法和迭代法的适用场景。4.1 列主元 Gauss 消元为什么必须选主元朴素 Gauss 消元在遇到小主元时会放大舍入误差。列主元每次选当前列绝对值最大的元素换到对角线上。function x gauss_pp(A, b) % 列主元 Gauss 消元 n size(A, 1); 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 if abs(Ab(k,k)) 1e-12 error(矩阵奇异或接近奇异); end for i k1:n m Ab(i,k) / Ab(k,k); Ab(i,k:n1) Ab(i,k:n1) - m * 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逻辑说明Ab是增广矩阵。第k步先在k到n行里找第k列绝对值最大的行换到第k行。消元时用乘子m更新当前行右侧部分。回代从最后一行往上。参数上A必须是方阵b是列向量。1e-12是奇异性阈值太小会漏判太大会误判。4.2 LU 分解与多次右端项复用如果同一个系数矩阵要解多个右端项LU 分解一次、回代多次更划算。function [L, U] lu_naive(A) n size(A,1); L eye(n); U A; for k 1:n-1 for i k1:n L(i,k) U(i,k) / U(k,k); U(i,k:n) U(i,k:n) - L(i,k) * U(k,k:n); end end end逻辑说明L是单位下三角U是上三角。分解后解Axb变成先解Lyb再解Uxy。参数上这个朴素版本没有选主元只适合对角占优或正定矩阵。实际用的时候我会加列主元或者直接调 Matlab 的lu。上机题目如果要求手写至少要在注释里说明没选主元的风险。4.3 Jacobi 与 Gauss-Seidel 迭代的收敛条件迭代法适合大型稀疏矩阵。Jacobi 每步用上一轮的全部值Gauss-Seidel 用本轮已更新的值。function [x, iter] gs_iter(A, b, x0, tol, maxit) % Gauss-Seidel 迭代 n length(b); x x0; for iter 1:maxit x_old x; for i 1:n s A(i,1:i-1) * x(1:i-1) A(i,i1:n) * x_old(i1:n); x(i) (b(i) - s) / A(i,i); end if norm(x - x_old, inf) tol return; end end warning(达到最大迭代次数未收敛); end逻辑说明内层循环里x(1:i-1)是本轮已更新的值x_old(i1:n)是上一轮的值。收敛判据用无穷范数。参数上tol一般取1e-6到1e-8maxit取几百到几千。收敛的充分条件是矩阵严格对角占优或者谱半径小于 1。不满足时迭代可能发散代码里的warning是后悔药别忽略。5. 非线性方程求根二分、Newton 与弦截法的避坑非线性方程求根上机通常要求实现二分法、Newton 法、弦截法并比较收敛速度和初值敏感性。5.1 二分法区间端点异号是硬前提二分法最稳但要求区间端点函数值异号。function [root, iter] bisection(f, a, b, tol, maxit) if f(a) * f(b) 0 error(端点函数值同号二分法不适用); end for iter 1:maxit 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; end逻辑说明每次取中点判断根在左半还是右半。参数上tol是区间半宽阈值不是函数值阈值。二分法线性收敛每步误差减半大约 10 步能到1e-320 步到1e-6。注意f(a)*f(b)可能溢出大数情况下改用sign判断。5.2 Newton 法导数写错是最常见的翻车点Newton 法收敛快但需要导数且对初值敏感。function [x, iter] newton_root(f, df, x0, tol, maxit) x x0; for iter 1:maxit fx f(x); dfx df(x); if abs(dfx) 1e-14 error(导数接近零Newton 法失效); end x_new x - fx / dfx; if abs(x_new - x) tol x x_new; return; end x x_new; end warning(达到最大迭代次数未收敛); end逻辑说明每步用切线零点更新。参数上df必须单独传入不能靠差分近似除非题目允许。1e-14是导数保护阈值。Newton 法是二阶收敛但初值离根太远可能发散或跳到别的根。我一般先用二分法缩小区间再切 Newton。5.3 弦截法不写导数但要用两个初值弦截法用差商代替导数收敛阶约 1.618。function [x, iter] secant_root(f, x0, x1, tol, maxit) for iter 1:maxit f0 f(x0); f1 f(x1); if abs(f1 - f0) 1e-14 error(差商接近零弦截法失效); end x2 x1 - f1 * (x1 - x0) / (f1 - f0); if abs(x2 - x1) tol x x2; return; end x0 x1; x1 x2; end x x1; warning(达到最大迭代次数未收敛); end逻辑说明用两点连线的零点作为新近似。参数上x0和x1是两个初值不能太近否则差商噪声大。弦截法不需要导数适合导数难求的函数但收敛不如 Newton 稳定。6. 上机代码的常见问题与排查清单这一章集中写我在数值分析上机里反复踩过的坑每条按现象、原因、解决来。6.1 矩阵维度不匹配点乘漏写是头号原因现象报错Matrix dimensions must agree或者结果是一个矩阵而不是向量。原因函数句柄里该用.*、./、.^的地方写了*、/、^。解决所有对向量逐元素操作的运算符都加点。检查f (x) exp(x) * sin(x)这种写法x是向量时必须改成exp(x) .* sin(x)。6.2 迭代不收敛先查对角占优再查初值现象迭代几百次还在振荡或者warning刷屏。原因Jacobi 和 Gauss-Seidel 要求矩阵对角占优或谱半径小于 1不满足就发散。Newton 法初值太远也会发散。解决先算A的对角占优条件不满足就换直接法。Newton 法先用二分法把区间缩到根附近再切换。6.3 差商表对角线取错节点顺序必须一致现象牛顿插值结果和拉格朗日对不上。原因差商表依赖节点顺序x没排序或者两个函数用的节点顺序不同。解决插值前统一sort并且x和y要同步排序。切比雪夫节点生成后默认不是递增的必须手动排。6.4 积分误差不降阶检查被积函数是否支持向量化现象n翻倍但误差几乎不变。原因f只支持标量输入f(x)对向量返回标量导致求和结果错误。解决用arrayfun包一层或者把函数改成向量化写法。测试方法f([1 2 3])应该返回三个值。6.5 奇异矩阵误判阈值设得太死现象明明可解的方程组报「矩阵奇异」。原因奇异性阈值1e-12对某些病态矩阵太宽松或太严格。解决用cond(A)看条件数条件数大于1e12时结果本身不可信不是代码问题。阈值可以调到1e-10试试但别为了跑通把阈值调到1e-6那是自欺欺人。7. 把四部分串成一个可复用的数值实验框架单独跑通每个实验只是第一步。真正省时间的是把插值、积分、线性求解、求根封装成统一入口用配置文件控制参数这样换题目只改配置不改代码。我一般建一个numlab目录下面分interp/、quad/、linsolve/、rootfind/每个目录一个主函数加一个测试脚本。% numlab/run_all.m addpath(genpath(pwd)); config.n 20; config.tol 1e-8; config.maxit 1000; % 插值实验 x linspace(0, 1, config.n); y sin(2*pi*x); xq linspace(0, 1, 200); yq lagrange_interp(x, y, xq); fprintf(插值最大误差: %.3e\n, max(abs(yq - sin(2*pi*xq)))); % 积分实验 f (x) exp(x) .* sin(x); I_trap trap_comp(f, 0, pi, config.n); I_simp simpson_comp(f, 0, pi, config.n); fprintf(梯形: %.6f Simpson: %.6f\n, I_trap, I_simp); % 线性方程组实验 A gallery(poisson, 5); % 25x25 稀疏矩阵 b ones(size(A,1), 1); x_direct gauss_pp(full(A), b); fprintf(直接法残差: %.3e\n, norm(A*x_direct - b)); % 求根实验 froot (x) x.^3 - 2*x - 5; dfroot (x) 3*x.^2 - 2; [root, it] newton_root(froot, dfroot, 2, config.tol, config.maxit); fprintf(Newton 根: %.10f 迭代: %d\n, root, it);逻辑说明addpath(genpath(pwd))把子目录都加进路径。config集中管理参数换实验只改这里。gallery(poisson, 5)生成一个经典稀疏矩阵用来测直接法和迭代法。残差norm(A*x - b)是验证解对不对的硬指标比看x本身可靠。参数上config.n同时控制插值节点数和积分区间数实际用的时候建议分开因为两者最优值不同。config.tol对求根和迭代法都适用但积分没有迭代终止条件靠n控制。config.maxit只对迭代法有效。验证方法上每个实验都要有一个已知解析解或已知性质的测试用例。插值用sin函数积分用exp(x)*sin(x)线性方程组用残差求根用x^3 - 2x - 5这个经典方程。没有解析解的时候用两种不同方法交叉验证比如 Newton 和二分法结果一致才放心。进阶技巧上我习惯把每次运行的参数和误差写进一个log.txt格式是时间戳 实验名 参数 误差。这样调参的时候能看出趋势不用反复翻命令窗口。另外Matlab 的profile能看每个函数耗时插值点超过几千时向量化写法和循环写法差距很明显值得测一次。最后说个我自己的教训早期做上机我总想一次把四个实验写在一个脚本里结果一个报错全挂排查半天。后来改成每个实验独立函数加独立测试脚本主脚本只负责调用和打印出问题定位快得多。数值分析上机的代码不需要多聪明需要的是每个函数都能单独验证、每个参数都有明确含义、每次失败都能定位到具体行。希望帮到你。本文还有配套的精品资源点击获取
