企业官网建设流程全解析

简介：这份数值分析插值算法实验报告及MATLAB代码，面向正在学习数值分析、计算方法或MATLAB编程的本科生与自学者，围绕拉格朗日插值与牛顿插值的编程实现、结果对比和误差分析展开。内容涵盖自定义拉格朗日与牛顿插值函数，针对 f(x)=2x³+x²+1、lnx、1/(1+25x²) 等函数，验证插值节点增多时曲线变化，展示牛顿差商表，并借助函数图像和函数值表格比较两种方法，同时检验拉格朗日余项定理。压缩包为1个PDF文件，约1.3MB，方便直接阅读和打印。目前已有2666人学习下载。读者可据此理清插值公式推导、MATLAB循环与矩阵运算、差商递推和误差估计的实现思路，也可用于课程实验报告撰写、期末复习与编程实践参考。

1. 一组实验数据摆上桌，插值算法才算真正开始

一组实测数据往往只有十几个点，采样成本高、实验不可重做，却要在第 7 个点附近读出中间值，这时候插值算法就是唯一能给出结果的工具。数值分析把插值放在课程中段，不是因为它公式难，而是它把节点、误差余项、病态和稳定性这几件事一次性摊开：拉格朗日给的是解析表达式，牛顿差商给的是可增量更新的系数，三次样条给的是工程上真正能用的光滑曲线。围绕插值算法的 MATLAB 实现，从基函数、差商表一直写到节点选择与样条边界，再落到实验报告里能直接复现的误差表和收敛阶验证，是一条完整的动手路线。适合正在做数值分析实验、需要同时交代码和报告的人，也适合手头有离散标定数据、想把它变成可查询曲线的工程师。

2. 拉格朗日插值与牛顿差商：两种写法、一套结果

插值的理论地基是唯一性：在 n+1 个互异节点上，次数不超过 n 的插值多项式存在且唯一。正因为唯一，拉格朗日形式和牛顿形式算出来的多项式本质是同一条曲线，差别只在系数怎么组织、求值怎么加速、后续加节点时要不要全部重算。理解这一点，才不会在实验报告里把两种方法写成两个不同的算法。

误差余项是插值实验里最该写进报告的一行：R(x) = f^(n+1)(ξ)/(n+1)! · ∏(x - x_i)，ξ 落在包含 x 与全部节点的区间内。它把误差拆成两半，前半段由函数本身的光滑度决定，后半段由节点分布决定，第 3 章要处理的就是后半段。

2.1 拉格朗日基函数的向量化实现

基函数 l_k(x) 是一个连乘，直接双重循环写出来最直观，但 MATLAB 里每次只对第 k 个基点做向量化内积，能避免对每个待求点单独循环，几十行数据下性能差不了多少，可读性却好很多。

function yy = lagrange_interp(x, y, xx) % x, y : 已知节点与函数值（同长度向量，节点互异） % xx : 待求点，可为向量 % yy : 插值多项式在 xx 处的取值 x = x(:); y = y(:); xx = xx(:).'; % 统一成行向量便于广播 n = numel(x); yy = zeros(size(xx)); for k = 1:n Lk = ones(size(xx)); % 第 k 个基函数 l_k(xx) for j = 1:n if j ~= k Lk = Lk .* (xx - x(j)) / (x(k) - x(j)); end end yy = yy + y(k) * Lk; % 加权求和 end end

逻辑说明：内层循环维护的是单个基函数在整条待求点向量上的取值，外层累加 y(k)·l_k，正好对应拉格朗日公式的线性组合。参数上，x 与 y 必须等长且 x 互异，否则分母会出现零；xx 允许任意长度，甚至可以超出节点区间，但超出区间属于外推，误差余项里的 ξ 不再有保证，实验报告里应单独说明。

调用方式是xx = linspace(min(x), max(x), 500); yy = lagrange_interp(x, y, xx);，画图时把节点用散点叠加在曲线上，plot(xx, yy); hold on; plot(x, y, 'o')这几行是 MATLAB 画图里最省事的对比手段。

2.2 牛顿差商表：增量计算与 Horner 求值

牛顿形式的价值在于差商可以列表递推，多来一个数据点只补一行，不必重算全部系数。做在线标定、分批读数据的场景时，这一点比拉格朗日方便得多。

function [c, D] = divided_diff(x, y) % 计算牛顿插值的差商系数 c 与完整差商表 D x = x(:); y = y(:); n = numel(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-j+1)); end end c = diag(D).'; % 对角线元素即牛顿系数 end function yy = newton_eval(x, c, xx) % Horner 式（秦九韶）求值：p = c1 + (x-x1)(c2 + (x-x2)(c3 + ...)) x = x(:); xx = xx(:).'; n = numel(c); yy = c(n) * ones(size(xx)); for k = n-1:-1:1 yy = c(k) + (xx - x(k)) .* yy; end end

