数值计算-多项式插值方法

Dec 1, 2020 · 数值计算 ·

分享到:

多项式插值方法

本笔记框架主要来自于知乎：https://www.zhihu.com/question/22320408,作者：最爱麦丽素

多项式插值综述
Lagrange插值法
重心拉格朗日插值法
Newton插值法
Newton插值法余项与泰勒公式
切比雪夫节点与勒贝格常数
切比雪夫节点
勒贝格常数
Hermite插值
插值法的变形
附录
克莱姆法则证明

多项式插值综述

给定（n+1）个数据点\((x_i,y_i),i\in\{0,...,n\}\)，且任意两个\(x_i\)都不相等。可以用一个p阶多项式\((0\leq p \leq n)\)进行插值，使： \[p(x_i)=y_i,0 \leq i \leq n\] \(p(x)\)具体可以写为： \[p(x)=a_nx^n+a_{n-1}x^{n-1}+\dotsb+a_1x+a_0\] 在\(p(x)\)中最多会有n阶多项式，n+1个未知参数，而给定的n+1个数据点\((x_i,y_i)\)，和这些未知数恰好可以构成n+1元一次方程组。如果用矩阵的形式可以写成： \[\begin{bmatrix}x_0^n&x_0^{n-1}&\dotsb&x_0&1\\ x_1^n&x_1^{n-1}&\dotsb&x_1&1\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ x_n^n&x_n^{n-1}&\dotsb&x_n&1\\ \end{bmatrix} \begin{bmatrix}a_n\\ a_{n-1}\\ \vdots\\ a_0\\ \end{bmatrix}= \begin{bmatrix}y_0\\ y_1\\ \vdots\\ y_n\\ \end{bmatrix} \] 简写为\(X\vec{a}=\vec{y}\)。我们所要求的就是多项式的系数向量\(\vec{a}\)，最直观的方式当然是求逆矩阵，从而得到\(\vec{a}=X^{-1}\vec{y}\)。由于X为范德蒙矩阵，且\(x_i,x_j\)两两不相等，这就是的范德蒙矩阵\(X\)必然可逆且唯一（行列式不为0）。

但是计算大矩阵的逆是一项艰巨的工作，因此我们可以使用克莱姆法则(Cramer‘s Rule) （证明见附录），直接求出每一项系数，即： \[a_i=\frac {\begin{vmatrix} x_0^n&x_0^{n-1}&\dotsb&y_0&\dotsb&x_0&1\\ x_1^n&x_1^{n-1}&\dotsb&y_1&\dotsb&x_1&1\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ x_n^n&x_n^{n-1}&\dotsb&y_n&\dotsb&x_n&1\\ \end{vmatrix}} {\begin{vmatrix} x_0^n&x_0^{n-1}&\dotsb&x_0^i&\dotsb&x_0&1\\ x_1^n&x_1^{n-1}&\dotsb&x_1^i&\dotsb&x_1&1\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ x_n^n&x_n^{n-1}&\dotsb&x_n^i&\dotsb&x_n&1\\ \end{vmatrix}}(0 \leq i \leq n)\] 分子中被替换的是范德蒙行列式的第i列。且由于范德蒙矩阵的性质，x只要两两不相同，则行列式不等于0。

从这个式子中，我们也可以发现无论是求矩阵的逆，还是用克莱姆法则，他们的结果是一样且唯一的，只是计算量复杂。因此后来改进的Newton法，Lagrange法等等都只是降低了计算量，求出来的多项式系数都是唯一的。

Lagrange插值法

Lagrange插值法利用了基函数构造的特点，用基函数的线性组合来组成插值多项式。从这个角度来看，Lagrange插值法可以用以下式子表示： \[L(x)=a_0l_0(x)+a_1l_1(x)+a_2l_2(x)+\dotsb+a_nl_n(x)\] 其中\(l_i(x)\)就是基函数，\(a_i\)为基函数线性组合的系数。

