I am trying to use quintic Hermite basis functions for FEM applications, could someone please direct me to the general formula that would help me generate quintic Hermite shape functions?
In natural coordinates from -1 to 1, for instance in cubic Hermite functions it is
$$\phi = 0.25 x^3 - 0.75x + 0.5\, .$$
I mentioned how to do it for cubic polynomials in a previous answer. And has an expanded version in my blog.
You can do the derivation in global coordinates and obtain the following global basis functions.
\begin{align} H_1(x) &= \frac{(x - b)^3 (10 a^2 - 5 a b - 15 a x + b^2 + 3 b x + 6 x^2)}{(b - a)^5}\, ,\\ H_2(x) &= -\frac{(x - a)^3 (a^2 - 5 a b + 3 a x + 10 b^2 - 15 b x + 6 x^2)}{(b - a)^5}\, ,\\ H_3(x) &= -\frac{(x - a) (x - b)^3 (3 x - 4 a + b)}{(b - a)^4}\, ,\\ H_4(x) &= - \frac{(x - a)^3 (x - b) (3x - 4 b + a)}{(b - a)^4}\, ,\\ H_5(x) &= \frac{(x - a)^2 (x - b)^3}{2 (b - a)^3}\, ,\\ H_6(x) &= - \frac{(x - a)^3 (x - b)^2}{2 (b - a )^3}\, , \end{align}
where $x \in [a, b]$. The interpolation for a function $f(x)$, can be written as
$$p_5(x) = H_1(x) f(a) + H_2(x) f(b) + H_3(x) f'(a) + H_4(x) f'(b) + H_5(x) f''(a) + H_6(x) f''(b)\, . $$
Nevertheless, I would suggest using the local basis functions and expressing the global interpolation in terms of them. I think
$$p_5(s) = N_1(s) f(a) + N_2(s) f(b) + |J| (N_3(s) f'(a) + N_4(s) f'(b)) + |J|^2 (N_5(s) f''(a) + N_6(s) f''(b))\, , $$
with $|J| = (b - a)/2$, and
\begin{align} N_1(x) &= - \frac{(x - 1)^3 (3 x^2 + 9 x + 8)}{16}\, ,\\ N_2(x) &=\frac{(x + 1)^3 (3 x^2 - 9 x + 8)}{16}\, ,\\ N_3(x) &=- \frac{(x - 1)^3 (x + 1) (3 x + 5)}{16}\, ,\\ N_4(x) &=- \frac{(x - 1) (x + 1)^3 (3 x - 5)}{16}\, ,\\ N_5(x) &=- \frac{(x - 1)^{3} (x + 1)^2}{16}\, ,\\ N_6(x) &=\frac{(x - 1)^2 (x + 1)^3}{16}\, . \end{align}
I tested it, and it seems to be working. The following snippet shows it
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(3x2 + 9x + 8)/16
N2 = (x + 1)3(3x2 - 9x + 8)/16
N3 = -(x - 1)3(x + 1)(3x + 5)/16
N4 = -(x - 1)(x + 1)3(3x - 5)/16
N5 = -(x - 1)3*(x + 1)2/16
N6 = (x - 1)2*(x + 1)3/16
interp = N1f1 + N2f2 + jaco(N3g1 + N4g2) + jaco2(N5h1 + N6h2)
return interp
def fun(x):
return np.sin(2np.pix)/(2np.pix)
def grad(x):
return np.cos(2np.pix)/x - np.sin(2np.pix)/(2np.pix**2)
def hess(x):
return -2np.pinp.sin(2np.pix)/x - 2np.cos(2np.pix)/x2
+ np.sin(2np.pix)/(np.pix**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()
with the following result
