仪器只在几个位置留下读数，你却需要知道两个位置之间的值。把点连起来是很自然的想法；问题在于，用直线连、用一个多项式连，还是分段连，会给出不同的答案。

本章把“穿过已知数据”与“接近未知函数”分开检查。插值条件能在节点上核验，节点之间的可靠性则需要误差分析。即使每个节点都一丝不差，整条曲线仍可能很不合适。

让每一个基函数只负责一个数据值

给定 $n+1$ 个互异节点 $x_0,\ldots,x_n$ 和数据 $y_0,\ldots,y_n$，多项式插值要求找到次数不超过 $n$ 的多项式 $p_n$，使

$p_n(x_j)=y_j,\qquad j=0,\ldots,n.$

这里的数据可以来自函数求值，也可以来自测量。若数据来自测量，就要确认测量误差是否可以忽略；若来自程序求值，也要留意计算误差。准确穿过记录下来的点，还不能自动保证准确复现真实函数在节点上的值。

怎样构造这个多项式？我们给第 $j$ 个数据配一个基函数 $L_j$，要求它在自己的节点取 1，在其他节点取 0。为了在其他节点为零，分子放入所有 $x-x_k$，但跳过 $k=j$；为了在 处变成 1，再除以这个分子在 处的值：

$L_j(x)=\prod_{k\ne j}\frac{x-x_k}{x_j-x_k}.$

节点互异保证分母非零。于是

$p_n(x)=\sum_{j=0}^n y_jL_j(x)$

满足全部插值条件：在 $x_i$ 处，除了 $y_iL_i(x_i)=y_i$，其余各项都消失了。这是 Lagrange 形式。

例：三个不等距节点。 取数据 $(0,1),(1,3),(3,1)$。相应基函数为

$$L_0=\frac{(x-1)(x-3)}3,\qquad L_1=-\frac{x(x-3)}2,\qquad L_2=\frac{x(x-1)}6.$$

在 $0,1,3$ 上逐个代入，可以验证每个基函数都是“自己的位置为 1，其余为 0”。组合后得到

$p_2=L_0+3L_1+L_2=1+3x-x^2.$

例如 $p_2(2)=3$。这一步是在推算没有给出的数据；代回三个已知点，只能验证插值条件，不能证明真实对象在 2 处一定等于 3。

在这组指定节点上，每个基函数只保留自己的数据；加权相加便能满足全部插值条件。表格没有规定节点之间的值必须介于 0 和 1。

这样的多项式是否还有另一个？若 $p_n,q_n$ 都满足条件，则 $p_n-q_n$ 次数至多为 ，却有 个互异零点。非零的 次多项式不可能有这么多零点，因此差只能恒为零。刚才的构造给出了存在性，这个零点论证给出了唯一性。

也可以把 $p_n=a_0+a_1x+\cdots+a_nx^n$ 代入所有节点，得到 Vandermonde 系统 ，其中 。它在互异节点下可逆，但这并不表示幂基系数求解总是数值良好。唯一性是精确算术的结论，计算敏感性仍要单独检查。

新增一个点，给旧多项式补一项

现在收到新的数据 $(4,5)$。原来 $p_2(4)=-3$，需要把这里抬高 8，又不破坏 $0,1,3$ 三处已经满足的条件。

乘积 $x(x-1)(x-3)$ 在旧节点全为零，所以可设

$p_3(x)=p_2(x)+c\,x(x-1)(x-3).$

在 4 处代入，$-3+12c=5$，得到 $c=2/3$。这个修正既照顾新点，又完整保留旧点。一般地，新增节点后只需补上

$$p_n(x)=p_{n-1}(x)+c_n\prod_{j=0}^{n-1}(x-x_j),\qquad c_n=\frac{y_n-p_{n-1}(x_n)}{\prod_{j=0}^{n-1}(x_n-x_j)}.$$

把每次补上的项排在一起，就得到 Newton 插值形式：

$p_n(x)=c_0+c_1(x-x_0)+c_2(x-x_0)(x-x_1)+\cdots+c_n\prod_{j=0}^{n-1}(x-x_j).$

系数可以通过差商表计算。约定零阶差商 $f[x_i]=y_i$，并递推定义

