In mathematics, the finite element method (FEM) is a numerical technique for finding approximate

solutions to boundary value problems for differential equations. It uses variational

methods (the calculus of variations) to minimize an error function and produce a stable solution.
Analogous to the idea that connecting many tiny straight lines can approximate a larger circle, FEM
encompasses all the methods for connecting many simple element equations over many small
subdomains, named finite elements, to approximate a more complex equation over a larger domain.
Contents
[hide]
1 History
2 Technical discussion
o 2.1 General principles
o 2.2 Illustrative problems P1 and P2
o 2.3 Weak formulation
2.3.1 The weak form of P1
2.3.2 The weak form of P2
2.3.3 A proof outline of existence and uniqueness of the solution
3 Discretization
o 3.1 For problem P1
o 3.2 For problem P2
o 3.3 Choosing a basis
o 3.4 Small support of the basis
o 3.5 Matrix form of the problem
o 3.6 General form of the finite element method
4 Various types of finite element methods
o 4.1 AEM
o 4.2 Generalized finite element method
o 4.3 Mixed finite element method
o 4.4 hp-FEM
o 4.5 hpk-FEM
o 4.6 XFEM
o 4.7 S-FEM
o 4.8 Fiber Beam Method
o 4.9 Spectral methods
o 4.10 Meshfree methods
o 4.11 Discontinuous Galerkin methods
o 4.12 Finite element limit analysis
o 4.13 Stretched grid method
5 Comparison to the finite difference method
6 Application
7 See also
8 References
9 Further reading
10 External links
History[edit]
While it is difficult to quote a date of the invention of the finite element method, the method originated
from the need to solve complex elasticity and structural analysis problems incivil and aeronautical
engineering. Its development can be traced back to the work by A. Hrennikoff and R. Courant. In
China, in the later 1950s and early 1960s, based on the computations of dam constructions, K.
Feng proposed a systematic numerical method for solving partial differential equations. The method
was called the finite difference method based on variation principle, which was another independent
invention of finite element method. Although the approaches used by these pioneers are different,
they share one essential characteristic: mesh discretization of a continuous domain into a set of
discrete sub-domains, usually called elements.
Hrennikoff's work discretizes the domain by using a lattice analogy, while Courant's approach divides
the domain into finite triangular subregions to solve second order ellipticpartial differential equations
(PDEs) that arise from the problem of torsion of a cylinder. Courant's contribution was evolutionary,
drawing on a large body of earlier results for PDEs developed by Rayleigh, Ritz, and Galerkin.
The finite element method obtained its real impetus in the 1960s and 1970s by the developments
of J. H. Argyris with co-workers at the University of Stuttgart, R. W. Clough with co-workers at UC
Berkeley, O. C. Zienkiewicz with co-workers at the University of Swansea, Philippe G. Ciarlet at the
University of Paris 6 and Richard Gallagher
[1]
with co-workers at Cornell University. Further impetus
was provided in these years by available open source finite element software programs. NASA
sponsored the original version of NASTRAN, and UC Berkeley made the finite element program
SAP IV
[2]
widely available. A rigorous mathematical basis to the finite element method was provided
in 1973 with the publication by Strang and Fix.
[3]
The method has since been generalized for
the numerical modeling of physical systems in a wide variety of engineering disciplines,
e.g., electromagnetism,heat transfer, and fluid dynamics.
[4][5]

Technical discussion[edit]
General principles[edit]
The subdivision of a whole domain into simpler parts has several advantages:
[6]

