Finite element method

The finite element method replaces an infinite-dimensional weak formulation by finitely many unknowns attached to a finite element mesh. Simple local functions, often piecewise polynomial functions, join to form a global approximation. The attraction is that each finite element can be integrated separately, while the assembled equations enforce a global balance. A piecewise-linear hat function only meets its immediate neighbors, so a large problem can have a sparse matrix rather than a dense one. Gilbert Strang's MIT lecture develops this connection from the one-dimensional weak formulation through local hats to the stiffness matrix and load integrals.^[ref]

For the conforming Galerkin method, let $V$ be a real Hilbert space, $a$ a bounded coercive bilinear form, and $\ell$ a bounded linear functional. Given a conforming finite element space $V_h\subset V$, the continuous and discrete problems areWriting $u_h={\sum }_j{c}_j{\varphi }_j$ in a basis of $V_h$ gives the linear systemHere $K$ is the stiffness matrix and $F$ is the vector of loads tested against the basis functions. The Galerkin method uses the same trial and test space in this construction. Finite element methods also use the Petrov-Galerkin method, with different trial and test spaces; testing against the trial space is therefore a property of this construction, rather than a definition of every finite element method.

The finite element solution of the one-dimensional Poisson equation below gives every entry explicitly. Refining the finite element mesh then improves the approximation by adding local degrees of freedom. For a symmetric coercive bilinear form, Galerkin orthogonality makes the discrete solution the best approximation in the energy norm; the geometry of orthogonal projection explains why this local construction yields a globally meaningful answer.

Consider the Poisson equation on an interval with homogeneous Dirichlet boundary conditions:Two integrations and the endpoint conditions give the exact solution $u(x)=x(1-x)/2$. To approximate it, multiply by a test function $v$ that vanishes at both endpoints. Integration by parts givesThe boundary term is zero. Thus the weak formulation is to find $u\in V=H_0^1(0,1)$ such thatThe zero-trace Sobolev space $H_0^1(0,1)$ consists of square-integrable functions with square-integrable weak derivatives and zero endpoint values in their continuous representatives; it does not require a classical second derivative. This is why continuous piecewise linear functions are admissible despite corners. The Poincaré inequality makes $a$ a coercive bilinear form on $V$. Strang derives this one-dimensional weak formulation, the piecewise-linear hat functions, and their load and derivative-product integrals in the cited lecture; its endpoint example uses a free left endpoint, while both endpoints here are fixed.^[ref]

Take the uniform finite element mesh $x_j=jh$, $j=0,\dots ,n$, with $h=1/n$. The continuous $P_1$ trial space isHere $P_1$ means polynomials of degree at most one on each interval. This is a conforming finite element space in $V$, with $n-1$ interior piecewise-linear hat functions ${\varphi }_i$ as a basis. Their endpoint values enforce the Dirichlet boundary conditions automatically. Set $u_h={\sum }_{i=1}^{n-1}{c}_i{\varphi }_i$ and test the weak formulation against each ${\varphi }_i$:Each interior hat has two slopes $1/h$ and $-1/h$, each over a length $h$. Its squared derivative therefore integrates to $2/h$. Adjacent hats have opposite slopes on their one shared interval, giving $-1/h$. All other derivative products vanish. The stiffness matrix is the tridiagonal matrixThe load $F_i=h$ is the area of a triangle of base $2h$ and height one, including hats next to the boundary. Boundary hats are omitted because the boundary values are prescribed. The piecewise-linear hat function article explains how the same entries arise by adding local two-node contributions.

With four mesh intervals, $h=1/4$, there are three unknowns:For example, the first row gives $8(3/32)-4(1/8)=1/4$ and the middle gives $-4(3/32)+8(1/8)-4(3/32)=1/4$. The coefficients are the nodal values $u_h(x_i)=c_i$. At $x=3/8$, linear interpolation gives $u_h(3/8)=7/64=0.109375$, whereas $u(3/8)=15/128=0.1171875$: exactness at the nodes still leaves a visible error between them.

