Hermite interpolation using quintic polynomials

hermite_quintic.ipynb

hermite_quintic.py

# -*- coding: utf-8 -*-
"""
Compute a (quintic) Hermite interpolation

@author: Nicolas Guarin-Zapata
@date: January 2024
"""
import numpy as np
import matplotlib.pyplot as plt


def hermite_interp(fun, grad, hess, x0=-1, x1=1, npts=101):
    jaco = (x1 - x0)/2
    x = np.linspace(-1, 1, npts)
    f1 = fun(x0)
    f2 = fun(x1)
    g1 = grad(x0)
    g2 = grad(x1)
    h1 = hess(x0)
    h2 = hess(x1)
    N1 = -(x - 1)**3*(3*x**2 + 9*x + 8)/16
    N2 = (x + 1)**3*(3*x**2 - 9*x + 8)/16
    N3 = -(x - 1)**3*(x + 1)*(3*x + 5)/16
    N4 = -(x - 1)*(x + 1)**3*(3*x - 5)/16
    N5 = -(x - 1)**3*(x + 1)**2/16
    N6 = (x - 1)**2*(x + 1)**3/16
    interp = N1*f1 + N2*f2 + jaco*(N3*g1 + N4*g2) + jaco**2*(N5*h1 + N6*h2)
    return interp


def fun(x):
    return np.sin(2*np.pi*x)/(2*np.pi*x)


def grad(x):
    return np.cos(2*np.pi*x)/x - np.sin(2*np.pi*x)/(2*np.pi*x**2)


def hess(x):
    return -2*np.pi*np.sin(2*np.pi*x)/x - 2*np.cos(2*np.pi*x)/x**2\
        + np.sin(2*np.pi*x)/(np.pi*x**3)

#%% Test
a = 2
b = 5
nels = 7
npts = 200
x = np.linspace(a, b, npts)
y = fun(x)
plt.plot(x, y, color="black")
xi = np.linspace(a, b, num=nels, endpoint=False)
dx = xi[1] - xi[0]


for x0 in xi:
    x1 = x0 + dx
    x = np.linspace(x0, x1, npts)
    y = hermite_interp(fun, grad, hess, x0=x0, x1=x1, npts=npts)
    plt.plot([x[0], x[-1]], [y[0], y[-1]], marker="o", color="#4daf4a",
         linewidth=0)
    plt.plot(x, y, linestyle="dashed", color="#4daf4a")
plt.xlabel("x")
plt.ylabel("y")
plt.legend(["Exact function", "Interpolation"])
plt.show()