逻辑说明：差商表的上三角第 j 列是 j-1 阶差商，对角线自左上到右下依次是一阶、二阶……系数，转置后即为从常数项到最高次项的系数。Horner 求值把嵌套乘法从内向外展开，把 n 次多项式的求值复杂度压到 O(n)，同时减少中间变量的舍入。参数上，D 建议保留输出，实验报告里贴上三角差商表能直观展示计算过程，也能在节点加密时观察高阶差商是否剧烈放大——放大明显就是等距节点不稳定的信号。

2.3 两种实现的结果对照与参数说明

在 f(x)=e^x、x∈[0,1] 等距取 6 个点时，两种方法在待求点上的差值应在 1e-14 量级，超过这个量级说明实现有问题。

提示：节点数超过 15 且等距分布时，两种形式都会出现条件数恶化，此时应该换节点分布或直接上样条，不要靠提高数值精度硬扛。

参数上还有两个容易忽略的细节：一是节点向量建议先做sort，差商公式依赖 x 的有序性；二是待求点用列向量还是行向量要统一，MATLAB 的隐式扩展在向量形状不匹配时会报维度错误，而不是默默给出错误结果。

3. 龙格现象、节点选择与误差估计：插值不是点越多越好

把节点加密到 20 个，曲线反而在区间两端翘得更高，这是插值实验里最反直觉的现象。原因在误差余项的节点乘积项 ∏(x - x_i)：等距节点在区间两端分布稀疏、在中间密集，端点的乘积因子随 n 增大而急剧变大，压过了 1/(n+1)! 的收缩速度。换句话说，加密等距节点是在拿阶乘的收益去补端点乘积的亏损，一旦亏损占优，插值就发散。

3.1 龙格现象复现：等距节点为什么越加密越糟

龙格函数 f(x)=1/(1+25x²) 在 [-1,1] 是经典反例，用它跑一遍就能把现象坐实。

f = @(x) 1./(1+25*x.^2); xx = linspace(-1, 1, 2001); for n = [6 10 14 18] x = linspace(-1, 1, n+1); % 等距节点 y = f(x); yy = lagrange_interp(x, y, xx); fprintf('n=%2d, max|err|=%.4e\n', n, max(abs(yy - f(xx)))); end

逻辑说明：n+1 是节点个数，n 是多项式次数，两者别混。运行后一般会看到最大误差从 1e-1 量级随 n 增大而不降反升，端点附近的振荡被放大。参数上，xx 取 2001 个点是经验值，太稀会漏掉端点处的尖峰，太密只是浪费内存，最大误差对采样密度并不敏感，真正敏感的是节点分布。

注意：判断是否出现龙格现象，要看端点 5% 区间内的最大偏差，而不是全区间误差。中间段的误差可能一直在下降，把全区间最大误差掩盖成缓慢收敛。

3.2 切比雪夫节点：把误差余项里的节点乘积项压到最小

切比雪夫零点的分布是把区间端点加密、中间放稀，正好抵消等距节点的缺陷。节点的显式写法是 x_k = cos((2k+1)π/(2n+2))，k=0…n，落在 [-1,1] 上，映射到任意区间 [a,b] 时用 x = (a+b)/2 + (b-a)/2·cos(…)。

m = 10; % 多项式次数 k = (0:m).'; xc = sort(cos((2*k+1)*pi/(2*m+2))); % 切比雪夫零点，升序排列 yc = f(xc); yyc = lagrange_interp(xc, yc, xx); fprintf('Chebyshev n=%d, max|err|=%.4e\n', m, max(abs(yyc - f(xx))));

逻辑说明：sort不是可选项，差商和基函数都要求节点有序，直接使用 cos 的输出是降序的。参数上，m 是多项式次数，节点数为 m+1；同一个 m 下切比雪夫节点的最大误差通常比等距节点低几个数量级，而且随 m 增大稳定下降，直到触及双精度舍入底噪（约 1e-15 量级）后不再改善。

3.3 误差估计与收敛表生成

实测误差和理论余项要能对上。实验报告里最好给出一张节点数对最大误差的表，一行等距、一行切比雪夫，趋势一眼可见。

表内数值只表示量级走向，实际结果随计算平台略有差异，报告里填自己跑出来的数才算数。生成表格的循环可以直接复用 3.1 的框架，把节点分布换成一个匿名函数参数，一行改动就能同时跑两套。收敛阶可以用最小二乘拟合 log(err) 对 n 的斜率来估计，但要注意误差触底后会失真，这时候把最大误差换成均方根误差，或者用更高精度重算。

4. 三次样条与分段插值：工程数据更该用的方案

多项式插值的次数一高就不稳，工程上的出路是把高次换成低次加分段。分段线性最稳但一阶导数不连续，曲线上会有折角；三次样条在每个小区间上是三次多项式，并强制一阶、二阶导数在内节点处连续，兼顾了光滑性和稳定性。代价是需要额外的边界条件把自由度配平，这也是样条实现里最常被写错的地方。

