【科学计算】插值理论

插值理论及其应用

插值(interpolation)是一种典型的化繁为简的方法，利用既定函数\(f(x)\)在某个区间内若干点（插值节点，往往还是等距点）的函数值，作出适当的较简单的特定函数\(\varphi(x)\)，在这些点上取已知值，在其他点上用特定函数\(\varphi(x)\)的值作为函数\(f(x)\)的近似值。往往特定函数\(\varphi(x)\)会取为最简单的多项式，称之为插值多项式。

插值主要应用于既定函数\(f(x)\)难以计算的时候，使用特定函数\(\varphi(x)\)将大大简化计算量，同时由于特定函数与目标函数在一定点上是相交的，所以误差可能不会太大。

常用的插值法有：多项式插值（又可细分为Lagrange插值和Newton插值）、Hermite插值、分段插值和样条插值。本文将对以上的插值法进行介绍。

目录

• 插值理论及其应用
  • 一、多项式插值
    • 多项式插值存在吗？
    • Lagrange插值法
    • Newton插值法
    • 误差分析
  • 二、分段插值
    • Runge现象与Faber定理
    • 分段线性插值
    • 两点三次Hermite插值
    • Hermite插值
    • 分段三次Hermite插值
  • 三、样条插值
    • 样条插值特点
    • 样条插值函数的一般求解
    • 自然边界求解

一、多项式插值

多项式插值，指的是在给定函数\(f(x)\)在\(n+1\)个互不相同的点\(a=x_0<x_1<\cdots<x_{n+1}<x_n\)上的函数值\(f(x_i)=y_i\)时，求一个\(n\)次多项式\(p_n(x)\)，使得\(p_n(x_i)=y_i\)。由于多项式是最简单的函数，所以多项式插值的操作也是最简洁的。

在多项式插值中，往往会涉及这样的一个函数：设\(x_0,\cdots,x_n\)为插值节点，

\[\omega_n(x)=\prod_{i=0}^{n}(x-x_i). \]

由于此函数比较常见，所以在前面预先定义。

多项式插值存在吗？

首先考虑多项式插值的存在性，也就是说我们要能够找到一个\(n\)次多项式满足条件，自然用待定系数法最简单，不妨设

\[p_n(x)=a_0+a_1x+\cdots+a_nx^n=\sum_{i=0}^{n}a_ix^n, \]

则结合\(p_n(x_i)=y_{i}\)这\(n+1\)个条件，我们有

\[\begin{array}{c} \\ a_0+a_1x_0+\cdots+a_nx_0^n=y_0, \\ a_0+a_1x_1+\cdots+a_nx_1^n=y_1, \\ a_0+a_1x_2+\cdots+a_nx_2^n=y_2, \\ \vdots \\ a_0+a_1x_n+\cdots+a_nx_n^n=y_n. \end{array} \]

考虑我们需要解决的问题，在实际应用时，\(x_1,x_2,\cdots,x_n\)以及\(y_1,y_2,\cdots,y_n\)都是已知量，我们需要求出插值多项式的具体表达式，也就是求出\(a_0,a_1,\cdots,a_n\)这\(n+1\)个数的值，显然这是一个线性方程组，因此可以改写为

\[\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{bmatrix}\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n} \end{bmatrix}, \]

系数矩阵是一个方阵（这也是我们要求插值多项式是\(n\)次的原因），因此可以用Cramer法则直接判断解的存在唯一性，又注意到系数矩阵的行列式恰是范德蒙行列式，结合插值节点各不相同的条件，有

\[|A|=\prod_{0\le i<j\le n}^n (x_j-x_i)\ne 0, \]

这就表明线性方程组有唯一解，即插值多项式有唯一的系数，因而是唯一存在的。这样，我们就解决了插值多项式的存在性问题，还顺便证明了其唯一性。

理论上，可以直接用Cramer法则求出每一个\(a_i\)了，但是由于这样做，除去范德蒙行列式还需要计算\(n+1\)个\(n+1\)阶行列式，计算量十分庞大，在插值节点较多的时候可行度低，因此需要寻找其他求解插值多项式的方法。典型的方法是Lagrange插值法和Newton插值法。

Lagrange插值法

Lagrange插值法求插值多项式的核心是“基函数”，即将目标插值多项式\(p_n(x)\)分解成若干个基础函数的加和。如果我们能够找到\(n+1\)个这样的\(n\)次多项式\(l_i(x)\)，它们在\(x=x_i\)时取值为\(1\)，在其他插值节点处取值为\(0\)，那么将它们所叠加得到的\(n\)次多项式

\[L_n(x)=\sum_{i=0}^{n} y_il_i(x), \]

