A crash course on the Finite Element Method

This document is meant as a quick start into the FEM.

The idea here is to give a motivation and picture of the overall working principle of the method itself.

The analysis of the method will require a more careful treatment (see e.g. NPDE 1). We concentrate on the conceptional aspects of the method.

Three examples (in NGSolve):

• Navier-Stokes equations (flow fields in fluids)
• Maxwell's equations (electromagnetism)
• Elasticity (mechanics, structure deformation)

What we saw (in short) are:

• comparibly simple models for complex physical phenomena
• the unknowns are fields (functions) defined on a (possibly complex) domain $\Omega$
• the equations are partial differential equations (PDEs)

Many practical problems are even more complex due to e.g.

• more complex geometries
• more complex physics
• non-linearities
• stochasticity
• coupling of several PDEs
• ...

We will however take a step back on consider rather simple problems:

• linear PDEs
• scalar unknowns
• 1D / 2D

But for these, we want to develop/understand solvers - more or less - from scratch.

Solving the Poisson equation

We search for the solution $u: \Omega \to \mathbb{R}$, so that $$ -\Delta u(x) = f(x) \quad \forall \, x \in \Omega $$

• $f: \Omega \to \mathbb{R}$ is a given function
• the domain $\Omega$ (open, bounded, typically Lipschitz boundary or smoother) is a subset of $\mathbb{R}^n$.
• The Poisson equation is a model for many physical phenomena:
  • f can be a heat source distribution, and u is the temperature
  • f can be a electrical charge distribution, and u is the electrostatic potential
• To select a unique solution $u$ we have to specify boundary conditions, for example homogeneous Dirichlet boundary conditions $$ u(x) = 0 \quad \forall \, x \in \partial \Omega. $$

Quick motivation of the Poisson equation

Let $V$ be a control volume in $\Omega$ and let $S$ be the surface of $V$. Then a simple model for the heat flow out of $V$ is given by the energy balance:

"The change of heat energy in $V$ is given by the heat flux flowing in through $S$ plus the energy transfered into $V$ due to a heat source $f: \Omega \to \mathbb{R}$."

If we further assume that the heat distribution is stationary, i.e. the heat energy in $V$ does not change in time, we obtain

"The change of heat energy in $V$ is given by the heat flux $\mathbf{q}: S \to \mathbb{R}^d$ flowing in through $S$ plus the energy transfered into $V$ due to a heat source $f: \Omega \to \mathbb{R}$ is zero."

$$ \int_{\partial V} -\mathbf{q} \cdot n \, ds + \int_V f \, dx = 0. \tag{B} $$

A classical model for the heat flux is given by Fourier's law: $$ \mathbf{q} = - \lambda \nabla u, $$ where $\lambda$ is the thermal conductivity and $u$ is the temperature (or $\lambda$ is the diffusivity and $u$ the concentration for other diffusion processes).

Then the divergence theorem allows us to transfer the surface integral into a volume integral: $$ \int_{\partial V} -\mathbf{q} \cdot n \, ds = \int_{\partial V} \lambda \nabla u \cdot n \, ds = \int_V \operatorname{div}( \lambda \nabla u ) \stackrel{\lambda = \text{const}}{=} \int_V \lambda \Delta u \, dx. $$

Plugging this into $(B)$ we obtain $$ \int_V \lambda \Delta u + f \, dx = 0. $$

As $V$ is arbitrary, we can (formally) localise the equation to obtain the PDE: Find $u: \Omega \to \mathbb{R}$ so that for all $x \in \Omega$ there holds $$ -\operatorname{div}(\lambda \nabla u) = - \lambda \Delta u= f. \tag{P} $$

On the domain boundary we further need to specify on each part of the boundary conditions, either

• the value of $u$, (Dirichlet boundary conditions)
• the flux $\mathbf{q}\cdot n$ (Neumann boundary conditions) or
• a more complex boundary condition

The simplest examples are homogeneous conditions, i.e. either $u=0$ or $\mathbf{q} \cdot n = 0$ and $\lambda = 1$.

Why/how to come up with discretization methods for PDEs?

