【数值分析——插值法】牛顿插值多项式及Matlab实现

牛顿多项式插值法（Newton polynomial interpolation）

1. 算法思想及公式定义

   对于要插值的数据点，如果是只存在一对数据点，即 $\left(x_0,y_0\right)$ ，那么插值多项式显然为：
$$f_0\left(x\right)=y_0$$
其中， $a_0=y_0$

  如果现在我们新增一个数据点 $\left(x_1,y_1\right)$ ，那么，插值多项式的图像就是经过这两个点 $\left(x_0,y_0\right)、\left(x_1,y_1\right)$ 的直线，不难得到插值函数为：
$$f_1\left(x\right)=y_0+a_1\cdot \left(x-x_0\right)$$
带入 $f_1\left(x_1\right)=y_1$ 可得， $a_1=\frac{y_1-y_0}{x_1-x_0}$

  现在，我们在新增一个数据点 $\left(x_2,y_2\right)$ ，还是写成之前的 $f_1\left(x\right)$ 的形式，即：
$$f_2\left(x\right)=y_0+\frac{y_1-y_0}{x_1-x_0}\cdot \left(x-x_0\right)+a_2\cdot \left(x-x_0\right)\cdot \left(x-x_1\right)$$
带入 $f_2\left(x_2\right)=y_2$ 可得， $a_2=\frac{\frac{y_2-y_0}{x_2-x_0}-\frac{y_1-y_0}{x_1-x_0}}{x_2-x_1}$

  我们可以看到，每增加一个数据点后，我们只需要对之前的插值多项式做一个修正即可，即增加一项。如果我们可以求出 $a_i$ ，那么我们就能求出新增数据点之后的插值多项式。

  根据上面的分析，我们定义牛顿插值多项式为：
$$f_n\left(x\right)=a_0+a_1\cdot \left(x-x_0\right)+a_2\cdot \left(x-x_0\right)\cdot \left(x-x_1\right)+\cdots +a_n\cdot \left(x-x_0\right)\cdot \left(x-x_1\right)\cdots \left(x-x_{n-1}\right)$$
其中，$a_n\left(n=0,1,2\cdots n\right)$ 为待定系数。

2. 商差（均差）的概念及性质

  自变量之差与因变量之差叫做商差，函数 $y=f\left(x\right)$ 在区间 $\left[x_i,x_{i+1}\right]$ 的平均变化率：
$$f\left[x_i,x_{i+1}\right]=\frac{f\left(x_{i+1}\right)-f\left(x_i\right)}{x_{i+1}-x_i}$$
称为 $f\left(x\right)$ 关于 $x_i,x_{i+1}$ 的一阶差商，并记为 $f\left[x_i,x_{i+1}\right]$ 。

  同理可得，二阶商差为：
$$f\left[x_i,x_{i+1},x_{i+2}\right]=\frac{f\left[x_{i+1},x_{i+2}\right]-f\left[x_i,x_{i+1}\right]}{x_{i+2}-x_i}$$
$$⋮$$
  m阶商差为：
$$f\left[x_0,x_1,x_2,\cdots ,x_m\right]=\frac{f\left[x_1,x_2,\cdots ,x_m\right]-f\left[x_0,x_1,\cdots ,x_{m-1}\right]}{x_m-x_0}$$

商差表

$x_i$	$f\left(x_i\right)$	$f\left[x_i,x_{i+1}\right]$	$f\left[x_i,x_{i+1},x_{i+2}\right]$	$f\left[x_i,x_{i+1},x_{i+2},x_{i+3}\right]$	$\cdots$
$x_0$	$f\left(x_0\right)$	-	-	-	$\cdots$
$x_1$	$f\left(x_1\right)$	$f\left[x_0,x_1\right]$	-	-	$\cdots$
$x_2$	$f\left(x_2\right)$	$f\left[x_1,x_2\right]$	$f\left[x_0,x_1,x_2\right]$	-	$\cdots$
$x_3$	$f\left(x_3\right)$	$f\left[x_2,x_3\right]$	$f\left[x_1,x_2,x_3\right]$	$f\left[x_0,x_1,x_2,x_3\right]$	$\cdots$
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$	$\ddots$

3. 牛顿插值多项式的余项

  定义：若在 $\left[a,b\right]$ 上用 $L_n\left(x\right)$ 近似 $f\left(x\right)$ ，则其截断误差为 $R_n\left(x\right)=f\left(x\right)-L_n\left(x\right)$ ，也称为插值多项式的余项。

  根据均差的定义：