显然可以满足\(L_n(x_i)=y_i\)的条件，再结合插值多项式的存在唯一性，它就是插值多项式，这种方法求得的多项式\(L_n(x)\)就称为Lagrange插值多项式。

关键在于如何求出这些基函数，应当从其性质分析它的结构。由于\(l_i(x_j)=\delta_{ij}\)，因此它必然含有以下因式：

\[(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})(x-x_n), \]

这保证代入每一个不等于\(x_i\)的插值节点时都能让基函数的值为\(0\)，再结合\(l_i(x_i)=1\)的条件，可知

\[l_i(x)=\prod_{j\ne i}\frac{(x-x_j)}{(x_i-x_j)}. \]

这样的基函数也是唯一的，因为我们相当于用同样的插值节点进行了一次更简单的插值（是否可以发现\(l_i(x)\)的形式与开篇所提到的\(\omega_n(x)\)有所类似？实际上可以用\(\omega_n(x)\)表示之）。由此，我们得到Lagrange插值法得到的插值多项式：

\[L_n(x)=\sum_{i=0}^n y_il_i(x)=\sum_{i=0}^{n}y_i\prod_{j\ne i,1\le j\le n}\left(\frac{x-x_j}{x_i-x_j}\right). \]

现在我们用Lagrange插值法求函数\(f(x)=x^4\)的插值多项式：

\[\begin{array}{|c|c|c|c|} \hline x & -1 & 0 & 1 & 2 \\ \hline f(x) & 1 & 0 & 1 & 16 \\ \hline \end{array} \]

插值多项式应当是三次的，且

\[l_0(x)=\frac{(x-0)(x-1)(x-2)}{(-1-0)(-1-1)(-1-2)}=-\frac{x(x-1)(x-2)}{6}, \\ l_1(x)=\frac{(x-(-1))(x-1)(x-2)}{(0-(-1))(0-1)(0-2)}=\frac{(x+1)(x-1)(x-2)}{2}, \\ l_2(x)=\frac{(x-(-1))(x-0)(x-2)}{(1-(-1))(1-0)(1-2)}=-\frac{(x+1)x(x-2)}{2}, \\ l_3(x)=\frac{(x-(-1))(x-0)(x-1)}{(2-(-1))(2-0)(2-1)}=\frac{(x+1)x(x-1)}{6}, \]

于是

\[L_3(x)=1l_0(x)+0l_1(x)+1l_2(x)+16l_3(x)=2x^3+x^2-2x. \]

此多项式的求解可以使用多项式计算器。

Newton插值法

要使用Newton插值法，就必须介绍一个重要概念：差商(difference quotient)。定义\(f[x_0,x_1,\cdots,x_k]\)为\(f(x)\)关于节点\(x_0,\cdots,x_k\)的\(k\)阶差商，递归定义为

\[f[x_i]=f(x_i), \\ f[x_i,x_j]=\frac{f[x_j]-f[x_i]}{x_j-x_i}, \\ f[x_i,x_j,x_l]=\frac{f[x_j,x_l]-f[x_i,x_j]}{x_l-x_i},\\ \vdots \\ f[x_0,x_1,\cdots,x_k]=\frac{f[x_1,\cdots,x_k]-f[x_0,\cdots,x_{k-1}]}{x_k-x_0}. \]

可以观察到，每一个差商都是一个分式，其中分母是首尾两个节点的差，分子是除去第一个节点的差商和除去最后一个节点的差商的差。差商具有两个重要性质，一是对称性，即\(f[x_0,\cdots,x_k]\)的值与其中节点的排列顺序无关；二是以下事实：

\[\begin{aligned} &\quad f[x_0,x_1,\cdots,x_k]\\ &= \sum_{i=0}^k\frac{f(x_i)}{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})(x_i-x_k)}\\ &= \sum_{i=0}^k\left(f(x_k)\cdot\prod_{j\ne i,0\le j\le k}\frac{1}{x_i-x_j} \right). \end{aligned} \]

这个事实可由数学归纳法证明，这里略去证明的细节，直接给出利用Newton法求插值多项式的推论：若插值节点为\(x_0<x_1<\cdots<x_n\)，则Newton法给出的插值多项式为

\[\begin{aligned} N(x)&=f[x_0] \\ &+ f[x_0,x_1](x-x_0) \\ &+ f[x_0,x_1,x_2](x-x_0)(x-x_1) \\ &+\cdots \\ &+ f[x_0,x_1,\cdots,x_{n+1}](x-x_0)(x-x_1)\cdots(x-x_n). \end{aligned} \]

