首页 文章

牛顿插值

提问于
浏览 次
0

嗨，我有这个函数来计算牛顿多项式的系数列表：

function p = polynom(x,y,c)
m = length(x);
p = c(m)*ones(size(y));
for k = m-1:-1:1
 p = p.*(y-x(k)) + c(k);
end

我已经有另一个程序可以正确找到划分的差异 . 对于 x=[3 1 5 6] ， y=[1 -3 2 4] 我得到 c=[1.0000 2.0000 -0.3750 0.1750] 这是正确的 .

但是当我使用上面的函数时，它会给出结果：

p =

   -3.0000  -53.6000   -0.1000    1.3500

但正确的答案应该是：

p = 

   0.1750   -1.9500    7.5250   -8.7500

我的功能出了什么问题？

matlab

1 回答

  • 0

    这是我用过的代码 . 我很久以前就用过了，所以我会发布整个代码 . 我希望这是你需要的 .

    table.m

    %test table
X = [1 1.05 1.10 1.15 1.2 1.25];
Y = [0.68269 0.70628 0.72867 0.74986 0.76986 0.78870];

    position.m

    %M-file position.m with function position(x)
%that for argument x returns 
%value 1 if x<X(2),
%value 2 if x>X(n-1)
%else 0 

function position=position(x)
table;
n=length(X);
if x<X(2)
    position=1;
else
if x>X(n-1)
    position=2;
else position=0;
end
end

    Newton.m

    function Newton=Newton(x)
table;
%position(x);
n=length(X);
A=zeros(n,n+1);
% table of final differences - upper triangle
for i=1:n
    A(i,1)=X(i);
    A(i,2)=Y(i);
end
for i=3:n+1
    for j=1:(n-i+2)
        A(j,i)=A(j+1,i-1)-A(j,i-1);
    end
end
%showing matrix of final differences, A
A
h=X(2)-X(1);
q1=(x-X(1))/h;
q2=(x-X(n))/h;
s1=q1;
s2=q2;
%making I Newton polynomial
if (position(x)==1)
    p1=A(1,2);
    f=1;
    for i=1:n-1
        f=f*i;
        p1=p1+(A(1,i+2)*s1)/f;
        s1=s1*(q1-i);
    end
    Newton=p1;
else
%making II Newton polynomial
if (position(x)==2)
    p2=A(n,2);
    f1=1;
    for i=1:n-1
        f1=f1*i;
        p2=p2+(A(n-i,i+2)*s2)/f1;
        s2=s2*(q2+i);
    end
    Newton=p2;
    %else, printing error
else Newton='There are more suitable methods than Newton(I,II)';

end
end
    回复于

相关问题