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

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

先来看三个概念



差分与均差的关系如下:



牛顿(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;
else
    disp('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);   %求△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
    end
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;
else
    disp('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);
    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
    end
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')


评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

风翼冰舟

额~~~CSDN还能打赏了

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值