因为\(L(x)\)需要经过每一个点\((x_i,y_i)\)，即\(y_i=L(x_i)\)。我们可以让\(a_i\)直接等于\(y_i\), 那我们就只需要构造对应的\(l_i(x_i)=1,l_i(x_j)=0\)就好了。Lagrange精巧第构造了这些基函数： \[l_i(x)=\prod_{j=0,j\neq i}^n\frac{(x-x_j)}{(x_i-x_j)}\] \[L(x)=\sum_{i=0}^n{y_i \prod_{j=0,j\neq i}^n\frac{(x-x_j)}{(x_i-x_j)}}\] 从表达式中，我们不难得到以下性质：

\(i \neq j\)，则\(l_i(x_j)=0\)，因为\(\exist x_j-x_j=0,\)使得分子为0
\(i=j\)，则\(l_i(x_j)=1\)，因为分子和分母完全相同。
\(L(x)\)的阶数最大为n。\(\prod_{0 \leq j \leq n,j\neq i}\frac{x-x_j}{x_i-x_j}\)一共有n个x相乘。
插值多项式是唯一的，和范德蒙矩阵求出来的是一样的。

拉格朗日插值法的公式结构整齐紧凑，在理论分析中十分方便，然而在计算中，当插值点增加或减少一个时，所对应的基本多项式就需要全部重新计算，于是整个公式都会变化，非常繁琐。这时可以用重心拉格朗日插值法或牛顿插值法来代替。此外，当插值点比较多的时候，拉格朗日插值多项式的次数可能会很高，因此具有数值不稳定的特点，也就是说尽管在已知的几个点取到给定的数值，但在附近却会和“实际上”的值之间有很大的偏差（如图1）。这类现象也被称为龙格现象，解决的办法是分段用较低次数的插值多项式。

图1 拉格朗日插值法的数值稳定性：如图，用于模拟一个十分平稳的函数时，插值多项式的取值可能会突然出现一个大的偏差（图中的14至15中间）

重心拉格朗日插值法

来源维基百科：https://zh.wikipedia.org/wiki/%E6%8B%89%E6%A0%BC%E6%9C%97%E6%97%A5%E6%8F%92%E5%80%BC%E6%B3%95

重心拉格朗日插值法是拉格朗日插值法的改进。在拉格朗日插值多项式中，运用多项式： \[l(x)=(x-x_0)(x-x_1)\dotsb (x-x_k)\] 可以将原拉格朗日多项式的项改写： \[l_j(x)=\frac{l(x)}{x-x_j}\frac{1}{\prod_{i=0,i\neq j}^k(x_j-x_i)}\] 如果我们将重心权定义为： \[w_j=\frac{1}{\prod_{i=0,i\neq j}^k(x_j-x_i)}\] 拉格朗日多项式的项可以变为： \[l_j(x)=l(x)\frac{w_j}{x-x_j}\] 同样，原拉格朗日多项式可以把\(l(x)\)提取出来： \[L(x)=l(x)\sum_{j=0}^k\frac{w_j}{x-x_j}y_j\] 即所谓的重心拉格朗日插值公式第一型或改进拉格朗日插值公式。它的优点是当插值点的个数增加一个时，将每个\(w_{j}\)都除以\((x_{j}-x_{k+1})\)，就可以得到新的重心权\(w_{k+1}\)，计算复杂度为\(O(n)\)，比重新计算每个基本多项式所需要的复杂度\(O(n^{2})\)降了一个量级。

将以上的拉格朗日插值多项式用来对函数\(g(x)\equiv 1\)插值，可以得到： \[\because y_i \equiv 1 \\ \therefore \forall x,g(x)=l(x)\sum_{j=0}^k\frac{w_j}{(x-x_j)} \] 我们用\(L(x)\)除以\(g(x) \equiv 1\)，可以得到： \[L(x)=\frac{\sum_{j=0}^k\frac{w_j}{x-x_j}y_j}{\sum_{j=0}^k\frac{w_j}{(x-x_j)}}\] 这个公式被称为重心拉格朗日插值公式（第二型）或真正的重心拉格朗日插值公式。它继承了第一型容易计算的特点，并且在代入\(x\)值计算\(L(x)\)的时候不必计算多项式\(l(x)\)。它的另一个优点是，结合切比雪夫节点进行插值的话，可以很好地模拟给定的函数，使得插值点个数趋于无穷时，最大偏差趋于零。同时，重心拉格朗日插值结合切比雪夫节点进行插值可以达到极佳的数值稳定性。第一型拉格朗日插值是向后稳定的，而第二型拉格朗日插值是向前稳定的，并且勒贝格常数很小。

