牛顿插值法在Matlab上的实现

---

一、均差及其性质

1.均差的定义

在开始介绍牛顿插值多项式之前，需要先引入均差的定义

定义 :称$f[x_0,x_k]=\frac{f(x_k)-f(x_0)}{x_k-x_0}$为函数$f(x)$关于点$x_0$,$x_k$的一阶均差.

$f[x_0,x_1,x_k]=\frac{f[x_0,x_k]-f[x_0,x_1]}{x_k-x_1}$称为$f(x)$的二阶均差. 一般地，称$$f[x_0,x_1,\cdots ,x_k]=\frac{f[x_0,\cdots ,x_{k-2},x_k]-f[x_0,x_1,\cdots ,x_{k-1}]}{x_k-x_{k-1}}$$为$f(x)$的$\mathbf{k}$阶均差（均差也称为差商）。

2.均差的性质

（1）$k$阶均差可表示为函数值$f(x_0)$，$f(x_1)$，$\cdots$，$f(x_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)}$$可用归纳法证明此性质。这个性质也表明均差与节点的排列次序无关，称为均差的对称性，即$$f[x_0,x_1,\cdots ,x_k]=f[x_1,x_0,x_2,\cdots ,x_k]=\cdots =f[x_1,\cdots ,x_k,x_0]$$

（2）由性质（1）及$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}$$

（3）若$f(x)$在$[a,b]$上存在$n$阶导数，且节点$x_0$，$x_1$，$\cdots$，$x_n\in [a,b]$，则$n$阶均差与导数的关系为$$f[x_0,x_1,\cdots ,x_k]=\frac{f^{(n)}(\xi )}{n!},\xi \in [a,b]$$这个公式可直接用罗尔定理证明。

3.均差表

均差计算可列为均差表如下：

$x_k$	$f(x_k)$	一阶均差	二阶均差	三阶均差	四阶均差	$\cdots$
$x_0$	$\underline{f(x_0)}$
$x_1$	$f(x_1)$	$\underline{f[x_0,x_1]}$
$x_2$	$f(x_2)$	$f[x_1,x_2]$	$\underline{f[x_0,x_1,x_2]}$
$x_3$	$f(x_3)$	$f[x_2,x_3]$	$f[x_1,x_2,x_3]$	$\underline{f[x_0,x_1,x_2,x_3]}$
$x_4$	$f(x_4)$	$f[x_3,x_4]$	$f[x_2,x_3,x_4]$	$f[x_1,x_2,x_3,x_4]$	$\underline{f[x_0,x_1,x_2,x_3,x_4]}$
$⋮$	$⋮$	$⋮$	$⋮$	$⋮$	$⋮$	$⋮$

二、牛顿插值多项式

1.牛顿插值多项式

设待插函数为$f(x)$，利用拉格朗日插值多项式和均差的概念可以得到$$f(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\cdots +f[x_0,x_1,\cdots ,x_n](x-x_0)\cdots (x-x_{n-1})+f[x,x_0,\cdots ,x_n]{\omega }_{n+1}(x)=P_n(x)+R_n(x)$$
其中牛顿均差插值多项式:
$$P_n(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)+\cdots +f[x_0,x_1,\cdots ,x_n](x-x_0)\cdots (x-x_{n-1})$$
牛顿插值余项:
$$R_n(x)=f[x,x_0,\cdots ,x_n]{\omega }_{n+1}(x)$$

三、牛顿插值多项式在Matlab上的实现

1.Matlab代码

function [f0,f1] = Newtown_f(x0,y0)
%Newtown_f x0,y0分别是牛顿插值法所用点的x坐标、y坐标向量
%f0 输出多项式的句柄函数形式
%f1 输出多项式的符号表达式
syms x y
N=length(x0);   %统计样本点个数
mean_difference=zeros(N);
mean_difference(:,1)=y0';
mean_difference_temp=diff(y0);  %生成用于储存均差的矩阵,并存储一阶均差
%计算i-1阶均差
for i=2:N
    delta=i-1;
    for j=delta+1:N
        mean_difference(j,i)=mean_difference_temp(j-delta)/(x0(j)-x0(j-delta));
    end
   mean_difference_temp=diff(mean_difference(i:end,i));
end
%计算多项式
for i=0:(N-1)
    if i==0
        alpha=mean_difference(i+1,i+1);
        polynomial=1;
        f1=alpha*polynomial;
    else
        alpha=mean_difference(i+1,i+1);
        polynomial=polynomial*(x-x0(i));
        f1=f1+alpha*polynomial;
    end
end
f1=simplify(f1);    %简化多项式
f0=matlabFunction(f1);

2.代码使用演示

在多项式函数 $f(x)=x^2$ 上取三个样本点进行插值

%% 样本点生成
N=3; 
x0=linspace(-1,1,N);
y0=x0.^2;

%% 使用插值法进行计算
[f1,f2]=Newtown_f(x0,y0);

%% 画图
xo_o=linspace(-1,1,500);    %生成画图用点的x值
yo_o=f1(xo_o);  %代入多项式计算画图用点的y值
plot(xo_o,yo_o);
f2 %在命令行窗口展示多项式符号表达式

上图分别为插值多项式在$[-1,1]$区间上的图像和得到插值多项式$f2=x^2$，可以看到插值多项式图像与原多项式图像保持一致。

---

备注

一般情况下，插值多项式所插的点越多，多项式的次数就越高，插值结果越接近原函数，但对于龙格函数，例如：$f(x)=\frac{1}{1+25x^2}$，会存在龙格现象，即插值次数越高，插值结果越偏离原函数的特殊情况。