In this model the Galerkin method solution equals the nodal approximation $I_{h}u$ obtained by linear interpolation on any finite element mesh. Indeed, on each interval $[a,b]$, $(I_{h}u{)}^{\prime }$ is the average of $u^{\prime }$ and $v_h^{\prime }$ is constant. Therefore ${\int }_a^b(u^{\prime }-(I_{h}u{)}^{\prime })v_h^{\prime }dx=0$. Summing proves Galerkin orthogonality for $I_{h}u$, and uniqueness gives $u_h=I_{h}u$. This argument relies on the constant diffusion coefficient, conformity and exact integration; nodal exactness is not a general promise for finite element methods with other operators or integration rules.

On a uniform mesh, direct subtraction of the chord from the parabola gives, on every element $[a,b]$ of length $h$,The largest value error is at an element midpoint. Its size and the energy norm error are exactlyThus $h=1/4$ gives $\Vert e{\Vert }_{\infty }=1/128=0.0078125$ and $\Vert e{\Vert }_a=1/(4\sqrt{12})\approx 0.0721688$. With eight elements, $h=1/8$, these become $1/512=0.001953125$ and $1/(8\sqrt{12})\approx 0.0360844$. Halving the mesh spacing quarters the maximum value error but only halves the energy norm error: values and slopes converge at different rates. The figure shows the small value error at its own scale, where the four-to-one reduction becomes clear.

The spatial Galerkin method for the heat equation gives a mass matrix $M$ and a diffusion stiffness matrix $S$. With zero forcing, $\frac{d}{dt}(c^{T}Mc)=-2c^{T}Sc\le 0$. The quadratic form $c^{T}Mc$ is the squared L2 norm of the represented finite-element function, proving an energy contraction. For nonzero forcing, the corresponding identity has the additional term $2c^{T}F$, which can be estimated using the Cauchy-Schwarz inequality.

For a uniform piecewise-linear Galerkin heat discretization, the consistent mass matrix is $M=(h/6)\mathrm{tridiag}(1,4,1)$ and the stiffness matrix is $K=h^{-1}\mathrm{tridiag}(-1,2,-1)$. The generalized eigenvalues are $6(1-\cos \theta )/[h^2(2+\cos \theta )]$. Applying Forward Euler method requires $k\lambda \le 2$, giving the displayed uniform bound for $\mu =k/h^2$. For $J$ intervals with Dirichlet endpoints the exact finite-grid upper bound is $(2-\cos (\pi /J))/[3(1+\cos (\pi /J))]$, which approaches $1/6$. Replacing the consistent mass matrix by a lumped one changes the stability restriction.

A finite element specifies a reference region, a finite-dimensional space of local trial functions, and degrees of freedom that uniquely determine each such function. Assembling compatible copies over a finite element mesh gives a global trial space. For a one-dimensional linear element, the region is an interval, the functions are affine functions, and the degrees of freedom are the endpoint values.

A conforming finite element space is a finite-dimensional vector subspace of the solution's variational space $V$. For a second-order problem with homogeneous Dirichlet boundary conditions, continuous piecewise-polynomial functions vanishing on the boundary form such a subspace of the zero-boundary Sobolev space $H_0^1(\Omega )$. Conformity allows exact Galerkin orthogonality and the Céa lemma.

A finite element mesh partitions a domain into simple regions called elements, such as intervals, triangles or tetrahedra. In a conforming partition, distinct element interiors are disjoint, their closures cover the domain closure, and any nonempty intersection is a common face, edge or vertex of both elements. Chen's author notes give this compatibility condition for simplicial meshes.^[ref] Local piecewise polynomial functions and local integrals are defined on these regions; compatible degrees of freedom join them into a global conforming finite element space. The assembled stiffness matrix and mass matrix then inherit the local overlap pattern.

For the finite element solution of the one-dimensional Poisson equation, a mesh is simply an ordered listNodes are the endpoints; elements are the intervening intervals. Four equal elements have five nodes but only three unknown nodal values after imposing the zero Dirichlet boundary conditions. A nonuniform partition such as $0,1/8,1/4,1/2,1$ has element lengths $1/8,1/8,1/4,1/2$ and $h=1/2$, so counting elements alone does not measure the largest unresolved scale.