前向稳定：算法的前向误差是结果与真解之间的差别，即\(\Delta y=y^{*}-y\)。

后向稳定：后向误差是满足\(f(x+\Delta x)=y^{*}\)的最小\(\Delta x\)，也就是说后向误差说明算法的所解决的问题。如果对于任意的输入\(x\)来说后向误差都很小，那么算法就是后向稳定的。

前向误差和后向误差通过条件数发生关系：前向误差的幅度最多是条件数乘以后向误差的幅度。

Newton插值法

回头看看Lagrange的基函数，如果我们能找到一组新的基函数，使得新增的基函数对原有函数进行修正，就能防止有新节点插入时Lagrange的重复计算。考虑以下基函数。 \[\begin{aligned} \phi_0(x)&=1\\ \phi_1(x)&=(x-x_0)\\ \phi_2(x)&=(x-x_0)(x-x_1)\\ \dotsb&=\dotsb\\ \phi_{n+1}(x)&=\prod_{i=0}^n(x-x_i)\\ \end{aligned} \] 根据线性无关的基函数我们可以用其线性组合得到待定系数的插值多项式： \[p(x)=a_0\phi_0(x)+a_1\phi_1(x)\dotsb+a_n\phi_n(x)\] 首先，对于1个点\((x_0,y_0),p_0(x)=a_0=y_0\)，显然。

对于2个点\((x_0,y_0),(x_1,y_1)\)。 \[p_1(x)=y_0+a_1(x-x_0)\\ a_1=\frac{y_1-y_0}{x_1-x_0}\]

对于三个点\((x_0,y_0),(x_1,y_1),(x_2,y_2)\)。 \[p_2(x)=y_0+\frac{y_1-y_0}{x_1-x_0}+a_2\phi_2(x)\\ a_2=\frac{\frac{y_2-y_1}{x_2-x_1}-\frac{y_1-y_0}{x_1-x_0}}{x_2-x_0}\]

根据系数的特定，我们定义均差（或差商）： >一阶均差： >\[f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}\] >二阶均差：一阶均差的均差。 >\[f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}\] >N阶均差：N-1阶均差的均差。

均差具有对称性：均差\(f[x_0,x_1,\dotsb,x_k]\)与插值节点的顺序无关，即 \[f[x_0,x_1,\dotsb,x_k]=f[x_{i_0},x_{i_1},\dotsb,x_{i_k}]\] 其中，\(x_{i_0},x_{i_1},\dotsb,x_{i_k}\)是\(x_0,x_1,\dotsb,x_k\)的任意一个排列。

根据上述定义和规律，我们知道对于n+1个节点，其插值多项式为： \[\begin{aligned} p_n(x)=&f(x_0)+\\ &f[x_1,x_0](x-x_0)+\\ &f[x_2,x_1,x_0](x-x_0)(x-x_1)+\\ &f[x_3,x_2,x_1,x_0](x-x_0)(x-x_1)+(x-x_2)+\\ &\dotsb \end{aligned}\] 每新增一个节点\((x_{n+1},y_{n+1})\)时，只需要增加一项新的基函数\(\phi_{n+1}=(x-x_0)(x-x_1)\dotsb(x-x_n)\)和新的系数\(f[x_0,x_1,\dotsb,x_{n+1}]\)即可。其中， \[f[x_0,x_1,\dotsb,x_{n+1}]=\frac{f[x_0,x_1,\dotsb,x_n]-f[x_1,x_2,\dotsb,x_{n+1}]}{x_0-x_{n+1}}\] 下面这个图显示了均差表的计算，由于实际运算中，这些都是数值，因此计算并不困难。也是O(n)的计算复杂度。

均差计算1

新增一个插值点，只需要计算相关的差分就可以了。

均差计算2

Newton插值法余项与泰勒公式

