Accelerating the pace of engineering and science

# Documentation

## Finite Element Method (FEM) Basics

The core Partial Differential Equation Toolbox™ algorithm is a PDE solver that uses the Finite Element Method (FEM) for problems defined on bounded domains in the plane.

The solutions of simple PDEs on complicated geometries can rarely be expressed in terms of elementary functions. You are confronted with two problems: First you need to describe a complicated geometry and generate a mesh on it. Then you need to discretize your PDE on the mesh and build an equation for the discrete approximation of the solution. The PDE app provides you with easy-to-use graphical tools to describe complicated domains and generate triangular meshes. It also discretizes PDEs, finds discrete solutions and plots results. You can access the mesh structures and the discretization functions directly from the command line (or from a file) and incorporate them into specialized applications.

Here is an overview of the Finite Element Method (FEM). The purpose of this presentation is to get you acquainted with the elementary FEM notions. Here you find the precise equations that are solved and the nature of the discrete solution. Different extensions of the basic equation implemented in Partial Differential Equation Toolbox software are presented. A more detailed description can be found in Elliptic Equations, with variants for specific types in Systems of PDEs, Parabolic Equations, Hyperbolic Equations, Eigenvalue Equations, and Nonlinear Equations.

You start by approximating the computational domain Ω with a union of simple geometric objects, in this case triangles. The triangles form a mesh and each vertex is called a node. You are in the situation of an architect designing a dome. He has to strike a balance between the ideal rounded forms of the original sketch and the limitations of his simple building-blocks, triangles or quadrilaterals. If the result does not look close enough to a perfect dome, the architect can always improve his work using smaller blocks.

Next you say that your solution should be simple on each triangle. Polynomials are a good choice: they are easy to evaluate and have good approximation properties on small domains. You can ask that the solutions in neighboring triangles connect to each other continuously across the edges. You can still decide how complicated the polynomials can be. Just like an architect, you want them as simple as possible. Constants are the simplest choice but you cannot match values on neighboring triangles. Linear functions come next. This is like using flat tiles to build a waterproof dome, which is perfectly possible.

A Triangular Mesh (left) and a Continuous Piecewise Linear Function on That Mesh

Now you use the basic elliptic equation (expressed in Ω)

$-\nabla \cdot \left(c\nabla u\right)+au=f.$

If uh is the piecewise linear approximation to u, it is not clear what the second derivative term means. Inside each triangle, ∇uh is a constant (because uh is flat) and thus the second-order term vanishes. At the edges of the triangles, cuh is in general discontinuous and a further derivative makes no sense.

What you are looking for is the best approximation of u in the class of continuous piecewise polynomials. Therefore you test the equation for uh against all possible functions v of that class. Testing means formally to multiply the residual against any function and then integrate, i.e., determine uh such that

$\underset{\Omega }{\int }\left(-\nabla ·\left(c\nabla {u}_{h}\right)+a{u}_{h}-f\right)v\text{\hspace{0.17em}}dx=0$

for all possible v. The functions v are usually called test functions.

