Lagrange Polynomial Interpolation#

../_images/lag_interp.png

Creating interpolation polynomial from basis polynomials.#

Only at their own points are the basis polynomials equal to 1, at the other points they are zero. The final Lagrange polynomial interpolates through all the data points. Point (0, 1) lies on the P1 lagrange polynomial and the interpolation polynomial.

With the interpolation seen so far the common theme is that the function used in each method requires one more reference point than the degree of polynomial. This can be generalised for the higher degree polynomials. One curve is fitted through all the reference points.

The final polynomial can be built from so called Lagrange basis polynomials that have the property of having a value of 1 at its own data point and 0 at the other points.

\[P_i(x) = \prod_{j=1,j\neq i}^n \frac {x - x_j}{x_i - x_j} L(x) = \sum_{i=1}^n y_i P_i(x) \cdot\]

So for 3 points the polynomials are:-

\[\begin{split}P_1(x) &= \frac {(x - x_2) \cdot (x - x_3)}{(x_1 - x_3) \cdot (x_1 - x_3)} \\ P_2(x) &= \frac {(x - x_1) \cdot (x - x_3)}{(x_2 - x_1) \cdot (x_2 - x_3)} \\ P_1(x) &= \frac {(x - x_1) \cdot (x - x_2)}{(x_3 - x_1) \cdot (x_3 - x_2)}\end{split}\]
\[L = x_1 \cdot P_1 + x_2 \cdot P_2 + x_3 \cdot P_3\]

Show/Hide Code lagrange_poly0.py

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

sns.set_theme()

def quad_interp(x0, y0, x1, y1, x2, y2, x):
    y = y0 * (x - x1)*(x - x2)/((x0 - x1)*(x0 - x2)) + \
            y1 * (x - x0)*(x - x2)/((x1 - x0)*(x1 - x2)) + \
            y2 * (x - x0)*(x - x1)/((x2 - x0)*(x2 - x1))
    return y

x = [0, 1, 2]
y = [1, 3, 2]

xs = np.arange(-1.0,3.1,0.1)
ys0 = quad_interp(x[0],1,x[1],0,x[2],0, xs)
ys1 = quad_interp(x[0],0,x[1],1,x[2],0, xs)
ys2 = quad_interp(x[0],0,x[1],0,x[2],1, xs)

fig = plt.figure(figsize = (10,8))
plt.plot(xs, ys0, 'b', label = 'P1 $basis\,poly$')
plt.plot(xs, ys1, 'y', label = 'P2 $basis\,poly$')
plt.plot(xs, ys2, 'g', label = 'P3 $basis\,poly$')

L = ys0*y[0] + ys1*y[1] + ys2*y[2]
#print(L)

plt.plot(x, np.ones(len(x)), 'ko', x, np.zeros(len(x)), 'ko', label='basis points')
plt.plot(xs, L, 'r', x, y, 'ro', label='Interpolation')
plt.scatter(0, 1, fc='r', ec='k', s=60, label='P1 and Interpolation')

plt.legend()

plt.title('Lagrange Basis and Interpolation Polynomials')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')

plt.show()
#plt.savefig('../../figures/lag_interp.png')

One could just as easily multiply the basis polynomials before they are summed to make the interpolation polynomial. When interpolating over more than three reference points the Lagrange basis polynomial has to be enlarged, rather than expand the quadratic polynomial by hand use a script that manages the process automatically.

../_images/lag_interp_multi.png

Lagrange basis and interpolation polynomials for 4 points#

The basis polynomials go through its own point, at the other points they have a value of zero. The interpolation polynomial goes through all the data points.


Show/Hide Code lagrange_plot_multi.py

# https://stackoverflow.com/questions/4003794/lagrange-interpolation-in-python

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme()

class LagrangePoly:

    def __init__(self, X, Y):
        self.n = len(X)
        self.X = np.array(X)
        self.Y = np.array(Y)

    def basis(self, x, j):
        b = [(x - self.X[m]) / (self.X[j] - self.X[m])
             for m in range(self.n) if m != j]
        return np.prod(b, axis=0) * self.Y[j]

    def interpolate(self, x):
        b = [self.basis(x, j) for j in range(self.n)]
        return np.sum(b, axis=0)

X  = [-9, -4, -1, 7]
Y  = [5, 2, -2, 9]

lp = LagrangePoly(X, Y)
xx = np.arange(-100, 100) / 10

fig = plt.figure(figsize = (10,8))
plt.plot(xx, lp.basis(xx, 0), 'b', label = 'P1 $basis\,poly$')
plt.plot(xx, lp.basis(xx, 1), 'y', label = 'P2 $basis\,poly$')
plt.plot(xx, lp.basis(xx, 2), 'g', label = 'P3 $basis\,poly$')
plt.plot(xx, lp.basis(xx, 3), 'm', label = 'P4 $basis\,poly$')
plt.scatter(X, np.zeros(len(X)), fc='none', ec='k', label='basis points')
plt.plot(xx, lp.interpolate(xx), 'r', X, Y, 'ro', label='Interpolation')

plt.ylim(-5, 15)
plt.legend()
plt.title('Lagrange Basis and Interpolation Polynomials $4\,points$')
plt.xlabel('$x$')
plt.ylabel('$f(x)$')
plt.show()
#plt.savefig('../../figures/lag_interp_multi.png')

This method works best with only a few data points, otherwise the polynomial can fluctuate widely at the end points.

../_images/lag_interp_knot_points.png

Lagrange interpolation polynomial for 11 points#

The polynomial is well behaved because of the large data spread


../_images/lag_interp_knot_points2.png

Lagrange interpolation polynomial for 11 points#

The polynomial is ill behaved but still passes through all the data points

This Page