Accurate representation of complex geometry
Inclusion of dissimilar material properties
Easy representation of the total solution
Capture of local effects.
A typical work out of the method involves (1) dividing the domain of the problem into a collection of
subdomains, with each subdomain represented by a set of element equations to the original
problem, followed by (2) systematically recombining all sets of element equations into a global
system of equations for the final calculation. The global system of equations has known solution
techniques, and can be calculated from the initial values of the original problem to obtain a numerical
answer.
In the first step above, the element equations are simple equations that locally approximate the
original complex equations to be studied, where the original equations are oftenpartial differential
equations (PDE). To explain the approximation in this process, FEM is commonly introduced as a
special case of Galerkin method. The process, in mathematics language, is to construct an integral
of the inner product of the residual and the weight functions and set the integral to zero. In simple
terms, it is a procedure that minimizes the error of approximation by fitting trial functions into the
PDE. The residual is the error caused by the trial functions, and the weight functions
are polynomial approximation functions that project the residual. The process eliminates all the
spatial derivatives from the PDE, thus approximating the PDE locally with
a set of algebraic equations for steady state problems,
a set of ordinary differential equations for transient problems.
These equation sets are the element equations. They are linear if the underlying PDE is linear, and
vice versa. Algebraic equation sets that arise in the steady state problems are solved
using numerical linear algebra methods, while ordinary differential equation sets that arise in the
transient problems are solved by numerical integration using standard techniques such as Euler's
method or the Runge-Kutta method.
In step (2) above, a global system of equations is generated from the element equations through a
transformation of coordinates from the subdomains' local nodes to the domain's global nodes. This
spatial transformation includes appropriate orientation adjustments as applied in relation to the
reference coordinate system. The process is often carried out by FEM software
using coordinate data generated from the subdomains.
FEM is best understood from its practical application, known as finite element analysis (FEA). FEA
as applied in engineering is a computational tool for performing engineering analysis. It includes the
use of mesh generation techniques for dividing a complex problem into small elements, as well as
the use of software program coded with FEM algorithm. In applying FEA, the complex problem is
usually a physical system with the underlying physics such as the Euler-Bernoulli beam equation,
the heat equation, or the Navier-Stokes equations expressed in either PDE or integral equations,
while the divided small elements of the complex problem represent different areas in the physical
system.
FEA is a good choice for analyzing problems over complicated domains (like cars and oil pipelines),
when the domain changes (as during a solid state reaction with a moving boundary), when the
desired precision varies over the entire domain, or when the solution lacks smoothness. For
instance, in a frontal crash simulation it is possible to increase prediction accuracy in "important"
areas like the front of the car and reduce it in its rear (thus reducing cost of the simulation). Another
example would be in numerical weather prediction, where it is more important to have accurate
predictions over developing highly nonlinear phenomena (such as tropical cyclones in the
atmosphere, or eddies in the ocean) rather than relatively calm areas.

FEM mesh created by an analyst prior to finding a solution to a magnetic problem using FEM software. Colours
indicate that the analyst has set material properties for each zone, in this case a conductingwire coil in orange;
a ferromagnetic component (perhaps iron) in light blue; and air in grey. Although the geometry may seem simple, it
would be very challenging to calculate the magnetic field for this setup without FEM software, using equations alone.

FEM solution to the problem at left, involving acylindrically shaped magnetic shield. Theferromagnetic cylindrical part
is shielding the area inside the cylinder by diverting the magnetic fieldcreated by the coil (rectangular area on the
right). The color represents the amplitude of the magnetic flux density, as indicated by the scale in the inset legend,
red being high amplitude. The area inside the cylinder is low amplitude (dark blue, with widely spaced lines of
magnetic flux), which suggests that the shield is performing as it was designed to.
Illustrative problems P1 and P2[edit]
We will illustrate the finite element method using two sample problems from which the general
method can be extrapolated. It is assumed that the reader is familiar with calculusand linear algebra.
P1 is a one-dimensional problem

where is given, is an unknown function of , and is the second derivative of with
respect to .
P2 is a two-dimensional problem (Dirichlet problem)

