hermite-and-lagrange.ipynb

Kernel: Unknown Kernel

Comparison of Lagrange and Hermite polynomials using divided differences

import numpy as np
%pylab inline
import matplotlib.pyplot as plt
# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation
from IPython.core.display import Image 

Populating the interactive namespace from numpy and matplotlib

Coefficients for Lagrange polynomials

#Construction of a kth divided difference (recursive code)
def D( i, k, Xn, Yn ):
    #If k+i>N
    if i+k>=len(Xn):
        return 0
    #Zeroth divided difference
    elif k == 0:
        return Yn[i]
    #If higher divided difference
    else:
        return (D(i+1, k-1, Xn, Yn)-D(i, k-1, Xn, Yn))/(Xn[i+k]-Xn[i])

Coefficients for Hermite polynomials

#Construction of a kth divided difference for Hermite polynomials (recursive code)
def Dh( j, k, Zn, Yn, Ypn ):
    #If k+j>N
    if j+k>=len(Zn):
        return 0
    #Zeroth divided difference
    elif k == 0:
        return Yn[j/2]
    #First order divided difference (even indexes)
    elif k == 1 and j%2 == 0:
        return Ypn[j/2]
    #If higher divided difference
    else:
        return (Dh(j+1, k-1, Zn, Yn, Ypn)-Dh(j, k-1, Zn, Yn, Ypn))/(Zn[j+k]-Zn[j])

Lagrange polynomials

#Lagrange Interpolating Function
def Lagrange( x, Xn, Yn ):
    n = len(Xn)
    
    #Detecting size of x
    try:
        Ninter = len(x)
    except:
        Ninter = 1
        x = np.array([x,])
    
    #First coeficient
    Lp = D(0, 0, Xn, Yn)*np.ones(Ninter)
    
    #High-order coeficients
    for l in xrange(Ninter):
        for k in xrange(1, n):
            Lp[l] += D(0, k, Xn, Yn)*np.prod(x[l]-Xn[:k])
        
    return Lp

Hermite polynomials

#Hermite Interpolating Function
def Hermite( x, Xn, Yn, Ypn ):
    n = len(Xn)
    #Auxiliar array
    Zn = np.zeros( 2*len(Xn) )
    #Even indexes
    Zn[::2] = Xn
    #Odd indexes
    Zn[1::2] = Xn
    
    #Detecting size of x
    try:
        Ninter = len(x)
    except:
        Ninter = 1
        x = np.array([x,])
    
    #First coeficient
    Hp = Dh(0, 0, Xn, Yn, Ypn)*np.ones(Ninter)
    
    #High-order coeficients
    for l in xrange(Ninter):
        for k in xrange(1, 2*n+2):
            Hp[l] += Dh(0, k, Zn, Yn, Ypn)*np.prod(x[l]-Zn[:k])
        
    return Hp

Comparison

#Generating arrays of sin(x)
def function(x):
    return np.sin(x)**2

#Generating arrays of derivative of sin(x)
def dfunction(x):
    return 2*np.cos(x)*np.sin(x)

#Number of intervals for data
Ndat = 4
#Limits
xmin = 0
xmax = 2*np.pi

Xn = np.linspace( xmin, xmax, Ndat )
Yn = function(Xn)
Ypn = dfunction(Xn)

#Obtaining Hermite and Lagrange interpolation
Ninter = 100
x = np.linspace( xmin, xmax, Ninter )
y = Hermite( x, Xn, Yn, Ypn )
yl = Lagrange( x, Xn, Yn )
f = function(x)

#Plotting
plt.figure( figsize=(12,6) )
plt.plot( x, y, color="green", linewidth=2, label="Hermite interpolation" )
plt.plot( x, yl, color="red", linewidth=2, label="Lagrange interpolation" )
plt.plot( x, f, color="blue", linewidth=2, label="real function" )
plt.plot( Xn, Yn, "o", color="red", label="data" )

#Formatting
plt.legend( loc = 'lower left' )
plt.grid()
plt.xlabel( "x" )
plt.xlabel( "y" )
plt.ylim( (-1,1.5) )
plt.title( "Linear interpolation of $\sin(x)$" )

<matplotlib.text.Text at 0x7fe51943aa50>