Finite Element Methods

Finite element methods approximate the solution to a differential equation on a small element of the computational domain.

Distributions

Denote by $\mathrm{H}^s(a, b)$, for $s \geq 1$, the space of the functions $f \in C^{s-1}(a, b)$ such that $f^{(s-1)}$ is continuous and piecewise differentiable, so that $f^{(s)}$ exists unless for a finite number of points and belongs to $\mathrm{L}^2(a, b)$. The space $\mathrm{H}^s(a, b)$ is known as the Sobolev function space of order $s$ and is endowed with the norm $\left\Vert \cdot \right\Vert_{H^s(a,b)}$ defined as

$$ \begin{equation*} \left\Vert f \right\Vert_s = \left( \sum_{k=0}^s \left\Vert f^{\left(k\right)} \right\Vert_{\mathrm{L}^2\left(a,b\right)}^2 \right)^{1/2}. \end{equation*} $$

Let

$$ \begin{align*} C_0^\infty & = \lbrace \varphi \in C^\infty \, | \, \exists \, a, b \in (0,1) \quad \text{such that} \quad \varphi(x)=0 \quad \text{for} \quad 0 \leq x < a \quad \text{or} \quad b < x \leq 1 \rbrace . \end{align*} $$

Then for a function $v \in \mathrm{L}^2(0,1)$ we say $v^\prime$ is the weak derivative (or distributional derivative) if

$$ \begin{equation*} \int \limits_0^1 v^\prime \varphi \, \mathrm{d}x = - \int \limits_0^1 v \varphi^\prime \, \mathrm{d}x \quad \forall \, \varphi \in C^\infty_0 \left(0,1\right). \end{equation*} $$

Of interest is

$$ \begin{equation*} H^1 (0, 1) = \lbrace v \in \mathrm{L}^2(0, 1) \, : \, v^\prime \in \mathrm{L}^2(0, 1) \rbrace \end{equation*} $$

where $v^\prime$ is the distributional derivative of $v$, and

$$ \begin{equation*} H^1_0 (0, 1) = \lbrace v \in \mathrm{L}^2(0, 1) \, : \, v^\prime \in \mathrm{L}^2(0, 1), \, v(0) = v(1) = 0 \rbrace \end{equation*} $$

On $H^1$ there is the semi-norm:

$$ \left| v \right|_{\mathrm{H}^1(0,1)} = \left( \int_0^1 \left\| v^\prime \left( x \right) \right\|^2 \mathrm{d}x \right)^{1/2} = \left\Vert v^\prime \right\Vert_{\mathrm{L}^2(0,1)}. $$

To see that it is a semi-norm and not a norm, consider $v$ a constant, so ${v^{\prime}=0}$ thus ${\left| v \right|_{\mathrm{H}^1(0,1)} =0}$ for ${v \neq 0}$ and thus by definition is a semi-norm, rather than a norm. Now consider the integral on functions in $H_0^1$, it is the case that if the integral is zero so the function is constant, but as it must be zero on the boundaries, so the function is zero and hence a norm.

Galerkin Method

Consider the elementary problem:

$$ \begin{equation*} -\left( \alpha u^{\prime} \right)^{\prime} + \beta u^{\prime} + \gamma u = f\left( x \right) \quad \text{on} \quad (0,1) \quad \text{with} \quad u(0) = u(1) = 0 \end{equation*} $$

where $\alpha$, $\beta$, $\gamma \in C^0 \left( \left[ 0, 1 \right] \right)$ and $\alpha(x) \ge \alpha_0 >0$ for all $x \in \left[ 0, 1 \right]$.

Next, on $\mathrm{L}^2(0,1)$, define the scalar product

$$ \begin{equation*} \left( f, v \right) = \int \limits_0^1 f v \, \mathrm{d}x \end{equation*} $$

and a bilinear form $a : \left( \cdot, \cdot \right)$ which maps $H_0^1 \times H^1_0 \rightarrow \mathbb{R}$

$$ \begin{equation*} a\left( u, v \right) = \int_0^1 \left( \alpha u^\prime v^\prime + \beta u^\prime v + \gamma u v \right) \, \mathrm{d}x \end{equation*} $$

and consider the weak form

$$ \begin{equation*} \text{Find} \quad u \in H^1_0 \quad \text{such that} \quad a\left(u,v \right) =\left( f, v\right) \quad \forall \, v \in H^1_0\left(0,1\right). \end{equation*} $$
Theorem:

The following hold:

  1. Let $u$ be a $C^2$ be a solution of the elementary problem, then ${u \in H^1_0}$ also solves the weak form
  2. Let $u \in H^1_0$ be a solution of the weak problem. If and only if ${u \in C^2\left( \left[0,1 \right] \right)}$ then $u$ also solves the elementary problem.
