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 (2-D geometry) or tetrahedra (3-D geometry). (This discussion applies to both triangles and tetrahedra, but speaks of triangles.) The triangles form a mesh and each vertex is called a node. You are in the situation of an architect designing a dome. The architect has to strike a balance between the ideal rounded forms of the original sketch and the limitations of the simple building-blocks, triangles or quadrilaterals. If the result does not look close enough to a perfect dome, the architect can always improve the result by 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 *u _{h}* is the piecewise
linear approximation to

What you are looking for is the best approximation of *u* in
the class of continuous piecewise polynomials. Therefore you test
the equation for *u _{h}* against
all possible functions

$$\underset{\Omega}{\int}\left(-\nabla \xb7\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 *u _{h}* should
satisfy

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

where ∂Ω is the boundary of Ω and $$\overrightarrow{n}$$ is the outward pointing normal
on ∂Ω. The integrals of this formulation are well-defined
even if *u _{h}* and

Boundary conditions are included in the following way. If *u _{h}* is
known at some boundary points (Dirichlet boundary conditions), we
restrict the test functions to

$${\int}_{\Omega}^{}\left(\left(c\nabla {u}_{h}\right)\cdot \nabla v+a{u}_{h}v\right)\text{\hspace{0.17em}}dx}+{\displaystyle {\int}_{\partial {\Omega}_{1}}^{}q{u}_{h}v\text{\hspace{0.17em}}ds}={\displaystyle {\int}_{\Omega}f}v\text{\hspace{0.17em}}dx+{\displaystyle {\int}_{\partial {\Omega}_{1}}^{}gv\text{\hspace{0.17em}}ds}\text{}\forall v,$$

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

Any continuous piecewise linear *u _{h}* is
represented as a combination

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

where *ϕ*_{i} are
some special piecewise linear basis functions and *U*_{i} 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 *U*_{i}.
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 *u _{h}*,
which are also the values of

If the exact solution *u* is smooth, then FEM
computes *u _{h}* 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

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 *u _{h}*.
The function

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}(x,t)={\displaystyle \sum _{i=1}^{N}{U}_{i}(t){\varphi}_{i}}(x).$$

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`

).$$-\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*=*λ*_{h}*MU*to find*u*and_{h}*λ*as approximations to_{h}*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

*K**U*=*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.

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

Was this topic helpful?