• The unknown is a function $u: \Omega \to \mathbb{R}$, i.e. an infinite dimensional object
• $(P)$ is formulated for all $x \in \Omega$, i.e. it is an infinite number of equations
• To solve $(P)$ on a computer we need to discretize the problem, i.e. replace the infinite dimensional problem by a finite dimensional one
• Three options are:
  • Finite difference methods. Take only a finite number of points in $\Omega$ and replace the derivatives by finite difference quotients
  • Finite volume methods. Take only a finite number of (discrete) control volumes $V$ and replace the integrals by discrete (finite volume) approximations
  • Finite element method: Replace the problem by a formultaion w.r.t. test function, i.e. formulate an equation for all test functions $\phi$ in a space $V$. Afterwards replace $V$ by a finite dimensional subspace $V_h$.

Special case: 1D, $\lambda =1$

On the interval $\Omega = (a,b)$ the Poisson equation reads: Find $u:(a,b) \to \mathbb{R}$ so that $$ - \frac{d^2}{dx^2} u(x) = - u''(x) = f(x) \quad \forall \, x \in (a,b) $$ with boundary conditions:

• $u(a) = u(b) = 0$ (homogeneous Dirichlet boundary conditions) or
• $\frac{\partial u}{\partial n}(a) = \frac{\partial u}{\partial n}(b) = 0$ (homogeneous Neumann boundary conditions)

Weak formulation

• The formulation above in $(P)$ is called the strong form of the Poisson PDE.
• We will take a look at the weak form as the starting point for the finite element discretization method.

• We multiply the Poisson equation by a so called test function. It is an essentially arbitrary function (some restriction will be given later as needed): $$

• \Delta u(x) v(x) = f(x) v(x) \qquad \forall x \in \Omega $$

• We integrate over the domain $\Omega$:

$$ - \int_\Omega \Delta u(x) v(x) dx = \int_\Omega f(x) v(x) dx $$

• From Gauss' Theorem applied to the vector field $\nabla u v$ we obtain $$ \int_{\partial \Omega} \nabla u \cdot n \, v \,dx= \int_\Omega \operatorname{div} (\nabla u v) \,dx = \int_{\Omega} \Delta u v + \nabla u \nabla v \,dx $$

This allows to rewrite the left hand side such that (for suff. smooth fcts.) $$ \int_\Omega \nabla u \nabla v - \int_{\partial \Omega} \frac{\partial u}{\partial n} v = \int_\Omega f v $$

In the case of Dirichlet boundary conditions we allow only test-functions $v$ such that $v(x) = 0$ on the boundary $\partial \Omega$.

We have derived the weak form: find $u$ such that $u = 0$ on $\partial \Omega$ and

$$ \int_\Omega \nabla u \nabla v = \int_\Omega f v \tag{W} $$

holds true for all test-functions $v$ with $v = 0$ on $\partial \Omega$.

• Note that the weak formulation needs only first order derivatives of $u$ and $v$ (and also only in $L^2$), in contrast to the strong form which requires second order derivatives of $u$.

What is the right regularity to ask for $u$ and $v$?

(Without further details) the right (Sobolev) space is: $$ V = H^1(\Omega) := \{ u \in L_2(\Omega) : \nabla u \in L_2(\Omega) \}, $$ or, depending on the boundary conditions: $$ V = H_0^1(\Omega) = \{ u \in H^1(\Omega) : u_{|\partial \Omega} = 0 \}. $$

Let us consider the term on the left hand side of the variational formulation:

$$ A(u,v) := \int_{\Omega} \nabla u \nabla v \,dx $$

• Given $u$ and $v$ from the Sobolev space, we can compute a number $\int \nabla u \nabla v$
• $A(.,.)$ is a function mapping from two elements from $V$ into $\mathbb{R}$:

$$ A(.,.) : V \times V \rightarrow \mathbb{R} $$

The function $A(.,.)$ is linear in both arguments, and thus we call it a bilinear form.

Similarly, the right hand side

$$ f(v) := \int_{\Omega} f v \,dx $$

is a linear function

$$ f(\cdot) : V \rightarrow \mathbb{R}, $$

which we call a linear form.

Having these objects defined, our weak formulation reads now