$$f[x_i,\ldots,x_{i+k}] =\frac{f[x_{i+1},\ldots,x_{i+k}]-f[x_i,\ldots,x_{i+k-1}]}{x_{i+k}-x_i}.$$

一阶差商就是两点间斜率；更高阶继续比较相邻差商，但分母必须使用这一项所跨越的首尾节点距离，不能总拿相邻间距来除。Newton 系数是 $c_k=f[x_0,\ldots,x_k]$。

为什么这个递推会给出系数？用 $P_-$ 插值 $x_0,\ldots,x_{n-1}$，用 $$ 插值 ，则

$$P(x)=\frac{(x-x_0)P_+(x)-(x-x_n)P_-(x)}{x_n-x_0}$$

在全部 $n+1$ 个节点都取正确值，次数也至多为 $n$，所以由唯一性，它就是所求插值多项式。其 $x^n$ 系数，是 $P_+$ 与 $P_-$ 的最高次系数之差除以 。从一次插值逐阶重复这个论证，恰好得到差商递推；而 Newton 形式最后一项的乘积首项系数为 1，所以该最高次系数就是新增的 。

对前面的三点，差商表如下。每一列放在该差商的最后一个节点所在行：

因此 $p_2=1+2x-x(x-1)$，展开后仍为 $1+3x-x^2$。不同形式描述的是同一个多项式。若只是求值，不必展开全部系数，可以从内到外计算

$p_2(x)=1+x\bigl(2-(x-1)\bigr).$

一般从 $v=c_n$ 开始，依次执行 $v\leftarrow c_k+(x-x_k)v$，，得到 。这是嵌套求值，每次查询只需与次数成正比的工作量。

右侧每一项依赖左侧相邻两项，示意差商表的递推结构。连线只表示依赖关系，实际计算还必须除以所跨首尾节点的距离；左侧曲线不是正文算例的精确图。

下面的实验把基函数、数据权重和差商放在一起。移动一个数据值之前，先预测哪些节点上的结果应保持不变；再增加第四个点，观察新增项怎样在旧节点全部归零。曲线变了，不等于旧的插值条件被破坏。

误差公式里，函数和节点各占一部分

设节点都在 $[a,b]$ 内，数据确实为 $y_j=f(x_j)$，且 $f\in C^{n+1}[a,$。记节点乘积

$\omega_{n+1}(x)=\prod_{j=0}^n(x-x_j).$

对每个 $x\in[a,b]$，存在一个依赖于 $x$ 的 $\xi$，使

$f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x).$

这个式子并非只凭“在节点上都是零”猜出来的。固定一个不是节点的待估位置 $t$，选择

$C=\frac{f(t)-p_n(t)}{\omega_{n+1}(t)},\qquad H(x)=f(x)-p_n(x)-C\omega_{n+1}(x).$

$H$ 在 $n+1$ 个节点及 $t$ 处都为零，共有 $n+2$ 个互异零点。相邻零点之间用一次 Rolle 定理，$H'$ 至少有 $$ 个零点；重复 次，得到某处 。由于 ，而 ，便有 。代回 的定义，就是所求余项。若 本来就是节点，两边均为零，无须作除法。

二次插值余项证明中，辅助函数有三个节点零点及一个待估位置零点。逐次使用 Rolle 定理，可在三阶导数中找到所需的零点。

若已知 $|f^{(n+1)}|\le M$，便能使用上界

$|f(x)-p_n(x)|\le\frac{M}{(n+1)!}|\omega_{n+1}(x)|.$

例：给出一个不靠“看起来很贴近”的界。 用 $0,1/2,1$ 三个节点插值 $e^x$。在 $[0,1]$，三阶导数仍为 $e^x$，所以 $$。在 ，节点乘积为

$\frac14\left(-\frac14\right)\left(-\frac34\right)=\frac3{64},$

因此误差不超过 $e/128\approx0.02124$。这是保证，不是对实际误差的精确计算；未知的 $\xi$ 不应随手设成查询点或区间中点。

