Finite Elements

The finite element method (FEM) is a numerical method for solving partial differential equations (PDE) that occur in problems of engineering and mathematical physics. The basic concept of FEM is to divide continuous bodies into a mesh of simple parts, the so-called finite elements. Functions are represented by their values at certain support points of the mesh, so that the differential equation is discretized and replaced by a linear algebraic equation system.

As a practical example, we consider a wrench whose mechanical strain is to be determined under an external stress. Instead of a continuous mass distribution, let us assume that the wrench consists of a number \(n\) of mass points \(m_i\) connected by mechanical springs of different spring constants \(k_{ij}\). Then the static deflection of the individual mass points is obtained from the equilibrium of all forces acting on the point, resulting in a coupled equation system. This approximation for a real physical system is an example for a finite element method, that can also be applied to other problems in a more general form.

Nature and technology are full of complex forms, so that mathematical models are often difficult to implement in reality. A famous example is the determination of the circumference and area of ​​a circle, i.e. the calculation of \(\pi\). To approximate the circle, it can be replaced by a regular polygon with \(n\) sides, and by its circumference \(C_n\) an approximation for \(\pi\) can be obtained. In the limit \(n \to \infty\), the circumference converges to the true value \(\lim_{n\to \infty} C_n = 2 \pi r \) for the circle:

Here we have replaced the complex figure--the circle--by simple triangles. When using finite elements, a two-dimensional area is also subdivided into triangles, so that a mesh is created (triangulation). In similar fashion, a one-dimensional line is seperated into line segments, whereas a three-dimensional volume can be divided into tetrahedrons:

In the simplest case, there is a regular grid with right angles and uniform edge lengths \(h_i\), so that the nodal points in 2 dimensions are given by

\[ \left( \begin{array}{c} x_k \\ y_k \end{array} \right) = \left( \begin{array}{c} x_0 \\ y_0 \end{array} \right) + \left( \begin{array}{c} h_x \cdot k_x \\ h_y \cdot k_y \end{array} \right) \]

with indices \(k_x \in \{0, \dots, n_x-1\} \), \(k_y \in \{ 0, \dots, n_y-1 \} \). There are a total number of \(n = n_x \cdot n_y\) grid points, that can be labeled with a multi-index \(k = k_x + n_x \cdot k_y \in \{ 0, \dots, n-1\} \). Of the continuous rectangular shape, only a discrete number of points is considered.

In a similar way, physical fields \(f(x, y)\)--such as the temperature distribution on a weather map--are reduced to their function values \(​​f_k = f (x_k, y_k) \) at the grid points. In the world of finite elements, a function \(f\) therefore represents an \(n\)-dimensional vector

\[ \underline{f} = \left(\begin{array}{n} f_0 \\ f_1 \\ \vdots \\ f_{n-1} \end{array} \right) = \left(\begin{array}{n} f(x_0, y_0) \\ f(x_1,y_1) \\ \vdots \\ f(x_{n-1},y_{n-1}) \end{array} \right) \]

and can be stored as an array in computer memory. (This representation of a function is reminiscent of mathematical tables used before the introduction of pocket calculators and smartphones, in order to calculate trigonometric functions, for example.)

In order to approximate the function value \(f(x,y)\) at any point within a finite element, it is necessary to interpolate the values \(f_k\) ​​at the node points. In the simplest case, a linear interpolation can be performed, which can be calculated in one dimension by

\[f(x) \approx f_k + \frac{f_{k+1} - f_{k}}{x_{k+1} - x_{k}} (x - x_k) \quad x \in [x_k,x_{k+1}].\]

In the following section, we will show how interpolation methods can be generalized by the help of base function to higher dimensions and higher polynomials.

