利用均差的牛顿插值法(Newton)

函数f的零阶均差定义为 ，一阶定义均差为：

一般地，函数f 的k阶均差定义为：

或者上面这个式子求的k+1阶均差

利用均差的牛顿插值法多项式为：

简单计算的时候可以观看下面的差商（均差）表：

怎么利用差商表计算，可以看下面这个例子：

正常的话还有一个余项，在本文中先不考虑了。

总和上面的计算方法可以归纳出算法的大致思想：先计算差商表，类似于乘法口诀的思路，两个for循环就可以计算出，然后对于每一次内for循环以后，计算出了第一列，接着把相对应的f(x)计算出来，接着进入第二列的计算，接着计算相应的f(x).......一直到计算完毕最后一个f(x)，把所有的f(x)相加，便是最终的插值。

为了便于写算法，以五个样本点为例，计算了一下差商表：

【注】这里的覆盖是就是把这一列计算的值覆盖到前一列的对应地方，这样便于找到规律，可以发现分子始终是y(i+1)-y(i)，而分母则与列号有关系，详看代码，更易于理解

Newton1.m

function f = Newton1(x,y,x0)
%求已知数据点的均差形式牛顿插值多项式
%已知数据点的x坐标向量：x
%已知数据点的y坐标向量：y
%插值点的x坐标：x0
%求得的均差形式牛顿插值多项式或x0处的插值：f

syms t;
%计算输入的x y是否长度相等
if(length(x)==length(y))
    n=length(x);
    c(1:n)=0.0;
else
    dis('x和y的维数不相等！');
    return;  
end

f = y(1);%第0列的f(x)就是y(1)本身
y1 = 0; %这个y1不是y(1)，存的是差商表后面的值
l  = 1; %l是用来算f(x)后面对应的(x-x1)(x-x2).....的
 
for(i=1:n-1)   
    for j=1:i
        y1(j)=0;
    end
    for(j=i+1:n)
        y1(j) = (y(j)-y(j-1))/(x(j)-x(j-i)); %利用前面的列计算后面列的值
    end
    c(i) = y1(i+1);     
    l = l*(t-x(i));  
    f = f + c(i)*l;
    simplify(f);
    y = y1;
    
    if(i==n-1)
        if(nargin == 3)
            f = subs(f,'t',x0);%替换函数，用后面的替换前面的，把t替换为x0
        else
            f = collect(f);                %将插值多项式展开
            f = vpa(f, 6);
        end
    end
end

Newton1Insert.m

% x=[1 1.2 1.8 2.5 4];
% y=[1 1.44 3.24 6.25 16];
% f=Newton1(x,y)
% f=Newton1(x,y,2.0)

% x=[0.40 0.55 0.65 0.80 0.90];
% y=[0.41075 0.57815 0.69675 0.88811 1.02652];
% f=Newton1(x,y,0.596)


x1=0:2*pi;
y1=sin(x1);
xx=0:0.2:2*pi;
yy=Newton1(x1,y1,xx);
plot(x1,y1,'o:',xx,yy,'r')