# 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.

## Discretization

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 errorConsider 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{aligned} \tilde b &= b \\ \tilde c &= \frac{a}{8} h^2 + c. \end{aligned}$

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 AdaptationConsider 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?

## Interpolation

## Partial Differential Equations

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.