Lagrange polynomial interpolation

This article was kindly contributed by Vlad Gladkikh

Assume we have data (xiyi), i = 1, …, n.

The data don’t have to be equally spaced.

We can pass a Lagrange polynomial P(x) of degree n−1 through these data points.

function y = lagrange_interp(x,data)
    
    for i = 1:length(x)
        y(i) = P(x(i),data);
    end
    
endfunction

The polynomial P(x) is a linear combination of polynomials Li(x), where each Li(x) is of degree n−1

P(x) = y1L1(x) + y2L2(x) + … + ynLn(x)

The polynomials Li(x) are defined as

Li(x) = (xx1) … (xxi-1)(xxi+1) … (xxn) / αi

where

αi = (xix1) … (xixi-1)(xixi+1) … (xixn)

It is obvious that Li(xj) = δij and therefore P(xi) = yi meaning that the polynomial P passes through all the data (xiyi), i = 1, …, n.

If the the data (xiyi) are organized in an n×2 matrix, then the polynomial P can be computed as follows

function y = P(x,data)

    n = size(data,1);

    xi = data(:,1);
    yi = data(:,2);

    L = cprod_e(xi,x) ./ cprod_i(xi); 
    y = yi' * L;

endfunction

where the function cprod_e calculates the numerators of Li(x) for each i, namely all the products

(xx1) … (xxi-1)(xxi+1) … (xxn)

function y = cprod_e(x,a)
    
    n = length(x);
    y(1) = prod(a-x(2:$));
    for i = 2:n-1
        y(i) = prod(a-x(1:i-1))*prod(a-x(i+1:$));
    end
    y(n) = prod(a-x(1:$-1));
    
endfunction

The function cprod_i calculates all αi for i = 1, …, n.

function y=cprod_i(x)
    
    n = length(x);
    y(1) = prod(x(1)-x(2:$));
    for i = 2:n-1
        y(i) = prod(x(i)-x(1:i-1))*prod(x(i)-x(i+1:$));
    end
    y(n) = prod(x(n)-x(1:$-1));
    
endfunction

Now, let’s test our code.

I take two examples from the book “Fundamentals of Engineering Numerical Analysis” by Prof. Parviz Moin.

Both examples use data obtained from the Runge’s function

y = 1/(1+25 x2)

The data in the first example are equally spaced:

data1=[...
-1.0 0.038;
-0.8 0.058;
-0.6 0.1;
-0.4 0.2;
-0.2 0.5;
0.0 1.0;
0.2 0.5;
0.4 0.2;
0.6 0.1;
0.8 0.058;
1.0 0.038];

The commands to perform Lagrange interpolation are:

x = linspace(-1,1,1000);
y = lagrange_interp(x,data1);

Here’s the plot of the data (red circles), the Lagrange polynomial (solid line) and the Runge’s function (dashed line):

plot(x,y,'-b',...
x, 1.0 ./ (1+25.*x.^2),'--g',...
data1(:,1),data1(:,2),'or')

The second example takes non-equally spaced data:

data2 = [...
-1.0 0.038;
-0.95 0.042;
-0.81 0.058;
-0.59 0.104;
-0.31 0.295;
0.0 1.0;
0.31 0.295;
0.59 0.104;
0.81 0.058;
0.95 0.042;
1.0 0.038];

Interpolating and plotting them as we did in the previous example produces the following picture.

y = lagrange_interp(x,data2);

plot(x,y,'-b',...
x, 1.0 ./ (1+25.*x.^2),'--g',...
data2(:,1),data2(:,2),'or')