Partial Differential Equation Toolbox

Basics of the Finite Element Method

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 pdetool graphical user interface 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 M-file) and incorporate them into specialized applications.

Below 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 the PDE Toolbox are presented. A more detailed description can be found in The Finite Element Method.

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 )

If uh is the piecewise linear approximation to u, it is not clear what the second derivative term means. Inside each triangle, is a constant (because uh is flat) and thus the second-order term vanishes. At the edges of the triangles, 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

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

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

where is the boundary of and is the outward pointing normal on . Note that 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., . The FEM formulation reads: Find uh such that

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

Any continuous piecewise linear uh is represented as a combination

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 = j, j = 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).

The PDE Toolbox provides functions that assemble K and F. This is done automatically in the graphical user interface, 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 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. Below are some generalizations that you can access both through the graphical user interface and with command-line functions.

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

• can be approximated by

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

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

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

• Eigenvalue problems: Solve

• 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, 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. The toolbox provides a nonlinear solver called pdenonlin using a damped Gauss-Newton method.
• 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.

Getting Started

Using the PDE Toolbox Graphical User Interface