$$ \text{find} \, u \in V : \quad A(u,v) = f(v) \quad \forall \, v \in V $$

This abstract formalism of Hilbert spaces, bilinear and linear forms apply for a large class of (elliptic, but not only) partial differential equations.

We now have a achieved an equation for all test functions $v \in V$.

The Finite Element Method

We cannot compute the solution in an infinite dimensional Hilbert space. But, we can define a finite dimensional sub-space (Galerkin method): $$ V_h \subset H^1_0 $$

and restrict the weak formulation to $V_h$:

$$ \text{find} \, u_h \in V_h : \quad A(u_h,v_h) = f(v_h) \quad \forall \, v_h \in V_h $$

The finite element solution $u_h$ is the Galerkin approximation to the true solution $u$.

One can analyze the $discretization error $| u - u_h |_{H^1}$ and in the case of the Poisson problem bound it by constant times the best approximation in $V_h$ (in the $H^1$ norm):

$$\| u - u_h \|_{H^1} \leq C \inf_{v_h \in V_h} \| u - v_h \|_{H^1} \leq c N^{-k/n}$$

For computing the solution $u_h$ we have to choose a basis for the function space $V_h$, where $N = \operatorname{dim} V_h$

$$ V_h = \operatorname{span} \{ \phi_1(x), \ldots \phi_N(x) \} $$

By means of this basis we can expand the solution $u_h$ as

$$ u_h(x) = \sum_{i=1}^N u_i \phi_i(x) $$

$u_i$ are the degrees of freedoms in the FE approximation.

The coefficients $u_i$ are combined to the coefficient vector $u = (u_1, \ldots, u_N) \in \mathbb{R}^N$

Instead of testing with all test-functions from $V_h$, by linearity of $A(\cdot,\cdot)$ and $f(\cdot)$, it is enough to test only with the basis functions $\phi_j(x), j = 1, \ldots, N$

Thus, the finite element probem can be rewritten as

$$ \text{find } u \in \mathbb{R}^N : \quad A(\sum_i u_i \phi_i(x), \phi_j(x)) = f(\phi_j(x)) \qquad \forall \, j = 1, \ldots N $$

By linearity of $A(\cdot,\cdot)$ in the first argument we can write

$$ \text{find } u \in \mathbb{R}^N : \quad \sum_{i=1}^N A(\phi_i, \phi_j) \, u_i = f(\phi_j) \qquad \forall \, j = 1, \ldots N $$

Since the basis functions are known, we can compute the matrix $A \in \mathbb{R}^{N\times N}$ with entries

$$ A_{j,i} = A(\phi_j,\phi_i) = \int_\Omega \nabla \phi_j(x) \cdot \nabla \phi_i(x) $$

and the vector $f \in \mathbb{R}^N$ as

$$ f_j = f(\phi_j) = \int_\Omega f(x) \phi_j(x) $$

Solving the finite element problem results in the linear system of equations for the coefficient vector $u = (u_1, \ldots, u_N)$:

$$ \text{find } u \in \mathbb{R}^N : \quad A u = f $$

By means of the coefficient vector, we have a representation of the finite element solution

$$ u_h(x) = \sum u_i \phi_i(x) $$

Questions that remain for later:

• How to setup the space $V_h$, respecitvely the basis functions $\phi_i$?
• How to then evaluate the integrals for $A_{j,i}$ and $f_j$?
• How to solve the linear system $A u = f$ (efficiently)?

$\leadsto$ quick FEM example

Task:

We are given the 1D PDE: Find $u: (0,1) \to \mathbb{R}$ such that \begin{align*} u(x) + \partial_x u(x) - \partial_x^2 u(x) = u(x) + u'(x) - u''(x) & = f(x) \quad \text{ in } (0,1) , \\ u(s) & = 0 \quad \text{ for } s \in \{0,1\}. \end{align*} for a given function $f: (0,1) \to \mathbb{R}$.

• Formulate the corresponding weak form, i.e. find $u \in H^1_0((0,1))$ such that $$ .. $$

• Extract the bilinear and linear forms $A(\cdot,\cdot)$ and $f(\cdot)$ from the weak form.