本篇为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
注意:
- 书中算法存储数据下标从0开始,但实际要从1开始
- 第三步对书上算法稍稍改动了一点儿,同阶差商计算顺序从低到高。
运行
运行示例为《数值计算方法》—丁丽娟 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)
输出如下:
手算例题
计算插值多项式方法和程序算法一模一样。
此处多了一个误差估计: