函式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
 

% 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')