次数变大时，一方面分母出现阶乘，另一方面高阶导数和节点乘积也在变化。不能只看见阶乘就断言插值必然越来越准。经典反例是 $f(x)=1/(1+25x^2)$ 在 $[-1,1]$ 上的等距插值：提高次数后，端部可能出现越来越大的振荡。函数很光滑，问题仍会发生；这称为 Runge 现象，即使精确算出插值多项式也不能消除。

Chebyshev 节点到底改善了什么

我们暂时固定次数和高阶导数界 $M$，只尝试压低 $\max|\omega_{n+1}|$。在 $[-1,1]$，一个有用的选择是

$x_j=\cos\frac{(2j+1)\pi}{2(n+1)},\qquad j=0,\ldots,n.$

这是 Chebyshev 多项式 $T_{n+1}$ 的根。这里采用的是不包含端点的根节点；另一种常见选择为 $\cos(j\pi/n)$，它包含端点，是极值节点。两者都向端部聚集，但公式和权重不能混用。

同样使用五个节点，Chebyshev 根节点在两侧的相邻间隔较小。这里列的是保留四位小数的位置和间隔，方框不按距离比例排列。

从 $T_m(\cos\theta)=\cos(m\theta)$ 出发，余弦恒等式给出

$T_{m+1}(x)=2xT_m(x)-T_{m-1}(x),\qquad T_0=1,\quad T_1=x.$

因此 $T_m$ 确实是多项式；当 $m\ge1$，其最高次项系数为 $2^{m-1}$。$2^{-n}T_{n+1}$ 是首项系数为 1、根恰为上述节点的多项式，故

$\omega_{n+1}(x)=2^{-n}T_{n+1}(x),\qquad \max_{[-1,1]}|\omega_{n+1}|=2^{-n}.$

这还是同次数、首项系数为 1 的多项式所能达到的最小最大值。证明只用符号交替：$2^{-n}T_{n+1}$ 在 $n+2$ 个极值位置交替取 $+2^{-n}$ 与 。若另一首项系数为 1 的多项式 在整个区间绝对值都小于 ，那么 在这些位置必须交替变号，从而至少有 个零点；但最高次项相消后，它的次数至多为 ，矛盾。

若区间为 $[a,b]$，把 $[-1,1]$ 平移缩放过去，节点改为

$x_j=\frac{a+b}{2}+\frac{b-a}{2}\cos\frac{(2j+1)\pi}{2(n+1)}.$

每个节点差都多出因子 $(b-a)/2$，所以节点乘积的最大值为 $((b-a)/2)^{n+1}2^{-n}$。

“最优”在这里有准确的对象：它最小化节点乘积的最大值，从而改善使用统一导数界的误差保证。对于某个指定函数，实际误差还取决于导数；不能据此宣称每一点的实际误差都小于任何其他布点。测量位置若已经固定，当然也不能把节点换掉却沿用原来的读数。

下面保持同样的节点数量，对比等距节点和 Chebyshev 根节点。选择 Runge 函数，逐渐增加次数，把注意力放在两个端部；再换一个低次多项式，检查是否能重现它。图旁显示的是密集采样所得的最大误差估计，它可能漏掉采样点之间更高的峰，不能称作严格的连续区间上界。

数据有一点噪声，曲线会被拉动多少

设节点固定，只把数据改成 $y_j+\delta y_j$。Lagrange 形式直接给出

$\delta p(x)=\sum_{j=0}^n\delta y_jL_j(x).$

如果每个数据误差满足 $|\delta y_j|\le\varepsilon$，那么

$|\delta p(x)|\le\varepsilon\sum_{j=0}^n|L_j(x)|.$

和式 $\Lambda(x)=\sum|L_j(x)|$ 描述这个查询位置对数据误差的最坏放大。它在节点上等于 1，在节点之间可能大于 1；因为某些 $L_j$ 可以为负，$\sum L_{}$ 并不等于 。

仍用节点 $0,1,3$。在 $x=2$，三个基函数值为 $-1/3,1,1/3$，所以 $\Lambda(2)=5/3$。若数据扰动分别为 ，它们乘上基函数后方向一致，恰好得到 。这个例子说明上界并非总是松得无法达到。

