作者:LogM
本文原载于 https://blog.csdn.net/qq_28739605/article/details/72854371,不允许转载~
文章难免有错误之处,请在原文评论处指出~
牛顿插值(牛顿差商法)的Matlab实现。
原本是一道作业题,发上来给大家参考。不确定是否存在bug。
函数的代码:
function newton_divided_difference( f, min_x, max_x, n )
%UNTITLED 使用牛顿差商法插值
% f为原函数,min_x max_x定义了区间,n为多项式次数
%步骤1:求x,f,f[1],f[2]...f[n]
x = min_x : (max_x-min_x)/n : max_x;
for i = 1:1:n+1
y(i) = f(x(i));
end
a=[x',y'];
for i = 1:1:n
for j = 1:1:(n-i+1)
a(j,i+2) = (a(j+1,i+1)-a(j,i+1))/(a(j+i,1)-a(j,1));
end
end
%步骤2:画原始函数图像
figure(1);
STEP = 0.001;
x = min_x : STEP : max_x;
for i = 1:1:((max_x-min_x)/STEP+1)
y_1(i) = f(x(i));
end
plot(x,y_1,'r')
%步骤3:画插值函数图像
hold on
for i = 1:1:((max_x-min_x)/STEP+1)
y_2(i) = a(1,n+2);
for j = n:-1:1
y_2(i) = y_2(i) * (x(i) - a(j,1)) + a(1,j+1);
end
end
plot(x,y_2,'b')
hold off
%步骤4:画误差图像
figure(2);
y_error = y_1 - y_2;
plot(x,y_error)
end
调用:
结果:
原函数图像和插值得到的函数图像对比:
误差图像: