Interpolation is a mathematical technique used to estimate values between two known data points. It is particularly useful in various scientific and engineering applications, such as function approximation, signal processing, and data analysis. One common method of interpolation is the Hermite polynomial interpolation, which allows us to approximate a function using both its function values and derivative values at specific points. In this blog, we will explore Hermite polynomial interpolation and implement it using the concept of divided differences in Python. The data used in this example corresponds to the values of the gamma function.

Understanding Hermite Polynomial Interpolation

Hermite polynomial interpolation is a powerful interpolation technique that provides an accurate approximation of a function by taking into account not only its function values at given points but also its derivative values. This results in a more flexible interpolation method that can accurately represent complex functions.

The divided differences method is a key concept in Hermite polynomial interpolation. It involves computing differences between function values at distinct points and using these differences to construct a polynomial. In the case of Hermite interpolation, we use the divided differences method to create a polynomial that matches both function values and derivative values.

Divided Differences Method

The divided differences method is a recursive algorithm used to compute coefficients for interpolating polynomials. It takes a set of points and their corresponding function values (and derivative values in the case of Hermite interpolation) and constructs a polynomial that passes through these points. The divided differences are computed in a step-by-step manner, and the final polynomial is built by combining these differences.

Let’s implement Hermite polynomial interpolation using the divided differences method in Python with NumPy and Matplotlib. We’ll break the code into blocks and explain each step.

Import Libraries

In this first code block, we import the necessary libraries, including NumPy for array operations and Matplotlib for data visualization. Additionally, we import scipy.special.gamma to compare our Hermite interpolation with the gamma function.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma

Compute Divided Differences