在查询位置 x=2，三个带符号基函数与所选数据扰动相乘后同向累加，输出变化达到单个数据误差上限的 5/3 倍。

这是插值问题对数据的敏感性，与程序怎样求值还不是一回事。为了高效计算，可以预先求重心权重

$w_j=\frac1{\prod_{k\ne j}(x_j-x_k)}.$

在非节点位置，$L_j(x)=\omega_{n+1}(x)w_j/(x-x_j)$。常数 1 自己就是其插值多项式，因此 。用这个等式消去共同因子，得到

$p_n(x)=\frac{\displaystyle\sum_{j=0}^n\frac{w_jy_j}{x-x_j}}{\displaystyle\sum_{j=0}^n\frac{w_j}{x-x_j}}.$

若 $x$ 恰为节点，就直接返回对应的 $y_j$，不要计算除零。所有权重乘以同一个非零常数，结果不变。对 $0,1,3$，权重为 $1/3,-1/2,1/6$，可以统一乘以 6，用 $$ 计算。

例如在 $x=2$ 查询数据 $(1,3,1)$，分子为 $2\times1/2-3\times3/1+1\times1/(-1)=-9$，分母为 ，商为 3，与前面两种形式一致。靠近节点时不要先把多项式展开成很大的幂基系数，再期待相互抵消恢复小结果。

预计算权重后，每次查询只需线性数量的运算，通常比反复计算全部 Lagrange 乘积更合适。改写公式可以改善求值过程，却不能消除原有的数据放大，也不能让等距 Runge 插值自动收敛。前两章对“问题条件”与“算法稳定性”的区分，在这里又派上了用场。

分段之后，怎样把接头接好

如果节点越来越多，一定要用一个越来越高次的多项式吗？可以换一种做法：在相邻节点之间只使用低次多项式。

最直接的是分片线性插值。在 $[x_j,x_{j+1}]$ 上，令 $h_j=x_{j+1}-x_j>0$，使用

$S(x)=y_j\frac{x_{j+1}-x}{h_j}+y_{j+1}\frac{x-x_j}{h_j}.$

它只需要本段两个数据。若数据来自二阶连续可微的 $f$，且本段 $|f''|\le M$，把一次插值余项用于本段，可得

$|f(x)-S(x)|\le\frac M2|(x-x_j)(x-x_{j+1})|\le\frac{Mh_j^2}{8}.$

最后一步用了两个非负距离之和为 $h_j$，它们的乘积最大为 $h_j^2/4$。这解释了为什么缩小每段宽度有效。代价是折线的斜率通常在节点跳变，若后续计算需要平滑导数，就还不够。

三次样条在每段使用次数至多为 3 的多项式，并要求连接处的函数值、一阶导数、二阶导数连续。设共有 $n$ 段，每段四个系数，共 $4n$ 个未知数。两端插值给 $2n$ 条条件，内部一阶、二阶导数匹配各给 $n-1$ 条，还差两条边界条件。

自然三次样条选择 $S''(x_0)=S''(x_n)=0$。它并不要求端点斜率为零。若端点斜率已知，也可改成指定两端的一阶导数，这称为夹持条件；不同边界条件一般会给出不同曲线。自然条件是一个建模选择，不代表被插值的真实函数恰好在端点具有零二阶导数。

这两段多项式在 x=1 的值、一阶导数与二阶导数分别相同；自然边界另外规定两端的二阶导数为零。

为了构造自然样条，把节点二阶导数记为 $M_j=S''(x_j)$。每段三次式的二阶导数是直线，因此

$S_j''(x)=M_j\frac{x_{j+1}-x}{h_j}+M_{j+1}\frac{x-x_j}{h_j}.$

积分两次，并让两端分别取 $y_j,y_{j+1}$，得到

$$\begin{aligned} S_j(x)={}&\frac{M_j(x_{j+1}-x)^3}{6h_j} +\frac{M_{j+1}(x-x_j)^3}{6h_j}\\ &+\left(y_j-\frac{M_jh_j^2}{6}\right)\frac{x_{j+1}-x}{h_j} +\left(y_{j+1}-\frac{M_{j+1}h_j^2}{6}\right)\frac{x-x_j}{h_j}. \end{aligned}$$

