函数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处的插值：fsyms t;
%计算输入的x y是否长度相等
if(length(x)==length(y))n=length(x);c(1:n)=0.0;
elsedis('x和y的维数不相等！');return;
endf = 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:iy1(j)=0;endfor(j=i+1:n)y1(j) = (y(j)-y(j-1))/(x(j)-x(j-i)); %利用前面的列计算后面列的值endc(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替换为x0elsef = collect(f); %将插值多项式展开f = vpa(f, 6);endend
end
% 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')
