0
$\begingroup$

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\, .$$

$\endgroup$
4
  • $\begingroup$ I assume that you are asking for 1D interpolation. The example that you present is a cubic instead of a quintic polynomial, though. $\endgroup$
    – nicoguaro
    Commented Jan 29 at 12:25
  • $\begingroup$ Yes, you are a correct, it is a cubic Hermite polynomial, I just wanted to provide an example in what kind of coordinate system, I wanted to generate the quintic form. Do you have any recommendations on what formula or method to use to generate the quintic form on the macro mesh? $\endgroup$
    – Nomad
    Commented Jan 29 at 15:32
  • 2
    $\begingroup$ Do you want your final approximation to be $C^1$ or $C^2$? For $C^2$, googling yields, for example, rose-hulman.edu/~finn/CCLI/Notes/day09.pdf . Have a look at the functions $H_i^5$, which are defined on $[0,1]$. Shift them to $[-1,1]$ and you're done. $\endgroup$
    – cos_theta
    Commented Jan 30 at 8:42
  • 1
    $\begingroup$ I was just about to post the same reply. Section 9.2 on page 9-3 shows the quintic polynomials. $\endgroup$
    – IPribec
    Commented Jan 30 at 8:44

1 Answer 1

5
$\begingroup$

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*(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()

with the following result

enter image description here

$\endgroup$
1
  • $\begingroup$ Thank you very much for such elaborated answer! $\endgroup$
    – Nomad
    Commented Feb 3 at 18:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.