Theorem: Fundamental Theorem of the Calculus of Variations

Suppose that $f$ is integrable on $(0,1)$ and

$$ \begin{equation*} \int_0^1 \phi f \, \mathrm{d}x = 0 \quad \forall \, \phi \in C^\infty_0\left(\left[ 0,1 \right]\right) \end{equation*} $$

then $f=0$.

Approximate $H_0^1$ by $V_h$. The discrete weak problem is then:

$$ \begin{equation*} \text{Find a} \quad u_h \in V_h \quad \text{such that} \quad a\left(u_h, v_h\right) = \left(f, v_h \right) \quad \forall \, v_h \in V_h \end{equation*} $$

Let $\lbrace \varphi_1, \varphi_2, \ldots, \varphi_N \rbrace$ be a basis of $V_h$, then, with $N=\text{dim}V_h$, so that

$$ \begin{equation*} u_h \left(x \right) = \sum\limits_{j=1}^N u_j \varphi_j\left(x\right). \end{equation*} $$

So the problem can be written as: Find ${\left(u_1, \ldots u_N \right) \in \mathbb{R}^N}$ such that

$$ \begin{equation*} \sum\limits_{j=1}^N u_j a\left( \varphi_j, \varphi_i \right) = \left( f, \varphi_i \right) \quad i=1, \ldots, N. \end{equation*} $$

Denote ${a_{ij}=a\left( \varphi_j, \varphi_i \right)}$ as the elements of the matrix $A$, let ${u=\left(u_1, \ldots, u_N\right)}$ and ${f=\left(f_1, \ldots, f_N\right)}$ be vectors where each entry is given by ${f_i = f \varphi_i}$, so that the problem is equivalent to solving the linear problem ${Au=f}$

Theorem: Poincaré-Friedrich Inequality

Let $\Omega \subset \mathbb{R}^n$ be contained in $n$-dimensional cube of length $s$, then

$$ \begin{equation*} \bigl\Vert v \bigr\Vert_{L^2\left(\Omega\right)} \le s \left| v \right|_{H_0^1\left( \Omega \right)}. \end{equation*} $$

For functions which are zero on the boundary a simplified form is

$$ \begin{equation*} \int_{a}^b \left| v(x) \right|^2 \, \mathrm{d}x \leq C_p \int_a^b \bigl\| v^\prime \left(x\right) \bigr\|^2 \mathrm{d}x \quad \forall \, v \in V_0 \end{equation*} $$
Theorem:

Let

$$ \begin{align*} C & = \dfrac{1}{\alpha_0} \Bigl( \bigl\Vert \alpha \bigr\Vert_\infty + C_p^2 \bigl\Vert \gamma \bigr\Vert_\infty \Bigr) \end{align*} $$

then

$$ \left| u - u_h \right|^{}_{H^1\left(0,1\right)} \le C \min\limits_{ w_h \in V_h} \left| u - w_h \right|_{H^1\left(0,1\right)}. $$
Definition: Coercivity and Continuity of Bilinear Forms

A bilinear form $a\left(\cdot,\cdot\right)$ on $V$, with a norm $\bigl\Vert \cdot \bigr\Vert_V$, then a bilinear form is coercive if there exists an ${\alpha_0 > 0}$ such that

$$ \begin{equation*} a(v, v) \geq \alpha_0 \bigl\Vert v \bigr\Vert^2 \quad \forall \, v \in V. \end{equation*} $$

A bilinear form is said to be continuous if there exists an ${M > 0}$ such that

$$ \begin{equation*} \left| a\left(u, v\right) \right| \leq M \bigl\Vert u \bigr\Vert_V \bigl\Vert u \bigr\Vert_V \quad \forall \, u, v \in V. \end{equation*} $$
Theorem: Lax-Milgram

If coercive and continuous, and the right hand side $(f, v)$ satisfies the following inequality

$$ \begin{equation*} \left| \left( f, v \right) \right| \leq K \bigl\Vert v\bigr\Vert_V \quad \forall \, v \in V. \end{equation*} $$

Then the weak and discrete weak form problems admit unique solutions which satisfy

$$ \begin{equation*} \bigl\Vert u \bigr\Vert_V \leq \dfrac{K}{\alpha_0} \quad \text{and} \quad \bigl\Vert u_h \bigr\Vert_V \leq \dfrac{K}{\alpha_0}. \end{equation*} $$
Lemma: Céa

It can be shown that

$$ \begin{equation*} \left\Vert u - u_h \right\Vert_V \le \dfrac{M}{\alpha_0} \min_{\substack{w_h \in V_h}} \left\Vert u - w_h \right\Vert_V . \tag{3} \end{equation*} $$

