This piece of code is a Matlab/GNU Octave function to perform Lagrange interpolation. Inputs are the data points, that is, an array xi which specifies the x coordinates, and another array yi which specifies the corresponding y coordinates. The function returns the array f which is actually the coefficient array of the Lagrange polynomial.

function f = lagrangepoly(xi, yi)

  n = length(xi);
  f = zeros(1, n);
  Ilog = logical(eye(n));
  for i = 1:n
    Pi = poly(xi(~Ilog(i,:)));
    Pi = Pi ./ polyval(Pi, xi(Ilog(i,:)));
    Pi = Pi .* yi(Ilog(i,:));
    f = f + Pi;
  end

In order to understand the code, let us trace out the execution of the code for a simple example. Say, we need to find the quadratic polynomial that passes throught points \( (1,2) \), \( (3,1) \), and \( (4,5) \). Therefore, xi = [1 3 4] and <em>yi = [2 1 5]</em>

The complete Lagrange interpolated polynomial is the sum of three quadratic polynomials each passing through one of the three points. Each iteration of the for loop computes each one of these polynomials Pi and adds it the overall polynomial <em>f</em>.

I iteration: polynomial passing through (1,2)

Pi = poly(xi(~Ilog(i,:)));

This creates the numerator for the first polynomial \( (x-3)(x-4) \)

Pi = Pi ./ polyval(Pi, xi(Ilog(i,:)));

This calculates the denominator of the first polynomial \( (1-3)(1-4) = (-2)*(-3) = 6 \) and divides the numerator by it. Now, we have \[ Pi = \frac{(x-3)(x-4)}{6} \]

Pi = Pi .* yi(Ilog(i,:));

Then we multiply by the corresponding value of y, in this case \( 2 \). Now, we have \[ Pi = 2 \frac{(x-3)(x-4)}{6} \]

f = f + Pi;

Finally, we add the first polynomial to the overall polynomial f.

II iteration: polynomial passing through (3,1)

Pi = poly(xi(~Ilog(i,:)));

This creates the numerator for the second polynomial \( (x-1)(x-4) \)

Pi = Pi ./ polyval(Pi, xi(Ilog(i,:)));

This calculates the denominator of the first polynomial \( (3-1)*(3-4) = 2*(-1) = -2 \) and divides the numerator by it. Now, we have \[ Pi = \frac{(x-3)(x-4)}{-2} \]

Pi = Pi .* yi(Ilog(i,:));

Then we multiply by the corresponding value of y, in this case \( 1 \). Now, we have \[ Pi = \frac{(x-3)(x-4)}{-2} \]

f = f + Pi;

Finally, we add the second polynomial to the overall polynomial f.

III iteration: polynomial passing through (4,5)

Pi = poly(xi(~Ilog(i,:)));

This creates the numerator for the third polynomial \( (x-1)(x-3) \)

Pi = Pi ./ polyval(Pi, xi(Ilog(i,:)));

This calculates the denominator of the first polynomial \( (4-1)*(4-3) = 3*1 = 3 \) and divides the numerator by it. Now, we have \[ Pi = \frac{(x-1)(x-3)}{3} \]

Pi = Pi .* yi(Ilog(i,:));

Then we multiply by the corresponding value of y, in this case \( 5 \). Now, we have \[ Pi = 5 \frac{(x-3)(x-4)}{3} \]

f = f + Pi;

Finally, we add the third polynomial to the overall polynomial f.

After the completion of the three iterations, the polynomial f is the required Lagrange interpolated polynomial passing through the three given points.

For a detailed explanation of Lagrange interpolation, please look at Lagrange Interpolating Polynomial in Wolfram Mathworld.