Refining a finite element mesh replaces some or all elements by smaller ones. For continuous piecewise linear functions in one dimension, retaining all old nodes gives nested conforming finite element spaces $V_h\subset V_{h/2}$ when every interval is bisected. The old approximation remains available, and the new interior nodes allow additional changes of slope. With a symmetric coercive bilinear form and exact integration, best approximation from Galerkin orthogonality means that the energy norm error cannot increase under such nested refinement.

For the finite element solution of the one-dimensional Poisson equation, uniform refinement from four to eight elements increases the unknown count from three to seven. The exact formulas give maximum errors $0.0078125$ and $0.001953125$, and energy norm errors approximately $0.0721688$ and $0.0360844$. Refinement near a region of strong curvature can spend unknowns where the linear interpolation error matters most; convergence in a general problem still needs approximation and regularity assumptions such as those in the finite element interpolation estimate.

A family of simplicial meshes is shape-regular when the ratio of each element's diameter to its inscribed-ball radius stays uniformly bounded. This excludes degenerating shapes and makes finite element interpolation estimates uniform in the mesh size. The element sizes may still vary substantially across the domain.

On a shape-regular mesh in dimensions at most three, nodal piecewise-linear interpolation of $u\in H^2(\Omega )$ satisfieswhere $h$ is the largest element diameter and $C$ is independent of $h$. For homogeneous Dirichlet boundary conditions, $I_{h}u$ belongs to the conforming finite element space. Scaling each element to a reference element gives the estimate; shape regularity bounds the scaling constants, while subtraction of an affine function leaves an error controlled by the second derivatives.

The mass matrix is the Gram matrix of a finite element basis in the L2 inner product. Its coefficient quadratic form is the squared L2 norm of the represented trial function, so it is a positive-definite matrix for a linearly independent basis.

A mass eigenstate diagonalizes the Hermitian mass operator or mass matrix and has a definite mass. In a flavor basis with off-diagonal mass terms, it is a linear combination of flavor states. For unstable particles, decay propagation instead involves an effective complex matrix; its eigenvectors need not inherit the orthogonality of a Hermitian mass problem.

On a one-dimensional uniform finite element mesh of spacing $h$, let $w$ be an affine function and ${\varphi }_i$ the interior piecewise-linear hat functions. The weighted mass matrix has entriesand zero entries for nonadjacent indices. The diagonal product ${\varphi }_i^2$ is symmetric about $x_i$ and integrates to $2h/3$; the neighbouring product is symmetric about $(x_i+x_{i+1})/2$ and integrates to $h/6$. The odd part of an affine function integrates to zero around each centre, giving the formulas. This evaluates a variable reaction term exactly without substituting a constant coefficient or a quadrature approximation.

The stiffness matrix represents a variational bilinear form $a$ in a finite element basis. A symmetric coercive bilinear form gives a positive-definite matrix. Local support of the piecewise-linear hat functions produces a sparse matrix, allowing efficient numerical solution.

For $-u^{\prime \prime }+u=f$ with homogeneous Dirichlet boundary conditions, the conforming piecewise-linear hat functions on a uniform grid give a diffusion stiffness matrix $S$ and a mass matrix $M$. Their nonzero entries are $S_{ii}=2/h$, $S_{i,i+1}=-1/h$, $M_{ii}=2h/3$, and $M_{i,i+1}=h/6$, with symmetric counterparts. The sum is a positive-definite matrix because its coefficient quadratic form is $\int ((u_h^{\prime }{)}^2+u_h^2)$. For forcing $f(x)=x$, the load at $x_i$ is $h{x}_i$.

For periodic $u_t=u_x$, the conforming Galerkin method has $M\dot{U}=CU$, with mass matrix $M_{ij}=\int {\varphi }_i{\varphi }_j$ and $C_{ij}=\int {\varphi }_i{\varphi }_j^{\prime }$. Integration by parts and periodicity give $C^T=-C$. ThereforeThe conserved quantity is the squared L2 norm of the finite element function. On a uniform one-dimensional mesh of spacing $h$, the piecewise-linear hat functions give $M_{ii}=2h/3$, $M_{i,i\pm 1}=h/6$, $C_{i,i+1}=1/2$ and $C_{i,i-1}=-1/2$.

A cubic Hermite finite element is piecewise cubic and joins with continuous first derivative. It is therefore conforming for fourth-order variational problems whose energy space lies in $H^2$.

