差分牛顿插值法要求是等距的。
先来看三个概念
差分与均差的关系如下：
牛顿(Newton)插值的基本公式为：
由于差分插值是等距的，所以可以设x=x0+nh
对于上式
再由差分和均差的关系，可以将上面的黄色部分也就是牛顿插值基本公式转换成：
这里只介绍前插法：
此外还有一个余项公式：
余项部分暂不考虑。现在可以将公式转换为代码：
NewtonForward.m
function f = Newtonforward(x,y,x0)
%求已知数据点的向前差分牛顿插值多项式
%已知数据点的x 坐标向量:x
%已知数据点的y 坐标向量:y
%为插值点的x坐标:x0
%求得的向前差分牛顿插值多项式或x0处的插值:f
syms t;
if(length(x) == length(y))n = length(x);c(1:n) = 0.0;
elsedisp('x和y的维数不相等！');return;
end
f = y(1);
y1 = 0;
xx =linspace(x(1),x(n),n); %吧x(1)到x(n)分成n-1份，比如1到100分成11份就是：1 10 20...100
if(xx ~= x)disp('节点之间不是等距的！');return;
end
for(i=1:n-1)for(j=1:n-i)y1(j) = y(j+1)-y(j); %求△fendc(i) = y1(1); l = t;for(k=1:i-1)l = l*(t-k); %每一项的∏部分end;f = f + c(i)*l/factorial(i); %得到前i项的牛顿插值的和simplify(f);y = y1;if(i==n-1)if(nargin == 3)f = subs(f,'t',(x0-x(1))/(x(2)-x(1))); %如果有输入x0，那么就把未知数t替换成(x0-x(1))/(x(2)-x(1))，表示x0的位置在哪里elsef = collect(f);f = vpa(f, 6);endend
end
NewtonForwardInsert.m
xx=0:2*pi;
yy=sin(xx);
x1=0:0.5:2*pi;
y1=Newtonforward(xx,yy,x1);
plot(xx,yy,'o',x1,y1,'r')
【附】前插和后插的公式对比
前插公式：
余项：
后插公式：
余项：
后插法的matlab实现：
Newtonback.m
function f = Newtonback(x,y,x0)
%求已知数据点的向后差分牛顿插值多项式
%x:已知数据点的x坐标向量
%y:已知数据点的y坐标向量
%x0:插值点的x坐标
%f:求得的向后差分牛顿插值多项式或x0处的插值
syms t;
if(length(x)==length(y))n=length(x);c(1:n)=0.0;
elsedisp('x和y的维数不相等！！');return;
end
f=y(n);
y1=0;
xx=linspace(x(1),x(n),n); %将x(1)到x(n)分成x(2)-x(1)等间距段
if(xx~=x)disp('节点之间不是等距的！');return;
end
for(i=1:n-1)for(j=i+1:n)y1(j)=y(j)-y(j-1);endc(i)=y1(n);l=t;for(k=1:i-1)l=l*(t+k);endf=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(n))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endend
end
NewtonBackInsert.m
xx=0:2*pi;
yy=sin(xx);
x1=0:0.5:2*pi;
y1=Newtonback(xx,yy,x1);
plot(xx,yy,'o',x1,y1,'r')
