拉格朗日、牛顿、样条这些名字刚接触数值计算的同学基本都会背但真到了自己动手算题、写代码、处理一组实验数据的时候往往还是不知道选哪个。我当年学这一章最深的感受就是公式能看懂作业能做完但完全不理解为什么要发明这么多种插值方法更不清楚它们之间到底是什么关系。这章的内容在整门数值计算课里属于承上启下的位置。前面学的误差分析、浮点数、泰勒展开在这一章会第一次真正派上用场后面要学的数值积分、数值微分、常微分方程数值解又全都要依赖“用多项式近似函数”这套思想。所以如果能把这一章吃透后面学什么都顺畅。这篇内容我从这章真正想解决的问题讲起把拉格朗日、牛顿、分段插值、样条、最小二乘这条线串起来再结合我实际用Python跑实验的记录把书上没写的那些坑和心得一并整理出来。1. 这一章到底在解决什么问题1.1 一个问题引发的思考先说一个现象。你手里有一组测量数据比如某个物理实验在不同时间点记录的温度值一共测了8个点。现在你想知道第2.5秒的时候温度大概是多少但你没有那个时刻的实测数据。怎么办最直接的做法是拿尺子把相邻两个点连起来在中间那段直线上找对应的温度值。这就是最朴素的分段线性插值。本质上做的事情是在已知离散数据点之间用一个简单的函数把空隙填上让这个函数经过所有已知点。但如果我想知道整段时间范围内温度的变化趋势而不只是某两个点之间的值分段直线就太粗糙了。于是自然想到能不能构造一个多项式让它恰好穿过所有已知点再用这个多项式来预测任意时刻的数值这就是插值问题的核心。数学上这么定义给定n1个互异的节点x0, x1, ..., xn以及这些节点上对应的函数值y0, y1, ..., yn我们要找一个多项式p(x)使得p(xi) yi 对每个i都成立。然后就用这个p(x)去计算其他位置的近似值。这类问题在工程里太常见了。查表的时候需要估计表格之外的中间值这叫内插需要根据已有数据推测范围之外的趋势这叫外推。无论是机械设计里的应力系数查表还是气象数据插值做等值线图本质上都是同一个数学问题。1.2 为什么偏偏要用多项式理解了要解决的问题下一个自然的问题是为什么逼近函数的时候首选多项式一个重要的原因是Weierstrass逼近定理如果f(x)在闭区间上连续那么对于任意给定的精度都存在一个多项式能以这个精度一致逼近f(x)。通俗地说连续函数都可以被多项式逼近到任意程度。这个定理给了我们底气——用多项式做逼近理论上没有什么不可逾越的障碍。说得更实际一点多项式的优势在于运算简单。加减乘除都能快速计算求导、积分有现成的公式计算量可控。相比之下如果用三角函数或者指数函数去拟合数据光是函数求值本身就要消耗大量计算资源。在一台没有浮点协处理器的老式设备上做气象数据插值多项式的计算效率优势会被放大得非常明显。还有个容易被忽略的点多项式的泰勒展开形式让我们可以很方便地估计误差。利用插值余项公式能明确知道当前逼近的误差有多大从而判断是否需要提高多项式次数或者加密节点。1.3 插值和拟合的边界要先分清学这一章之前很多人容易把“插值”和“拟合”搞混。当年我也是纠结了好久才完全理清这两个概念的区别。插值的关键约束是p(xi) yi多项式必须严格穿过每一个数据点。这适用于数据本身比较可靠、没有明显噪声的情况。最小二乘拟合则允许多项式不穿过数据点只要求整体偏离最小适用于数据本身带有随机误差的工程测量场景。这个区分很重要。如果把带噪声的测量数据硬做插值得到的函数会因为这些噪声而剧烈波动反而失真。后面的章节我会专门讲最小二乘但在理解插值方法的时候心里要先有这根弦。2. 两类核心插值方法拉格朗日与牛顿2.1 拉格朗日插值思路直接适合理论分析拉格朗日插值的思路在不熟悉数学推导的人看来其实是比较直观的。既然要求多项式p(xi) yi那就构造一组特殊的基函数每个基函数只在一个节点上取值为1在其他所有节点上取值为0。比如给定三个点(x0, y0), (x1, y1), (x2, y2)拉格朗日基函数L0(x)要满足L0(x0) 1L0(x1) 0L0(x2) 0。可以这样构造L0(x) ((x - x1)(x - x2)) / ((x0 - x1)(x0 - x2))代入x0、x1、x2验证结果分别为1、0、0正好符合要求。于是插值多项式写成p(x) y0·L0(x) y1·L1(x) y2·L2(x)这个构造逻辑非常清晰每个基函数负责“接管”一个节点的贡献。好处是结构一目了然程序实现也简单每个基函数都可以独立计算彼此之间不干扰。我当时用Python写拉格朗日插值核心代码就十几行先写一个函数计算单个基函数再在主函数里对每个节点累加yi乘以对应的基函数值。实现过程很顺几乎没遇到什么障碍。拉格朗日形式的最大问题在于如果节点数量增加了哪怕只增加一个节点所有基函数都得重新算一遍。因为每个基函数的分母和分子都发生了变化。这在动态调整节点数的场景下非常不经济。2.2 牛顿插值差商递推适合数值计算牛顿插值形式走的是完全不同的路线。它引入差商的概念把插值多项式写成一个逐次累加的形式p(x) f[x0] f[x0,x1](x - x0) f[x0,x1,x2](x - x0)(x - x1) ...其中f[x0,x1]是f在x0和x1上的一阶差商f[x0,x1,x2]是二阶差商它们都可以用递推公式从函数值计算出来。写成这种形式的巧妙之处在于当增加一个新节点时之前已经算好的系数全部保留只需要在原有多项式后面追加一项即可。这个性质在逐步增加插值节点的场景下价值极大。我实际写牛顿插值的时候发现差商表的构建是个递推填表的过程。把函数值放在第一列然后用公式f[xi, xi1] (f[xi1] - f[xi]) / (xi1 - xi)算一阶差商再用一阶差商算二阶差商逐层往上填。这个过程用二维数组存非常直观代码结构也简单。条件数问题上牛顿插值算法本身的运算量大约需要n²/2次除法。在节点数量不大时两种方法差距不明显但n比较大时牛顿差商的递推计算会明显更节省而且代码的自然结构也方便增量计算。2.3 殊途同归两种形式同一个多项式这一点在学的时候一定要想清楚给定同样的n1个节点拉格朗日插值和牛顿插值得到的是同一个多项式只是表达形式不同。原因是插值多项式在“次数不超过n的多项式空间”里是唯一的。假设有两个n次多项式p和q都经过n1个互异节点那么它们的差r p - q也是一个次数不超过n的多项式而且至少有n1个不同的零点。一个n次多项式最多只能有n个零点除非它是零多项式。所以r必须恒等于0p和q完全相同。这个唯一性定理的价值在于不同形式的插值多项式只是同一事物的不同表达我们可以根据使用场景选择最方便的形式而不用担心计算结果不一致。拉格朗日形式适合理论推导和证明牛顿形式适合实际数值计算两者不矛盾。误把两个不同形式的插值多项式当成不同的方法、做了重复的习题是很多初学者绕的弯路。我当时也一度困惑后来用同一个插值函数分别用两种方法编程验证数值结果完全一致才彻底放心。2.4 插值余项逼近到什么程度只知道怎么构造插值多项式是不够的还得知道误差有多大。教科书上给出的是带拉格朗日余项的误差公式R(x) f(x) - p(x) f^(n1)(ξ) / (n1)! · Π(x - xi)i从0到n其中ξ是包含所有节点和x的区间内某个数。这个式子在实际使用中有一个难点ξ的具体位置很难确定所以通常用f^(n1)在区间上的上界M来代替得到误差估计|R(x)| ≤ M / (n1)! · |Π(x - xi)|这个形式的误差估计揭示了几件事。等式里的Π(x - xi)是一个和节点分布有关的量它的大小和x所在位置直接相关。如果x离所有节点都比较近这个乘积就小如果x在节点区间之外这个乘积可以非常大意味着外插的误差会迅速膨胀。所以实际应用里很少对外推结果抱有太高的信任。另一个方面分母上的(n1)!增长得非常快这让我们觉得提高次数可以快速减小误差。但不要忘了分子上的高阶导数上界M也可能随着n增大而剧烈增长尤其是对被逼近函数本身性质不好的情况。这两股力量互相竞争造成了后面要讲的龙格现象。3. 高次插值并不总是更好3.1 龙格现象等距节点的隐患学插值这一章几乎所有人都会撞上龙格现象。教科书里最经典的反例是f(x) 1 / (1 25x²)在区间[-1, 1]上取等距节点做高次多项式插值。随着插值次数提高在端点附近反而出现了越来越大的振荡误差甚至在端点区域发散。我一开始只是把龙格现象当作一个需要记住的结论并没有真正理解为什么会这样。直到自己写了程序、画出了不同次数下的插值曲线才直观感受到这种振荡有多剧烈。n10的时候曲线在端点附近已经出现明显波浪n20的时候端点附近的波峰几乎要把整条曲线掀翻而中间区域的逼近效果却很好。问题出在哪儿从误差公式来看等距节点在高次时使得|Π(x - xi)|在区间端部变得非常大增长速度超过了(n1)!的衰减作用。而f^(n1)(ξ)上界也在增长于是误差被完全压不住地放大了。这个现象的工程启示很明确不要盲目追求高次多项式插值尤其是在节点等距分布的情况下。高次并不意味着精确反而可能更不可靠。3.2 防止振荡的节点选择策略那是不是高次插值就完全不可取也不是。关键问题出在等距节点的分布方式上。如果改变节点分布让节点在区间端点附近更密集、中间更稀疏龙格现象就能被显著抑制。这里要用到切比雪夫节点。在[-1, 1]区间上n1个切比雪夫节点的位置是xk cos((2k1)π / (2n2))k 0, 1, ..., n这些节点的分布方式恰好是两端密集、中间稀疏从理论分析可知在这个节点集上的多项式插值误差可以更均匀地分布在整个区间上从而避免端点处误差的集中放大。我在实验里用f(x) 1 / (1 25x²)做对比同样是n15等距节点的插值在端点附近振荡严重、误差约在10量级改用切比雪夫节点后误差降到10^-3量级整条拟合曲线几乎贴合原函数。这个差距不是一点点是几个数量级的差异。如果实际问题的区间是[a, b]而不是[-1, 1]用线性变换把切比雪夫节点映射过去即可x (a b) / 2 (b - a) / 2 · xk这一小节的内容让我意识到在数值计算里问题不在于“要不要用高次插值”而在于“用什么节点做高次插值”。选择合适的节点分布往往比单纯降低次数更有效。3.3 高次插值的使用前提即便如此我现在的代码里仍然很少直接用高次多项式插值。原因很简单切比雪夫节点需要事先确定节点位置而这个条件在很多实际场景下并不成立。比如数据是从传感器采样得到的节点的位置是由采样时间决定的没有办法自由调整。再有就算节点分布优化得很好多项式次数一高数值稳定性也会变差。计算过程中的舍入误差会被放大高次多项式的条件数本身就很大甚至可能掩盖插值方法本身的理论精度。所以我的经验是高次多项式插值更适合用于理论推导、构造数值积分公式、以及那些节点位置可以自主调整的特定场景。对于一般的工程数据分段低次插值和样条插值是更稳妥的选择。4. 分段低次插值与样条稳定的工程选择4.1 分段线性插值简单到极致分段线性插值的思想一句话就能说清把每两个相邻节点连成一条直线段整体得到一条折线。这个方法在计算上极其简单每个区间只需要一次线性插值而且局部性很好——修改一个节点的值只会影响相邻两个区间不会波及全局。我最早接触分段线性插值是做温度数据曲线图的时候。把几百个采样点直接连成折线结果就是一条看起来很“硬”的曲线转折处棱角分明。对很多工程场景来说这种粗糙度可以接受因为数据本身精度有限曲线光滑不光滑并不重要。分段线性插值的另一个好处是稳定性。它的误差是O(h²)量级其中h是相邻节点的最大间距。想要提高精度只需要加密节点而不会出现高次插值那种越加密反而越振荡的情况。这一点在实际应用中非常关键。缺点是函数在节点处不可导导数不连续。如果后续需要计算导数或者要求曲线足够光滑分段线性插值就撑不住了。4.2 三次样条插值光滑与稳定的平衡点既想要分段插值的稳定性又想要整体足够光滑三次样条就是一个很好的平衡方案。三次样条在每个小区间上用一个三次多项式但在相邻区间的连接点上要求函数值、一阶导数、二阶导数都连续。这样就得到了一个整体上二阶光滑的分段函数既避免了高次插值的振荡又克服了分段线性的不光滑问题。我第一次真正理解三次样条的价值是在处理一条数据量不大的实验曲线时。用三次样条插值得到的曲线不仅通过了所有数据点而且曲线形态非常自然流畅没有分段折线那种生硬的转折感。三次样条的构造思路是每个区间上的三次多项式有4个未知系数n个区间共4n个未知数。每个区间内部要求插值条件给2个方程区间连接处的函数值连续给n-1个方程一阶导连续给n-1个方程二阶导连续给n-1个方程。总共加起来还差2个方程于是需要加上边界条件常用的有自然边界条件两端二阶导为零或者固定边界条件两端一阶导给定。求解三次样条最终归结为求解一个三对角线性方程组。这个方程组只有约3n个未知数而且三对角结构可以用追赶法非常高效地求解计算量是O(n)量级对于几千个数据点来说完全不是负担。4.3 实际计算中直接调库还是自己实现在实际工程里做三次样条插值通常不必自己从头实现。各类科学计算库都有现成模块Python的scipy.interpolate.CubicSpline和相关样条接口可以直接用MATLAB的csape接口更灵活还能配合不同边界条件生成样条。但学这一章的时候我还是强烈建议自己至少手写一遍三对角方程组的解法。不是为了“造轮子”而是为了真正理解样条的内部结构。十二年之后我遇到一个需要定制边界条件的工业应用场景就是因为在脑海里对样条方程的结构非常清楚才能快速写出适合自己需求的变体版本而不必被现成库的接口限制。如果只是做数据处理直接用现成库也完全没问题。我的建议是学习阶段手写验证工作阶段用库抽查两头都要有数。5. 当插值失灵时最小二乘逼近兜底5.1 数据有噪声时插值是错的插值要求函数严格穿过每个数据点。这个要求本身隐含了一个前提这些数据点是精确的。而在真实的工程测量中数据总是带有噪声的。传感器测量温度有误差测量位置有误差实验记录有误差。如果硬让一条多项式曲线穿过所有带噪声的点曲线就会被噪声的形状带着走产生大量无意义的起伏。这种现象我在处理实验数据时遇到过好多次。有一组标定数据总共20个点直接做插值得到的曲线波动大得离谱后来分析发现是因为个别测量点的误差比较大。改用最小二乘拟合一个低次多项式后曲线变得合理得多虽然不经过每个数据点但整体趋势和物理规律吻合得更好。这就是为什么数值计算课在讲完插值之后紧接着要讲最小二乘逼近。它代表的是另一种思路当数据本身不可靠时放弃“严格穿过每个点”的约束改为“让整体偏离最小”。5.2 最小二乘多项式拟合的思路最小二乘拟合的核心思想是对于一个给定的m次多项式计算它在所有数据点上与真实函数值的误差平方和然后通过选择多项式的系数使这个平方和最小。目标函数可以写成S(a0, a1, ..., am) Σ(yj - (a0 a1xj ... am·xjᵐ))²对每个系数ak求偏导并令其为零得到一组关于系数的线性方程叫作正规方程。这个方程组规模不大m次拟合就是m1个未知数的线性方程组求解起来并不困难。选择拟合次数是个需要权衡的问题。次数太低拟合曲线太粗糙无法反映数据的真实变化趋势次数太高又可能出现过拟合把噪声也拟合进去了。我一般是从1次开始尝试逐步升高观察误差平方和的变化。当误差平方和不再明显下降时就说明已经达到合适的复杂度继续升高次数只会拟合噪声。还有一个实践细节数据点x的数值范围如果很大直接构造正规方程时可能出现严重的条件数问题。我习惯先把x做归一化或者中心化处理再进行拟合。这能让数值稳定性好很多也能避免系数取值过大致使矩阵病态。5.3 插值、拟合如何选择总结我的判断标准数据可靠、节点数量不多且确定时用插值数据量大且包含噪声或者需要寻找整体趋势时用拟合既要求严格穿过关键点又要求曲线光滑时用样条插值。举个最简单的例子要生成一副插值曲线直接对数据点做三次样条插值再在曲线上取点结果又光滑又稳定。但如果要对一组带有流量噪声的测量数据做趋势分析三次样条会把噪声完全复现出来这时候应该改用低次最小二乘拟合。在数值计算里没有一个方法能通吃所有场景。每个方法都有它的适用边界知道什么时候不用这个方法比知道什么时候用它更重要。6. Python实操记录把理论跑一遍6.1 最简实现拉格朗日基函数拉格朗日插值的代码非常精简。教我数值计算的某位老师说过一句话我到现在还记得只要你能写出基函数整个插值就完成了一半。确实如此。核心代码如下def lagrange_basis(x_nodes, k, x): result 1.0 xk x_nodes[k] for j, xj in enumerate(x_nodes): if j ! k: result * (x - xj) / (xk - xj) return result def lagrange_interp(x_nodes, y_nodes, x): total 0.0 for k in range(len(x_nodes)): total y_nodes[k] * lagrange_basis(x_nodes, k, x) return total这段实现的时间复杂度是O(n²)n是节点个数。对于课程实验规模的数据来说完全够用。但如果要做频繁的大量求值建议用牛顿形式把系数先算好再用霍纳法则求多项式值可以把单次求值降到O(n)。用这段代码跑一遍龙格现象观察不同次数下的曲线形态是我强烈建议每个学习者都亲手做一次的事。这个实验只要十几行代码但比看十遍教科书都管用。6.2 牛顿差商表的建立与增量更新牛顿插值的核心数据结构是差商表。实现时用二维数组逐列递推代码逻辑非常清晰def divided_difference_table(x_nodes, y_nodes): n len(x_nodes) table [[0.0] * (n1) for _ in range(n)] for i in range(n): table[i][0] y_nodes[i] for j in range(1, n): for i in range(n-j): table[i][j] (table[i1][j-1] - table[i][j-1]) / (x_nodes[ij] - x_nodes[i]) return table表建好后再构造多项式并求值。系数就是表中第一行的第0~n-1列。求值的时候从后往前累加用类似霍纳法的技巧def newton_interp(x_nodes, y_nodes, x): table divided_difference_table(x_nodes, y_nodes) n len(x_nodes) result table[0][n-1] for k in range(n-2, -1, -1): result result * (x - x_nodes[k]) table[0][k] return result这段代码和拉格朗日版本跑出来的结果在数学上是完全一致的。但仔细观察运算过程会发现牛顿形式把所有基函数计算和线性组合合并成了更少的乘加操作而且差商表只需要构建一次后续任意点的求值都很快。做实验的时候我每次新加一个节点只需要在差商表最后追加一行已有的行全部保留。这个体验和拉格朗日需要全部重算是完全不同的用一次就能明白牛顿形式的工程价值。6.3 龙格现象与切比雪夫节点的对比实验我在实验脚本里分别用等距节点和切比雪夫节点构造插值多项式对f(x) 1 / (1 25x²)做了对比核心验证逻辑用下面几行就可以表示import numpy as np from numpy.polynomial import chebyshev def f(x): return 1 / (1 25 * x**2) # 等距节点 n 16 x_uniform np.linspace(-1, 1, n) y_uniform f(x_uniform) # 切比雪夫节点 k np.arange(n) x_cheb np.cos((2*k 1) * np.pi / (2*n)) y_cheb f(x_cheb) # 后续用牛顿插值分别构造并采样比较误差实验结果显示等距节点在端点附近产生了非常明显的振荡最大误差在1左右切比雪夫节点插值在各点上的误差都小得多曲线和真实函数几乎重合。用一张误差对比曲线图从视觉上就能把两种策略的差距看得明明白白。这个实验给我最大的启发是误差分布比误差大小更值得关注。等距节点插值的问题不只是误差大更关键的是误差分布极不均匀——中间区域很好端点区域糟糕透顶。切比雪夫节点的本质是让误差在整个区间上尽量趋于均匀。6.4 三次样条的快速上手处理三次样条我推荐直接用现成库但前提是理解边界条件的含义。scipy.interpolate.CubicSpline最常见的用法如下from scipy.interpolate import CubicSpline import numpy as np x_data np.linspace(0, 10, 11) y_data np.sin(x_data) cs CubicSpline(x_data, y_data, bc_typenatural) x_fine np.linspace(0, 10, 200) y_fine cs(x_fine)bc_type可以选择natural、clamped等边界条件。默认的自然样条通常足够用但如果数据本身有导数信息clamped条件能利用这些信息让曲线更贴合物理背景。我处理实验曲线时经常根据物理过程的初末状态把端点导数信息填进去效果比自然条件更好。6.5 实操中踩过的坑这里整理几个我实际遇到过的问题都是教科书不会专门提醒的。第一节点必须严格互异。计算凸性时如果两个节点取值完全相同分母就为零程序直接崩溃。我在用算法自动生成节点时因为浮点舍入导致两个节点非常接近虽然没有完全相等但条件数已经很大计算出来的插值结果偏差很大。解决方法是让程序自动检查节点间距低于某个阈值就报错或者重新生成。第二浮点精度对于高次插值影响巨大。用n20的等距节点做插值时由于龙格现象和舍入误差叠加曲线端点附近的振荡加剧得吓人。我当时以为是公式写错了排查了很久最后确认是数值稳定性问题。后来用numpy提升精度以及改用切比雪夫节点问题才缓解。第三外插永远要谨慎。误差公式中的Π(x - xi)在x超出节点区间后会迅速增大外插距离越远误差越大。我在一次应用里试图用插值多项式外推一小段范围结果推测出的数值完全背离物理规律。从那时起凡是顺着数据的整体趋势做预测我一律用最小二乘拟合而不是插值。第四scipy的CubicSpline接口默认不允许外推求值。如果传入的x超出数据范围默认会抛出异常。如果需要外推需要显式设置extrapolateTrue但外推结果的可靠性依然要自己负责。6.6 实现层面的细节建议代码实现插值算法有几个可以提升性能的细节值得专门提一下。用numpy的数组操作批量计算拉格朗日基函数比纯Python循环快得多。比如要计算y_pred把所有基函数一次性用数组表示然后做一次点积。节点数量几百个、采样点几千个的场景下这种向量化可以节省几十倍的时间。牛顿差商表构建时默认用list of list或二维numpy数组都可以。关键在于列的索引顺序要对应“第几阶差商”不要搞混。我写代码时习惯把差商表横着存每一行是同一个节点的0阶、1阶、2阶差商这样索引更直觉出错率也低。对于工程上几十万、上百万的数据点三次样条的线性方程组用追赶法求解是O(n)的完全吃得消。如果数据点再多要考虑样条的空间局部性可能需要对数据分段处理减少单次构建样条的规模。7. 这一章的后续走向插值与多项式逼近学完之后后面几章的内容基本都会和这一章产生联系。数值积分中常使用的复合求积公式本质上是对被积函数做分段多项式插值再积分数值微分同样可以通过构造插值多项式后求导来得到近似公式常微分方程数值解里的多步方法也依赖插值思想来构造。我在学完这一章之后形成的一个习惯是遇到任何新的数值方法先想清楚它背后的逼近思想是什么误差主要从哪里来适用范围是什么。掌握了这些具体公式反而不需要死记用到的时候查一下就能快速上手。如果想把这一章的内容再往外扩展可以看看Bézier曲线和B样条它们是对样条思想的进一步延伸在计算机辅助几何设计中有着非常广泛的应用。也可以接触一下径向基函数插值在多维散乱数据插值里很常用。这些内容都是在多项式插值的基础上发展起来的底子打好了再往上走会顺畅很多。