The Ritz method minimizes a quadratic energy over a finite-dimensional trial space. For basis functions ${\varphi }_j$ and bilinear form $a$, its linear system is $Ka=F$ with $K_{jk}=a({\varphi }_k,{\varphi }_j)$ and $F_j=⟨f,{\varphi }_j⟩$.

For a self-adjoint Sturm-Liouville operator $L_{ϵ}=L_0+ϵL_1$ with common boundary conditions and an $L^2$-normalized lowest eigenfunction $U_0$ of $L_0$, evaluating the Rayleigh quotient on $U_0$ gives the displayed upper bound. More generally the lowest eigenvalue is the infimum of affine functions of $ϵ$, so it is concave in that parameter. The trial state must satisfy the boundary requirements of the variational domain; changing the operator without changing the trial function's normalization can still be legitimate when the weight is fixed.

For $-(p{y}^{\prime }{)}^{\prime }=\lambda wy$ on $[0,1]$ with $p=(1+x)e^{2x}$, $w=e^{2x}$, $y(0)=0$ and $y^{\prime }(1)=0$, the Rayleigh quotient is $\int p{y}^{\prime 2}/\int w{y}^2$. The trial $y_1=x{e}^{-x}$ gives numerator $5/12$ and denominator $1/3$, hence $0<{\lambda }_{\min }\le 5/4$. Every displayed $y_n$ obeys the same mixed boundary conditions; so does $\sin [(m+1/2)\pi x]$. Adding any of these to the trial space cannot worsen the optimized upper bound. Matrix elements $A_{ij}=\int p{y}_i^{\prime }{y}_j^{\prime }$ and $M_{ij}=\int w{y}_i{y}_j$ give a finite-dimensional generalized eigenvalue problem $Av=\lambda Mv$. The case $n=1$ repeats the original trial and adds no variational freedom.

For a bounded symmetric coercive bilinear form on a Hilbert space, a solution of $a(u,w)=\ell (w)$ minimizes $J(v)=a(v,v)/2-\ell (v)$. The displayed identity proves strict minimality and uniqueness. Restricting to a conforming finite element space makes the discrete Ritz method equivalent to the Galerkin method. Galerkin orthogonality gives exact best approximation in the energy norm; the Céa lemma translates it into a reference-norm estimate.

A continuous symmetric coercive bilinear form $a$ on a real Hilbert space and a continuous linear functional $\ell$ give a unique minimum of $I$. The parallelogram identity makes any minimizing sequence Cauchy in the energy norm; coercivity and completeness give its limit. Differentiation gives $a(u,w)=\ell (w)$, and $I(v)-I(u)=a(v-u,v-u)$ proves uniqueness. A nonsymmetric operator instead differentiates to its symmetric part. Mere pointwise positivity without a lower bound in the chosen norm need not yield a minimizing vector.

A symmetric coercive bilinear form $a$ defines the energy norm $\Vert v{\Vert }_a=\sqrt{a(v,v)}$. In this norm the Ritz method is an orthogonal projection onto its finite-dimensional trial space.

Let $V_h\subset V$ be a conforming finite element space. Suppose the exact solution satisfies $a(u,v)=\ell (v)$ for every $v\in V$ and its Galerkin method approximation satisfies $a(u_h,v_h)=\ell (v_h)$ for every $v_h\in V_h$, using the same bilinear form and linear functional. Subtracting the two identities gives Galerkin orthogonality:This identity itself does not require symmetry. Its interpretation as an orthogonal projection in an energy norm does require a bounded symmetric coercive bilinear form: then $a$ is an inner product, $e=u-u_h$ is perpendicular to every available discrete correction, andConsequently $u_h$ is the unique closest point in $V_h$ in the energy norm. This Pythagorean theorem in a function space explains the beauty of the Galerkin method: integrated local equations compute a global best approximation. The Céa lemma supplies the corresponding estimate in a chosen reference norm.

For the finite element solution of the one-dimensional Poisson equation, $a(w,v)={\int }_0^1{w}^{\prime }{v}^{\prime }dx$ and $\Vert w{\Vert }_a=\Vert w^{\prime }{\Vert }_{L^2(0,1)}$. Testing the error against the $i$th piecewise-linear hat function on a uniform mesh givesIn this model $e$ vanishes at every node, so each of these integrals is individually zero, although $e$ is positive inside every element. Galerkin orthogonality controls the derivative-weighted balance; it does not say that the pointwise error vanishes.

