# Solving PDEs using FEM using cubic Hermite polynomials

everyone. I am a beginner in Numerical mathematics, I have some idea of how to use Galerkin method to solve PDEs numerically, but so far I had no luck finding an example of how to solve a simple PDE using cubic Hermite polynomials. If you could help me out and provide an example of how to solve some poisson equation using cubic Hermite polynomials or direct me to some resources that could be helpful, that would be greatly appreciated. Thank you.

• It's a similar approach to when you use polynomial functions expressed wrt the monomial basis. E.g. if you solve a heat diffusion equation you still compute the stiffness matrix $W_{ij} = \int_{\Omega} \nabla \phi_i \cdot \nabla \phi_j$, the mass matrix $M_{ij} = \int_{\Omega} \phi_i \phi_j$, and a potential inhomogeneity $b_i = \int_{\Omega} \phi_i f$. The only thing that is different is that $\phi_i$ are expressed wrt the Hermite basis and support first derivatives as well as pointwise interpolation at nodes. This easily allows for enforcing $C^1$ continuity. Commented Nov 29, 2023 at 19:49
• My answer below is effectively what @lightxbulb said. Btw, Hermite shape functions are also discussed in 1) scicomp.stackexchange.com/questions/3620/…, and 2) scicomp.stackexchange.com/questions/36958/… Commented Nov 29, 2023 at 23:41

This method was studied in the context of boundary-value problems in,

Wheeler, M. F. (1973). An optimal $$L_\infty$$ error estimate for Galerkin approximations to solutions of two-point boundary value problems. SIAM Journal on Numerical Analysis, 10(5), 914-917.

It likely appeared already in Mary F. Wheeler's master's thesis, and also several of her other works.

The piecewise-cubic Galerkin FEM of Wheeler is mentioned in

Douglas, J., & Dupont, T. (1973). A finite element collocation method for quasilinear parabolic equations. Mathematics of computation, 27(121), 17-28. https://doi.org/10.2307/2005243

Personally, I encountered the method in the chemical engineering textbook: Nonlinear Analysis in Chemical Engineering (Finlayson, 1980), McGraw-Hill Inc.

Below is a complete example of the method, in Python, solving the linear problem, $$u''(x) - \mathit{Th}^2 u(x) = 0, \quad x \in [0,1]$$ with boundary conditions $$u'(0) = 0$$ and $$u(1) = 1$$. The constant $$\mathit{Th}$$ is known as the Thiele modulus. This is a classic chemical engineering model representing reaction and diffusion in a catalyst slab.

Within an element $$I_k = [x_{k-1},x_k]$$, the function is represented by the cubic function: $$u_k(x) = \sum_{i=1}^4 a_{k,i} H_i(\xi(x)), \quad x \in I_k$$ where $$H_i$$ are the polynomials: \begin{align*} H_1(\xi) &= (1 - \xi)^2(1 + 2\xi) \\ H_2(\xi) &= \xi (1 - \xi)^2 h_k \\ H_3(\xi) &= \xi^2 (3 - 2 \xi) \\ H_4(\xi) &= \xi^2 (\xi - 1) h_k \end{align*} specified in terms of local element coordinates, $$\xi = (x - x_{k-1})/h_k$$, and $$h_k = x_k - x_{k-1}$$. A graph of the polynomials is found below. The degrees of freedom $$a_{k,1}$$ and $$a_{k,3}$$ are the values at the breakpoints, whereas $$a_{k,2}$$ and $$a_{k,4}$$ are associated with the first derivatives, hence making this a $$C^1$$-continuous basis as already mentioned in the comments.

A few things to note:

