# Piecewise Hermite interpolation

In this post, we find the Hermite interpolation functions for the domain [-1, 1]. And then, we use it for a pieciwise interpolation. Notice that this interpolation has $C^1$ continuity compared to the $C^0$ continuity that is common in Lagrange interpolation.

To compute the polynomials explicitly we use sympy.

from __future__ import division
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt


We want to find a set of basis function that satisfy the following

\begin{align*} N_1(x_1) &= 1\\ N_1(x_2) &= 0\\ N_2(x_1) &= 1\\ N_2(x_2) &= 0\\ N'_3(x_1) &= 1\\ N'_3(x_2) &= 0\\ N'_4(x_1) &= 1\\ N'_4(x_2) &= 0 \end{align*}

These can be written as

\begin{equation*} \begin{bmatrix} 1 &x_1 &x_1^2 &x_1^3\\ 1 &x_2 &x_2^2 &x_2^3\\ 0 &1 &2x_1 &3x_1^2\\ 0 &1 &2x_2 &3x_2^2 \end{bmatrix} \begin{bmatrix} a_{11} &a_{12} &a_{13} &a_{14}\\ a_{21} &a_{22} &a_{23} &a_{24}\\ a_{31} &a_{32} &a_{33} &a_{34}\\ a_{41} &a_{42} &a_{43} &a_{44} \end{bmatrix} = \begin{bmatrix} 1 &0 &0 &0\\ 0 &1 &0 &0\\ 0 &0 &1 &0\\ 0 &0 &0 &1 \end{bmatrix} \end{equation*}

The formula for interpolation would be

\begin{equation*} f(x) \approx N_1(x) u_1 + N_2(x) u_2 + |J|(N_3(x) u'_3 + N_4(x) u'_4)\quad \forall x\in [a, b] \end{equation*}

with $|J| = (b - a)/2$ the Jacobian determinant of the transformation.

x1, x2, x = sym.symbols("x1 x2 x")
V = sym.Matrix([
[1, x1, x1**2, x1**3],
[1, x2, x2**2, x2**3],
[0, 1, 2*x1, 3*x1**2],
[0, 1, 2*x2, 3*x2**2]])
V

\begin{equation*} \left[\begin{matrix}1 & x_{1} & x_{1}^{2} & x_{1}^{3}\\ 1 & x_{2} & x_{2}^{2} & x_{2}^{3}\\ 0 & 1 & 2 x_{1} & 3 x_{1}^{2}\\ 0 & 1 & 2 x_{2} & 3 x_{2}^{2}\end{matrix}\right] \end{equation*}

We can see that the previous matrix is a confluent Vandermonder matrix [1]. It is similar to a Vandermonde matrix for the first $n$ nodes and the derivatives of each row for the following ones.

The coefficients for our polynomials are given by the inverse of this matrix.

sym.simplify(V.inv())

\begin{equation*} \left[\begin{matrix}\frac{x_{2}^{2} \left(3 x_{1} - x_{2}\right)}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{x_{1}^{2} \left(x_{1} - 3 x_{2}\right)}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{x_{1} x_{2}^{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & - \frac{x_{1}^{2} x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\- \frac{6 x_{1} x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{6 x_{1} x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{x_{2} \left(2 x_{1} + x_{2}\right)}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & \frac{x_{1} \left(x_{1} + 2 x_{2}\right)}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\\frac{3 x_{1} + 3 x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{3 x_{1} + 3 x_{2}}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & - \frac{x_{1} + 2 x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & - \frac{2 x_{1} + x_{2}}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\\- \frac{2}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{2}{x_{1}^{3} - 3 x_{1}^{2} x_{2} + 3 x_{1} x_{2}^{2} - x_{2}^{3}} & \frac{1}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}} & \frac{1}{x_{1}^{2} - 2 x_{1} x_{2} + x_{2}^{2}}\end{matrix}\right] \end{equation*}

And we replace the values for the extreme nodes, -1 and 1.

V_inv = sym.simplify(V.subs({x1:-1, x2:1}).inv())
V_inv

\begin{equation*} \left[\begin{matrix}\frac{1}{2} & \frac{1}{2} & \frac{1}{4} & - \frac{1}{4}\\ - \frac{3}{4} & \frac{3}{4} & - \frac{1}{4} & - \frac{1}{4}\\ 0 & 0 & - \frac{1}{4} & \frac{1}{4}\\ \frac{1}{4} & - \frac{1}{4} & \frac{1}{4} & \frac{1}{4}\end{matrix}\right] \end{equation*}

The polynomials are the product of the coefficients and the monomials

sym.factor(V_inv.T * sym.Matrix([1, x, x**2, x**3]))

\begin{equation*} \left[\begin{matrix}\frac{1}{4} \left(x - 1\right)^{2} \left(x + 2\right)\\- \frac{1}{4} \left(x - 2\right) \left(x + 1\right)^{2}\\\frac{1}{4} \left(x - 1\right)^{2} \left(x + 1\right)\\\frac{1}{4} \left(x - 1\right) \left(x + 1\right)^{2}\end{matrix}\right] \end{equation*}

The interpolation basis would be

\begin{align*} N_1 (x) &= \frac{1}{4} (x - 1)^2 (2 + x)\\ N_2 (x) &= \frac{1}{4} (x + 1)^2 (2 - x)\\ N_3 (x) &= \frac{1}{4} (x - 1)^2 (x + 1)\\ N_4 (x) &= \frac{1}{4} (x + 1)^2 (x - 1)\, , \end{align*}

and the following function computes the interpolation for a given function and derivative

def hermite_interp(fun, grad, x0=-1, x1=1, npts=101):
jaco = (x1 - x0)/2
x = np.linspace(-1, 1, npts)
f1 = fun(x0)
f2 = fun(x1)
N1 = 1/4*(x - 1)**2 * (2 + x)
N2 = 1/4*(x + 1)**2 * (2 - x)
N3 = 1/4*(x - 1)**2 * (x + 1)
N4 = 1/4*(x + 1)**2 * (x - 1)
interp = N1*f1 + N2*f2 + jaco*(N3*g1 + N4*g2)
return interp


In this case, we interpolate the sinc function

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

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


The following snippet compute the interpolation and plot it.

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, x0=x0, x1=x1, npts=npts)
plt.plot(x, y, linestyle="dashed", color="#4daf4a")
plt.plot([x[0], x[-1]], [y[0], y[-1]], marker="o", color="#4daf4a",
linewidth=0)
plt.xlabel("x")
plt.ylabel("y")
plt.legend(["Exact function", "Interpolation"])
plt.show()