Weird behavior of Matlab poly() function

Question

I'm currently dealing with matlab code about finding Lagrange interpolation of given datapoints.

Here's my code:

X=1998:1:2008;
Y=[21300 23057 24441 25917 27204 28564 29847 31200 32994 34800 36030];
plot(X,Y,'-o'); %intuitive interpolation
hold on

I=zeros(1,length(X));

for k=1:1:length(X)
    r=X;
    r(k)=[];
    p=poly(r);
    m=polyval(p,X(k));
    m=1/m;
    p=p*m;
    p=p*Y(k);
    I=I+p;
end
t=1997:0.001:2010;
funcval=polyval(I,t);
plot(t,funcval)
hold on

As you can see, points of (X,Y) are designated datapoints and I is the polynomial coefficient vector of interpolating polynomial. This code works well until there are 5 points when

X=1998:1:2002
Y=[21300 23057 24441 25917 27204]

But when more datapoints are added, p=poly(r) suddenly produces weird output which leads us to polynomial with complex roots. For example, when there are 6 datapoints

X=1998:1:2003;
Y=[21300 23057 24441 25917 27204 28564];
plot(X,Y,'-o'); %intuitive interpolation
hold on

I=zeros(1,length(X));

for k=1:1:length(X)
    r=X;
    r(k)=[];
    p=poly(r);
    m=polyval(p,X(k));
    m=1/m;
    p=p*m;
    p=p*Y(k);
    I=I+p;
end
n=roots(I);

n must be root vector, and it is

2006.55014379998 + 0.00000000000000i
2002.53696731502 + 5.24636347412357i
2002.53696731502 - 5.24636347412357i
1995.63481881447 + 3.41460653963521i
1995.63481881447 - 3.41460653963521i

which is surely not a desired result of

[1998 1999 2000 2001 2002 2003]

I suspect this might be related with overflow error, but not sure.
Has anybody confronted with similar case? Any advice would be greatly helpful.

Thanks.

Weird behavior of Matlab poly() function

0 Answers