4.1 spline 与 pchip 的接口差异

MATLAB 里最省事的两个入口是 interp1 的 'spline' 和 'pchip' 方法，前者默认 not-a-knot 边界，后者是保形分段三次埃尔米特插值。

x = [0 1 2 3 4 5]; y = [0 0.8 0.9 0.2 -0.6 -1.0]; xx = linspace(0, 5, 501); y_spline = interp1(x, y, xx, 'spline'); % 三次样条，C2 连续 y_pchip = interp1(x, y, xx, 'pchip'); % 保形，C1 连续 y_linear = interp1(x, y, xx, 'linear'); % 分段线性，对照用 plot(xx, y_spline, xx, y_pchip, xx, y_linear, x, y, 'o'); legend('spline','pchip','linear','data');

逻辑说明：三个方法共用同一套调用签名，替换字符串即可切换，便于在同一张图上比较。参数上，x 必须严格递增，重复节点会直接报错；xx 超出 x 的范围时三者都在外推，'spline' 的外推是三阶多项式的自然延伸，可能给出夸张的值，'pchip' 相对保守，但都不建议把外推结果写进正式结论。

4.2 边界条件与 pp 结构的使用

默认的 not-a-knot 边界意味着首尾两个内节点处的三阶导数连续，在多数均匀数据上表现良好。但如果你知道端点的斜率，就应当把边界条件传进去，否则样条在端点附近会略微偏离物理预期。用 spline 返回 pp 结构，可以避免重复求解。

pp = spline(x, y); % 返回分段多项式结构体 d1 = ppval(pp, 2.5); % 求值 dp = fnder(pp, 1); % 一阶导数（需要曲线拟合工具箱） fprintf('x=2.5 处值 %.4f\n', d1);

逻辑说明：pp 结构保存了断点与各段的系数，反复求值不必重解三对角方程组，做参数扫描或实时查询时收益明显。参数上，fnder 依赖 Curve Fitting Toolbox，若环境里没有这个工具箱，可以自行对每段三次多项式求导系数，或者改用diff(ppval(pp,xx))./diff(xx)做数值差分近似。

4.3 非均匀采样、单调数据与外推的边界

传感器标定表通常要求单调，用 spline 有可能在相邻两点之间产生小的非单调起伏，导致反查表出现多个解；这类数据换 pchip 更稳妥。非均匀采样时，样条公式本身不要求等距，但节点间距比值超过 10 会让三对角方程组的条件数变差，实验里可以先用diff(x)看一眼间距分布，必要时对数据做分段处理或重采样。

5. 实验报告的验证闭环：误差表、收敛阶与可复现脚本

实验报告最容易失分的地方不是代码写错，而是结果无法复现、误差没有量化。把验证做成脚本化的闭环，比贴一堆截图有用得多。核心思路是固定解析函数、固定评测网格，只改节点数，然后观察最大误差随节点数的变化趋势，斜率就是经验收敛阶。

f = @(x) 1./(1+25*x.^2); xx = linspace(-1, 1, 2001); ns = [4 8 16 32 64]; errE = zeros(size(ns)); errC = zeros(size(ns)); for i = 1:numel(ns) n = ns(i); xe = linspace(-1, 1, n+1); % 等距 ye = f(xe); errE(i) = max(abs(lagrange_interp(xe, ye, xx) - f(xx))); k = (0:n).'; xc = sort(cos((2*k+1)*pi/(2*n+2))); % 切比雪夫零点 yc = f(xc); errC(i) = max(abs(lagrange_interp(xc, yc, xx) - f(xx))); end p = polyfit(log(ns), log(errC), 1); % 斜率即经验收敛阶 fprintf('empirical order = %.3f\n', p(1));

逻辑说明：polyfit拟合的是误差对数对节点数的斜率，切比雪夫节点在龙格函数上会给出明显负斜率，等距节点则常常给出正斜率，两条曲线的走向差异就是报告里最有说服力的一张图。参数上，ns 取 2 的幂次是为了在双对数图上等距，便于肉眼读数；评测网格 xx 全程固定，不能随节点数变化，否则误差表不可比。

有三点验证细节值得单独写进报告。第一，误差触到 1e-15 后停止下降是舍入误差底噪，不是算法失效，此时应改用更高精度或换更难的测试函数，而不是继续加密节点。第二，用已知多项式做测试时插值应当精确到舍入误差，这是代码正确性的第一道自检，比如取 f(x)=3x⁴-2x+1，任意节点下最大误差都应在 1e-12 以内。第三，把max换成norm(err, 2)/sqrt(numel(err))得到均方根误差，能避免单一异常点主导结论，对含尖峰的函数尤其有用。最后，脚本开头固定随机种子与数据来源，末尾输出机器精度eps和 MATLAB 版本，别人拿到脚本后一次运行就能复现同一张表。

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