lagrangian polynomial function giving the incorrect answer

I have to write a function that gives the lagrangian interpolating polynomial estimation for a point, given a set of points.

Here is what I wrote:

function polyatc = lagrange(Xlist, Ylist, c)
    polyatc = 0;
    xprod=1;
    for i=1:size(Ylist)
        for j= setdiff(1:size(Xlist), i)
            xprod=xprod*((c-Xlist(j))/(Xlist(i)-Xlist(j)));
        end
        polyatc = polyatc+Ylist(i)*xprod;
    end
    disp(polyatc);
end

where polyatc should be the polynomial estimation for the function value at point x=c and xprod should be the section to the right of in the interpolating polynomial:

This all looks right to me and looks very similar to ones I have found online.

However, when I input: lagrange([1,5],[3,11],2) I get 3 when I should be getting 5 as the subsequent interpolating polynomial is 2x+1. I've tried following the logic of the code but I always come out with 5 and can't figure out where it is going wrong.

Answer 1

[John D'Errico]

A good try at some code, with three errors in it. Well, essentially two errors, but you made one of those errors twice so I'm not sure where the error count really comes out. ;-)

Your first mistake is in this line:

for i=1:size(Ylist)

Lets look at what happens.

Ylist
Ylist =
     3    11
>> size(Ylist)
ans =
     1     2

The size of Ylist is NOT the number of elements in Ylist. It is a vector, of length 2. That is because Ylist is a 1x2 vector. Now what happens when you use 1:size(Ylist)?

1:size(Ylist)
ans =
     1

It returns the number 1, so a bad idea, because that confuses the for loop. It goes only from 1 to 1, not 1 to the number of elements in Ylist as you want. Use numel instead. Use numel to learn how many elements are in an array of any size, although length will also tell you how many elements are in a vector. You actually made that mistake TWICE, in two for loops.

Next, you defined xprod as 1 BEFORE the outer loop. But the second time through the loop, it will use the last value from the previous iteration of the loop. Again, BAD. xprod needs to be initialized inside the main loop, so it will be properly initialized for the inner loop.

As a quick verification, what should the line be that passes through those points?

polyfit(Xlist,Ylist,1)
ans =
            2            1

So, y = 2*x + 1, as you suggested it should be. (Sometimes I'm just too lazy to think.) So the prediction at x==2 should indeed be 5.

Here is my modified version of your code. As you can see, quite close, but with those spots fixed.

function polyatc = lagrange(Xlist, Ylist, c)
  polyatc = 0;
  for i=1:numel(Ylist)
    xprod=1;
    for j= setdiff(1:numel(Xlist), i)
      xprod=xprod*((c-Xlist(j))/(Xlist(i)-Xlist(j)));
    end
    polyatc = polyatc+Ylist(i)*xprod;
  end
  disp(polyatc);
end

When I run THAT code, I get, for c=2, I see this:

5

In fact, as a nice test, if I do this:

syms c

and then run my version of your code, I get this:

2*c + 1

And that is what we should get. For any input c, the interpolant should produce 2*c+1.

Answer 2

[the cyclist]

I don't have a solution for you, but I can say that the line

xprod=xprod*((c-Xlist(j))/(Xlist(i)-Xlist(j)));

is never reached for your inputs.

I don't think you want

for i=1:size(Ylist)

because size(Ylist) is [1 2]. Maybe you meant

for i=1:length(Ylist)

?