本篇为Newton插值法，构造插值多项式
拉格朗日（Lagrange）插值法链接如下：
链接: Lagrange插值法的matlab实现.

算法

关于Newton插值多项式的具体内容和算法见《数值计算方法》—丁丽娟，P130-134

该算法核心是差商表的建立

上图来自《数值计算方法》—丁丽娟，P133-134

程序

Matlab代码如下：

//  输入量
xi: 离散样点的横坐标值
yi: 离散样点的纵坐标值
x: 插值多项式中自变量符号

// 输出量：
y: Newton插值多项式

// An highlighted block
function p= Newton_fun(x,xi,yi)
n=length(xi);
f=zeros(n,n);

% 对差商表第一列赋值
for k=1:n      
    f(k)=yi(k);
end
% 求差商表
for i=2:n       % 差商表从0阶开始；但是矩阵是从1维开始存储！！！！！！
    for k=i:n
        f(k,i)=(f(k,i-1)-f(k-1,i-1))/(xi(k)-xi(k+1-i));  
    end
end
disp('差商表如下：');
disp(f);

%求插值多项式
p=0;          
for k=2:n
    t=1;
    for j=1:k-1
        t=t*(x-xi(j));
        disp(t)
    end
    p=f(k,k)*t+p;
    disp(p)
end
p=f(1,1)+p;

end

注意：

1. 书中算法存储数据下标从0开始，但实际要从1开始
2. 第三步对书上算法稍稍改动了一点儿，同阶差商计算顺序从低到高。

运行

运行示例为《数值计算方法》—丁丽娟 P-134 例三
例一：

// Command Window 输入
>> xi=[11 12 13];
>> yi=[2.3979 2.4849 2.5649];
>> x=sym('x');
>> p= Newton_fun(x,xi,yi)

运行结果如下：

实际上上式与 p = (87*(x - 11))/1000 - (7*(x - 11)*(x - 12))/2000 + 23979/10000 一模一样。
例二：

// Command Window 输入
>> xi=[10 11 12 13 14];
>> yi=[2.3026 2.3979 2.4849 2.5649 2.6391];
>> x=sym('x');
>> p= Newton_fun(x,xi,yi)
>> a=10:0.1:14;
>> p=subs(p,x,a);
>> plot(a,p)

运行结果如下：

插值多项式：
p = (953x)/10000 - (18689938453587(x - 10)(x - 11))/4503599627370496 + (7993589098605227(x - 10)(x - 11)(x - 12))/36893488147419103232 - (614891469122219*(x - 10)(x - 11)(x - 12)*(x - 13))/147573952589676412928 + 1687/1250

对已知点计算

运行示例为《数值计算方法》—丁丽娟 P-134 例三

Command Window 输入
>> p= Newton_fun(11.5,xi,yi)

输出如下：

手算例题

计算插值多项式方法和程序算法一模一样。
此处多了一个误差估计：