更有matlab、python版本的代码在仓库中可以对照使用。仓库地址见文章开头

*拉格朗日插值余项

既然是拟合，总归是有误差的，接下来我们研究一下最终的误差。对于拉格朗日插值法，误差实际上是一个阶段误差\(R(x) = f(x) - p(x)\)，根据之前我们推导出的拉格朗日插值函数，误差余项记为\(R(x) = f(x) - p(x)\)。

那么我们作一辅助函数如下

容易发现\(g^{n+1}(t) = (n+1)!\)，而插值函数只有n阶，故\(p^{n+1}(t) = 0\)，因此我们构造的辅助函数是连续且具有n+1阶导数。当\(t=x\)时，则会有\(\varphi(t)=0\)，根据Rolle定理，则必有以下：

• 当\(\varphi(t)\)在插值区间内有n+2个互异零点时，则一阶导数至少有n+1个互异零点
• 当\(\varphi^{(1)}(t)\)在插值区间内有n+1个互异零点时，则\(\varphi^{(2)}(t)\)至少有n个互异零点

以此类推，容易知道\(\varphi^{(n)}(t)\)有2个零点时，则\(\varphi^{(n+1)}(t)\)至少有一零点，设此零点为\(t=\xi\)，对\(\varphi(t)\)求n+1次导数后，有

故我们插值的截断误差为\(f(x)-p(x) = \frac{f^{n+1}(\xi)}{(n+1)!}\prod_{i=0}^{n}(x-x_i)\)，其中\(\xi\)取决于我们需要预测的值。自此，我们得到了拉格朗日插值法的所有内容。

Newton差商插值

拉格朗日插值多项式的得出并不容易，至少在计算机眼中，他是复杂的，需要计算多个n-1阶的多项式。同时，这里我们给出一个定理：所有的插值多项式和拉格朗日插值多项式的结果应当是相同的，只是以不同的方式计算和书写。。

差商

首先我们介绍以下差商的概念，实际差商是一个很简单的东西，定义如下：假设顺序排列的节点\(x_0,x_1,\dotsb,x_n\)对应函数值为\(f(x_0),f(x_1),\dotsb,f(x_n)\)，那么一阶差商为

二阶差商为

那么K阶差商为：

k阶差商亦可写成

实际上就是个递推的公式，看起来似乎很复杂，但实际上自己尝试去计算的时候会发现，是一个很简单的东西，同时差商是函数的一个线性组合，这又可以说到线性代数了，我们假设组合出的多项式是一个类似坐标的玩意，而每一阶的差商实际上就是向量的基，因此通过差商，我们可以表示出多项式空间中任一多项式。

牛顿差商公式

牛顿差商公式就是将我们的差商做出组合：

随便给出几个插值采样点，进行计算拉格朗日和牛顿插值公式，我们在计算出各阶差商以后，会发现最终得到的插值多项式和拉格朗日多项式是一致的。那我们为什么还需要多此一举引入牛顿法进行插值而不是直接使用拉格朗日插值呢？

答案很简单，拉格朗日法在构建多项式时，我们必须先知道采样点的数量，如果后续有了新的采样点，我们就需要重新计算整个多项式，这是对于性能浪费极大的。如果采用牛顿法，那么在有新的采样点时，我们只需要多计算一次差商（因为之前的差商已经计算好），可以很容易的将新的点加入我们的多项式中。

牛顿法的代码如下（更多代码见仓库）

double Newton(double[] x, double[] y,double xi)
{
    if (x.Length != y.Length)
    {
        throw new ArgumentException("dim x not equal dim y");
    }
    int n = x.Length;
    double[,] args = new double[n,n];
    for (int i = 0; i < n; i++)
    {
        args[i,0] = y[i];
    }
    for (int i = 1; i <= n; i++)
    {
        for (int j = i; j <= n; j++)
        {
            args[j,i] = (args[j - 1,i - 1] - args[j,i - 1]) / (x[j - i] - x[j]);
        }
    }
    double result = y[0];
    for (int i = 1; i < x.Length; i++)
    {
        double mulValue = 1.0;
        for (int j = 0; j < i - 1; j++)
        {
            mulValue = mulValue * (xi - x[j]);
        }
        result += args[i, i] * mulValue;
    }
    return result;
}

插值余项

实际上，牛顿法因为得到的结果和拉格朗日插值一样，他们的误差也应该是一样的，不过既然我们选择了不同的计算方法，那么误差也可以使用新的手段进行计算。如果我们把我们的插值点看作一固定点，带入我们的差商公式，依次递推，会得到最后的插值余项。

高次插值的Runge(龙格)现象