Finite Element Method

The finite element method (FEM) is a special technique for constructing a subspace $V_h$ based on piecewise polynomial interpolation.

Introduce a partition $\mathcal{T}_h$ of $[0,1]$ into $n$ subintervals $I_j = [x_j, x_{j+1}]$, $n \geq 2$, of width ${h_j~=~x_{j+1}~-~x_j}$, $j = 0, \ldots, n − 1$, with $0 = x_0 < x_1 < \ldots < x_{n−1} < x_n = 1$ and let $h~=~\max\limits_{\mathcal{T}_{h}}\left(h_{j}\right)$ .

Since functions in $\mathrm{H}_{0}^{1}(0,1)$ are continuous it makes sense to consider for ${k \geq 1}$ the family of piecewise polynomials $X_{h}^{k}$ introduced in $(2)$, but where now $[a, b]$ must be replaced by $[0,1]$. Any function $v_{h} \in X_{h}^{k}$ is a continuous piecewise polynomial over $[0,1]$ and its restriction over each interval $I_{j} \in \mathcal{T}_{h}$ is a polynomial of degree less than or equal to $k$.

Considering the cases $k=1$ and $k=2$, set

$$ \begin{equation*} V_{h} = X_{h}^{k, 0}=\lbrace v_{h} \in X_{h}^{k} \, : \, v_{h}(0) = v_{h}(1) = 0 \rbrace . \end{equation*} $$

The dimension $N$ of the finite element space $V_{h}$ is equal to ${n k -1}$. In the following the two cases $k=1$ and $k=2$ will be examined.

To assess the accuracy of the Galerkin FEM, first notice that, due to Céa’s lemma

$$ \begin{equation*} \min_{w_{h} \in V_{h}}\left\|u-w_{h}\right\|_{\mathrm{H}_{0}^{1}(0,1)} \leq\left\|u-\Pi_{h}^{k} u\right\|_{\mathrm{H}_{0}^{1}(0,1)} \tag{4} \end{equation*} $$

where $\Pi_{h}^{k} u$ is the interpolant of the exact solution $u \in V$ from the weak form of the governing equation. From inequality \eqref{eq:cea_estimate} estimating the Galerkin approximation error $\left\| u - u_{h} \right\|_{\mathrm{H}_{0}^{1}(0,1)}$ is then equivalent to estimating the interpolation error $\left\| u -\Pi_{h}^{k} u \right\|_{H_{0}^{1}(0,1)}$ . When ${k=1}$, using \eqref{eq:cea} and the bounds on the interpolation errors $(1)$

$$ \left\| u - u_{h} \right\|^{}_{\mathrm{H}_{0}^{1}(0,1)} \leq \frac{M}{\alpha_{0}} C h \|u\|_{\mathrm{H}^{2}(0,1)} $$

provided that $u \in \mathrm{H}^{2}(0,1)$. This estimate can be extended to the case ${k>1}$ as stated in the following convergence result.

$$ \begin{equation*} \left\| u - u_{h} \right\|_{ \mathrm{H}_{0}^{1}(0,1) } \end{equation*} $$
Theorem: Lax-Milgram

Let $u \in H_0^1\left( 0,1 \right)$ be the exact solution of

$$ \begin{equation*} a\left( u, v \right) = f\left(v\right) \quad \forall \, v \in H_0^1\left(0\right) \end{equation*} $$

and let $u_h \in V_h$ be it finite element approximation using a continuous piecewise polynomial of degree less than or equal to $k$, where $k \ge 1$. Furthermore, assume that $u \in \mathrm{H}^{s}(0,1)$ for some $s \geq 2$. Then the error is bounded as

$$ \begin{equation*} \bigl\Vert u - u_{h}\bigr\Vert_{\mathrm{H}_{0}^{1}(0,1)} \le \dfrac{M}{\alpha_0} C h^l \bigl\Vert u \bigr\Vert_{H^{l+1}(0,1)} \end{equation*} $$

where $l = \min\left( k, s-1 \right)$. Additionally, under the same assumptions it is possible to show that

$$ \begin{equation*} \bigl\Vert u - u_{h}\bigr\Vert_{\mathrm{L}^{2}(0,1)} \leq C h^{l+1} \bigl\Vert u \bigr\Vert_{\mathrm{H}^{l+1}(0,1)} . \end{equation*} $$

The error estimate shows that the Galerkin method is convergent, that is the approximation error tends to zero as $h \rightarrow 0$. The order of convergence is $k$.

Example:

An .ipynb notebook, with an example of the finite element method in one-dimension can be viewed here, or can be downloaded from here.

An example of the finite element method in two-dimensions can be viewed here, or can be downloaded from here.