Assume we have data (*x*_{i}, *y*_{i}), *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 *L*_{i}(*x*), where each *L*_{i}(*x*) is of degree *n*−1

*P*(*x*) = *y*_{1}*L*_{1}(*x*) + *y*_{2}*L*_{2}(*x*) + ... + *y*_{n}*L*_{n}(*x*)

The polynomials *L*_{i}(*x*) are defined as

*L*_{i}(*x*) = (*x*−*x*_{1}) ... (*x*−*x*_{i-1})(*x*−*x*_{i+1}) ... (*x*−*x*_{n}) / α_{i}

where

α_{i} = (*x*_{i}−*x*_{1}) ... (*x*_{i}−*x*_{i-1})(*x*_{i}−*x*_{i+1}) ... (*x*_{i}−*x*_{n})

It is obvious that *L*_{i}(*x*_{j}) = δ_{ij} and therefore *P*(*x*_{i}) = *y*_{i} meaning that the polynomial *P* passes through all the data (*x*_{i}, *y*_{i}), *i* = 1, ..., *n*.

If the the data (*x*_{i}, *y*_{i}) 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 *L*_{i}(*x*) for each *i*, namely all the products

(*x*−*x*_{1}) ... (*x*−*x*_{i-1})(*x*−*x*_{i+1}) ... (*x*−*x*_{n})

```
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 *x*^{2})

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')
```