与Lagrange插值法相比，Newton插值法似乎显得有些“不讲道理”，但它显得更为简洁美观，因此也更适合我们用于求插值多项式。剩下的问题，就是我们如何求上述表达式中出现的差商，差商表是常用的求差商工具，下表中，我们用\(\partial^n\)表示\(n\)阶差商。

\[\begin{array}{|c|} \hline k & f(x_k) & \partial^1 & \partial^2 & \partial^3 & \cdots \\ \hline x_0& \underline{f(x_0)} \\ \hline & & \underline{f[x_0,x_1]} \\ \hline x_1& f(x_1)&& \underline{f[x_0,x_1,x_2]} \\ \hline & & f[x_1,x_2] && \underline{f[x_0,x_1,x_2,x_3]} \\ \hline x_2& f(x_2)&& f[x_1,x_2,x_3] && \cdots \\ \hline & & f[x_2,x_3] && f[x_1,x_2,x_3,x_4] \\ \hline x_3& f(x_3)&& f[x_2,x_3,x_4] \\ \hline & & f[x_3,x_4] \\ \hline x_4& f(x_4) \\ \hline \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \end{array} \]

有了差商表，就可以很快地计算各阶差商，同时下划线部分就是Newton多项式中，各项前面的系数。

现在我们用Newton插值法求函数\(f(x)=x^4\)的插值多项式：

\[\begin{array}{|c|c|c|c|} \hline x & -1 & 0 & 1 & 2 \\ \hline f(x) & 1 & 0 & 1 & 16 \\ \hline \end{array} \]

构造差商表如下：

\[\begin{array}{|c|} \hline k & f(x_k) & \partial^1 & \partial^2 & \partial^3\\ \hline -1& \underline{1} \\ \hline & & \underline{-1} \\ \hline 0& 0&& \underline{1} \\ \hline & & 1 && \underline{2} \\ \hline 1& 1&& 7 &\\ \hline & & 15 && \\ \hline 2& 16 && \\ \hline \end{array} \]

因此Newton插值多项式为

\[\begin{aligned} N_3(x)&=1 \\ &+(-1)(x-(-1)) \\ &+1(x-(-1))(x-0) \\ &+2(x-(-1))(x-0)(x-1) \\ &=2x^3+x^2-2x \end{aligned} \]

误差分析

需要注意的是，由于Lagrange插值法和Newton插值法得到的都是\(n\)次多项式，结合我们在第一部分论证的存在唯一性，它们得到的结果一定是一致的。

我们需要用插值函数作为被插函数的估计，所以误差分析很有必要，如果误差较大，那么插值函数就不适合作为估计。

记\(R(x)=f(x)-p_n(x)\)为用插值函数估计被插函数的截断误差，定理表明，取任一包含所有插值节点的闭区间\([a,b]\)（常常就是\([x_0,x_n]\)），，如果\(f(x)\in C^{n+1}[a,b]\)（即\(f^{(n+1)}(x)\)在\([a,b]\)上连续），则\(\forall x\in [a,b]\)，\(\exists \xi \in (a,b)\)，使得

\[R(x)=f(x)-p_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_n(x). \]

这个截断误差和Taylor展开的Lagrange余项很像，其中的\(\xi\)都是与\(x\)相关的、一般不可能求得的数。

对任意固定的\(x\in [a,b]\setminus\{x_0,\cdots,x_n\}\)（这时候\(x\)就视为常数），考虑辅助函数

\[\phi(t)=R(t)-\frac{R(x)}{\omega_n(x)}\omega(t), \]

由于\(f(t)\)、\(p_n(t)\)和\(\omega_n(t)\)在\(t\in \{x_0,\cdots,x_n\}\)时都为\(0\)，且特别当\(t=x\)时，\(\phi(x)=0\)，所以\(\phi(t)\)在\([a,b]\)上至少具有\(n+2\)个零点：\(x,x_0,\cdots,x_n\)。