The discretization of the field \(f (x, y)\) is, of course, more accurate the smaller the edge length \(h\) of the triangles are. The value of the parameter \(h\) depends on the desired accuracy in the interpolation and on the exact function course. If the function \(f (x, y)\) already represents a linear plane by itself, the interpolation exactly matches the function. The interpolation error therefore depends on the curvature of the function, which results in one dimension by the second derivative \(f''(x)\) and in two dimensions by the eigenvalues ​​of the Hessian matrix

\[ \mathbf{H} = \left(\begin{array}{cc} \frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \partial y} \\ \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial y^2} \end{array} \right). \]

For an accurate and efficient representation of a function, the local node density of the mesh must be adapted to the function profile (mesh adaptation). While areas with low curvature can be interpolated with a coarse grid with sufficient accuracy, areas with high curvature in the function require a high node density.

Interpolation error

Consider the function \(f: \big[-\frac{h}{2}, \frac{h}{2}\big] \to \mathbb{R}\) with

\[ f(x) = \frac{a}{2} x^2 + b x + c \]

and its linear interpolation

\[ \tilde f(x) = \tilde b x + \tilde c \]

with \(\tilde f\big(-\frac{h}{2}\big) = f\big(-\frac{h}{2}\big)\) and \(\tilde f\big(\frac{h}{2}\big) = f\big(\frac{h}{2}\big)\). What is the maximal value of the interpolation error

\[ e = \big|f(x) - \tilde f(x)\big| \]

on the interval \(\big[-\frac{h}{2}, \frac{h}{2}\big]?\)

---

The equations \(\tilde f\big(-\frac{h}{2}\big) = f\big(-\frac{h}{2}\big)\) and \(\tilde f\big(\frac{h}{2}\big) = f\big(\frac{h}{2}\big)\) are solved by

\[ \begin{align*} \tilde b &= b \\ \tilde c &= \frac{a}{8} h^2 + c. \end{align*}\]

The interpolation error reads

\[e(x) = \frac{|a|}{2} \left(x - \frac{h}{2} \right)^2\]

with a maximal value

\[\max_x e(x) = e(0) = \frac{|a|}{8} h^2.\]

The error depends linearly on the absolute curvature \(|a| = \big|f''\big|\) and quadratically on the edge length \(h\).

Mesh Adaptation

Consider the ''hat-function'' \(f(x, y) = \exp\left(-50 \cdot \big(x^2+y^2\big)^4\right)\) on the unit circle \(x^2 + y^2 \leq 1,\) as shown below:

Which of the meshes shown below is optimally suited for the linear interpolation of the function?

Probably the best-known example for a PDE is the so-called Poisson equation, which is in general

\[ - \vec \nabla \cdot \big(k \vec \nabla u\big) = f \]

with the desired solution function \(u\), a material parameter \(k,\) and a density function \(f\). For a homogeneous medium \((k = \text{const})\) in \(d = 3\) dimensions, the PDE results in

\[ - \left[\frac{\partial^2}{\partial x^2} + \frac{\partial^2 }{\partial y^2} + \frac{\partial^2 }{\partial z^2} \right] u(x,y,z) = \frac{1}{k} f(x,y,z). \]

For a given density \(f\), the solution \(u: \Omega \to \mathbb{R}\) of the Poisson equation on a certain domain \(\Omega \subset \mathbb{R}^d\) is determined by the boundary conditions on its surface \(\partial \Omega\), so that this problem is called a boundary value problem. The Poisson equation appears in different physical contexts; for example,

• electrostatics with an electric potential \(u = \phi\), dielectric constant \(k = \varepsilon,\) and charge density \(f = \rho\)
• heat diffusion with a temperature field \(u = T\), heat conductivity \(k,\) and heat source \(f = \dot q\)
• elasticity with a displacement vector \(u = \vec u\), tensor of elastic constants \(k = \, \stackrel{=}{C},\) and a force density \(f = \vec f\)

and many more. The Poisson equation determines, for example, the electric field in the interior of a capacitor due the fixed electric potential on the metal plates as boundary condition.