拉格朗日插值多项式，公式简洁，理论分析方便。但缺点是当改变节点或增加节点时，必须全部重新计算，这在实际计算中不方便。为此，介绍另一种形式的插值多项式 — 牛顿(Newton) 插值多项式。

一、差商及其性质

 设函数 $f(x)$ 在互异点 $x_0,x_1,\cdots$ 上的值为 $f(x_0),f(x_1),\cdots$ ，定义

$f(x)$ 在 $x_i,x_j$ 上的一阶差商为

$$f[x_i,x_j] = \frac{f(x_j)-f(x_i)}{x_j-x_i},$$

$f(x)$ 在$x_i,x_j,x_k$ 上的二阶差商为

$$f[x_i,x_j,x_k] = \frac{f[x_j,x_k]-f[x_i,x_j]}{x_k-x_i},$$

递推下去，$f(x)$ 在 $x_0,x_1,\cdots,x_k$ 上的 $k$ 阶差商为

$$f[x_0,x_1,\cdots,x_k] = \frac{f[x_1,x_2,\cdots,x_k]-f[x_0,x_1,\cdots,x_{k-1}]}{x_k-x_0},$$

注：规定 $f(x)$ 在 $x_i$ 的零阶差商为$f[x_i] = f(x_i)$。

差商性质 1、$k$ 阶差商可以表示成 $k+1$ 个节点上函数值的线性组合，即

$$f[x_0,x_1,\cdots,x_k] = \sum_{j=0}^{k}\dfrac{1}{\overline{\omega}_{k+1}^\prime} f(x_j);$$

2、 差商对其节点具有对称性，即

$$f[x_{i0},x_{i1},\cdots,x_{ik}] = f[x_0,x_1,\cdots,x_k],$$

其中 $x_{i0},x_{i1},\cdots,x_{ik}$ 是~$x_0,x_1,\cdots,x_k$ 的任一排列. 这一性质由性质1直接得到.

3、若 $f(x)$ 具有 $k$ 阶连续导数，则$$f[x_0,x_1,\cdots,x_k]=\frac{f^{(k)}(\xi)}{k!},$$

其中 $\xi$ 在 $k+1$ 个节点之间.

二、牛顿插值多项式及其余项

由差商定义，把 $x$ 看成~$[a,\,b]$ 内一点，可得

$$f(x) =f(x_0)+(x-x_0)f[x_0,x],$$

由性质(2)中节点对称性, $f[x_0,x_1,x] = f[x_1,x_0,x]=\dfrac{f[x_0,x]-f[x_1,x_0]}{x-x_1}=\dfrac{f[x_0,x]-f[x_0,x_1]}{x-x_1},$

$$f[x_0,x] = f[x_0,x_1]+(x-x_1) f[x_0,x_1,x],$$

$$f[x_0,x_1,x] =f[x_0,x_1,x_2]+(x-x_2) f[x_0,x_1,x_2,x],$$

$$\cdots \cdots$$

$$f[x_0,x_1,\cdots,x_{n-1},x]   =f[x_0,x_1,\cdots,x_{n-1},x_n]  +(x-x_n)f[x_0,x_1,\cdots,x_{n},x],$$

只要把后一式代入前一式可得

$$f(x) = f(x_0)+(x-x_0)f[x_0,x_1]+ (x-x_0)(x-x_1)f[x_0,x_1,x_2]+\cdots+(x-x_0)(x-x_1)\cdots(x-x_{n-1})f[x_0,x_1,\cdots,x_n]+ +(x-x_0)(x-x_1)\cdots(x-x_{n})f[x_0,x_1,\cdots,x_{n},x],$$

若记

$$N_n(x)  = f(x_0)+(x-x_0)f[x_0,x_1]+ (x-x_0)(x-x_1)f[x_0,x_1,x_2]+\cdots+(x-x_0)(x-x_1)\cdots(x-x_{n-1})f[x_0,x_1,\cdots,x_n],$$

$$R_n(x) =(x-x_0)(x-x_1)\cdots(x-x_{n})f[x_0,x_1,\cdots,x_{n},x],$$

则

$$f(x) = N_n(x)+R_n(x),$$

下面给出Newton基本插值的MATLAB程序。

function  yi = NewtonInterp(x,y,xi)

%  Newton 基本插值公式，调用格式为

%   yi = NewtonInterp(x,y,xi)

%  其中

%  x 为插值节点，y为节点处函数值，

%  xi 为为估计函数自变量，yi 为xi处函数估计值

%

%计算差商表

n = length(x);  Y= zeros(n); Y(:,1) = y';

for k = 1 : n-1

   for i = 1 : n - k

      if abs(x(i+k) - x(i))<eps

         error('% 输入的插值节点必须互异！');

      end

    Y(i,k+1) = (Y(i+1,k) - Y(i,k))/(x(i+k) - x(i));

   end

end

m = length(xi); yi = zeros(1,m);

for i = 1 : n   % 计算Newton插值公式

   z = ones(1,m);

   for k = 1 :  i - 1

       z = z .*(xi - x(k));

   end

   yi = yi + Y(1,i) * z;

end

测试代码

x = [1 4 9]; y = [1 2 3]; xi = [5 6];

%x = [-2 , 0,1 , 2]; y = [17,1,2,17];

yi = NewtonInterp(x,y,xi)

>> NewtonInterpmain

yi =

2.2667    2.5000