可以当场检查：代入左右端点，三次项与线性修正恰好抵消，分别留下 $y_j$、$y_{j+1}$；二阶求导又回到前一条直线。值与二阶导数的匹配已由这种写法保证，只剩一阶导数要接好。

设 $d_j=(y_{j+1}-y_j)/h_j$。上式在本段左右两端的导数分别为

$$S_j'(x_j)=d_j-\frac{h_j}{6}(2M_j+M_{j+1}),\qquad S_j'(x_{j+1})=d_j+\frac{h_j}{6}(M_j+2M_{j+1}).$$

让第 $j-1$ 段的右端导数等于第 $j$ 段的左端导数，整理后得到

$$h_{j-1}M_{j-1}+2(h_{j-1}+h_j)M_j+h_jM_{j+1} =6(d_j-d_{j-1}),\qquad j=1,\ldots,n-1.$$

再加上自然条件 $M_0=M_n=0$，便可求全部 $M_j$。每条方程只涉及相邻三个未知量，系数矩阵是三对角的；它可以用线性数量的运算求解，不必当成稠密大矩阵处理。

也能从这些方程看出唯一性。假设两组解不同，它们的差满足右端为零的同一系统。选差的绝对值最大的一个内部位置，记最大值为 $E>0$。该行主对角项的绝对值为 $2(h_{j-1}+h_j)E$，两侧项的绝对值之和至多为 $$，不可能抵消。于是只能 ，解唯一。这一步使用了所有间距严格为正。

例：从方程算出真正的分段曲线。 对 $(0,1),(1,3),(2,1)$，两个间距都为 1，斜率分别为 2 与 $-2$。自然条件给 $M_0=M_2=0$，中间方程为 ，所以 。

代入段式，得到

$$S(x)=\begin{cases} 1+3x-x^3,&0\le x\le1,\\ 3-3t^2+t^3,\quad t=x-1,&1\le x\le2. \end{cases}$$

在接头 $x=1$，两侧的值都是 3，一阶导数都是 0，二阶导数都是 $-6$；在 0 与 2，二阶导数都为零。因此它确实满足所有条件。$S(1/2)=2.375$，而穿过同样三点的全局二次式 $1+4x-2x^2$ 在这里等于 2.5。两者都满足数据，采用的模型条件不同；仅凭三个数据无法断定哪一个就是未知真值。

分片计算不等于数据影响只在局部

自然样条的一个查询只用所在段的表达式，但这些表达式的系数来自同一组联立方程。改变一个数据，通常会影响所有段，只是影响大小不同。分片线性插值的系数确实只依赖相邻两点，不能把这个性质直接移给自然三次样条。

可以用一个很小的反例亲手核对。节点为 $0,1,2,3$，原数据全为零，原曲线自然也为零。只把 $x=1$ 的数据改为 $\delta$。两条内部矩方程为

$4M_1+M_2=-12\delta,\qquad M_1+4M_2=6\delta,$

解得 $M_1=-18\delta/5$、$M_2=12\delta/5$。在远处 $[2,3]$ 段，两端的数据仍为零，但 已不为零；代入段式，

$S(x)=\frac{2\delta}{5}\bigl((3-x)^3-(3-x)\bigr),\qquad 2\le x\le3.$

特别地，$S(2.5)=-3\delta/20$，只要 $\delta\ne0$，这里就被改变了。相同数据的分片线性插值在 $[2,3]$ 仍然是零。

下面拖动一个节点，对照分片线性与自然样条。查看接头两侧的导数，再启用“只有一个数据改变”的例子，观察远处曲线与原曲线的差。自然边界约束二阶导数，不能从图中端点看起来平缓，就误认成斜率为零。

数据若含明显噪声，强迫曲线准确通过全部读数，可能把噪声一起画进结果。此时需要依据任务选择低次拟合、平滑或其他模型；最小二乘将在第 7 章展开。当前至少要检查三个问题：读数是否需要被精确满足，关心的区间是否真的被数据覆盖，后续是否需要可靠的导数。出了数据区间的推算属于外推，前面的区间内误差保证不能原封不动搬过去。

把公式变成可以核对的判断