利用Rolle定理，\(\phi'(t)\)在\((a,b)\)上至少具有\(n+1\)个零点，\(\phi''(t)\)在\((a,b)\)上至少具有\(n+2\)个零点，以此类推得到\(\phi^{(n+1)}(t)\)在\((a,b)\)上至少具有一个零点，记作\(\xi\)，也就是

\[\phi^{(n+1)}(\xi)=R^{(n+1)}(\xi)-\frac{R(x)}{\omega_n(x)}\omega_n^{(n+1)}(\xi)=0, \]

由于\(R^{(n+1)}(t)=f^{(n+1)}(t)-p_n^{(n+1)}(t)\)，且\(p_n(t)\)是不超过\(n\)次的多项式，故\(n+1\)阶导数为\(0\)；\(\omega_n(t)\)是最高次系数为\(1\)的多项式，所以\(\omega_n^{(n+1)}(t)=(n+1)!\)，于是

\[f^{(n+1)}(\xi)=\frac{R(x)}{\omega_n(x)}(n+1)!, \]

从中解得

\[R(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_n(x). \]

如果\(f^{(n+1)}(x)\)在\([a,b]\)上有界，即\(|f^{(n+1)}(\xi)|\le M\)恒成立，那么截断误差也有上界为

\[|R(x)|\le \frac{M}{(n+1)!}\omega_n(x). \]

作为误差估计的应用，我们举出以下例子来看看如何使用误差估计来求插值多项式。依然以\(f(x)=x^4\)为例，求其在\(-1,0,1,2\)处的插值多项式\(p_3(x)\)，我们有

\[f(x)-p_3(x)=\frac{f^{(4)}(\xi)}{4!}(x-(-1))(x-0)(x-1)(x-2). \]

但特别的是，此时\(f^{(4)}(x)=4!\)是一个常数，因此不管\(\xi\)取\([-1,2]\)中的哪一个值，都有

\[f(x)-p_3(x)=(x+1)x(x-1)(x-2)=x^4-2x^3-x^2+2x, \]

这就说明\(p_3(x)=2x^3+x^2-2x\)，与Lagrange法和Newton法计算出来的一致，并且大大简化了计算。不过，这是一个特殊的例子，因为恰好有\(f^{(n+1)}(x)=C\)，倘若\(f(x)\)是更高次的多项式或者甚至不是多项式，则由于\(\xi\)的未知性，无法用这种方式得出插值多项式。

二、分段插值

分段插值指的是在插值节点分成的每一个小区间内使用某种方法插值，将所有插值函数拼接在一起形成插值函数。

Runge现象与Faber定理

在统计中，许多情况下都是样本容量\(n\)越大，对未知的估计就越准确，因此我们自然会想，在一个区间上，是否插值节点越多，插值多项式就会越来越逼近被插函数。但事实并非如此，Runge现象给出了一个反例：对于函数

\[f(x)=\frac{1}{1+25x^2},\quad -1\le x\le 1. \]

在\([-1,1]\)上取等距节点，节点个数为\(n+1\)，当\(n\to \infty\)时，\(\varphi_n(x)\)依然不会收敛于\(f(x)\)，而只是在某一区间上收敛。这也说明了即使节点不断加密，插值多项式次数不断升高，高次插值多项式也不一定收敛到相应的被插函数。

更一般地有如下的Faber定理：对\([a,b]\)上任给的三角矩阵

\[\begin{matrix} x_0^{(0)} \\ x_0^{(1)} & x_1^{(1)} \\ x_0^{(2)} & x_{1}^{(2)} & x_2^{(2)} \\ \cdots \\ x_0^{(n)} & x_1^{(n)} & x_2^{(n)} & \cdots & x_n^{(n)} \end{matrix} \]

总存在\(f(x)\in C[a,b]\)，使得由三角形矩阵中的任一行元素为插值节点所生成的插值多项式\(\varphi_n(x)\)，当\(n\to \infty\)时不能一致收敛到\(f(x)\)。

因此，需要对原有的多项式插值作出修改，以便插值函数能更好地应用在实际中，分段插值就是一种解决Runge现象的方法。

分段线性插值

分段线性插值是最简单的分段插值，设有如下的\(n+1\)个插值节点：\(x_0<x_1<\cdots<x_n\)，它将区间\([x_0,x_n]=[a,b]\)划分为\(n\)个小区间。在每一个小区间上使用多项式插值，由于每一个小区间上只有两个节点，因此多项式只能是一次函数。直观上看，在\([a,b]\)上的分段线性插值结果，就是将所有的节点用一条折线连接起来。

分段线性插值的构造依然使用基函数，此时对\([x_i,x_{i+1}]\)上的基函数，就应该是

\[l_i(x)=\left[y_{i}+\frac{y_{i+1}-y_i}{x_{i+1}-x_i}(x-x_i)\right]I_{[x_i,x_{i+1}]}, \]

这里\(I_{[x_{i},x_{i+1}]}\)是示性函数，当\(x\in [x_i,x_{i+1}]\)时为\(1\)，否则为\(0\)。总的分段线性插值函数就是将各个\(l_i(x)\)相加。

由于分段线性插值是多项式插值的组合，那么在每一个小区间上，就自然有相应的误差估计。在\([x_{i},x_{i+1}]\)上，我们已经得到了误差估计为

\[R(x)=\frac{f''(\xi)}{2}(x-x_i)(x-x_{i+1}), \]

注意，此时\((x-x_i)(x-x_{i+1})\)的最值在区间的中点处取到。如果在此区间上\(f''(x)\)有界，即\(|f''(x)|\le M\)，那么在此小区间上，就有

\[|R(x)|\le \frac{\max |f''(x)|}{2}\left(\frac{x_i+x_{i+1}}{2}-x_i\right)\left(x_{i+1}-\frac{x_i+x_{i+1}}{2} \right)\le \frac{M}{8}(x_{i+1}-x_i)^2. \]

特别当\(f''(x)\)在整个区间\([a,b]\)上都有界\(M\)，记\(h=\max_i(x_{i+1}-x_i)\)，那么分段线性插值的误差函数就是

\[|R(x)|\le \frac{M}{8}h^2. \]

在许多时候，结点是等距的，故\(h=x_{i+1}-x_i\)，可以用此式给出误差估计。

尽管分段线性插值克服了Runge现象，但是由于函数是折线，自然也就带来了不光滑的缺陷，即\(\varphi_n(x)\)在\([a,b]\)上不是处处可微的，这样就缺少了许多连续函数具有的优良性质。

两点三次Hermite插值

为了克服分段线性插值的不可微缺陷，我们需要一种带微分的插值方法（Hermite插值），具体表现为在某一个小区间\([x_0,x_1]\)的两个端点处，不但指定函数值，还指定导函数的值，并用一个多项式来拟合。这样拟合的好处就是，在分段插值时，如果限定了在端点处的导数，就可以让小区间的左导数等于右导数，从而让分段插值函数在\([a,b]\)上整体可微，即\(\varphi_n(x)\in C^1[a,b]\)。

由于现在有\(f(x_0)=y_0\)，\(f(x_1)=y_1\)，\(f'(x_0)=y_0'\)，\(f'(x_1)=y_1'\)四个约束条件，所以如果使用待定系数法，要获得唯一的插值多项式，应当是三次的。这里，我们不使用待定系数法，而是用Lagrange法和Newton法的思想来解决带微商的插值问题。

Lagrange法插值的基本思想是基函数，对于这个问题，对应于四个约束条件，也需要四个基函数，满足

\[\begin{array}{|c|} \hline & h(x_0) & h(x_1) & h'(x_0) & h'(x_1) \\ \hline h_0(x) & 1 & 0 & 0 & 0 \\ \hline h_1(x) & 0 & 1 & 0 & 0 \\ \hline H_0(x) & 0 & 0 & 1 & 0 \\ \hline H_1(x) & 0 & 0 & 0 & 1 \\ \hline \end{array} \]

对于\(h_0(x)\)，为使\(h(x_1)=h'(x_1)=0\)，因式中应当含有\((x-x_1)^2\)，同时为满足\(h(x_0)=1\)，可以设其形式为

\[h_0(x)=(1+b(x-x_0))\left(\frac{x-x_1}{x_0-x_1} \right)^2, \]

再通过\(h_0'(x_0)=b+2=0\)得到\(b=-2/(x_0-x_1)\)，由此可以确定\(h_0(x)\)的形式，\(h_1(x)\)也类似得到，即

\[h_0(x)=\left(1+2\frac{x-x_0}{x_1-x_0} \right)\left(\frac{x-x_1}{x_0-x_1} \right)^2, \\ h_1(x)=\left(1+2\frac{x-x_1}{x_0-x_1} \right)\left(\frac{x-x_0}{x_1-x_0} \right)^2. \]

对于\(H_0(x)\)，为使\(h(x_1)=h'(x_1)=0\)，因式中应含有\((x-x_1)^2\)，同时为满足\(h(x_0)=0\)，因式中应含有\((x-x_0)\)，可以设其形式为

\[H_0(x)=a(x-x_0)\left(\frac{x-x_1}{x_0-x_1} \right)^2, \]

结合\(H_0'(x_0)=1\)可推出\(a=1\)，由此可以确定\(H_0(x_)\)的形式，\(H_1(x)\)也类似得到，即

\[H_0(x)=(x-x_0)\left(\frac{x-x_1}{x_0-x_1} \right)^2,\\ H_1(x)=(x-x_1)\left(\frac{x-x_0}{x_1-x_0} \right)^2. \]

最后得到Hermite插值多项式为

\[H(x)=y_0h_0(x)+y_1h_1(x)+y_0'H_0(x)+y_1'H_1(x). \]

利用基函数求过\(0,1\)两点的三次Hermite插值多项式，要求

\[f(0)=1,\quad f'(0)=\frac{1}{2},\\ f(1)=2,\quad f'(1)=\frac{1}{2}. \]

套用以上形式，容易算出基函数为

\[h_0(x)=(1+2x)(x-1)^2,\\ h_1(x)=(3-2x)x^2,\\ H_0(x)=(x-1)^2x,\\ H_1(x)=x^2(x-1), \]

因此三次Hermite插值多项式为

\[\varphi(x)=1h_0(x)+2h_1(x)+\frac{1}{2}H_0(x)+\frac{1}{2}H_1(x)=-x^3+\frac{3}{2}x^2+\frac{1}{2}x+1. \]

基函数的思想虽然直白，但运算量上还是不小，而利用Newton法的核心——差商，则显得巧妙很多。在此，我们要先注意到一阶差商的形式为

\[f[x_0,x_1]=\frac{f(x_1)-f(x_0)}{x_1-x_0}, \]

而导数恰好就是差商的极限，也就是说，当\(x_1\to x_0\)时，就有

\[f'(x_0)=\lim_{x_1\to x_0}\frac{f(x_1)-f(x_0)}{x_1-x_0}=\lim_{x_1\to x_0}f[x_1,x_0]. \]

运用这个关系，我们可以在端点处，构造两个极限节点，即\(x_0'\to x_0\)，\(x_1'\to x_1\)，这样，二阶差商的值就通过导数给出，从而我们有

\[\begin{aligned} \varphi(x)&= f(x_0)\\ &+ f[x_0,x_0'](x-x_0)\\ &+f[x_0,x_0',x_1](x-x_0)(x-x_0')\\ &+ f[x_0,x_0',x_0,x_1](x-x_0)(x-x_0')(x-x_1)\\ &=f(x_0)+f'(x_0)(x-x_0)+\frac{f[x_0,x_1]-f'(x_0)}{x_1-x_0}(x-x_0)^2 \\ &+\frac{f'(x_0)+f'(x_1)-2f[x_0,x_1]}{(x_1-x_0)^2}(x-x_0)^2(x-x_1). \end{aligned} \]

尽管看起来这个式子很繁琐，但实际应用起来是很简便的。

使用差商求过\(0,1\)两点的三次Hermite插值多项式，依然要求

\[f(0)=1,\quad f'(0)=\frac{1}{2},\\ f(1)=2,\quad f'(1)=\frac{1}{2}. \]

先构造差商表为

\[\begin{array}{|c|} \hline k & f(x_k) & \partial^1 & \partial^2 & \partial^3\\ \hline 0 & \underline{1} \\ \hline & & \underline{\frac{1}{2}} \\ \hline 0^+ & 1 && \underline{\frac{1}{2}} \\ \hline & & 1 && \underline{-1} \\ \hline 1 & 2 && -\frac{1}{2} &\\ \hline & & \frac{1}{2} && \\ \hline 1^+ & 2 && \\ \hline \end{array} \]

于是Hermite插值多项式为

\[\begin{aligned} \varphi(x)&= 1 \\ &+ \frac{1}{2}(x-0) \\ &+ \frac{1}{2}(x-0)^2 \\ &+ (-1)(x-0)^2(x-1)\\ &= -x^3+\frac{3}{2}x^2+\frac{1}{2}x+1. \end{aligned} \]

在计算过程上，Newton法求Hermite插值多项式要更快捷。

Hermite插值

两点三次Hermite插值是Hermite插值的特殊情况，它只在两个点处给出了被插函数的值和导数，实际上Hermite插值可以推广到\(n+1\)个节点\(x_0,\cdots,x_n\)上，此时需给定这些节点的函数值和导数值：

\[\begin{array}{|c|} \hline x & x_0 & x_1 & x_2 & \cdots & x_n \\ \hline y & y_0 & y_1 & y_2 & \cdots & y_n \\ \hline y' & y_0' & y_1' & y_2' & \cdots & y_n' \\ \hline \end{array} \]

由于有\(2n+2\)个约束条件，故构造出的Hermite插值多项式应当是\(2n+1\)次的（为了含有\(2n+2\)个未知数）。可以证明，插值基函数是

\[h_i(x)=\left(1-\frac{\omega_n''(x_i)}{\omega_n'(x_i)}(x-x_i) \right)^2\prod_{j\ne i}\frac{x-x_j}{x_i-x_j},\\ H_i(x)=\left(\prod_{j\ne i}\frac{x-x_j}{x_i-x_j} \right)^2(x-x_i). \]

于是Hermite插值多项式是

\[H(x)=\sum_{i=0}^{n}(y_ih_i(x)+y_i'H_i(x)). \]

接下来给出其误差估计定理：如果\(f^{(2n+2)}(x)\)在\([a,b]\)存在，则\(\forall x\in [a,b]\)，总存在\(\xi\in (a,b)\)，使得

\[R(x)=f(x)-H(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}\omega_n^2(x). \]

证明过程与多项式插值是类似的。

引入辅助函数为

\[\psi(t)=R(t)-\frac{R(x)}{\omega_n^2(x)}\omega^2_n(t), \]

由于此时\(x_0,\cdots,x_n\)都是\(\psi(t)\)的二重零点，故它们依然是\(\psi'(t)\)的零点。此外，\(\psi(t)\)在\([a,b]\)上拥有\(x,x_0,\cdots,x_n\)一共\(n+2\)个零点，所以由Rolle定理，\(\psi'(t)\)在这些区间内还拥有\(n+1\)个异于\(x_0,\cdots,x_n\)的零点，故\(\psi'(t)\)在\([a,b]\)内至少应拥有\(2n+2\)个零点。反复使用Rolle定理，可知\(\psi^{(2n+2)}(t)\)在\([a,b]\)内拥有一个零点，记作\(\xi\)，即

\[R^{(2n+2)}(\xi)=\frac{R(x)}{\omega_n^2(x)}[\omega_n^2(\xi)]^{(2n+2)}. \]

由于\(H(t)\)是\(2n+1\)次多项式，\(\omega_n^2(t)\)是最高次项为\(1\)的\(2n+2\)次多项式，故

\[f^{(2n+2)}(\xi)=(2n+2)!\frac{R(x)}{\omega_n^2(x)}, \]

即

\[R(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}\omega_n^2(x). \]

分段三次Hermite插值

如果我们知道了\(x_0,\cdots,x_n\)处的函数值以及导数值，除了能直接构造一个\(2n+1\)次多项式进行Hermite插值，也可以类似分段线性插值一样，构造分段Hermite插值。此时，在每一段上的插值函数是一个三次多项式，类似于已知\(x_0,x_1\)处的函数值\(y_0,y_1\)和\(y_0',y_1'\)一样。

由于此过程类似分段线性插值，只是使用了不同的插值多项式，因此不过多展开。要注意的是，分段线性插值是一条折线，因此不是\([a,b]\)上可导的。但是分段Hermite插值保证了每一个结点处的微商值，因此在每一个小区间上左导数等于右导数，所以分段三次Hermite插值是\(C^1[a,b]\)的。

现在来看在每一个小区间上的误差估计，前面我们已经得出了一般情况下Hermite插值的误差估计，特别在两个点\(x_0,x_1\)处，有

\[R(x)=\frac{f^{(4)}(\xi)}{4!}(x-x_0)^2(x-x_1)^2, \]

如果\(f^{(4)}(x)\)在\([x_0,x_1]\)上有界\(M\)，则

\[|R(x)|\le\frac{M}{4!}(x-x_0)^2(x-x_1)^2\le \frac{M}{4!}\left(\frac{x_0+x_1}{2}-x_0 \right)^2\left(\frac{x_0+x_1}{2}-x_1 \right)^2=\frac{M(x_1-x_0)^4}{384}. \]

因此，如果\(f^{(4)}(x)\)在\([x_0,x_1]\)上有界\(M\)，且每一个小区间的长度最大值为\(h\)，那么分段三次Hermite插值的误差估计上限为

\[|R(x)|\le \frac{Mh^4}{384}. \]

三、样条插值

最后简要介绍一下样条插值。分段Hermite插值虽然解决了分段线性插值不可导的问题，但也存在一些问题，一是实际生活中，插值点处的微商\(y_0',y_1',\cdots,y_n'\)难以获得；二是此函数是\(C^1\)的而不是\(C^2\)的。样条插值便是在每一个区间\([x_i,x_{i+1}],i=0,1,\cdots,n-1\)上用三次函数进行插值，使得插值函数在连接点具有连续曲率。

样条插值特点

具体说来，三次样条插值函数\(s(x)\)满足下列条件：

1. \(s(x_i)=y_i\)，\(i=0,1,\cdots,n\)；
2. 在每一个小区间\([x_i,x_{i+1}]\)上，\(s(x)\)是一个不高于三次的多项式；
3. \(s(x)\in C^2[a,b]\)。

不过，这三个条件不能唯一确定\(s(x)\)，因为约束条件个数不足。如果要使\(s(x)\)是唯一确定的，还需要两个边界条件（下一部分将解释为什么是两个）。常用的边界条件有：

1. 自然边界条件：\(s''(x_0)=s''(x_n)=0\)。
2. 固定边界条件：\(s'(x_0)=y_0'\)，\(s'(x_n)=y_n'\)已知。
3. 周期条件：\(f(x)\)是周期函数，基本周期为\(x_n-x_0\)，从而\(s(x)\)是周期函数。

样条插值函数的一般求解

设每一个小区间\([x_i,x_{i+1}]\)上的插值函数是\(s_i(x)\)，\(i=0,1,\cdots,n-1\)，其区间长度为\(x_{i+1}-x_i=h_i\)，此时\(y_i,h_i\)都是已知量。设

\[s_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i)^2+d_i(x-x_i)^2, \]

则由\(s_i(x_i)=y_i\)可得\(a_i=y_i\)。另一边，需要满足\(s_i(x_{i+1})=y_{i+1}\)，也就是

\[y_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_{i+1}.\quad (1) \]

接下来考虑微商之间的关系，有\(s_i'(x)=b_i+2c_i(x-x_i)+3d_i(x-x_i)^2\)，由\(s_i'(x_{i+1})=s'_{i+1}(x_{i+1})\)，有

\[b_i+2c_ih_i+3d_ih_i^2=b_{i+1}\quad (2), \]

由于样条插值是\(C^2\)的，所以考虑二阶微商，即\(s_i''(x)=2c_i+6d_i(x-x_i)\)，由\(s_{i}''(x_{i+1})=s''_{i+1}(x_{i+1})\)，有

\[2c_i+6d_ih_i=2c_{i+1}\quad (3). \]

现设\(m_i=2c_i\)，则由\((3)\)，有

\[d_i=\frac{2c_{i+1}-2c_{i}}{6h_i}=\frac{m_{i+1}-m_i}{6h_i}, \]

再由\((1)\)，有

\[b_i=\frac{y_{i+1}-y_{i}-\frac{h_i^2(m_{i+1}-m_{i})}{6}-\frac{m_ih_i^2}{2}}{h_i}=\frac{y_{i+1}-y_i}{h_i}-\frac{h_i(m_{i+1}+2m_{i})}{6}, \]

这样，所有未知数都是关于\(m_i\)的函数，最后回到关系式\((2)\)，得到

\[\frac{y_{i+1}-y_{i}}{h_i}-\frac{h_i(m_{i+1}+2m_i)}{6}+m_ih_i+\frac{h_i(m_{i+1}-m_i)}{2}=\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{h_{i+1}(m_{i+2}+2m_{i+1})}{6}, \]

稍加整理，得到

\[h_im_i+2(h_i+h_{i+1})m_{i+1}+h_{i+1}m_{i+2}=6\left[\frac{y_{i+2}-y_{i+1}}{h_{i+1}}-\frac{y_{i+1}-y_i}{h_i} \right],\quad i=0,1,\cdots,n-1. \]

这是一个关于\(m_0,m_1,\cdots,m_{n}\)的线性方程组，含有\(n+1\)个变量，却只有\(n-1\)个等式，因此还缺少两个约束条件。

自然边界求解

现在对\(s''(x_0)=s''(x_n)=0\)的情况给出求解。此时\(s_0''(x_0)=2c_0=0\)，故\(m_0=0\)；\(s_n''(x_n)=0\)有\(m_n=0\)（尽管我们只使用了\(s_0,\cdots,s_{n-1}\)个等式，但\(s_n\)可以类似定义），这样，我们可以得到

现在对\(s''(x_0)=s''(x_n)=0\)的情况给出求解。此时\(s_0''(x_0)=2c_0=0\)，故\(m_0=0\)；\(s_n''(x_n)=0\)有\(m_n=0\)（尽管我们只使用了\(s_0,\cdots,s_{n-1}\)个等式，但\(s_n\)可以类似定义），这样，我们可以得到一个线性方程组：

\[\begin{bmatrix} 1 & 0 & 0 & 0 & \cdots & 0 &0 \\ h_0 & 2(h_0+h_1) & h_1 & 0 & \cdots & 0 & 0 \\ 0 & h_1 & 2(h_1+h_2) & h_2 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & \cdots & h_{n-2} & 2(h_{n-2}+h_{n-1}) & h_{n-1} \\ 0 & 0 & 0 & 0 & \cdots &0 & 1 \end{bmatrix}\begin{bmatrix} m_0 \\ m_1 \\ m_2 \\ \vdots \\ m_{n-1} \\ m_n \end{bmatrix}\\ =6\begin{bmatrix} 0 \\ \frac{y_2-y_1}{h_1}-\frac{y_1-y_0}{h_0} \\ \frac{y_3-y_2}{h_2}-\frac{y_2-y_1}{h_1} \\ \vdots \\ \frac{y_{n}-y_{n-1}}{h_{n-1}}-\frac{y_{n-1}-y_{n-2}}{h_{n-2}} \\ 0 \end{bmatrix} \]

线性方程组的求解是容易的，这样我们解出了\(m_0,\cdots,m_n\)之后，就可以确定各个\(b_i,c_i,d_i\)的值，从而得到每一区间上的三次Hermite插值多项式。