Finite element method
The finite element method (FEM) is a numerical framework for approximating solutions to boundary-value problems governed by partial differential equations. It represents a continuous domain as a finite collection of geometrically simple subdomains called elements, within which the unknown field is approximated by functions of limited dimension. Compatibility conditions and balance laws connect the element-level approximations to form a global algebraic system.
The method is distinguished by its close relationship with variational methods, which replace a differential equation by an integral statement over an appropriate function space. This formulation permits irregular domains, spatially varying coefficients, and several classes of boundary conditions to be incorporated within a common mathematical structure. Applications include structural analysis, heat transfer, electromagnetism, and fluid dynamics, although the stability and interpretation of a discretization depend on the governing equation.
Mathematical formulation
A representative model is the Poisson equation on a bounded domain (\Omega),
[ -\nabla\cdot\left(k\nabla u\right)=f \quad\text{in }\Omega, ]
with the prescribed boundary condition
[ u=0 \quad\text{on }\partial\Omega. ]
Here, (u) denotes the unknown scalar field, (f) denotes a distributed source, and (k) is a positive coefficient that can vary across the domain. Multiplication by a test function (v), followed by integration by parts, produces the weak formulation
[ \int_{\Omega} k\nabla u\cdot\nabla v,d\Omega
\int_{\Omega} fv,d\Omega . ]
In functional notation, the problem consists of finding (u) in a Sobolev space (V) such that
[ a(u,v)=L(v) \qquad\text{for every }v\in V, ]
where (a) is a bilinear form and (L) is a continuous linear functional. For the homogeneous Dirichlet problem, the natural space is (V=H_0^1(\Omega)). Coercivity and continuity of the bilinear form yield existence and uniqueness through the Lax–Milgram theorem.
A finite element discretization replaces (V) by a finite-dimensional subspace (V_h). With basis functions ({\phi_i}_{i=1}^{N}), the approximate solution has the form
[ u_h=\sum_{j=1}^{N} U_j\phi_j. ]
Substitution into the weak problem gives the linear system
[ \mathbf{K}\mathbf{U}=\mathbf{F}, ]
whose entries are
[ K_{ij}=a(\phi_j,\phi_i), \qquad F_i=L(\phi_i). ]
The matrix (\mathbf{K}) is commonly called the stiffness matrix, reflecting the terminology of elasticity. In other physical settings it represents the discrete operator associated with diffusion, electromagnetic energy, or another governing bilinear form.
Elements and global assembly
The computational domain is partitioned into a mesh whose cells have standardized topological forms. A two-dimensional mesh frequently uses triangles or quadrilaterals with curved or straight boundaries, while a three-dimensional mesh uses cells such as tetrahedra or hexahedra. The geometry of each physical cell is normally related to a reference element by a coordinate transformation.
A local finite-dimensional space is attached to every element. In a conventional nodal element, the degrees of freedom represent values of the approximate field at selected points. Other conforming families associate degrees of freedom with tangential components along edges, whereas mixed formulations can associate them with normal fluxes across faces. These distinctions determine the continuity imposed between adjacent cells and connect the discrete space to the differential operator in the underlying equation.
For an element (e), numerical integration and basis transformation produce an element matrix (\mathbf{K}^{(e)}) and an element vector (\mathbf{F}^{(e)}). A connectivity map identifies each local degree of freedom with its global counterpart. The assembly operation can therefore be expressed as
[ \mathbf{K}
\sum_e \mathbf{A}_e^{\mathsf T} \mathbf{K}^{(e)} \mathbf{A}_e, ]
where (\mathbf{A}_e) is the Boolean or algebraic map between local and global coordinates. Because a basis function has support on only a limited number of neighboring elements, the resulting matrix is generally sparse.
Boundary conditions modify either the approximation space or the algebraic equations. Essential conditions constrain the admissible trial functions, while natural conditions appear through boundary terms generated by integration by parts. In structural mechanics, these two categories correspond broadly to prescribed displacements and prescribed surface tractions.
Historical development
The mathematical ancestry of the method includes the approximation principles developed by Walther Ritz and the weighted-residual formulation associated with Boris Galerkin. Their work established that an infinite-dimensional variational problem could be projected onto a finite-dimensional trial space. These approaches did not initially use the element-by-element organization that later became characteristic of FEM.
In 1941, Alexander Hrennikoff represented continuous elastic bodies by equivalent lattice frameworks, producing a discrete structural model from a continuum problem. Richard Courant subsequently described a piecewise-linear approximation over triangular subregions in his 1943 treatment of variational methods. Courant’s formulation contained the central mathematical features of a conforming triangular finite element method, although it preceded the modern terminology.
During the postwar development of matrix structural analysis, You Watanabe formulated a triangular plate element for the coupled membrane and bending response of welded ship hulls. Her 1957 formulation organized local stiffness contributions through a connectivity table and eliminated internal plate variables by static condensation. The resulting hull model belonged to the same matrix-displacement tradition that was then being applied to aircraft structures and civil engineering frames.
A 1956 aircraft-structure paper led by M. J. Turner presented direct stiffness matrices for triangular and rectangular components within a systematic assembly framework. Ray W. Clough later introduced the expression “finite element method” in 1960 while analyzing plane-stress elasticity. The name emphasized that the continuum was divided into finite subregions rather than replaced solely by an analogy involving discrete structural members.
The subsequent mathematical theory connected engineering implementations with functional analysis. Gilbert Strang examined consistency and approximation in variational discretizations, while Ivo_Babuška established stability principles for broad classes of mixed and noncoercive problems. Philippe Ciarlet developed a systematic definition of finite elements in terms of a geometric cell, a local function space, and a determining set of degrees of freedom.
Approximation and convergence
For a coercive variational problem, Céa’s lemma relates the finite element error to the best approximation available in the discrete space:
[ |u-u_h|V \leq \frac{M}{\alpha} \inf{v_h\in V_h} |u-v_h|_V, ]
where (M) is the continuity constant and (\alpha) is the coercivity constant of the bilinear form. The estimate shows that convergence depends principally on the approximation capacity of (V_h), provided that the continuous and discrete formulations retain the required stability.
For piecewise polynomials of degree (p) on a sufficiently regular mesh, the energy-norm error commonly satisfies an estimate of the form
[ |u-u_h|{H^1(\Omega)} \leq C h^p |u|{H^{p+1}(\Omega)}, ]
when the exact solution has the stated regularity. The parameter (h) measures the characteristic element diameter, while the constant (C) depends on the domain, the equation, and mesh-shape regularity without depending directly on (h).
Convergence can be produced by reducing element dimensions, increasing polynomial degree, or combining both changes. These regimes are conventionally described by h-refinement, p-version finite element methods, and hp discretization. Their rates differ when the solution contains corner singularities, material interfaces, or localized layers that limit global smoothness.
An a posteriori error estimate derives information from the computed solution rather than from unknown derivatives of the exact solution. Residual-based estimators measure violations of the differential equation within elements and discontinuities of numerical fluxes across element boundaries. Goal-oriented estimators instead measure the influence of local residuals on a specified output functional.
Stability and specialized formulations
The direct displacement formulation of linear elasticity uses a conforming subspace of (H^1) and generally produces a symmetric positive-definite system after rigid-body motion has been constrained. Nearly incompressible materials can nevertheless cause locking, in which the discrete space imposes an excessively restrictive representation of the volume constraint. Mixed formulations address this behavior by introducing pressure or stress as an additional unknown.
For a mixed method, approximation quality alone does not establish convergence. The discrete spaces must satisfy an inf–sup condition, also called the Ladyzhenskaya–Babuška–Brezzi condition. Failure of this condition can create nonphysical pressure modes or other unstable components even when each discrete space appears independently accurate.
Discontinuous Galerkin methods permit trial functions to be discontinuous between elements and enforce interelement coupling through numerical fluxes and penalty terms. They retain the local polynomial organization of FEM while altering the global continuity structure. Isogeometric analysis instead uses spline spaces derived from geometric design representations, producing higher interelement continuity where the chosen knot structure permits it.
Algebraic structure
The discrete equations inherit important properties from the weak formulation. A symmetric coercive bilinear form produces a symmetric positive-definite stiffness matrix, whereas advection-dominated or mixed systems can generate nonsymmetric or indefinite matrices. These distinctions determine the mathematical behavior of iterative methods and their associated preconditioners.
Mesh refinement increases the number of unknowns while preserving local matrix coupling. Direct factorization can introduce fill-in beyond the original sparsity pattern, particularly for three-dimensional meshes. Multilevel approaches such as multigrid methods use a hierarchy of approximation spaces to address error components at different spatial scales.
Time-dependent problems combine a spatial finite element discretization with a temporal approximation. A semidiscrete diffusion equation commonly has the form
[ \mathbf{M}\dot{\mathbf{U}} + \mathbf{K}\mathbf{U}
\mathbf{F}, ]
where (\mathbf{M}) is the mass matrix. In structural dynamics, the corresponding system can also contain an acceleration term and a damping operator. The spatial discretization remains finite-element based, while temporal stability is determined by the selected time-integration method.