where is a connected open region in the plane whose boundary is "nice"
(e.g., a smooth manifold or a polygon), and and denote the second derivatives
with respect to and , respectively.
The problem P1 can be solved "directly" by computing antiderivatives. However, this method
of solving the boundary value problem works only when there is one spatial dimension and
does not generalize to higher-dimensional problems or to problems like . For
this reason, we will develop the finite element method for P1 and outline its generalization to
P2.
Our explanation will proceed in two steps, which mirror two essential steps one must take to
solve a boundary value problem (BVP) using the FEM.
In the first step, one rephrases the original BVP in its weak form. Little to no computation
is usually required for this step. The transformation is done by hand on paper.
The second step is the discretization, where the weak form is discretized in a finite-
dimensional space.
After this second step, we have concrete formulae for a large but finite-dimensional linear
problem whose solution will approximately solve the original BVP. This finite-dimensional
problem is then implemented on a computer.
Weak formulation[edit]
The first step is to convert P1 and P2 into their equivalent weak formulations.
The weak form of P1[edit]
If solves P1, then for any smooth function that satisfies the displacement boundary
conditions, i.e. at and , we have
(1)
Conversely, if with satisfies (1) for every smooth
function then one may show that this will solve P1. The proof is easier for twice
continuously differentiable (mean value theorem), but may be proved in
a distributional sense as well.
By using integration by parts on the right-hand-side of (1), we obtain
(2)
where we have used the assumption that .
The weak form of P2[edit]
If we integrate by parts using a form of Green's identities, we see that if solves P2, then
for any ,

where denotes the gradient and denotes the dot product in the two-dimensional
plane. Once more can be turned into an inner product on a suitable space of
"once differentiable" functions of that are zero on . We have also assumed
that (see Sobolev spaces). Existence and uniqueness of the solution can
also be shown.
A proof outline of existence and uniqueness of the solution[edit]
We can loosely think of to be the absolutely continuous functions
of that are at and (see Sobolev spaces). Such functions are
(weakly) "once differentiable" and it turns out that the symmetric bilinear map then
defines an inner product which turns into a Hilbert space (a detailed proof is
nontrivial). On the other hand, the left-hand-side is also an inner
product, this time on the Lp space . An application of the Riesz representation
theorem for Hilbert spaces shows that there is a unique solving (2) and therefore P1.
This solution is a-priori only a member of , but using elliptic regularity, will be
smooth if is.
Discretization[edit]


A function in with zero values at the endpoints (blue), and a piecewise linear approximation
(red)
P1 and P2 are ready to be discretized which leads to a common sub-problem (3). The
basic idea is to replace the infinite-dimensional linear problem:
Find such that

with a finite-dimensional version:
(3) Find such that

where is a finite-dimensional subspace of . There are many
possible choices for (one possibility leads to the spectral method).
However, for the finite element method we take to be a space of
piecewise polynomial functions.
For problem P1[edit]
We take the interval , choose values
of with and we
define by:

where we define and . Observe that functions
in are not differentiable according to the elementary definition of
calculus. Indeed, if then the derivative is typically not
defined at any , . However, the derivative
exists at every other value of and one can use this derivative for
the purpose ofintegration by parts.


A piecewise linear function in two dimensions
For problem P2[edit]
We need to be a set of functions of . In the figure on the right,
we have illustrated a triangulation of a 15
sided polygonal region in the plane (below), and a piecewise
linear function (above, in color) of this polygon which is linear on
each triangle of the triangulation; the space would consist of
functions that are linear on each triangle of the chosen
triangulation.
One often reads instead of in the literature. The reason is
that one hopes that as the underlying triangular grid becomes finer
and finer, the solution of the discrete problem (3) will in some sense
converge to the solution of the original boundary value problem P2.
The triangulation is then indexed by a real valued
parameter which one takes to be very small. This
parameter will be related to the size of the largest or average
triangle in the triangulation. As we refine the triangulation, the
space of piecewise linear functions must also change with ,
hence the notation . Since we do not perform such an analysis,
we will not use this notation.
Choosing a basis[edit]


Basis functions vk (blue) and a linear combination of them, which is
piecewise linear (red)
To complete the discretization, we must select a basis of . In the
one-dimensional case, for each control point we will choose the
piecewise linear function in whose value is at and zero
at every , i.e.,