Partial integration (Green's formula) yields that uh should satisfy

$\underset{\Omega }{\int }\left(\left(c\nabla {u}_{h}\right)\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\nabla v+a{u}_{h}v\right)dx-\underset{\partial \Omega }{\int }\stackrel{\to }{n}\text{\hspace{0.17em}}·\text{\hspace{0.17em}}\left(c\nabla {u}_{h}\right)v\text{\hspace{0.17em}}ds=\underset{\Omega }{\int }fv\text{\hspace{0.17em}}dx\text{ }\forall v,$

where ∂Ω is the boundary of Ω and $\stackrel{\to }{n}$ is the outward pointing normal on ∂Ω. The integrals of this formulation are well-defined even if uh and v are piecewise linear functions.

Boundary conditions are included in the following way. If uh is known at some boundary points (Dirichlet boundary conditions), we restrict the test functions to v = 0 at those points, and require uh to attain the desired value at that point. At all the other points we ask for Neumann boundary conditions, i.e., $\left(c\nabla {u}_{h}\right)\cdot \stackrel{\to }{n}+q{u}_{h}=g$. The FEM formulation reads: Find uh such that

where ∂Ω1 is the part of the boundary with Neumann conditions. The test functions v must be zero on ∂Ω – ∂Ω1.

Any continuous piecewise linear uh is represented as a combination

${u}_{h}\left(x\right)=\sum _{i=1}^{N}{U}_{i}{\varphi }_{i}\left(x\right),$

where ϕi are some special piecewise linear basis functions and Ui are scalar coefficients. Choose ϕi like a tent, such that it has the "height" 1 at the node i and the height 0 at all other nodes. For any fixed v, the FEM formulation yields an algebraic equation in the unknowns Ui. You want to determine N unknowns, so you need N different instances of v. What better candidates than v = ϕi, i = 1, 2, ... , N? You find a linear system KU = F where the matrix K and the right side F contain integrals in terms of the test functions ϕi, ϕj, and the coefficients defining the problem: c, a, f, q, and g. The solution vector U contains the expansion coefficients of uh, which are also the values of uh at each node xi, since uh(xi) = Ui.

If the exact solution u is smooth, then FEM computes uh with an error of the same size as that of the linear interpolation. It is possible to estimate the error on each triangle using only uh and the PDE coefficients (but not the exact solution u, which in general is unknown).

There are Partial Differential Equation Toolbox functions that assemble K and F. This is done automatically in the PDE app, but you also have direct access to the FEM matrices from the command-line function assempde.

To summarize, the FEM approach is to approximate the PDE solution u by a piecewise linear function uh. The function uh is expanded in a basis of test-functions ϕi, and the residual is tested against all the basis functions. This procedure yields a linear system KU = F. The components of U are the values of uh at the nodes. For x inside a triangle, uh(x) is found by linear interpolation from the nodal values.

FEM techniques are also used to solve more general problems. The following are some generalizations that you can access both through the PDE app and with command-line functions.

• Time-dependent problems are easy to implement in the FEM context. The solution u(x,t) of the equation

$d\frac{\partial u}{\partial t}-\nabla \cdot \left(c\nabla u\right)+au=f,$

can be approximated by

${u}_{h}\left(x,t\right)=\sum _{i=1}^{N}{U}_{i}\left(t\right){\varphi }_{i}\left(x\right).$

• This yields a system of ordinary differential equations (ODE)

$M\frac{dU}{dt}+KU=F,$

which you integrate using ODE solvers. Two time derivatives yield a second order ODE

$M\frac{{d}^{2}U}{d{t}^{2}}+KU=F,$

etc. The toolbox supports problems with one or two time derivatives (the functions parabolic and hyperbolic).

• Eigenvalue problems: Solve

$-\nabla \cdot \left(c\nabla u\right)+au=\lambda du,$

for the unknowns u and λ (λ is a complex number). Using the FEM discretization, you solve the algebraic eigenvalue problem KU = λhMU to find uh and λh as approximations to u and λ. A robust eigenvalue solver is implemented in pdeeig.

• If the coefficients c, a, f, q, or g are functions of u or ∇u, the PDE is called nonlinear and FEM yields a nonlinear system K(U)U = F(U). You can use iterative methods for solving the nonlinear system. For elliptic equations, the toolbox provides a nonlinear solver called pdenonlin using a damped Gauss-Newton method. The parabolic and hyperbolic functions call the nonlinear solver automatically.

• Small triangles are needed only in those parts of the computational domain where the error is large. In many cases the errors are large in a small region and making all triangles small is a waste of computational effort. Making small triangles only where needed is called adapting the mesh refinement to the solution. An iterative adaptive strategy is the following: For a given mesh, form and solve the linear system KU = F. Then estimate the error and refine the triangles in which the error is large. The iteration is controlled by adaptmesh and the error is estimated by pdejmps.

Although the basic equation is scalar, systems of equations are also handled by the toolbox. The interactive environment accepts u as a scalar or 2-vector function. In command-line mode, systems of arbitrary size are accepted.

If cδ > 0 and a ≥ 0, under rather general assumptions on the domain Ω and the boundary conditions, the solution u exists and is unique. The FEM linear system has a unique solution which converges to u as the triangles become smaller. The matrix K and the right side F make sense even when u does not exist or is not unique. It is advisable that you devise checks to problems with questionable solutions.

## References

[1] Cook, Robert D., David S. Malkus, and Michael E. Plesha, Concepts and Applications of Finite Element Analysis, 3rd edition, John Wiley & Sons, New York, 1989.