泰勒把牛顿插值法做了一些改造。首先，设\(f(x)\)是一个函数，它在\(x_0,x_0+\Delta x,x_0+2\Delta x,x_0+3\Delta x,\cdots,x_0+n\Delta x\)的值已知（和之前的相比，相当于每个点都是等距离间隔的，间隔\(\Delta x\)），令： \[ \Delta f(x_0)=f(x_0+\Delta x)-f(x_0)，也称为一阶差分，\] \[\Delta f(x_0+\Delta x)=f(x_0+2\Delta x)-f(x_0+\Delta x), \] 进一步令： \[\Delta^2 f(x_0)=\Delta f(x_0+\Delta x)-\Delta f(x_0)，也称为二阶差分（为一阶差分的差分）\\ \Delta^3 f(x_0)=\Delta^2 f(x_0+\Delta x)-\Delta^2 f(x_0)，也称为三阶差分。\] 做了这些假设之后我们来看看之前提到的\(a_1\)会变成什么样子： \[a_1=\frac{f(x_1)-f(x_0)}{x_1-x_0}\implies a_1=\frac{\Delta f(x_0)}{\Delta x}。\] 而\(f_1(x)\)会变成： \[f_1(x)=f(x_0)+\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0)\\ \implies f_1(x)=f(x_0)+\frac{\Delta f(x_0)}{\Delta x}(x-x_0)。\] 同样的\(f_2(x)\)就变成了： \[f_2(x)=f(x_0)+\frac{\Delta f(x_0)}{\Delta x}(x-x_0)+\frac{\Delta^2 f(x_0)}{2\Delta x}(x-x_0)(x-x_1)。\] 泰勒断言，当\(\Delta x=0\)时： \[f_1(x)=f(x_0)+f'(x_0)(x-x_0)\\ f_1(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2\\ （\Delta x=0时有x-x_1=x-x_0）\] 以此类推泰勒就得到了大名鼎鼎的泰勒公式： \[f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\cdots\]

切比雪夫节点与勒贝格常数

https://zh.wikipedia.org/wiki/%E5%A4%9A%E9%A1%B9%E5%BC%8F%E6%8F%92%E5%80%BC

切比雪夫节点与勒贝格常数都是在没有给定插值节点时，自己选择插值节点的准则。

切比雪夫节点

对于一个插值区间\([a, b]\) 如果要在这个插值区间上选取\(n\)个点作为插值基点，使得上面的最大误差最小，则基点的选法如下: \[x_i=\frac{b+a}{2}+\frac{b-a}{2}\cos\frac{(2i-1)\pi}{2n}\\ (i=1,2,\dotsb,n)\] 这些节点称为切比雪夫（插值）节点。

勒贝格常数

在插值节点\(x_0、\dotsb、x_n\)以及包含所有插值节点的区间\([a, b]\)确定下来，插值的过程就是将函数\(f\)映射到多项式\(p\)。这定义了一个\([a, b]\)区间上所有连续函数从空间\(C([a, b])\)到其自身的映射。映射\(X\)是线性的，并且它是函数\(f\)在子空间\(Π_n\)上的投影。勒贝格常数\(L\)定义为\(X\)的算子范数，它满足勒贝格引理： \[\|f-X(f)\|\leq (L+1)\|f-p^{*}\|.\] 换句话说，插值多项式的误差最多是最优逼近的\(L+1\)倍，这表明我们要找使\(L\)最小的插值节点。尤其是选择切比雪夫节点时： \[L \ge \frac{2}{\pi}log(n+1)+C \quad for\quad some\quad costant\quad C \] 我们再次证明切比雪夫节点是多项式插值中比较好的选择，但是这些节点并不是最优的。

Hermite插值

埃尔米特插值是另一类插值问题，这类插值在给定的节点处，不但要求插值多项式的函数值与原函数值相同。同时还要求在节点处，插值多项式的一阶直至指定阶的导数值，也与被插函数的相应阶导数值相等，这样的插值称为埃尔米特(Hermite)插值。 Hermite插值在不同的节点，提出的插值条件个数可以不同，若在某节点\(x_{i}\)，要求插值函数多项式的函数值，一阶导数值，直至\(m_{i}-1\)阶导数值均与被插函数的函数值相同及相应的导数值相等。我们称\(x_{i}\)为\(m_{i}\)重插值点节,因此，Hermite插值应给出两组数，一组为插值点\(\{x_{i}\}_{i=0}^{n}\)节点，另一组为相应的重数标号\(\{m_{i}\}_{i=0}^{n}\)。

插值法的变形

讲到这里基本上所有插值都是上面这些方法的结合体或者扩展了，根据不同的应用我们会使用不同的算法来解决我们的问题，像已知很多点的一次导数和他的位置，一般只有在应用物体位置和速度都知道的情况下，而这也只有在方程次数很高的情况下才会应用高次Hermite插值。通常我们接触的插值应用也只在3次方程进行。比如当我们仅仅知道点的位置而想让曲线更加平滑，不考虑误差，用最普通的Newton来算就好。为了增加视觉上的平滑或者说和周围点的连贯行，我们可以根据已知相邻点的位置来计算当前点的导数，进而使用Hermite算法。

如果我们将所有的点分段进行计算，则变为样条插值。

对于样条插值的连续性，如果我们选择3次样条插值来满足一阶导数连续的话，实际上是满足不了2阶导数连续的，举个例子，\(x^3\)在y轴左侧，\(x^2\)在y轴右侧构成的曲线，在原点处的一阶导数都为0，而\(f_+''(0)=2,f_-''(0)=0\)就会导致看起来略微有些别扭的情况，但是实际应用中，这种不合适经常被忽略。

在多段插值中，已知点的位置，对于点上的导数的计算的选择的不同则变为不同的插值。插值法的变形

附录

克莱姆法则证明

来源：维基百科https://zh.wikipedia.org/wiki/%E5%85%8B%E8%90%8A%E5%A7%86%E6%B3%95%E5%89%87

克莱姆法则：一个线性方程组可以用矩阵与向量的方程来表示： \[Ax=c \quad(1)\] 其中的\(A\)是一个\(n\times n\)的方块矩阵，而向量\(x=( x_1, x_2, \cdots x_n )^T\)是一个长度为n的列向量\(c=( c_1, c_2, \cdots c_n )^T\)也一样。

克莱姆法则说明：如果\(A\)是一个可逆矩阵（\(\det{A} \neq 0\)），那么方程(1)有解\(x=( x_1, x_2, \cdots x_n )^T\)，其中

\[x_i = { \det(A_i) \over \det(A)} \quad(1)\]

当中\(A_{i}\)是被列向量\(c\)取代了\(A\)的第i列的列向量后得到的矩阵。为了方便，我们通常使用\(\Delta\)来表示\(\det(A)\)，用\(\Delta_i\)来表示\(\det(A_i)\)。所以等式(1)可以写成为： \[x_i = { \Delta_i \over \Delta }。\]

证明：

对于n元线性方程组\(Ax=c\)

把系数矩阵\(A\)表示成列向量的形式 \[A=\left(u_{1},u_{2},\cdots ,u_{n}\right)\] 由于系数矩阵可逆，故方程组一定有解\(x^{*}=A^{{-1}}c\).

设\(x^{*}=(x_{1},x_{2},\cdots ,x_{n})^{T}\)，即 \[Ax^{*}=\sum _{{k=1}}^{n}x_{k}u_{k}=c\]

考虑\(\Delta_i\)的值，利用行列式的线性和交替性质，有

\[{\begin{aligned} \Delta _{i}&=det\left(\cdots ,u_{{i-1}},c,u_{{i+1}},\cdots \right) \\&=det\left(\cdots ,u_{{i-1}},\sum _{{k=1}}^{n}x_{k}u_{k},u_{{i+1}},\cdots \right) \\&=\sum _{{k=1}}^{n}x_{k}\cdot det\left(\cdots ,u_{{i-1}},u_{k},u_{{i+1}},\cdots \right) \\&=x_{i}\cdot det\left(\cdots ,u_{{i-1}},u_{i},u_{{i+1}},\cdots \right) \\&=x_{i}\Delta \end{aligned}}\]

于是\(x_{i}=\frac{\Delta_{i}}{\Delta}\)。