for ; this basis is a shifted and scaled tent
function. For the two-dimensional case, we choose again one
basis function per vertex of the triangulation of the
planar region . The function is the unique function
of whose value is at and zero at every .
Depending on the author, the word "element" in "finite element
method" refers either to the triangles in the domain, the
piecewise linear basis function, or both. So for instance, an
author interested in curved domains might replace the triangles
with curved primitives, and so might describe the elements as
being curvilinear. On the other hand, some authors replace
"piecewise linear" by "piecewise quadratic" or even "piecewise
polynomial". The author might then say "higher order element"
instead of "higher degree polynomial". Finite element method is
not restricted to triangles (or tetrahedra in 3-d, or higher order
simplexes in multidimensional spaces), but can be defined on
quadrilateral subdomains (hexahedra, prisms, or pyramids in 3-
d, and so on). Higher order shapes (curvilinear elements) can
be defined with polynomial and even non-polynomial shapes
(e.g. ellipse or circle).
Examples of methods that use higher degree piecewise
polynomial basis functions are the hp-FEM and spectral FEM.
More advanced implementations (adaptive finite element
methods) utilize a method to assess the quality of the results
(based on error estimation theory) and modify the mesh during
the solution aiming to achieve approximate solution within
some bounds from the 'exact' solution of the continuum
problem. Mesh adaptivity may utilize various techniques, the
most popular are:
moving nodes (r-adaptivity)
refining (and unrefining) elements (h-adaptivity)
changing order of base functions (p-adaptivity)
combinations of the above (hp-adaptivity).
Small support of the basis[edit]


Solving the two-dimensional problem in the
disk centered at the origin and radius 1, with zero boundary
conditions.
(a) The triangulation.


(b) The sparse matrix L of the discretized linear system


(c) The computed solution,
The primary advantage of this choice of basis is that the inner
products

and

will be zero for almost all . (The matrix
containing in the location is known as
the Gramian matrix.) In the one dimensional case,
the support of is the interval . Hence,
the integrands of and are
identically zero whenever .
Similarly, in the planar case, if and do not share
an edge of the triangulation, then the integrals

and

are both zero.
Matrix form of the problem[edit]
If we
write and
then problem (3),
taking for ,
becomes
for (4)
If we denote by and the column
vectors and
, and if we let

and

be matrices whose entries are

and

