牛顿插值多项式.m

来自「拉格朗日插值多项式拟合,牛顿插值多项式,欧拉方程解偏微分方程,使用极限微分求解导」· M 代码 · 共 26 行

%_____________________________________________________
%函数调用方法：输入的是两个对应的向量x和y
%              输出的是插值多项式C，和差分数表D     
%____________________________________________________
function [C,D]=newtonpoly(X,Y)
 N=length(X);
 D=zeros(N,N);
 D(:,1)=Y';
 for j=2:N
     for k=j:N
         D(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));
        %相当于造差商表，其结构如下：
        % D(1,1)     0           0
        % D(2,1)  D(2,2)  ....   0
        %  ............ ......   0
        % D(N,1)  D(N,2)  ....D(N,N)
        %说明：只有下三角才有数值，每一列是有前一列得出的结果即D(2,2)=D(2,1)-D(1,1),由此得出上面的的推公式   
     end 
end
 C=D(N,N);
 %一下没搞明白
 for k=(N-1):-1:1
     C=conv(C,poly(X(k)));
     m=length(C);
     C(m)=C(m)+D(k,k);
 end