在多项式插值以及基于多项式插值的拉格朗日插值方法和含有高次幂的牛顿插值中，我们提到过龙格现象，那么什么是龙格现象呢？首先我们先提出这样一个问题，当我们在拉格朗日多项式之中提高项数和多项式的次数，是否可以提高我们数据的精度？这个答案是否定的，我们需要减少误差，也就是说针对任意的\(\epsilon>0\)，总是存在正整数N，当n大于N时，所有区间内的\(x\)满足:

我们可以看到，我们的误差，随着n的增大，(n+1)!逐渐趋于无穷，而后向的乘法，并不会增大太多会变小太多，因此误差很大程度取决于n+1阶导数的范围。对于类似\(sinx,e^x\)，这类导数的值域范围恒定的函数，误差会随着n的增加逐步趋于0。但是如果对于某些函数而言，随着项数次数增加，误差反而会变大。龙格现象是一种出现在多项式次数很高是，在区间端点出现拟合曲线剧烈振荡的一种现象。下面通过一个例子来进行说明。

假设给定\(f(x) = \frac{1}{x^2+1},x\in [-5,5]\)，等距离取5，10个点，构建多项式拟合，观察图像。

我们可以很清楚的发现，当靠近端点的时候。拟合曲线出现了剧烈的波动，误差到了一个不太可以接受的程度。并且次数越高的时候，抖动出现的更频繁、更剧烈。龙格现象的出现是由于插值节点是等距所造成的。同时，针对导函数不稳定的函数而言，盲目增加插值节点数列，则更加容易出现龙格现象。

一般而言，为了避免龙格现象，我们可以有以下的手段：

• 避免使用高次的插值多项式，一般而言，即使为了提高插值精度，也尽可能保证\(n<7\)。
• 分段进行插值，在每一小段使用一次插值，从而降低阶数。
• 对于性态不稳定的函数，选用其他插值手段，如Hermite插值等。

Hermite插值

在之前的插值方法中，我们的要求仅仅是在满足\(f(x)=p(x)\)或者两者足够接近作为一个条件，但是这并不足以做出完美的拟合，在许多时候，我们往往还会要求导数及高阶导数尽可能相同或接近，这种拓展了条件的插值方法，就是我们的Hermite插值方法。

最常见的是我们通常会要求曲线的一阶导数接近或相同，这样可以保证我们曲线的一个平滑性能。那么假设曲线\(y=f(x)\)，一阶导数为\(y^{'} = f^{'}(x)\)，那么Hermite插值需要保证：

不难发现，我们不仅要求了相等的条件，更是要求了相切的条件，因此上述公式具有2(n+1)个采样点（曲线和曲线导数各一半）,那么我们仿造拉格朗日插值条件。可以构造一个这样的插值多项式：

其中，需要满足：

接下来我们用比较长的篇幅证明推导以下我们的Hermite插值的基函数\(a(x),b(x)\)，以及证明他的唯一性。

首先我们来求\(b(x)\)，它需要满足

易知满足上述条件的该点为二重根，所以必有因子\((x-x_i)^2 ,i \neq k\)，

易知满足上述条件的该点为单根，因此还会有另一个因子，\((x-x_i) ,i=k\)。又因为b(x)是次数不超过2n+1的多项式，根据拉格朗日多项式的定义，可以假定