then we may rephrase (4)
as
(5)
It is not necessary to
assume
. For a general
function ,
problem (3)
with
for
becomes actually
simpler, since no
matrix is used,
, (6)
where
and
for
.
As we have
discussed before,
most of the
entries
of and are
zero because the
basis
functions hav
e small support.
So we now have
to solve a linear
system in the
unknown wher
e most of the
entries of the
matrix , which
we need to invert,
are zero.
Such matrices
are known
as sparse
matrices, and
there are efficient
solvers for such
problems (much
more efficient
than actually
inverting the
matrix.) In
addition, is
symmetric and
positive definite,
so a technique
such as
the conjugate
gradient
method is
favored. For
problems that are
not too large,
sparse LU
decompositions a
nd Cholesky
decompositions s
till work well. For
instance, MATLA
B's backslash
operator (which
uses sparse LU,
sparse Cholesky,
and other
factorization
methods) can be
sufficient for
meshes with a
hundred
thousand
vertices.
The matrix is
usually referred
to as the stiffness
matrix, while the
matrix is
dubbed the mass
matrix.
General
form of the
finite
element
method[edit]
In general, the
finite element
method is
characterized by
the following
process.
One chooses
a grid for .
In the
preceding
treatment, the
grid consisted
of triangles,
but one can
also use
squares or
curvilinear
polygons.
Then, one
chooses
basis
functions. In
our
discussion,
we used
piecewise
linear basis
functions, but
it is also
common to
use
piecewise
polynomial
basis
functions.
A separate
consideration is
the smoothness
of the basis
functions. For
second
order elliptic
boundary value
problems,
piecewise
polynomial basis
function that are
merely
continuous
suffice (i.e., the
derivatives are
discontinuous.)
For higher order
partial differential
equations, one
must use
smoother basis
functions. For
instance, for a
fourth order
problem such
as
, one may use
piecewise
quadratic basis
functions that
are .
Another
consideration is
the relation of the
finite-dimensional
space to its
infinite-
dimensional
counterpart, in
the examples
above .
A conforming
element
method is one in
which the
space is a
subspace of the
element space for
the continuous
problem. The
example above is
such a method. If
this condition is
not satisfied, we
obtain
anonconforming
element method,
an example of
which is the
space of
piecewise linear
functions over the
mesh which are
continuous at
each edge
midpoint. Since
these functions
are in general
discontinuous
along the edges,
this finite-
dimensional
space is not a
subspace of the
original .
Typically, one
has an algorithm
for taking a given
mesh and
subdividing it. If
the main method
for increasing
precision is to
subdivide the
mesh, one has
an h-method (h is
customarily the
diameter of the
largest element in
the mesh.) In this
manner, if one
shows that the
error with a
grid is bounded
above by ,
for
some a
nd , then
one has an
order p method.
Under certain
hypotheses (for
instance, if the
domain is
convex), a
piecewise
polynomial of
order method
will have an error
of
order
.
If instead of
making h smaller,
one increases the
degree of the
polynomials used
in the basis
function, one has
a p-method. If
one combines
these two
refinement types,
one obtains
an hp-method
(hp-FEM). In the
hp-FEM, the
polynomial
degrees can vary
from element to
element. High
order methods
with large
uniform p are
called spectral
finite element
methods (SFEM).
These are not to
be confused
with spectral
methods.
For vector partial
differential
equations, the
basis functions
may take values
in .
Various
types of
finite
element
methods[e
dit]
AEM[edit]
The Applied
Element Method,
or AEM combines
features of both
FEM
and Discrete
element method,
or (DEM).
Main
article: Applied
element method
Generalized
finite
element
method[edit]
The Generalized
Finite Element
Method (GFEM)
uses local spaces
consisting of
functions, not
necessarily
polynomials, that
reflect the
available
information on
the unknown
solution and thus
ensure good local
approximation.
Then a partition
of unity is used to
bond these
spaces together
to form the
approximating
subspace. The
effectiveness of
GFEM has been
shown when
applied to
problems with
domains having
complicated
boundaries,
problems with
micro-scales, and
problems with
boundary
layers.
[7]

Mixed finite
element
method[edit]
Main
article: Mixed
finite element
method
hp-FEM[edit]
The hp-
FEM combines
adaptively,
elements with
variable
size h and
polynomial
degree p in order
to achieve
exceptionally fast,
exponential
convergence
rates.
[8]

hpk-
FEM[edit]
The hpk-
FEM combines
adaptively,
elements with
variable size h,
polynomial
degree of the
local
approximations p
and global
differentiability of
the local
approximations (k
-1) in order to
achieve best
convergence
rates.
XFEM[edit]
Main
article: Extended
finite element
method
S-FEM[edit]
Main
article: Smoothed
finite element
method
Fiber Beam
Method[edit]
Main article: Fiber
beam element
Spectral
methods[edit
]
Main
article: Spectral
method
Meshfree
methods[edit
]
Main
article: Meshfree
methods
Discontinuo
us Galerkin
methods[edit
]
Main
article: Discontinu
ous Galerkin
method
Finite
element
limit
analysis[edit]
Main
article: Finite
element limit
analysis
Stretched
grid
method[edit]
Main
article: Stretched
grid method
Comparis
on to the
finite
difference
method[e
dit]

T
h
i
s

s
e
c
t
i
o
n

d
o
e
s

n
o
t

c
i
t
e

a
n
y

r
e
f
e
r
e
n
c
e
s

o
r

s
o
u
r
c
e
s
.

P
l
e
a
s
e

h
e
l
p

i
m
p
r
o
v
e

t
h
i
s

s
e
c
t
i
o
n

b
y

a
d
d
i
n
g

