7 函数插值
7.1 引言
已知函数\(y=f(x)\)在\(n+1\)个点\(x_0,x_1,\cdots,x_n\)上的函数值\(y_i=f(x_i),i=0,1,\cdots,n\)，但是不知道\(f(x)\)的具体表达式。若要估计\(f(x)\)在某点\(x'\)的值的时候，就需要找到一个函数\(P(x)\)来近似\(f(x)\)，然后用\(P(x')\)估计\(f(x')\)的值。
并且对于函数\(P(x)\)有如下要求：
\[P(x_i)=y_i,\;i=0,1,\cdots,n\tag{7-1}
\]
我们将\(f(x)\)称为被插值函数，\(P(x)\)称为插值函数，\(x_0,x_1,\cdots,x_n\)称为插值节点，求解插值函数的方法称为插值法。如果需要估计的点\(x'\)位于包含插值节点的最小闭区间内，将其称之为内插，否则称为外插。最后还应有一个误差函数
\[R(x)=f(x)-P(x)
\]
称为插值余项。
综上可知，函数插值实质上就是寻找一个函数\(P(x)\)来近似\(f(x)\)的方法。常用的插值函数是多项式函数，因此，在这一部分将主要介绍几种基本的多项式插值方法。
最简单的插值函数是多项式函数，这时插值问题就变成了，寻找一个\(n\)次多项式
\[P_n(x)=a_0+a_1x+a_2x^2+\cdots+a_nx^n
\]
并使其满足式子\((7-1)\)。最简单直接的方法就是求解一个线性方程组：
\[\left\{\begin{aligned}
a_0+a_1x_0+a_2x^2_0+\cdots+a_nx^n_0&=y_0\\
a_0+a_1x_1+a_2x^2_1+\cdots+a_nx^n_1&=y_1\\
\cdots\qquad&\\
a_0+a_1x_n+a_2x^2_n+\cdots+a_nx^n_n&=y_n\\
\end{aligned}\right.
\]
其系数行列式为
\[\left|\begin{matrix}
1&x_0&x_0^2&\cdots&x_0^n\\
1&x_1&x_1^2&\cdots&x_1^n\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&x_n&x_n^2&\cdots&x_n^n\\
\end{matrix}\right|=
\prod_{i=1}^n\prod_{j=0}^{i-1}(x_i-x_j)\ne0
\]
根据范德蒙行列式的性质很容易就可以知道上述行列式的值不为0，因此方程组有唯一解，即\(n+1\)个插值节点能够唯一确定一个\(n\)次多项式。但阶数高的时候，直接求解方程组的计算量太大，不实用。因此，我们通常寻求其他的方法来求得插值多项式。
此外根据这个结论，我们可以有以下推论：
- 如果被插值函数本身就是一个不超过\(n\)次的多项式，那么插值函数就是该多项式本身。
- 无论采用什么方法得到的插值多项式，本质上都是相同的，只是表达形式上的不同。
7.2 拉格朗日（Lagrange）插值
插值方法：
设满足要求的\(n\)次多项式可以表示为
\[L_n(x)=y_0l_0(x)+y_1l_1(x)+\cdots+y_nl_n(x)=\sum_{i=0}^ny_il_i(x)
\]
其中\(l_i(x)\)是\(n\)次多项式，并且要求\(l_i(x_i)=1,l_i(x_j)=0(j\ne i)\)。我们将\(l(x)\)称为插值基函数。
以构造\(l_1(x)\)为例：
-
首先，要满足\(l_1(x_j)=0,j\ne1\)，这说明\(l_1(x)=0\)具\(n\)个根\(x_j(j\ne1)\)。因此可以被因式分解为
\[l_1(x)=a(x-x_0)(x-x_2)\cdots(x-x_n)
\]
为了满足\(l_1(x_1)=1\)，即\(a(x_1-x_0)(x_1-x_2)\cdots(x_1-x_n)=1\)，因此
\[a=\frac{1}{(x_1-x_0)(x_1-x_2)\cdots(x_1-x_n)}
\]
于是
\[l_1(x)=\frac{(x-x_0)(x-x_2)\cdots(x-x_n)}{(x_1-x_0)(x_1-x_2)\cdots(x_1-x_n)}
\]
-
同理可以构造出其他的\(l_i(x)\)。
因此
\[L_n(x)=\sum_{i=0}^ny_i\frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}
{(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}
=\sum_{i=0}^ny_i\prod_{j=0\\j\ne i}^n\frac{x-x_j}{x_i-x_j}
\]
令\(\omega_{n+1}(x)=(x-x_0)(x-x_1)\cdots(x-x_n)\)，则可以将插值基函数表示成
\[l_i(x)=\frac{\omega_{n+1}(x)}{(x-x_i)\omega'_{n+1}(x_i)}
\]
因此
\[L_n(x)=\sum_{i=0}^ny_il_i(x)=\sum_{i=0}^ny_i\frac{\omega_{n+1}(x)}{(x-x_i)\omega'_{n+1}(x_i)}
\]
特殊地，当\(n=1\)时称为为线性插值，\(n=2\)时称为抛物线型插值。
插值余项：
假设\(a=x_0<x_1<\cdots<x_n=b\)，设函数\(f(x)\)在区间\([a,b]\)上具有\(n\)阶连续导数，并且在区间\((a,b)\)上存在\(n+1\)阶导数，那么拉格朗日插值法的插值余项为
\[R_n(x)=f(x)-L_n(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x),\;\xi\in(a,b)
\]
证明：显然当\(x=x_i(i=0,1,\cdots,n)\)时，上述式子成立。现将\(x\in[a,b]\)看成是一个不同于插值节点的固定点，构造一个辅助函数
\[\varphi(t)=f(t)-L_n(t)-k(x)\omega_{n+1}(t)
\]
其中\(k(x)\)是一个关于\(x\)的参数，目的是为了使得\(\varphi(x)=0\)。于是，\(\varphi(t)\)一共有\(x,x_0,x_1,\cdots,x_n\)这\(n+2\)个零点。接着，根据罗尔（Rolle）定理可知，\(\varphi'(t)\)一定在\((a,b)\)内至少有\(n+1\)个不同的零点。按照此思路，最终可以得到\(\varphi^{(n+1)}(t)\)在区间\((a,b)\)内至少有一个零点。因此
\[\varphi^{(n+1)}(\xi)=f^{(n+1)}(\xi)-0-(n+1)!k(x)=0\quad
\Rightarrow\quad k(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!},\;\xi\in(a,b)
\]
这也意味着
\[\varphi(x)=0\quad\Rightarrow\quad R_n(x)=f(x)-L_n(x)=k(x)\omega_{n+1}(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)
\]
有了上述插值余项理论上就能够估计插值的误差了，但是我们必须知道\(|f^{(n+1)}(x)|\)的一个上界，否则插值余项只有理论意义，无法用于实际估计。假设我们知道了\(|f^{(n+1)}(x)|\)的一个上界\(M\)，那么
\[|R_n(x)|\le\frac{M}{(n+1)!}|\omega_{n+1}(x)|\le\frac{M}{(n+1)!}(b-a)^{n+1}
\]
拉格朗日插值存在的问题： 在实际应用中，为了提高插值精度，有时需要增加插值节点的个数，那么原本已经求出的\(L_n(x)\)将作废，不得不重新计算所有的插值基函数。
7.3 牛顿（Newton）插值
为了解决拉格朗日插值的缺点，这一小节将介绍牛顿插值，类似于泰勒展开式求近似值
\[f(x)=f(x_0)+\frac{f'(x_0)}{1!}(x-x_0)+\cdots+\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n+R_n(x)
\]
为了增加精度只需要在之前的数据基础上，增加新的项即可。
7.3.1 差商及其性质
下面将给出差商的定义
一阶差商：
\[f[x_0,x_1]=\frac{f(x_0)-f(x_1)}{x_0-x_1}
\]
二阶差商（差商的差商）：
\[f[x_0,x_1,x_2]=\frac{f[x_0,x_1]-f[x_1,x_2]}{x_0-x_2}
\]
n阶差商：
\[f[x_0,x_1,\cdots,x_n]=\frac{f[x_0,x_1,\cdots,x_{n-1}]-f[x_1,x_2,\cdots,x_n]}{x_0-x_n}
\]
特殊地，我们称\(f[x_i]=f(x_i)\)是零阶差商。
介绍完差商的基本定义后，给出差商的几个基本性质。
性质1 \(k\)阶差商可以表示为相应节点上的函数值的线性组合，即
\[f[x_0,x_1,\cdots,x_k]=\sum_{j=0}^k\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_k)}
\]
证明：利用归纳法来证明，容易验证\(k=1\)时成立
\[f[x_0,x_1]=\frac{f(x_0)}{x_0-x_1}+\frac{f(x_1)}{x_1-x_0}
\]
假设\(k=l\)的时候也成立，那么
\[f[x_0,x_1,\cdots,x_l]=\sum_{j=0}^l\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_l)}\\
f[x_1,x_2,\cdots,x_{l+1}]=\sum_{j=1}^{l+1}\frac{f(x_j)}{(x_j-x_1)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}
\]
因为
\[\begin{aligned}
\frac{1}{x_0-x_{l+1}}f[x_0,x_1,\cdots,x_l]&=\frac{f(x_0)}{(x_0-x_1)\cdots(x_0-x_{l+1})}\\&+\sum_{j=1}^l\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}\cdot\frac{x_j-x_{l+1}}{x_0-x_{l+1}}\\
\frac{1}{x_{l+1}-x_0}f[x_1,x_2,\cdots,x_{l+1}]&=\frac{f(x_{l+1})}{(x_{l+1}-x_0)\cdots(x_{l+1}-x_l)}\\&+\sum_{j=1}^l\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}\cdot\frac{x_j-x_0}{x_{l+1}-x_0}
\end{aligned}
\]
于是
\[\begin{aligned}
&\sum_{j=1}^l\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}\cdot\frac{x_j-x_{l+1}}{x_0-x_{l+1}}\\
+&\sum_{j=1}^l\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}\cdot\frac{x_j-x_0}{x_{l+1}-x_0}\\
=&\sum_{j=1}^l\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}\cdot\left(\frac{x_j-x_0}{x_{l+1}-x_0}+\frac{x_j-x_{l+1}}{x_0-x_{l+1}}\right)\\
=&\sum_{j=1}^l\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}
\end{aligned}
\]
综上可知
\[\begin{aligned}
f[x_0,x_1,\cdots,x_{l+1}]&=\frac{1}{x_0-x_{l+1}}f[x_0,x_1,\cdots,x_l]+\frac{1}{x_{l+1}-x_0}f[x_1,x_2,\cdots,x_{l+1}]\\
&=\sum_{j=0}^{l+1}\frac{f(x_j)}{(x_j-x_0)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{l+1})}
\end{aligned}
\]
因此\(k=l+1\)也成立，所以对于任意的\(k\in\mathbb{N}\)都是成立的。
性质2 差商与节点的排列顺序无关，即
\[f[x_{i_0},x_{i_1},\cdots,x_{i_n}]=f[x_0,x_1,\cdots,x_n]
\]
其中\(i_0,i_1,\cdots,i_n\)是\(0,1,2,\cdots,n\)的任意一种排列。
性质2是性质1的推论，将差商表示成插值节点处函数值的线性组合很容易就可以发现性质2成立。
性质3 若\(f[x,x_0,\cdots,x_k]\)是关于\(x\)的\(m\)次多项式，则\(f[x,x_0,\cdots,x_{k+1}]\)是关于\(x\)的\(m-1\)次多项式。
证明：
\[f[x,x_0,\cdots,x_{k+1}]=\frac{f[x,x_0,\cdots,x_k]-f[x_0,\cdots,x_k,x_{k+1}]}{x-x_{k+1}}
\]
等式右侧的分子是关于\(x\)的\(m\)次多项式，并且当\(x=x_{k+1}\)时，分子为0。这说明分子含有因子\((x-x_{k+1})\)，该因子可以和分母相约分，因此\(f[x,x_0,\cdots,x_{k+1}]\)是关于\(x\)的\(m-1\)次多项式。
性质4 若\(f(x)\)是关于\(x\)的\(n\)次多项式，则\(f[x,x_0,\cdots,x_n]=0\)。
证明：
根据性质3可知，\(f[x,x_0]\)是\(n-1\)次多项式，\(f[x,x_0,x_1]\)是\(n-2\)次多项式。同理可得\(f[x,x_0,\cdots,x_{n-1}]\)是\(0\)次多项式（是个常数）。所以性质4成立。
7.3.2 插值方法
差商的计算： 通过列差商表的方式计算差商，例子如下
| \(x_i\) |
\(f(x_i)\) |
一阶差商 |
二阶差商 |
三阶差商 |
| 1 |
0 |
|
|
|
| 3 |
2 |
1 |
|
|
| 4 |
15 |
13 |
4 |
|
| 7 |
12 |
-1 |
-3.5 |
-1.25 |
通过列表的形式计算差商表达起来比较清晰，并且如果需要增加新的节点，那么只需要在原有的差商表上增加新的行，就可以计算更高阶的差商。
牛顿插值多项式：
根据差商的定义，设\(x\in[a,b]\)，那么
\[f[x,x_0]=\frac{f(x)-f(x_0)}{x-x_0}\;\Rightarrow\;f(x)=f(x_0)+f[x,x_0](x-x_0)
\]
同理
\[f[x,x_0,x_1]=\frac{f[x,x_0]-f[x_0,x_1]}{x-x_1}\;\Rightarrow\;
f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x,x_0,x_1](x-x_0)(x-x_1)
\]
最终我们可以得到\(f(x)=N_n(x)+R_n(x)\)，其中
\[\begin{aligned}
&N_n(x)=f(x_0)+f[x_0,x_1](x-x_0)+\cdots+f[x_0,x_1,\cdots,x_n](x-x_0)\cdots(x-x_{n-1})\\
&R_n(x)=f[x,x_0,x_1,\cdots,x_n](x-x_0)\cdots(x-x_n)=f[x,x_0,x_1,\cdots,x_n]\omega_{n+1}(x)
\end{aligned}
\]
在\(N_n(x)\)中的差商部分都不含有\(x\)，因此\(N_n(x)\)就是关于\(x\)的\(n\)次多项式，就是我们所要求得牛顿插值多项式。所以\(R_n(x)\)就是插值余项。
特殊地，当\(n=1\)时
\[N_1(x)=f(x_0)+\frac{f(x_0)-f(x_1)}{x_0-x_1}(x-x_0)
\]
显然就是直线得点斜式方程。
插值余项（被插值函数可导）：
根据插值多项式的唯一性，拉格朗日插值多项式\(L_n(x)\)与牛顿插值多项式\(N_n(x)\)是完全相同的，因此它们的余项也是相同的，因此
\[f[x,x_0,x_1,\cdots,x_n]\omega_{n+1}(x)=\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}(x)
\]
这就意味着\(n\)阶差商和被插值函数的导数具有如下关系：
\[f[x_0,x_1,\cdots,x_n]=\frac{f^{(n)}(\xi)}{n!},\;\xi\in[a,b]
\]
牛顿插值的优点：
- 对于牛顿插值而言，当增加插值节点的个数之后，只需要在原有的差商表上添加新的行，计算更高阶的差商就可以得到更高阶的插值多项式。而拉格朗日插值则必须重新计算所有的基函数。
- 牛顿插值的余项公式对于被插值函数\(f(x)\)是离散的或导数不存在时均适用，而拉格朗日插值余项则要求被插值函数可导。
7.4 埃尔米特（Hermite）插值
不少实际问题不但要求在节点上的函数值相等，而且还有要求其导数值也要相等（即插值节点上具有一阶光滑度），甚至要求更高阶的导数值也要相等，满足这种要求的插值多项式就是埃尔米特插值。显然要求的条件更多了，所需要的多项式次数也就更高了。在这一部分将介绍函数值和导数值个数相等的情况。
目标： 设节点\(a=x_0<x_1<\cdots<x_n=b\)上有\(f(x_i)=y_i,f'(x_i)=m_j(i=0,1,\cdots,n)\)。要求插值多项式\(H(x)\)满足条件：\(H(x_i)=y_i,H'(x_i)=m_i\)。
插值方法： 这里给出了\(2n+2\)个条件可以唯一确定一个次数不超过\(2n+1\)的多项式
\[H_{2n+1}(x)=a_0+a_1x+\cdots+a_{2n+1}x^{2n+1}
\]
如果直接解方程组求系数，计算量太大。这里仍然基于求拉格朗日插值基函数的想法来构造埃尔米特插值多项式，即埃尔米特插值应满足下述形式
\[H_{2n+1}(x)=\sum_{i=0}^n[y_i\alpha_i(x)+m_i\beta_i(x)]
\]
其中\(\alpha_i(x),\beta_i(x)\)都是\(2n+1\)次多项式，且满足下列条件：
\[\begin{aligned}
&(1)\;\alpha_i(x_i)=1\\
&(2)\;\alpha_i(x_j)=0(j\ne i)\\
&(3)\;\alpha_i'(x_k)=0(k=0,1,\cdots,n)\\
&(4)\;\beta_i'(x_i)=1\\
&(5)\;\beta_i'(x_j)=0(j\ne i)\\
&(6)\;\beta_i(x_k)=0(k=0,1,\cdots,n)
\end{aligned}
\]
令\(\alpha_i(x)=(ax+b)l_i^2(x)\)，其中\(l_i(x)\)是拉格朗日插值基函数，可以验证此时\(x\ne{x_i}\)时，条件\((2),(3)\)均是满足的，于是只需令
\[\left\{\begin{aligned}
\alpha_i(x_i)&=1\\
\alpha_i'(x_i)&=0
\end{aligned}\right.\Rightarrow
\left\{\begin{aligned}
&ax_i+b=1\\
&a+2(ax_i+b)l'_i(x_i)=0
\end{aligned}\right.\Rightarrow
\left\{\begin{aligned}
a&=-2l'_i(x_i)\\
b&=1+2x_il'_i(x_i)
\end{aligned}\right.
\]
因此\(\alpha_i(x)=[1-2(x-x_i)l_i'(x_i)]l_i^2(x)\)，现在只需求得\(l'_i(x_i)\)即可。
对\(l_i(x)\)，取对数可以得到
\[\ln{l_i(x)}=\ln{C}+\sum_{j=0\\j\ne i}^n(x-x_j)
\]
其中\(C\)只是一个系数，两边求导即可得到
\[\frac{l_i'(x)}{l_i(x)}=\sum_{j=0\\j\ne i}^n\frac{1}{x-x_j}\;\Rightarrow\;
l'_i(x)=l_i(x)\sum_{j=0\\j\ne i}^n\frac{1}{x-x_j}\;\Rightarrow\;
l'_i(x_i)=\sum_{j=0\\j\ne i}^n\frac{1}{x_i-x_j}
\]
于是
\[\alpha_i(x)=\left[1-2(x-x_i)\sum_{j=0\\j\ne i}^n\frac{1}{x_i-x_j}\right]l_i^2(x)
\]
令\(\beta_i(x)=(x-x_i)l_i^2(x)\)，可以验证\(\beta_i(x)\)已经满足条件\((4),(5),(6)\)。
综上可知，埃尔米特插值公式如下：
\[\left\{\begin{aligned}
&H_{2n+1}(x)=\sum_{i=0}^n[y_i\alpha_i(x)+m_i\beta_i(x)]\\
&\alpha_i(x)=\left[1-2(x-x_i)\sum_{j=0\\j\ne i}^n\frac{1}{x_i-x_j}\right]l_i^2(x)\\
&\beta_i(x)=(x-x_i)l_i^2(x)
\end{aligned}\right.\quad i=0,1,2,\cdots,n
\]
插值余项：
仿照证明拉格朗日插值余项的方法可以证明，埃尔米特插值余项为
\[R_{2n+1}(x)=f(x)-H_{2n+1}(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}\omega_{n+1}^2(x)
\]
其中\(\xi\in(a,b)\)并且与\(x\)的值相关。
证明：构造辅助函数
\[\varphi(t)=f(t)-H_{2n+1}(t)-k(x)\omega_{n+1}^2(t)
\]
于是\(x,x_0,x_1,\cdots,x_n\)均是\(\varphi(t)\)的零点，根据罗尔定理可知\(\varphi'(t)\)会有异于\(x,x_0,x_1,\cdots,x_n\)的另外\(n+1\)个零点，由于\(x_0,x_1,\cdots,x_n\)仍然是\(\varphi'(t)\)的零点，因此\(\varphi'(t)\)共有\(2n+2\)个零点。反复利用罗尔定理最终可以得到\(\varphi^{(2n+2)}(t)\)在区间\((a,b)\)上有一个零点\(\xi\)，于是
\[0=\varphi^{(2n+2)}(\xi)=f^{(2n+2)}(\xi)-0-(2n+2)!k(x)\;\Rightarrow\;
k(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!},\xi\in(a,b)
\]
最终可以得到
\[\varphi(x)=0\;\Rightarrow\;
R_{2n+1}(x)=f(x)-H_{2n+1}(x)=\frac{f^{(2n+2)}(\xi)}{(2n+2)!}\omega_{n+1}^2(x)
\]
唯一性：
通过反证法可以证明其唯一性，假设\(H_{2n+1}(x)\)和\(\bar{H}_{2n+1}(x)\)都是满足埃尔米特插值条件的函数，令\(\phi(x)=H_{2n+1}(x)-\bar{H}_{2n+1}(x)\)，那么
\[\begin{aligned}
&\phi(x_k)=H_{2n+1}(x_k)-\bar{H}_{2n+1}(x_k)=y_k-y_k=0\\
&\phi'(x_k)=H_{2n+1}'(x_k)-\bar{H}_{2n+1}'(x_k)=m_k-m_k=0
\end{aligned}\quad k=0,1,\cdots,n
\]
这说明\(\phi(x)\)在每一个插值节点\(x_k\)上都有2重根，即\(\phi(x)\)至少有\(2n+2\)个根（包含重根），但是\(\phi(x)\)是次数不超过\(2n+1\)的多项式，因此只能有\(\phi(x)=0\)。至此，埃尔米特插值多项式的唯一性得证。
7.5 分段低次插值
高次插值计算量大，然而精度却不一定好，根本原因在于收敛性不能够得到保证，即下述极限式
\[\forall x\in[a,b],\lim_{n\rightarrow\infty}P_n(x)=f(x)
\]
对于某些函数是不成立的。
例如，龙格（Runge）就给出了一个不收敛的实例：对于函数\(f(x)=1/(x^2+1)\)在\([-5,5]\)上的各阶导数均存在。如果在该区间上取\(n+1\)个等间距点\(x_k=-5+10k/n\)进行插值，所构造的拉格朗日插值多项式当\(n\rightarrow\infty\)时，仅在\(|x\le3.63|\)内收敛，而在其他位置发散。并且如果\(n\)越大，在区间端点处函数值抖动越厉害，这种现象称为龙格现象。
因此随着插值节点数的增加，插值多项式的次数也相应增加，而如果出现不收敛的情况，那么高次插值将带来函数值的剧烈震荡，从而使得数值不稳定。为了增加插值节点的个数，增加插值精度，又想避免由于高次插值带来的数值不稳定，就可以采用分段插值的方法。
定义 设\(f(x)\)是定义在\([a,b]\)上的函数，在\([a,b]\)上的节点
\[a=x_0\le x_1\le x_2\le\cdots\le x_n=b
\]
的函数值分别为\(y_0,y_1,\cdots,y_n\)，若函数\(\varphi(x)\)满足
- 在区间\([a,b]\)上连续
- 在每个子区间\([x_i,x_{i+1}]\;(i=0,1,\cdots,n-1)\)上是次数为\(m\)的多项式。
则称\(\varphi(x)\)是\(f(x)\)在\([a,b]\)上的分段\(m\)次插值多项式。
特殊地，当\(m=1\)时称为分段线性插值；当\(m=2\)时，称为分段抛物线插值。
分段线性插值：
分段线性插值就是每一个区间内都用一条线段来代替，因此实际上就是一个折线图。如果被插值函数\(f(x)\)本身是连续函数，那么分段线性插值是能够保证收敛性的。
如果\(f(x)\)具有连续二阶导数，那么每一个小区间\([x_i,x_{i+1}]\)上的插值余项为
\[R_i(x)=\frac{f''(\xi_i)}{2}(x-x_i)(x-x_{i+1}),\xi_i\in[x_i,x_{i+1}]
\]
令\(h(x)=(x-x_i)(x-x_{i+1})\)，因为在区间\([x_i,x_{i+1}]\)上有
\[-\frac{(x_{i+1}-x_i)^2}{4}\le h(x)\le0
\]
假设\(M=\max_{a\le x\le b}{f''(x)},h=\max_{0\le i\le n-1}{(x_{i+1}-x_i)}\)，那么总的插值余项
\[R(x)\le\frac{M}{8}h^2
\]
分段低次插值的优点是计算简便，收敛性有保证，但是却不能保证曲线的光滑性（即导数的连续性）。例如，分段线性插值在插值节点的一阶导数通常是不连续的。
7.6 三次样条插值
一方面为了保留低次插值的优点，又要保证插值曲线的光滑性，于是就发展出来了样条插值（spline）方法。样条一词来自工程船体和汽车等外形设计中，给出外形曲线上的一组离散点（样点），将具有弹性细长钢条（样条）在样点固定住，使其他地方自由弯曲，这样的所表示的曲线称为样条曲线（函数）。
在数学上，它近似于一条分段的三次多项式，并且在每一个插值节点上都有一阶和二阶连续导数。
定义 给定\([a,b]\)上的\(n+1\)个点：\(a=x_0<x_1<\cdots<x_n=b\)，及其对应的函数值\(y_i=f(x_i)\)，若\(S(x)\)满足
- \(S(x_i)=y_i,\;(i=0,1,\cdots,n)\)
- \(S(x)\)在\([x_i,x_{i+1}]\)上是三次多项式
- \(S(x)\in C^2[a,b]\)
则称\(S(x)\)是\(f(x)\)关于节点\(i=0,1,\cdots,n\)的三次样条插值函数，称\(x_i\)为样条节点。
\(S(x)\)在区间\([x_i,x_{i+1}]\)上的表达式为
\[S_i(x)=a_ix^3+b_ix^2+c_ix+d_i,\;i=0,1,2,\cdots,n-1
\]
一共有\(4n\)个未知数，因此需要有\(4n\)个约束条件。
- 定义1给出了\(n+1\)个约束。
- 定义3给出了\(3n-3\)个约束，即在点\(x_1,x_2,\cdots,x_{n-1}\)处\(S(x),S'(x),S''(x)\)都是连续的。
因此，定义只给出了\(4n-2\)个约束条件，如果要唯一确定一个三次样条函数，还需要2个额外的约束才行。我们把这2个额外的约束称为边界条件，常用的边界条件有三种：
-
第一种边界条件（固支边界条件）：\(S'(x_0)=f'_0,S'(x_n)=f'_n\)。
-
第二种边界条件：\(S''(x_0)=f''_0,S''(x)=f''_n\)。特殊地，当\(f''_0=f''_n=0\)称为自然边界条件。
-
第三种边界条件（周期边界条件）：当\(f(x)\)是以\(x_n-x_0\)为周期的周期函数时，要求\(S(x)\)也是周期函数，那么
\[S'(x_0+0)=S'(x_n-0)\\
S''(x_0+0)=S''(x_n-0)
\]
这种情况下要求\(y_0=y_n\)。
三弯矩法求三次样条插值函数：
三弯矩法是先假设样条节点上的二阶导数值，然后推导出它们所满足的线性方程组，最后求得每个小区间上的表达式。
在\([x_i,x_{i+1}](i=0,1,\cdots,n-1)\)上进行分段线性插值：
\[\begin{aligned}
S_i''(x)&=M_i\frac{x_{i+1}-x}{x_{i+1}-x_i}+M_{i+1}\frac{x-x_i}{x_{i+1}-x_i}\\
&=M_i\frac{x_{i+1}-x}{h_i}+M_{i+1}\frac{x-x_i}{h_i}
\end{aligned}
\]
其中\(h_i=x_{i+1}-x_i,S''(x_i)=M_i\)。上述式子实际上是保证了样条函数具有连续的二阶导数。
接着将上式积分可以得到
\[S'_i(x)=-M_i\frac{(x_{i+1}-x)^2}{2h_i}+M_{i+1}\frac{(x-x_i)^2}{2h_i}+C_1-C_2
\]
再次积分可以得到
\[S_i(x)=M_i\frac{(x_{i+1}-x)^3}{6h_i}+M_{i+1}\frac{(x-x_i)^3}{6h_i}+C_1(x-x_i)+C_2(x_{i+1}-x)
\]
利用\(S(x_i)=y_i,S(x_{i+1})=y_{i+1}\)，可以得到
\[C_1=\frac{1}{h_i}(y_{i+1}-\frac{M_{i+1}h_i^2}{6}),\;
C_2=\frac{1}{h_i}(y_{i}-\frac{M_{i}h_i^2}{6})
\]
此时可以保证样条函数本身的连续性和样点上的函数值满足要求，因此上述式子还需要满足的条件是具有一阶连续导数。
\[S'_i(x)=-M_i\frac{(x_{i+1}-x)^2}{2h_i}+M_{i+1}\frac{(x-x_i)^2}{2h_i}+
\frac{y_{i+1}-y_i}{h_i}-\frac{h_i}{6}(M_{i+1}-M_i)
\]
利用一阶导数连续的条件：\(S'_{i-1}(x_i-0)=S'_i(x_i+0)\)，注意这里要求\(i=1,2,\cdots,n-1\)，于是
\[M_{i}\frac{(x_i-x_{i-1})^2}{2h_{i-1}}+\frac{y_i-y_{i-1}}{h_{i-1}}-\frac{h_{i-1}}{6}(M_{i}-M_{i-1})=
-M_i\frac{(x_{i+1}-x_i)^2}{2h_i}+\frac{y_{i+1}-y_i}{h_i}-\frac{h_i}{6}(M_{i+1}-M_i)
\]
注意到\(h_i=x_{i+1}-x_i\)整理可得
\[\begin{aligned}
\frac{h_{i-1}}{6}M_{i-1}+\frac{h_{i-1}+h_i}{3}M_i+\frac{h_i}{6}M_{i+1}
&=\frac{y_{i+1}-y_i}{h_i}-\frac{y_{i}-y_{i-1}}{h_{i-1}}\\
&=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-\frac{y_{i}-y_{i-1}}{x_i-x_{i-1}}\\
&=f[x_i,x_{i+1}]-f[x_{i-1},x_i]
\end{aligned}
\]
等式两侧同乘以\(6/(h_i+h_{i-1})\)得到
\[\begin{aligned}
\frac{h_{i-1}}{h_i+h_{i-1}}M_{i-1}+2M_i+\frac{h_i}{h_i+h_{i-1}}M_{i+1}
&=6\frac{f[x_i,x_{i+1}]-f[x_{i-1},x_i]}{h_i+h_{i-1}}\\
&=6\frac{f[x_{i-1},x_i]-f[x_i,x_{i+1}]}{x_{i-1}-x_{i+1}}\\
&=6f[x_{i-1},x_i,x_{i+1}]
\end{aligned}
\]
令
\[\left\{\begin{aligned}
\mu_i&=\frac{h_{i-1}}{h_i+h_{i-1}}\\
\lambda_i&=1-\mu_i=\frac{h_{i}}{h_i+h_{i-1}}\\
d_i&=6f[x_{i-1},x_i,x_{i+1}]
\end{aligned}\right.
\]
于是可以得到
\[\mu_iM_{i-1}+2M_i+\lambda_iM_{i+1}=d_i,\quad i=1,2,\cdots,n-1\tag{7-6-1}
\]
上面\(n-1\)个方程共有\(n+1\)个未知数：\(M_0,\cdots,M_n\)，因此最后我们只需要引入边界条件即可解出所有的未知数。
第一种边界条件：
根据第一种边界条件可以得到
\[S'_0(x_0)=-\frac{h_0}{3}M_0-\frac{h_0}{6}M_1+\frac{y_1-y_0}{h_0}=f'_0\\
S'_{n-1}(x_n)=\frac{h_{n-1}}{3}M_n-\frac{h_{n-1}}{6}M_{n-1}+\frac{y_n-y_{n-1}}{h_{n-1}}=f'_n
\]
于是
\[2M_0+M_1=\frac{6}{h_0}(f[x_0,x_1]-f'_0)\\
M_{n-1}+2M_2=\frac{6}{h_{n-1}}(f'_n-f[x_{n-1},x_n])\tag{7-6-2}
\]
令
\[\left\{\begin{aligned}
\lambda_0&=1\\\mu_n&=1\\
d_0&=\frac{6}{h_0}(f[x_0,x_1]-f'_0)\\
d_n&=\frac{6}{h_{n-1}}(f'_n-f[x_{n-1},x_n])
\end{aligned}\right.
\]
将式(7-6-1)，(7-6-2)共\(n+1\)个方程改写成矩阵的形式，则有
\[\left[\begin{matrix}
2&\lambda_0&&&\\
\mu_1&2&\lambda_1&&\\
&\ddots&\ddots&\ddots&\\
&&\mu_{n-1}&2&\lambda_{n-1}\\
&&&\mu_n&2
\end{matrix}\right]
\left[\begin{matrix}
M_0\\M_1\\\vdots\\M_{n-1}\\M_n
\end{matrix}\right]=
\left[\begin{matrix}
d_0\\d_1\\\vdots\\d_{n-1}\\d_n
\end{matrix}\right]
\]
这是一个三对角线方程组，可以用追赶法求解。求出\(M_0,\cdots,M_n\)，即可求出样条函数\(S(x)\)的分段表达式。
第二种边界条件：
根据第二种边界条件可知：\(S_0''(x_0)=M_0=f''_0,S_{n-1}''(x_n)=M_n=f''_n\)。于是仅剩下\(M_1,\cdots,M_{n-1}\)未知，将式(7-6-1)改写成矩阵形式，则有
\[\left[\begin{matrix}
2&\lambda_1&&&\\
\mu_2&2&\lambda_2&&\\
&\ddots&\ddots&\ddots&\\
&&\mu_{n-2}&2&\lambda_{n-2}\\
&&&\mu_{n-1}&2
\end{matrix}\right]
\left[\begin{matrix}
M_1\\M_2\\\vdots\\M_{n-2}\\M_{n-1}
\end{matrix}\right]=
\left[\begin{matrix}
d_1-\mu_1f''_0\\d_2\\\vdots\\d_{n-2}\\d_{n-1}-\lambda_{n-1}f''_n
\end{matrix}\right]
\]
同样也可以通过追赶法求解。
第三种边界条件：
\[\begin{aligned}
&S''_0(x_0+0)=S''_{n-1}(x_n-0)\Rightarrow M_0=M_n\\
&S'_0(x_0+0)=S'_{n-1}(x_n-0)\Rightarrow \mu_nM_{n-1}+2M_n+\lambda_nM_1=d_n
\end{aligned}\tag{7-6-3}
\]
其中
\[\lambda_n=\frac{h_0}{h_{n-1}+h_0},\mu_n=1-\lambda_n=\frac{h_{n-1}}{h_{n-1}+h_0}\\
d_n=6\frac{f[x_0,x_1]-f[x_{n-1},x_n]}{h_{n-1}+h_0}
\]
将式(7-6-1)与(7-6-3)联立写成矩阵形式
\[\left[\begin{matrix}
2&\lambda_1&&&\mu_1\\
\mu_2&2&\lambda_2&&\\
&\ddots&\ddots&\ddots&\\
&&\mu_{n-1}&2&\lambda_{n-1}\\
\lambda_n&&&\mu_{n}&2
\end{matrix}\right]
\left[\begin{matrix}
M_1\\M_2\\\vdots\\M_{n-1}\\M_{n}
\end{matrix}\right]=
\left[\begin{matrix}
d_1\\d_2\\\vdots\\d_{n-1}\\d_{n}
\end{matrix}\right]
\]
综上可知，求解三次样条函数就需要求解对应的线性方程组，并且三种不同边界条件的方程组的系数矩阵都是严格对角占优的。
为什么叫三弯矩法：
在力学上\(M_i\)表示细梁在\(x_i\)处的截面弯矩，可以看到上面三种边界条件的方程组中每一个方程最多出现三个弯矩，因此各个方程称为三弯矩方程。因此，这种方法也称为三弯矩法。
浙公网安备 33010602011771号