实现牛顿插值的matlab代码

function [p,z]=myNewton(x,y,t);
%输入参数中x,y为元素个数相等的向量t为待估计的点可以为数字或向量。 
%输出参数中p为所求得的牛顿插值多项式z为利用多项式所得的t的函数值。
%x=[0.4 0.55 0.65 0.80 0.90 1.05]; 
%    y=[0.41075 0.57815 0.69675 0.88811 1.02652 1.25386]; 
%   t=0.4:0.1:1.05
p=[];%对向量p初始化为空
vectorLength=length(x);%计算向量x长度


%求n阶newton插值多项式
%计算每一阶的阶差
for i=1:vectorLength 
    %i=1，0阶均差为f(xi)
    if i==1
        p=[p,y(i)];
        
    %i~=1，计算f[x1...xi]
    else
        sum=0;%临时变量，存放累加结果
        %计算f(xj)/w'n+1(xj)(j~(1,i))
        for j=1:i
            multiple=1;%临时变量，存放累乘结果
            %计算(x(j)-x(1))...(x(j)-x(j-1))(x(j)-x(j+1))(x(j)-x(i))
            for k=1:i
                if j~=k
                multiple=multiple*(x(j)-x(k));
                end
            end
            sum=sum+y(j)/multiple;
        end
        p=[p,sum];%将计算结果放入p向量
    end
end


%根据阶差获得插值函数多项式


syms g mutiple f;%定义符号变量g mutiple f
f=0;


%对公式进行整合，形式为f=p(i)*(g-x(1))...(g-x(i-1))   i~(1,length(p))
for i=1:length(p)
    mutiple=1;
    for j=1:i-1
        mutiple=mutiple*(g-x(j));
    end
    f=f+mutiple*p(i);
end


f=expand(f);%将公式展开化简
disp('结果多项式为:');
p=sym2poly(f)%提取公式系数多项式并赋值给p






%求待估计点结果
z=[];%初始化待估计点结果向量为空
vectorLength=length(t);%计算待估计点数量
%求解待估计点结果向量
for i=1:vectorLength
    z=[z,polyval(p,t(i))];
end
disp('待估计值结果为');
t


%绘图，区间为给定x值的区间
ezplot(f,[x(1),x(end)]);%绘制插值函数图形
hold on,
plot(x,y,'r*','markersize',10);%绘制给定值图形
hold on,
plot(t,z,'ro','markersize',5);%绘制待估计值图形
hold on,
title('NewTon插值结果图');
grid on
legend('插值函数','给定值','待估计值');