c
i
t
a
t
i
o
n
s

t
o

r
e
l
i
a
b
l
e

s
o
u
r
c
e
s
.

U
n
s
o
u
r
c
e
d

m
a
t
e
r
i
a
l

m
a
y

b
e

c
h
a
l
l
e
n
g
e
d

a
n
d

r
e
m
o
v
e
d
.

(
N
o
v
e
m
b
e
r

2
0
1
0
)
The finite
difference
method (FDM) is
an alternative
way of
approximating
solutions of
PDEs. The
differences
between FEM
and FDM are:
The most
attractive
feature of the
FEM is its
ability to
handle
complicated
geometries
(and
boundaries)
with relative
ease. While
FDM in its
basic form is
restricted to
handle
rectangular
shapes and
simple
alterations
thereof, the
handling of
geometries in
FEM is
theoretically
straightforwar
d.
The most
attractive
feature of
finite
differences is
that it can be
very easy to
implement.
There are
several ways
one could
consider the
FDM a
special case
of the FEM
approach.
E.g., first
order FEM is
identical to
FDM for
Poisson's
equation, if
the problem
is discretized
by a regular
rectangular
mesh with
each
rectangle
divided into
two triangles.
There are
reasons to
consider the
mathematical
foundation of
the finite
element
approximatio
n more
sound, for
instance,
because the
quality of the
approximatio
n between
grid points is
poor in FDM.
The quality of
a FEM
approximatio
n is often
higher than in
the
correspondin
g FDM
approach, but
this is
extremely
problem-
dependent
and several
examples to
the contrary
can be
provided.
Generally, FEM is
the method of
choice in all types
of analysis in
structural
mechanics (i.e.
solving for
deformation and
stresses in solid
bodies or
dynamics of
structures)
whilecomputation
al fluid
dynamics (CFD)
tends to use FDM
or other methods
like finite volume
method (FVM).
CFD problems
usually require
discretization of
the problem into
a large number of
cells/gridpoints
(millions and
more), therefore
cost of the
solution favors
simpler, lower
order
approximation
within each cell.
This is especially
true for 'external
flow' problems,
like air flow
around the car or
airplane, or
weather
simulation.
Applicatio
n[edit]


Visualization of
how a car
deforms in an
asymmetrical
crash using
finite element
analysis.[1]
A variety of
specializations
under the
umbrella of the
mechanical
engineering
discipline (such
as aeronautical,
biomechanical,
and automotive
industries)
commonly use
integrated FEM in
design and
development of
their products.
Several modern
FEM packages
include specific
components such
as thermal,
electromagnetic,
fluid, and
structural working
environments. In
a structural
simulation, FEM
helps
tremendously in
producing
stiffness and
strength
visualizations and
also in minimizing
weight, materials,
and costs.
FEM allows
detailed
visualization of
where structures
bend or twist, and
indicates the
distribution of
stresses and
displacements.
FEM software
provides a wide
range of
simulation
options for
controlling the
complexity of
both modeling
and analysis of a
system. Similarly,
the desired level
of accuracy
required and
associated
computational
time
requirements can
be managed
simultaneously to
address most
engineering
applications. FEM
allows entire
designs to be
constructed,
refined, and
optimized before
the design is
manufactured.
This powerful
design tool has
significantly
improved both
the standard of
engineering
designs and the
methodology of
the design
process in many
industrial
applications.
[9]
Th
e introduction of
FEM has
substantially
decreased the
time to take
products from
concept to the
production
line.
[9]
It is
primarily through
improved initial
prototype designs
using FEM that
testing and
development
have been
accelerated.
[10]
In
summary,
benefits of FEM
include increased
accuracy,
enhanced design
and better insight
into critical design
parameters,
virtual
prototyping, fewer
hardware
prototypes, a
faster and less
expensive design
cycle, increased
productivity, and
increased
revenue.
[9]

FEA has also
been proposed to
use in stochastic
modelling, for
numerically
solving probability
models. See the
references
list.
[11][12]