The core Partial Differential Equation Toolbox™ algorithm uses the Finite Element Method (FEM) for problems defined on bounded domains in 2-D or 3-D space. In most cases, elementary functions cannot express the solutions of even simple PDEs on complicated geometries. The finite element method describes a complicated geometry as a collection of subdomains by generating a mesh on the geometry. For example, you can approximate the computational domain Ω with a union of triangles (2-D geometry) or tetrahedra (3-D geometry). The subdomains form a mesh, and each vertex is called a node. The next step is to approximate the original PDE problem on each subdomain by using simpler equations.

For example, consider the basic elliptic equation.

$$-\nabla \cdot \left(c\nabla u\right)+au=f\text{ondomain}\Omega $$

Suppose that this equation is a subject to the Dirichlet boundary condition $$u=r$$ on $$\partial {\Omega}_{D}$$ and Neumann boundary conditions on $$\partial {\Omega}_{N}$$. Here, $$\partial \Omega =\partial {\Omega}_{D}\cup \partial {\Omega}_{N}$$ is the boundary of Ω.

The first step in FEM is to convert the original differential
(*strong*) form of the PDE into an integral (*weak*)
form by multiplying with test function $$v$$ and integrating over the domain
Ω.

$$\underset{\Omega}{\int}\left(-\nabla \xb7\left(c\nabla u\right)+au-f\right)v\text{\hspace{0.17em}}d\Omega}=0\text{\hspace{1em}}\forall v$$

The test functions are chosen from a collection of functions
(functional space) that vanish on the Dirichlet portion of the boundary, $$v=0$$ on $$\partial {\Omega}_{D}$$. Above equation can be thought
of as weighted averaging of the residue using all possible weighting
functions $$v$$. The collection of functions
that are admissible solutions, *u*, of the weak form
of PDE are chosen so that they satisfy the Dirichlet BC, *u* = *r* on $$\partial {\Omega}_{D}$$.

Integrating by parts (Green's formula) the second-order term results in:

$$\underset{\Omega}{\int}\left(c\nabla u\text{\hspace{0.17em}}\nabla v+auv\right)d\Omega}-{\displaystyle \underset{\partial {\Omega}_{N}}{\int}\overrightarrow{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\nabla u\right)}\text{\hspace{0.17em}}v\text{\hspace{0.17em}}d\partial {\Omega}_{N}+{\displaystyle \underset{\partial {\Omega}_{D}}{\int}\overrightarrow{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\nabla u\right)}\text{\hspace{0.17em}}v\text{\hspace{0.17em}}d\partial {\Omega}_{D}={\displaystyle \underset{\Omega}{\int}fv\text{\hspace{0.17em}}d\Omega}\text{\hspace{1em}}\forall v$$

Use the Neumann boundary condition to substitute for second term on the left side of the equation. Also, note that $$v=0$$ on $$\partial {\Omega}_{D}$$ nullifies the third term. The resulting equation is:

$$\underset{\Omega}{\int}\left(c\nabla u\text{\hspace{0.17em}}\nabla v+auv\right)d\Omega}+{\displaystyle \underset{\partial {\Omega}_{N}}{\int}quv}\text{\hspace{0.17em}}d\partial {\Omega}_{N}={\displaystyle \underset{\partial {\Omega}_{N}}{\int}gv}\text{\hspace{0.17em}}d\partial {\Omega}_{N}+{\displaystyle \underset{\Omega}{\int}fv\text{\hspace{0.17em}}d\Omega}\text{\hspace{1em}}\forall v$$

Note that all manipulations up to this stage are performed on continuum Ω, the global domain of the problem. Therefore, the collection of admissible functions and trial functions span infinite-dimensional functional spaces. Next step is to discretize the weak form by subdividing Ω into smaller subdomains or elements $${\Omega}^{e}$$, where $$\Omega =\cup {\Omega}^{e}$$. This step is equivalent to projection of the weak form of PDEs onto a finite-dimensional subspace. Using the notations $${u}_{h}$$ and $${v}_{h}$$ to represent the finite-dimensional equivalent of admissible and trial functions defined on $${\Omega}^{e}$$, you can write the discretized weak form of the PDE as:

$$\underset{{\Omega}^{e}}{\int}\left(c\nabla {u}_{h}\nabla {v}_{h}+a{u}_{h}{v}_{h}\right)\text{\hspace{0.17em}}d{\Omega}^{e}}+{\displaystyle \underset{\partial {\Omega}_{N}^{e}}{\int}q{u}_{h}v{\text{\hspace{0.17em}}}_{h}d\partial {\Omega}_{N}^{e}}={\displaystyle \underset{\partial {\Omega}_{N}^{e}}{\int}gv{\text{\hspace{0.17em}}}_{h}d\partial {\Omega}_{N}^{e}}+{\displaystyle \underset{{\Omega}^{e}}{\int}f{v}_{h}d{\Omega}^{e}}\text{\hspace{1em}}\forall {v}_{h$$

Next, let *ϕ*_{i},
with *i* = 1, 2, ... ,* N*_{p},
be the piecewise polynomial basis functions for the subspace containing
the collections $${u}_{h}$$ and $${v}_{h}$$, then any particular
can be expressed
as a linear combination of basis functions:

$${u}_{h}={\displaystyle \sum _{1}^{{N}_{p}}{U}_{i}{\varphi}_{i}}$$

Here *U*_{i}
are yet undetermined scalar coefficients. Substituting
into to the discretized
weak form of PDE and using each $${v}_{h}={\phi}_{i}$$ as
test functions and performing integration over element yields a system
of *N*_{p} equations
in terms of *N*_{p} unknowns *U*_{i}.

Note that finite element method approximates a solution by minimizing
the associated error function. The minimizing process automatically
finds the linear combination of basis functions which is closest to
the solution *u*.

FEM yields a 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 *c*, *a*, *f*, *q*,
and *g* defining the problem. The solution vector *U* contains
the expansion coefficients of *u _{h}*,
which are also the values of

FEM techniques are also used to solve more general problems, such as:

Time-dependent problems. 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)$$

The result is a system of ordinary differential equations (ODEs)

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

Two time derivatives result in a second-order ODE

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

Eigenvalue problems. Solve

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

for the unknowns

*u*and*λ*, where*λ*is a complex number. Using the FEM discretization, you solve the algebraic eigenvalue problem*KU*=*λ**MU*to find*u*as an approximation to_{h}*u*. To solve eigenvalue problems, use`solvepdeeig`

.Nonlinear problems. 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*).

To summarize, the FEM approach:

Represents the original domain of the problem as a collection of elements.

For each element, substitutes the original PDE problem by a set of simple equations that locally approximate the original equations. Applies boundary conditions for boundaries of each element. For stationary linear problems where the coefficients do not depend on the solution or its gradient, the result is a linear system of equations. For stationary problems where the coefficients depend on the solution or its gradient, the result is a system of nonlinear equations. For time-dependent problems, the result is a set of ODEs.

Assembles the resulting equations and boundary conditions into a global system of equations that models the entire problem.

Solves the resulting system of algebraic equations or ODEs using linear solvers or numerical integration, respectively. The toolbox internally calls appropriate MATLAB

^{®}solvers for this task.

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

[2] Gilbert Strang and George Fix *An Analysis
of the Finite Element Method*. 2nd edition. Wellesley,
MA: Wellesley-Cambridge Press, 2008.

Was this topic helpful?