企业官网建设流程全解析

简介：这份资源是哈工大数值分析上机实验前四部分的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_{i=0}^{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, i+1: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-j+1)); end end coeff = table(1,:)'; % 对角线元素即为系数 end

table是一个下三角矩阵，第一列是函数值，第二列是一阶差商，以此类推。最终系数取第一行（或者对角线，取决于存储方式）。参数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), m+1); for k = 0:m A(:, k+1) = 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)); end

n是子区间数，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 开始），系数为 4；y(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); end

tol是收敛容差，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 = 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(abs(Ab(k:n, k)))找当前列从第k行到最后一行的最大绝对值，idx是相对下标，加k-1转成绝对下标。交换行用 Matlab 的矩阵索引一行搞定。消去步骤里factor是乘子，Ab(i, k:n+1)一次性更新整行，比逐元素循环快。

参数说明：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 = k+1:n LU(i,k) = LU(i,k) / LU(k,k); % 乘子存入下三角 LU(i,k+1:n) = LU(i,k+1:n) - LU(i,k) * LU(k,k+1:n); end end L = eye(n) + tril(LU, -1); U = triu(LU); end

tril(LU, -1)取严格下三角部分，加上单位矩阵得到L；triu(LU)取上三角包括对角线，得到U。这个实现没有选主元，所以只适用于顺序主子式全不为零的矩阵。代码包里一般会用一个对称正定矩阵做例子，保证分解不会失败。

提示：如果你要处理一般矩阵，把列主元的逻辑加进去，同时记录置换矩阵P，分解变成 ( PA = LU )。代码包里可能没有这个版本，但你可以基于上面的 Gauss 消去法自己改。

5. 避坑与排查：跑这份代码时最容易翻车的五个地方

5.1 现象：插值曲线在端点剧烈震荡，误差大得离谱

原因：节点数取得太多，用了高次多项式插值。等距节点下的高次插值会出现 Runge 现象，这是理论上的必然，不是代码写错了。

解决：把节点数降到 10 以内，或者改用分段插值（比如三次样条）。如果题目要求必须用高次，考虑改用 Chebyshev 节点，代码包里可能没有现成实现，但你可以把x = linspace(a, b, n)换成x = (a+b)/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 接口自动插入表格和图，但这属于另一个话题了。

最后说一个我自己的习惯：每次跑完一组参数，我会把n、h、最大误差、收敛阶记在一个表格里，跑三组不同n就能算出实际收敛阶。比如复化 Simpson 的理论收敛阶是 4，如果实测只有 2，那说明代码里有 bug 或者函数不光滑。这个验证步骤花不了两分钟，但能帮你避开“看起来对但其实错”的坑。从那以后我每次交实验报告前，都强制走一遍收敛阶验证，再也没被助教打回来过。希望帮到你。

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