• dense storage is used for simplicity, however banded or block-diagonal storage could be used instead
• in the solution array, odd elements represent the values $$u(x_k)$$, and even elements are the (scaled) first derivatives
• the same is true for the right-hand side vector, meaning that for a Dirichlet boundary condition at $$x = 1$$, the second-to-last element has an assigned value. For mixed-boundary conditions it may be necessary to adjust two array elements, and not just one as usual with Lagrange elements
• a smooth curve can be drawn with help of the interpolant function (not shown below)
• since the basis is $$C^1$$ continuous, material discontinuities have to be treated as internal boundary conditions (see this recent answer)
# hermite_fem.py

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

# Local mass matrix
mass_m = np.array([[13./35.,11./210,9./70.,-13./420],
[11./210.,1./105.,13./420.,-1./140.],
[9./70.,13./420,13./35.,-11./210],
[-13./420,-1./140,-11./210,1./105]])

# Local stiffness matrix
stiff_m = np.array([[6./5.,1./10.,-6./5.,1./10],
[1./10,2./15.,-1./10.,-1./30.],
[-6./5.,-1./10.,6./5.,-1./10],
[1./10,-1./30,-1./10,2./15]])

def _hermite_interp(U,x):
return U[0]*(2*x**3 - 3*x**2 + 1) + U[1]*(x**3 - 2*x**2 + x) + \
U[2]*(-2*x**3 + 3*x**2) + U[3]*(x**3 - x**2)

def hermite_fem(n,Th=1.0,plot=False):
"""Solves the equation  $$u'' - {Th}^2 u = 0$$ using a
Galerkin FEM with cubic Hermite spline shape functions

The boundary conditions are u' at x = 0, and u = 1 at x = 1.

Arguments:
n (int) ... number of breakpoints
Th (float) ... Thiele modulus
plot (boolean) ... plot the result or not
"""

nel = n-1
ndof = n*2

h = 1.0/nel

M = np.zeros((ndof,ndof))
rhs = np.zeros(ndof)

# Assembly loop
for k in range(nel):

kl, ku = 2*k, 2*k + 4
M[kl:ku,kl:ku] += (1.0/h)*stiff_m + (Th**2)*h*mass_m

# Dirichlet condition at x = 1
idr = ndof - 2
M[idr,:] = 0
M[idr,idr] = 1.0

rhs[idr] = 1.0

U = sp.linalg.solve(M,rhs)

if plot:

fig, ax = plt.subplots(nrows=1,ncols=1)

xx = np.linspace(0,1,61)
ax.plot(xx,np.cosh(Th*xx)/np.cosh(Th),'--',label="Analytic solution")

xu = np.linspace(0.,1.,n)
ax.plot(xu,U[0::2],'o',label="Hermite FEM")

ax.set_xlabel(r"$$x$$")
ax.set_ylabel(r"$$u(x)$$")
ax.set_title(r"Solution of $$u''(x) - \mathit{{Th}}\,u(x) = 0$$")
ax.legend()
plt.show()

if __name__ == '__main__':

hermite_fem(6,plot=True)


It is also possible to formulate an orthogonal collocation (OC) method using Hermite cubics. In this case, the residuals are forced to zero at the two knots corresponding to the shifted roots of a Gauss-Legendre polynomial. Continuity between elements is satisfied implicitly by the basis functions. Each element therefore contributes two equations to the global system. Finally, two equations are provided by the boundary conditions. A full description can be found in the book of Finlayson (1980).

The code is similar to the Galerkin FEM version:

sqrt3 = np.sqrt(3)

# Interpolation matrix
H_m = np.array([[(2*sqrt3/9+1./2.),(1./12. + sqrt3/36.),0.5-2*sqrt3/9,-(1./12.-sqrt3/36.)],
[(0.5-2*sqrt3/9),1./12.-sqrt3/36.,2*sqrt3/9+0.5,-(1./12.+sqrt3/36.)]])

# Second derivative matrix
B_m = np.array([[-2*sqrt3,-(1+sqrt3),2*sqrt3,(1-sqrt3)],
[2*sqrt3,-(1-sqrt3),-2*sqrt3,(1+sqrt3)]])

