MATLAB数值分析：Newton（牛顿）插值法详解

function [y,A]=NewtonInterpolation(X, Y, x)
%% 牛顿插值计算
%输入的量:  X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量;
%                Y是纵坐标向量;
%                x是以向量形式输入的m个插值点;
%输出的量:  y可以为m个插值构成的向量.
%                A为查商表
%% 范例
% x=0:0.5:3;
% y=[0 0.4794 0.8415 0.9975 0.9093 0.5985 0.1411];
% x0=0.1:0.5:2.9;
% NewtonInterpolation(x, y, x0)
%% 检错
if(length(X) == length(Y))
    n = length(X);
    m = length(x);
    y=x;
else
    error('▲▼▲▼▲▼▲▼警报: X和Y的维数不相等▲▼▲▼▲▼▲▼▲▼');
end
%% 计算
A=zeros(n,n);       %为差商矩阵预分配内存
A(:,1)=Y';              %赋值差商表第一列
for t=1:m             %使用循环来计算多个插值点
    for j=2:n           %循环嵌套计算差商表
        for i=j:n
            A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
        end
    end
    s=0;
    for k=1:n           %求解多项式
        p=1.0;
        for j=1:k-1
            p=p*(x(t)-X(j));
        end
        s=s+A(k,k)*p;
    end
    y(t)=s;
end
end