1. 用差商构造通过 $(-1,2),(0,1),(2,5)$、次数至多为 2 的插值多项式，分别用节点回代与唯一性说明结果。若同一位置 0 同时给出 1 和 2，会怎样？

查看解答

一阶差商为 $(1-2)/(0-(-1))=-1$ 和 $(5-1)/(2-0)=2$，二阶差商为 。Newton 系数是 ，所以 ，也可嵌套写为 。回代三点均成立，唯一性由两个候选之差至多二次却有三个互异零点推出。同一点被指定两个不同值时，连普通函数都无法满足；若重复的是同一个值，则只是重复条件，不能当作一个新的独立节点。给定导数的重复节点问题需要另一种插值设定。

2. 给第 1 题新增数据 $(3,13)$。保留旧多项式的表达，只补一项；计算新的 $p(1)$，并说明为何旧节点不受破坏。

查看解答

设 $p_3=x^2+1+c(x+1)x(x-2)$。在 3 处原值为 10，需要补 3，而乘积为 12，所以 。在 1 处，新值是 。新增项在 都为零，因此旧节点保持原值；节点之间的值没有同样的保护。

3. 在 $[0,\pi/2]$ 两端对 $\sin x$ 作线性插值，给出整个区间上的误差上界，并核对中点实际误差是否满足它。

查看推导

$|f''|=|\sin x|\le1$，区间宽 $h=\pi/2$。线性插值界为 $h^2/8=\pi^2/32\approx0.30843$。中点插值值为 ，实际函数值为 ，误差约为 ，小于上界。界用了整个区间的导数上限和距离乘积上限，通常不要求恰好等于实际误差。

4. 在 $[-1,1]$ 使用三个 Chebyshev 根节点，写出节点乘积及其最大绝对值。对照等距节点 $-1,0,1$，比较节点乘积的最大值，能否直接推出任意函数的实际误差都更小？

查看解答

根节点为 $-\sqrt3/2,0,\sqrt3/2$，乘积为 ，最大绝对值为 。等距节点的乘积为 ，在 达到绝对值 ，比 大。这个比较改善统一导数界下的保证；实际余项还含一个依赖位置和节点的导数值，不能作题目最后的无条件推断。

5. 对节点 $-1,0,1$，求 $x=1/2$ 处三个 Lagrange 基函数的值。若每个数据扰动绝对值至多为 $\varepsilon$，这里最大的输出变化是多少？怎样选扰动能达到它？

再对数据 $(1,0,1)$ 使用重心公式，验证同一位置的插值值。

查看解答

基函数为 $x(x-1)/2,1-x^2,x(x+1)/2$，在 $1/2$ 处分别为 $$。绝对值之和为 ，所以变化不超过 。选择数据扰动 ，三项都向正方向贡献，恰好达到该界。若把所有基函数只看成非负权重，就会漏掉这种放大。

6. 构造通过 $(0,0),(1,1),(3,1)$ 的自然三次样条。给出节点二阶导数，算出第一段及 $S(1/2)$，再核对右端斜率是否一定为零。

查看计算

间距为 $h_0=1,h_1=2$，两段斜率为 1 与 0。$M_0=M_2=0$，内部方程 给 。第一段为 ，所以 。第二段以 写为 ；在 ，导数为 ，并不为零，二阶导数才是 。两段在接头的斜率均为 ，二阶导数均为 。

7. 只有节点 $0,1,2,3$ 上的观测，且 $x=1$ 的读数可能有 $0.02$ 的误差。有人说“选自然样条，误差不会传到 $[2,3]$”。用正文的扰动例子反驳；再说明为什么换成分片线性，也不能保证它准确描述未知真函数。

查看分析

自然样条在 $x=2.5$ 的改变为 $-3\delta/20$，若 $\delta=0.02$，改变就是 $-0.003$，并非零。分片线性在远段只用两端读数，所以这一个数据扰动不会传过去；但真函数在两端之间仍可能弯曲。要控制后一种近似误差，需要二阶导数界等额外信息。数据误差传播与用直线近似曲线的误差不能混为一谈。

8

使用重心公式计算等距节点插值，就能消除 Runge 现象。

正确错误