def hermite_ocm(n,Th=1.,plot=False):
"""Orthogonal collocation on finite elements with cubic Hermite splines"""

nel = n - 1
ndof = 2*n

he = 1./nel

M = np.zeros((ndof,ndof))
rhs = np.zeros(ndof)

# Symmetry condition at the cneter
M[0,:] = 0
M[0,1] = 1

# Collocation
# The residual is forced to zero at the collocation points
# which are chosen as the Gauss quadrature knots

for k in range(nel):

rl = 2*k+1
ru = rl + 2  # We set two rows at a time,

cl = 2*k
cu = 2*k + 4  # and four columns

M[rl:ru,cl:cu] += (1./he**2)*B_m - Th**2*H_m

# Dirichlet condition
idr = ndof - 2
M[ndof-1,:] = 0
M[ndof-1,ndof-2] = 1.0
rhs[ndof-1] = 1.0

U = sp.linalg.solve(M,rhs)

if plot:

fig, ax = plt.subplots(nrows=1,ncols=1)

xx = np.linspace(0,1,61)
ax.plot(xx,np.cosh(Th*xx)/np.cosh(Th),'--',label="Analytic solution")

xu = np.linspace(0.,1.,n)
ax.plot(xu,U[0::2],'o',label="Hermite Orthogonal Collocation")

ax.set_xlabel(r"$$x$$")
ax.set_ylabel(r"$$u(x)$$")
ax.set_title(r"Solution of $$u''(x) - \mathit{{Th}}\,u(x) = 0$$")
ax.legend()
plt.show()


A plot of the discrete $$L_\infty$$-norm at the breakpoints, i.e.

    UT = lambda x, Th: np.cosh(Th*x)/np.cosh(Th)
xu = np.linspace(0.,1.,n)
linf = np.linalg.norm(U[0::2] - UT(xu,Th),np.inf)


shows the accuracy is $$O(h^4)$$, as described in the original works.

Plotting the results for grids of size $$N = 5,10,15,...,2000$$, makes the random nature of the round-off errors more visible. An explanation can be found at the Mathemathics Stack Exchange.

A more extensive comparison of the two approaches, FEM and OC, also in terms of work performed, is given in the work

Dyksen, W. R., Houstis, E. N., Lynch, R. E., & Rice, J. R. (1984). The performance of the collocation and Galerkin methods with Hermite bi-cubics. SIAM journal on numerical analysis, 21(4), 695-715.

• I'm confused if the round-off errors are expected to grow as $1/h$ or $1/h^2$. Or is it the former for FEM, and the latter for orthogonal collocation? Commented Nov 29, 2023 at 22:59
• Thank you for the answer, may I ask how you computed the local stiffness matrix and does not seem like you are using the _hermite_interp() function at all, why did you define it? Commented Nov 30, 2023 at 3:32
• I was using it for plotting in an earlier version of the code. The local matrices where computed symbolically, for an element shifted to the local element coordinates $\xi \in [0,1]$. Commented Nov 30, 2023 at 8:28
• if it is not too much trouble for you, would it be alright with you to give me an example of how to compute those, stiffness matrix and load vector, I wanted to train using pen and paper and then use it solve beam equations? I understand that you integrate the derivative of Hi*derivative of Hj, I tried computing it using your H functions, but I used [-1,1] as integration bounds, maybe that was the error. Also the idea of Mass Matrix is still unclear, could you clarify it, please? Commented Nov 30, 2023 at 10:55
• Integrating the product of two cubics by hand is tedious, although doable in principle. In the stiffness matrix you will have the product of two quadratics. I used a symbolic algebra tool to do it. If you write the equations out for 2 elements, you will have 6 degrees of freedom in total. The local element matrices will be 4-by-4. In practice, if your equation has non-linear properties, you would just use quadrature rule. A 3-point Gauss-Legendre quadrature will integrate the stiffness matrix exactly, but it will still miss the true mass matrix by one degree. Commented Nov 30, 2023 at 11:39