上述式子满足我们之前的规定，同时因为\(b^{'}_i(x) = 1\)，可以轻松的确定常系数为1。故\(b_i(x) = (x-x_i)l_i^{2}(x)\)。接下来可以确定一下\(a(x)\)。

同样的，因为\(a(x)=0,a^{'}_i(x)=0,i\neq k\)说明也具有二重根，同样基函数会含有\(（x-x_i)^2,i\neq k\)，但略微不同的是，\(k=i\)的时候，P我们不能在获得任何因子凑够2n+2项。因此我们利用待定系数法假定基函数为

同样代入我们的条件，最终可以求得：

最后我们证明一下我们得到的插值多项式是唯一的，我们利用反证法：

假设存在另一个\(h(x)\)是满足Hermite插值条件，那么设：

由于有n+1个i，因此\(\varphi(x)\)有2n+2个零点，但是\(\varphi(x) = H(x) - h(x) = 0\)次数不超过2n+1，所以\(\varphi(x) = 0\)，因此\(H(x) = h(x)\)，该多项式唯一。

Hermite插值余项

若我们的函数\(f(x)\)在插值区间上2n+2阶导数，那么我们假设余项为\(R(x) = f(x)-H(x)\)，由Hermite插值定义可知，\(x_i\)是\(R(x)\)的二重根，因此\(R(x)\)含有\((x-x_i)\)，那么设

做出辅助函数

根据定义易知

因此辅助函数至少有n+1个二重零点x_i，以及零点x，根据罗尔定理，容易知道在每两个二重零点之间都有一个零点，同时\(\varphi^{'}(x_i) = 0\)，\(\varphi^{'}(t)\)至少有2
n+2个零点，以此类推，那么有

所以余项公式为：

Hermite的代码如下（更多代码见仓库）

Chebyshev插值

我们之前谈论的插值函数通常采样点是一个平均分布的采样点，比如说我们一个实际的例子，我们需要采集子弹发射速度，我们往往会使用等间隔时间发射一颗子弹。但是往往等间隔这件事情本身就是有难度的，也就是说我们选取的点对结果是有影响的，如下图，在端点的时候，误差可能会达到一个无穷大的地步。切比雪夫插值就是采用最优点间距选取的插值方法。

切比雪夫插值的出发点就是在插值区间上，控制插值误差的上下界，即

将这个误差控制在绝对值为1以内，实际上这是一个关于x的多项式，那么我们的问题就是，这个多项式是否存在极值，是否有特定的\(x_i\)使得最大值足够的小呢？那么我们试着列出公式进行推导：

我们再考虑一个问题，我们既然要构造多项式，并且控制误差在1以内，那么考虑一下三角函数，因为三角有两个优秀的性质：

• 三角函数的值域在\([-1,1]\)
• 根据傅里叶展开，我们知道曲线可以用三角函数进行拟合

因此如果区间在\([-1,1],\)我们选择\(x_i = cos\frac{(2i-1)\pi}{2n}\)（如果区间为[a,b]，那么选取\(x_i = \frac{b+a}{2} + \frac{b-a}{2}cos\frac{(2i-1)\pi}{2n}\)），这样\(x_i\)最小值就是\(1/2^{n-1}\)。我们可以看看切比雪夫根构造的插值采样点和普通的插值采样点构造的函数有何区别：

我们利用切比雪夫根（建议学习切比雪夫微分方程），构造的多项式就是切比雪夫多项式，在此之前，我们看一看三角函数的展开：

通过观察，那么我们假设\(g(cos\theta) = cosn\theta\)，\(g(x)\)是一个多项式，系数为\(2^{n-1}\)。实际上这也是成立的，证明过程留给读者。那么我们令切比雪夫多项式公式为：\(g_n(cos\theta) = cosn\theta\)
，于是：

可以发现

我们就得到了切比雪夫多项式的递推公式。切比雪夫插值的代码如下：

两种曲线的拟合与预测

我们之前往往要求插值曲线严格的经过各个插值采样点，这个方式往往具有不可操控的地方。例如我们的采样点往往也只是一个粗略的估计，如果我们强行使插值曲线通过这些点，那么数据本身的误差也会反映在数据之中，例如下图中，使用一条直线去拟合获得的误差理论上要比经过这些点的曲线更加优秀。因此我们需要使用某种方式进行限制，这里我们介绍三种拟合方式。

直线拟合法

直线拟合是最为简单的拟合方法，也就是构造一套直线\(y=ax+b\)，使得我们的采样点近似的接近这条直线，我们有三种方式进行衡量拟合的优劣

• 准则一：使得误差最大绝对值最小
• 准则二：使得总误差的绝对值最小
• 准则三：使得误差平方和最小

其中准则三就是我们介绍的最小二乘法，我们就以准则三构造拟合的直线，那么设采样点有\((x_i,y_i),i\in n\)，那么总误差为：

容易发现我们的误差Q是关于a,b的函数，因此这里我们可以使用二元函数求导法则求得参数值。

联立上述式子可以解出a与b：

采用直线拟合法，则需要保证采样点大致像直线的分布，或者说均匀分布在直线两侧，如果不符合这个则尽量不使用直线拟合。

指数拟合

对于某些增长型的曲线，我们可以采用指数函数\(y=ae^{bx}\)进行拟合，虽然是指数函数，我们可以将它化为直线方程，如下：

利用这个方程，我们可以利用直线拟合的手段，确定系数。不过检验我们的拟合精度的指标略微不同，指标为：

• 误差总平方和尽可能小
• 标准差\(\sigma = Q/(n+2)\)
• 相关系数\(R = 1- \frac{\sum(\tilde{y}-y_i)^2}{\sum(y_i-y)^2}\)，R尽可能接近1

习题

Reference

《数值分析》（原书第二版） —— Timothy Sauer

《数值计算方法》（清华大学出版社） —— 吕同富等

About Me

作　　者：WarrenRyan
出　　处：https://www.cnblogs.com/WarrenRyan/
本文对应视频：BiliBili
关于作者：热爱数学、热爱机器学习，喜欢弹钢琴的不知名小菜鸡。
版权声明：本文版权归作者所有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文链接。若需商用，则必须联系作者获得授权。
特此声明：所有评论和私信都会在第一时间回复。也欢迎园子的大大们指正错误，共同进步。或者直接私信我
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【BiBili——小陈的学习记录
Github——StevenEco
BiBili——记录学习的小陈（计算机考研纪实）
掘金——小陈的学习记录
知乎——小陈的学习记录

联系方式

电子邮件：cxtionch@live.com





社交媒体联系二维码：