利用差分的牛顿插值法(Newton)

差分牛顿插值法要求是等距的。

先来看三个概念

差分与均差的关系如下：

牛顿(Newton)插值的基本公式为：

由于差分插值是等距的，所以可以设x=x0+nh

对于上式

再由差分和均差的关系，可以将上面的黄色部分也就是牛顿插值基本公式转换成：

这里只介绍前插法：

此外还有一个余项公式：

余项部分暂不考虑。现在可以将公式转换为代码：

NewtonForward.m

function f = Newtonforward(x,y,x0)%求已知数据点的向前差分牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsyms t;if(length(x) == length(y))    n = length(x);    c(1:n) = 0.0;else    disp('x和y的维数不相等！');    return;endf = y(1);y1 = 0;xx =linspace(x(1),x(n),n); %吧x(1)到x(n)分成n-1份，比如1到100分成11份就是：1 10 20...100if(xx ~= x)    disp('节点之间不是等距的！');    return;endfor(i=1:n-1)    for(j=1:n-i)        y1(j) = y(j+1)-y(j);   %求△f    end    c(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的位置在哪里        else            f = collect(f);            f = vpa(f, 6);        end    endend

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;else    disp('x和y的维数不相等！！');    return;endf=y(n);y1=0;xx=linspace(x(1),x(n),n);  %将x(1)到x(n)分成x(2)-x(1)等间距段if(xx~=x)    disp('节点之间不是等距的！');    return;endfor(i=1:n-1)    for(j=i+1:n)        y1(j)=y(j)-y(j-1);    end    c(i)=y1(n);    l=t;    for(k=1:i-1)        l=l*(t+k);    end        f=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)));        else            f=collect(f);            f=vpa(f,6);        end    endend

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')