The exact identity assumes conformity and exact evaluation of the continuous forms, including any numerical integration that is used. If assembly instead uses $a_h$ and ${\ell }_h$, subtraction gives the consistency defectThus inexact numerical integration need not preserve the projection. In the Petrov-Galerkin method the error is tested against a different test space, and for a nonsymmetric bilinear form the displayed identity need not describe perpendicularity in any energy norm.

For a continuous coercive bilinear form with continuity constant $M$ and coercivity constant $\alpha$, a conforming finite-element solution satisfiesFor a symmetric problem measured in its energy norm, the constant is one.

For any $v_h\in V_h$, Galerkin orthogonality gives $a(e,e)=a(e,u-v_h)$, where $e=u-u_h$. Coercivity and continuity yield $\alpha \Vert e{\Vert }_V^2\le M\Vert e{\Vert }_V\Vert u-v_h{\Vert }_V$. Divide by $\Vert e{\Vert }_V$ when it is nonzero and take the infimum. In the symmetric case, the Cauchy-Schwarz inequality in the energy norm gives the constant one.

For a symmetric elliptic variational problem with Galerkin orthogonality, let $e=u-u_h$ and solve the dual problem $a(v,z)=(e,v{)}_{L^2}$ for every $v\in H_0^1(\Omega )$. If elliptic regularity gives $\Vert z{\Vert }_{H^2}\le C\Vert e{\Vert }_{L^2}$, then a finite element interpolation estimate givesThus $\Vert e{\Vert }_{L^2}\le Ch\Vert e{\Vert }_{H^1}$. Combining this with the Céa lemma and an $O(h)$ energy norm error gives an $O(h^2)$ L2 norm error. The additional order depends on the dual regularity assumption, which can fail on unsuitable domains.

A piecewise-linear hat function is a continuous piecewise linear function that selects one node of a finite element mesh: it is one there and zero at every other node. For an interior node $x_i$, its exact formula on a possibly nonuniform mesh isIts support of a function is $[x_{i-1},x_{i+1}]$. It rises on the left interval and falls on the right, producing its hat-shaped graph. The interior hat functions form a basis of the continuous $P_1$ conforming finite element space with homogeneous Dirichlet boundary conditions. Chen's author notes define this nodal basis and illustrate it in one and two dimensions.^[ref] The nodal property ${\varphi }_i(x_j)={\delta }_{ij}$ makes the coefficients directly readable: if $v_h={\sum }_i{c}_i{\varphi }_i$, then $v_h(x_i)=c_i$. Unlike the global Lagrange cardinal polynomials, each of these functions is confined to two neighboring elements.

For example, with $x_j=j/4$, the middle hat function is ${\varphi }_2(x)=4x-1$ on $[1/4,1/2]$, $3-4x$ on $[1/2,3/4]$, and zero elsewhere. Hence ${\varphi }_2(3/8)=1/2$, while ${\varphi }_1(3/8)=1/2$ and ${\varphi }_3(3/8)=0$. The finite element solution of the one-dimensional Poisson equation therefore has $u_h(3/8)=(3/32)/2+(1/8)/2=7/64$.

Local assembly makes this locality computational. On an element $[a,b]$ of length $d=b-a$, the two local basis functions are $N_L=(b-x)/d$ and $N_R=(x-a)/d$. For the Poisson equation with unit forcing, their local stiffness matrix and load vector areAdd these four matrix entries and two loads to the rows and columns belonging to the element's endpoints, then impose the prescribed boundary values. An interior node receives contributions from its two neighboring elements. Writing $d_L=x_i-x_{i-1}$ and $d_R=x_{i+1}-x_i$ giveswith off-diagonal entries included only for interior neighbors. On a uniform mesh these reduce to $2/h$, $-1/h$ and $h$. Nonadjacent hat functions have no derivative overlap of positive length, so all distant entries vanish. This produces a sparse matrix with only three diagonals: the global coupling is built entirely from local two-node interactions.