Newton插值法 | matlab实现

function y_p = Newton(x, y, x_p)
% x为插值节点的横坐标
% y为插值节点的纵坐标
% x_p 为需计算的插值函数值的横坐标
% y_p 为插值函数值

i = length(x) - 1;
iter = 1;
num = length(x_p);
y_i = y(1);
y_p = [];

%
for n = 1: num
    dif_quo = []; % 存放差商
    while i > 0
        if length(dif_quo) ~= 0
                y = dif_quo;
        end
        dif_quo = []; 
        for j = 1: i 
            dq = (y(j) - y(j+1)) / (x(j) - x(j+iter))
            dif_quo = [dif_quo, dq];
        end  
        dif_quo % 显示每一次的差商
        arg = 1;
        for h = 1: num-i    
            arg = arg * (x_p(n) - x(h)); % 求出(x-x0)...(x-x_n-1)累乘
        end        
        y_i = y_i + dif_quo(1) * arg;   
        i = i - 1;
        iter = iter + 1;
    end
    y_p = [y_p, y_i];
end