One-dimensional Lagrange Interpolation

The problem of interpolation or finding the value of a function at an arbitrary point $X$ inside a given domain, provided we have discrete known values of the function inside the same domain is at the heart of the finite element method. In this notebooke we use Lagrange interpolation where the approximation $\hat f(x)$ to the function $f(x)$ is built like:

\begin{equation} \hat f(x)={L^I}(x)f^I \end{equation}

In the expression above $L^I$ represents the $I$ Lagrange Polynomial of order $n-1$ and $f^1, f^2,,...,f^n$ are the $n$ known values of the function. Here we are using the summation convention over the repeated superscripts.

The $I$ Lagrange polynomial is given by the recursive expression:

\begin{equation} {L^I}(x)=\prod_{J=1, J \ne I}^{n}{\frac{{\left( {x - {x^J}} \right)}}{{\left( {{x^I} - {x^J}} \right)}}} \end{equation}

in the domain $x\in[-1.0,1.0]$.

We wish to interpolate the function $ f(x)=x^3+4x^2-10 $ assuming we know its value at points $x=-1.0$, $x=1.0$ and $x=0.0$.



In [1]:

    
from __future__ import division
import numpy as np
from scipy import interpolate
import sympy as sym
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D



In [2]:

    
%matplotlib notebook
sym.init_printing()

First we use a function to generate the Lagrage polynomial of order $order$ at point $i$



In [3]:

    
def basis_lagrange(x_data, var, cont):
    """Find the basis for the Lagrange interpolant"""
    prod = sym.prod((var - x_data[i])/(x_data[cont] - x_data[i])
                    for i in range(len(x_data)) if i != cont)
    return sym.simplify(prod)

we now define the function $ f(x)=x^3+4x^2-10 $:



In [4]:

    
fun = lambda x: x**3 + 4*x**2 - 10



In [6]:

    
x = sym.symbols('x')
x_data = np.array([-1, 1, 0])
f_data = fun(x_data)

And obtain the Lagrange polynomials using:



In [10]:

    
basis = []
for cont in range(len(x_data)):
    basis.append(basis_lagrange(x_data, x, cont))
    sym.pprint(basis[cont])









    



x⋅(x - 1)
─────────
    2    
x⋅(x + 1)
─────────
    2    
   2    
- x  + 1

which are shown in the following plots/



In [31]:

    
npts = 101
x_eval = np.linspace(-1, 1, npts)
basis_num = sym.lambdify((x), basis, "numpy")  # Create a lambda function for the polynomials



In [37]:

    
plt.figure(figsize=(6, 4))
for k in range(3):    
    y_eval = basis_num(x_eval)[k]
    plt.plot(x_eval, y_eval)



In [41]:

    
y_interp = sym.simplify(sum(f_data[k]*basis[k] for k in range(3)))
y_interp









    Out[41]:





$$4 x^{2} + x - 10$$

Now we plot the complete approximating polynomial, the actual function and the points where the function was known.



In [36]:

    
y_interp = sum(f_data[k]*basis_num(x_eval)[k] for k in range(3))
y_original = fun(x_eval)

plt.figure(figsize=(6, 4))
plt.plot(x_eval, y_original)
plt.plot(x_eval, y_interp)
plt.plot([-1, 1, 0], f_data, 'ko')
plt.show()

The next cell change the format of the Notebook.



In [ ]:

    
from IPython.core.display import HTML
def css_styling():
    styles = open('../styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()