牛顿插值法

function y_h = newton_interpolation(x_vals,y_vals,n,x_h)
    %{
        牛顿插值法
        param x_vals：x轴数据点
        param y_vals：y轴数据点
        param n：插值阶数
        param x_h：x轴待插值点
        return y_h：对应x_h的插值
    %}
    len=min(length(x_vals),length(y_vals));
    if n>len-1
        error('当x和y的维度均为n时，牛顿插值的阶数不能超过n-1！') 
    end
%% 计算各阶差商
    quotient_diffs=cell(n,1);
    quotient_diffs{1}=diff(y_vals)./diff(x_vals);

    for ii=2:n
        temp=diff(quotient_diffs{ii-1});
        for jj=1:len-ii
            temp(jj)=temp(jj)/(x_vals(jj+ii)-x_vals(jj));
        end
        quotient_diffs{ii}=temp;
    end
%% 构造插值函数
    syms x real;
    N=y_vals(1);
    
    for ii=1:n
        temp=quotient_diffs{ii}(1);
        for jj=1:ii
            temp=temp*(x-x_vals(jj));
        end
        N=N+temp;
    end
%% 计算插值点
    y_h=double(subs(N,x,x_h));
end