\(z\) \(f[z]\) First Second
\(z_0 = x_0\) \(f [z_0] = f(x_0)\)    
\(z_1 = x_0\) \(f [z_1] = f(x_0)\) \(f[z_0, z_1] =f'(x_0)\)  
\(z_2 = x_1\) \(f [z_2] = f(x_1)\) \(f[z_1, z_2] =\frac{ f[z_2] - f[z_1]}{z_2 - z_1}\) \(f[z_0,z_1, z_2]= \frac{f[z_1, z_2] - f[z_0, z_1]}{z_2 - z_0}\)
\(z_3 = x_1\) \(f [z_3] = f(x_1)\) \(f[z_2, z_3] =f'(x_1)\) \(f[z_1, z_2, z_3] =\frac{f[z_2, z_3] - f[z_1, z_2]}{z_3 - z_1}\)
\(z_4 = x_2\) \(f [z_4] = f(x_2)\) \(f[z_3, z_4] = \frac{f[z_4] - f[z_3]}{z_4 - z_3}\) \(f[z_2, z_3, z_4] =\frac{f[z_3, z_4] - f[z_2, z_3]}{z_4 - z_2}\)
\(z_5 = x_2\) \(f [z_5] = f(x_2)\) \(f[z_4, z_5] = f'(x_2)\) \(f[z_3, z_4, z_5] =\frac{f[z_4, z_5] - f[z_3, z_4]}{z_5 - z_3}\)

The following code block defines the hermitedivdiff function, which calculates the divided differences required for Hermite interpolation. It takes three arrays as input: x_values (data points), y_values (function values), and y_prime_values (derivative values).

def hermitedivdiff(x_values, y_values, y_prime_values):
    m = len x_values
    l = 2 * m
    z = np.zeros(l)
    a = np.zeros(l)
    Q = []

    for i in range(m):
        z[2 * i:2 * i + 2] = x_values[i]
        a[2 * i:2 * i + 2] = y_values[i]
    Q.append(a.tolist())

    for i in np.flip(range(2, l, 2)):
        a[i] = (a[i] - a[i - 1]) / (z[i] - z[i - 1])
    for i in range(0, m):
        a[2 * i + 1] = y_prime_values[i]
    Q.append(a[1::].tolist())

    for j in range(2, l):
        for i in np.flip(range(j, l)):
            a[i] = (a[i] - a[i - 1]) / (z[i] - z[i - j])
        Q.append(a[j::].tolist())

    return a

Given the data:

xs = np.linspace(1, 5, 120)
ys = gamma(xs)
dx = xs[1] - xs[0]
y_prime_values = np.gradient(ys, dx)
sel=list(range(0,120,int(120/6)))
x_values = [xs[s] for s in sel]
y_values = [ys[s] for s in sel]
y_prime_values = [y_prime_values[s] for s in sel]
x_values y_values y_prime_values
1 1 -0.544959341
1.672268908 0.903676547 0.168339662
2.344537815 1.198951973 0.748269084
3.016806723 2.031372793 1.888625776
3.68907563 4.117888084 4.794227246
4.361344538 9.617465496 13.02618544

The divided difference table (\(Q\)) for the given data is as follows:

\(f[z]\) First Second Third 4th 5th 6th 7th 8th 9th 10th 11th
1                      
-0.54496 1                    
0.59750 -0.14328 0.90368                  
-0.19927 0.46354 0.16834 0.90368                
0.11468 -0.04507 0.40294 0.43922 1.19895              
-0.01365 0.09633 0.08444 0.45971 0.74827 1.19895            
0.00424 -0.00509 0.08605 0.20015 0.72881 1.23823 2.03137          
0.00447 0.01327 0.02166 0.11518 0.35501 0.96747 1.88863 2.03137        
-0.00208 -0.00112 0.01026 0.04235 0.20059 0.62471 1.80741 3.10369 4.11789      
0.00175 0.00262 0.00594 0.02224 0.08720 0.31784 1.05205 2.51467 4.79423 4.11789    
-0.00086 -0.00115 -0.00123 0.00263 0.02931 0.14632 0.61294 1.87618 5.03726 8.18062 9.61747  
0.00063 0.00125 0.00305 0.00696 0.02134 0.07236 0.29226 1.00590 3.22864 7.20778 13.02619 9.61747

The Hermite Interpolation Polynomial

At the core of Hermite interpolation is the construction of a polynomial, denoted as H(x), that represents the approximation of a given function. The Hermite polynomial H(x) is expressed as:

\[H(x) = Q_{0,0} + Q_{1,1}(x - x_0) + Q_{2,2}(x - x_0)^2 + Q_{3,3}(x - x_0)^2(x - x_1) + Q_{4,4}(x - x_0)^2(x - x_1)^2 + \ldots + Q_{2n+1,2n+1}(x - x_0)^2(x - x_1)^2 \ldots (x - x_{n-1})^2(x - x_n).\]

In this equation, H(x) is the Hermite interpolation polynomial, and it is defined by the coefficients \(Q_{i,i}\), where i ranges from 0 to \(2n+1\), where \(n\) is the degree of the polynomial. Each \(Q_{i,i}\) represents an entry in the divided-difference table.

Hermite Polynomial Approximation

In this code block, we define the hermite_poly_approx function, which uses the divided differences computed in the previous step to approximate the Hermite polynomial at a given point x.

def hermite_poly_approx(x_values, y_values, y_prime_values, x):
    m = len(x_values)
    Q = hermitedivdiff(x_values, y_values, y_prime_values)
    z = np.zeros(2 * m)
    for i in range(m):
        z[2 * i:2 * i + 2] = x_values[i]
    Hx = Q[0]
    pr = 1
    for j in range(2 * m - 1):
        pr *= x - z[j]
        Hx += Q[j + 1] * pr
    return Hx

Input Data and Visualization

In this final code block, we provide sample input data (x_values, y_values, and y_prime_values) and generate the Hermite polynomial interpolation over a specified range using the hermite_poly_approx function. We also visualize the results using Matplotlib.

xaxis=np.linspace(x_values[0],x_values[-1])
interp = hermite_poly_approx(x_values, y_values, y_prime_values, xaxis)
plt.plot(xaxis, interp, label='Hermite interpolation',linewidth=8.0,color='pink')
plt.plot(xaxis, gamma(xaxis), label="Gamma",linestyle='--',color='red')
plt.plot(x_values, y_values, 'o', label='data')
plt.legend(loc='upper left');

You can use this code to perform Hermite polynomial interpolation and visualize the results for your specific data and application. Hermite interpolation, with its ability to account for both function values and derivative values, is a versatile tool for approximating complex functions accurately.