$$f\left[x_0,x_1,x_2,\cdots ,x_m\right]=\frac{f\left[x_1,x_2,\cdots ,x_m\right]-f\left[x_0,x_1,\cdots ,x_{m-1}\right]}{x_m-x_0}$$
有
$$f\left[x,x_0,x_1,\cdots ,x_{m-1}\right]=f\left[x_0,x_1,x_2,\cdots ,x_m\right]+f\left[x,x_0,x_1,\cdots ,x_m\right]\cdot \left(x-x_m\right)$$
$$⋮$$
$$f\left[x,x_0\right]=f\left[x_0,x_1\right]+f\left[x,x_0,x_1\right]\cdot \left(x-x_1\right)$$
$$f\left(x\right)=f\left(x_0\right)+f\left[x,x_0\right]\cdot \left(x-x_0\right)$$
把前一式代入后一式中，就可以得到：
$$f\left(x\right)=f\left(x_0\right)+f\left[x,x_0\right]\cdot \left(x-x_0\right)+f\left[x,x_0,x_1\right]\cdot \left(x-x_0\right)\cdot \left(x-x_1\right)+\cdots +f\left[x_0,x_1,\cdots ,x_n\right]\cdot \left(x-x_0\right)\cdots \left(x-x_{n-1}\right)+f\left[x,x_0,x_1,\cdots ,x_n\right]\cdot w_{n+1}\left(x\right)$$
其中，
$$P_n\left[x\right]=f\left(x_0\right)+f\left[x,x_0\right]\cdot \left(x-x_0\right)+f\left[x,x_0,x_1\right]\cdot \left(x-x_0\right)\cdot \left(x-x_1\right)+\cdots +f\left[x_0,x_1,\cdots ,x_n\right]\cdot \left(x-x_0\right)\cdots \left(x-x_{n-1}\right)$$
$$w_{n+1}\left(x\right)=\left(x-x_0\right)\cdot \left(x-x_1\right)\cdots \left(x-x_n\right)$$
$P_n\left(x\right)$ 是近似 $f\left(x\right)$ 的插值多项式，前面也推导出了
$$f\left(x\right)=P_n\left(x\right)+R_n\left(x\right)$$
所以有
$$f\left(x\right)-P_n\left(x\right)=R_n\left(x\right)$$
那么插值余项为：
$$R_n\left(x\right)=f\left[x,x_0,x_1,\cdots ,x_n\right]\cdot w_{n+1}\left(x\right)$$

4. Matlab实现（只实现了插值函数，测试函数同前面的拉格朗日插值的测试）

function y=newton(x0,y0,x)
%   Newton插值函数
%   x0为已知数据点的x坐标
%   y0为已知数据点的y坐标
%   x为插值点的x坐标
%   y为各插值点函数值

% 获取数据点的个数
n=length(x0);
% 判断输入信息x0和y0的信息正确
if n~=length(y0)
    error('维数不等');
else
% 获取要预测插值点的个数
m=length(x);
% 初始化y值，并幅值为0，向量维度为m
y=zeros(1,m);

% 初始化A-->代表均差表，可参考均差表的形式来理解代码的实现
A=zeros(n,n);
A(:,1)=y0;
for j=2:n
    for i=j:n
        A(i,j)=(A(i,j-1)-A(i-1,j-1))/(x0(i)-x0(i-j+1));
    end
end

% 求整个插值函数
for i=1:n
    % 定义中间变量p
    p=1.0;

    for j=1:i-1
        p=p.*(x-x0(j));
    end
    y=y+A(i,i)*p;
end
end

pic.1 关于代码中A矩阵实现的由来

5. 结果分析

  根据原理分析，牛顿插值多项式是插值多项式 $P\left(x\right)$ 的另一种表示形式，与Lagrange多项式相比，它不仅克服了增加一个节点时整个计算工作重新开始的缺点, 且可以节省乘除法运算次数。

pic.2 Matlab对于sin函数的插值结果（未缩放）

pic.3 Matlab对于sin函数的插值结果（newton插值）

pic.4 Matlab对于sin函数的插值结果（lagrange插值）

  但是根据Matlab实现的计算结果来看，lagrange插值和newton插值的结果并不一样，这和之前推导出来的唯一性相违背，那么到底哪儿应该是正确的呢？为了寻找答案，我翻遍了网上的资料，以及看了我的Matlab代码实现，发现好像都没有什么问题，那么为什么会出现这种问题呢，难道真的是前面的多项式插值的唯一性存在错误吗？

  直到我一步步的将前面的实现过程在matlab的命令行敲出来，然后观察每个变量的具体数值的时候发现了如下图所示：

pic.5 具体值

  发现了什么，随着不断地计算，后面的项目全部趋近于0，直到全为0。看到这儿，我豁然开朗，这不就是前面提到的（数值分析第一章讲的）截断误差吗！就是这个截断误差，导致了两种插值方法的结果相差这么大！
  所以，虽然newton插值每增加一个节点，只需要修正插值多项值得一项，但是，在计算机计算的过程中，其截断误差也会随着项得不断增多而不断得增大。