The Finite Element Method

11 May 2026 · Download PDF · Download code

The aim of the Finite Element Method (FEM) is to solve (or more precisely approximate the solution of) a Partial Differential Equation (PDE) using a computer. In general the solution to a PDE is a function. Unfortunately it is computationally intractable to task a computer with finding this underlying function. The idea to resolve this issue is to parametrize a family of functions with the aim of computing the parameters through solving a linear system of equations. There is a plethora of ways to get to such a linear system of equations, with the most common being the so-called Galerkin approach.

1. Getting the computer to understand functions

In the most basic terms, computers are only capable of understanding numbers, which means that they are inherently unable to solve equations where the unknowns are functions. This is a problem when trying to solve differential equations. To solve this problem, the functions are instead approximated using parameterizations for which the coefficients (numbers) are unknown. Intuitively, these numbers are exactly what the computer is asked to find. For element methods, this parameterization is chosen to be the simplest possible: a linear combination of functions

𝑝 (𝐱) \approx 𝐓 (𝐱) 𝐩 = [\begin{matrix} 𝑇_{1} (𝐱) & 𝑇_{2} (𝐱) & \dots & 𝑇_{𝑛 (𝐱)} \end{matrix}] [\begin{matrix} 𝑝_{1} \\ 𝑝_{2} \\ ⋮ \\ 𝑝_{𝑛} \end{matrix}],

(1)

where $𝑝$ is the unknown function being approximated. Note that the linearity is with respect to the unknown parameters $𝐩$ , but not necessarily in the known basis functions $𝐓 (𝐱)$ .

One might ask: How does the above relate to the Finite Element Method? The answer is that the functions $𝑇_{𝑖}$ are chosen to be simpler functions with support equal to only a few subdomains of the original domain. These subdomains are commonly referred to as elements.

Interactive figure 1: Piecewise polynomial approximation

Interactive figure 1 (Piecewise polynomial approximation) In one dimension an approximation of the form (1) is easy to visualize. The graph shows the function $𝑝$ (dashed) and the approximation $𝐓 (𝑥) 𝐩$ (solid) on $𝑁$ elements (every other one shaded). The coefficients $𝑝_{𝑖} = 𝑝 (𝑥_{𝑖})$ are the large dots, and each basis function $𝑇_{𝑖}$ , scaled by its $𝑝_{𝑖}$ , is drawn faint. There are two ways to improve the approximation: more elements¹¹ Called $ℎ$ -refinement, as it decreases the element size $ℎ = 1 / 𝑁$ . or basis functions of a higher polynomial order $𝑝$ .²² Called $𝑝$ -refinement. Here $𝑝$ is the order, not the function being approximated.

Elements

𝑁

2. The Finite Element Method for the Helmholtz equation

In acoustics we are often interested in solving the Helmholtz equation, given as

Δ 𝑝 (𝐱) + 𝑘^{2} 𝑝 (𝐱) = 0, 𝐱 \in Ω,

(2)

where $𝑘$ is most often referred to as the wavenumber. In its current form the basis functions $𝐓 (𝐱)$ must be twice differentiable; piecewise linear functions, for example, have a vanishing second derivative inside each element and cannot be used. This constraint can be relieved using the so-called weak formulation. In order to do so we multiply the Helmholtz equation with a so-called test function $𝜑$ followed by an integration over the domain $Ω$

\int_{Ω} 𝜑 (𝐱) (Δ 𝑝 (𝐱) + 𝑘^{2} 𝑝 (𝐱)) d 𝐱 = 0 .

(3)

Utilizing integration by parts on the term including the Laplacian we can move derivatives from $𝑝 (𝐱)$ onto $𝜑 (𝐱)$ as

\int_{Ω} {(\nabla 𝜑 (𝐱))}^{⊤} \nabla 𝑝 (𝐱) d 𝐱 - 𝑘^{2} \int_{Ω} 𝜑 (𝐱) 𝑝 (𝐱) d 𝐱 - \int_{𝜕 Ω} 𝜑 (𝐱) \frac{𝜕 𝑝}{𝜕 𝑛} (𝐱) d 𝐱 = 0 .

(4)

In the above the constraint on the smoothness of $𝑝 (𝐱)$ has been reduced to first order. Furthermore the last term

\int_{𝜕 Ω} 𝜑 (𝐱) \frac{𝜕 𝑝}{𝜕 𝑛} (𝐱) d 𝐱

(5)

can be used to include the boundary conditions into the equation. For simplification purposes we will in the following assume that the above is equal to zero. As such we have

0 = \int_{Ω} {(\nabla 𝜑 (𝐱))}^{⊤} \nabla 𝑝 (𝐱) d 𝐱 - 𝑘^{2} \int_{Ω} 𝜑 (𝐱) 𝑝 (𝐱) d 𝐱 .

(6)

We will now use the Galerkin approach to discretize the above equation. Note that “Galerkin” simply refers to the approach where $𝜑 (𝐱)$ and $𝑝 (𝐱)$ are discretized using the same basis functions ( $𝐓 (𝐱)$ ). In short this means that we introduce

𝑝 (𝐱) \approx 𝐓 (𝐱) 𝐩, 𝜑 (𝐱) \approx 𝐚^{⊤} {𝐓 (𝐱)}^{⊤},

(7)

where $𝐚 \in ℂ^{𝑛}$ is a set of arbitrary coefficients that parametrize a whole family of test functions and $𝐩 \in ℂ^{𝑛}$ are the coefficients that we aim to find. Inserting the approximations it follows that

\begin{aligned} 0 & = \int_{Ω} {(\nabla 𝜑 (𝐱))}^{⊤} \nabla 𝑝 (𝐱) d 𝐱 - 𝑘^{2} \int_{Ω} 𝜑 (𝐱) 𝑝 (𝐱) d 𝐱 \\ \approx \int_{Ω} 𝐚^{⊤} {(\nabla 𝐓 (𝐱))}^{⊤} \nabla 𝐓 (𝐱) 𝐩 d 𝐱 - 𝑘^{2} \int_{Ω} 𝐚^{⊤} {𝐓 (𝐱)}^{⊤} 𝐓 (𝐱) 𝐩 d 𝐱 \\ = 𝐚^{⊤} (\int_{Ω} {(\nabla 𝐓 (𝐱))}^{⊤} \nabla 𝐓 (𝐱) d 𝐱) 𝐩 - 𝑘^{2} 𝐚^{⊤} (\int_{Ω} {𝐓 (𝐱)}^{⊤} 𝐓 (𝐱) d 𝐱) 𝐩 \\ = 𝐚^{⊤} (\int_{Ω} {(\nabla 𝐓 (𝐱))}^{⊤} \nabla 𝐓 (𝐱) d 𝐱 - 𝑘^{2} \int_{Ω} {𝐓 (𝐱)}^{⊤} 𝐓 (𝐱) d 𝐱) 𝐩 . \end{aligned}

(8)

We can make the above equality hold for all possible test functions, i.e. for all $𝐚 \in ℂ^{𝑛}$ , if we find $𝐩$ such that

(\int_{Ω} {(\nabla 𝐓 (𝐱))}^{⊤} \nabla 𝐓 (𝐱) d 𝐱 - 𝑘^{2} \int_{Ω} {𝐓 (𝐱)}^{⊤} 𝐓 (𝐱) d 𝐱) 𝐩 = 𝟎 .

(9)

This is exactly what we do in Finite Element computations. As a shorthand the above can be written as

(𝐊 - 𝑘^{2} 𝐌) 𝐩 = 𝟎,

(10)

where $𝐊$ is referred to as the stiffness matrix and $𝐌$ as the mass matrix. In the literature the equation can also be found as

(𝑐^{2} 𝐊 - 𝜔^{2} 𝐌) 𝐩 = 𝟎,

(11)

where $𝑐$ is the propagation speed and $𝜔$ is the angular frequency.

3. What is an element?

A key insight is that the element serves two purposes: it represents a subdomain of the original domain (also referred to as the geometry) while also describing parts of the unknown function(s) of interest. In general this subdomain is described by a parameterization, i.e. the element, as

𝐱^{𝑒} (𝐮) = 𝐗^{𝑒} 𝐍^{𝑒} (𝐮) \in Ω^{𝑒}, \forall 𝐮 \in {ℒ︀}^{𝑒},

(12)

where the superscript $𝑒$ denotes the element number, $𝐗^{𝑒}$ is a matrix with columns equal to the interpolation nodes of the geometry, $𝐍^{𝑒} (𝐮)$ are the so-called shape functions, $Ω^{𝑒}$ is the element in global coordinates and ${ℒ︀}^{𝑒}$ are the local coordinates. In addition to the geometric interpolation of each element, we need to further define interpolations of the unknown functions, which in acoustics is usually taken as the pressure $𝑝 (𝐱)$ . On element $𝑒$ this interpolation can be done as

𝑝 (𝐱^{𝑒} (𝐮)) = 𝐓 (𝐱^{𝑒} (𝐮)) 𝐩 = \underset{𝐓^{𝑒} (𝐮)}{\underset{⏟}{𝐓 (𝐱^{𝑒} (𝐮)) {(𝐋^{𝑒})}^{⊤}}} \underset{𝐩^{𝑒}}{\underset{⏟}{𝐋^{𝑒} 𝐩}} = 𝐓^{𝑒} (𝐮) 𝐩^{𝑒}, 𝐮 \in {ℒ︀}^{𝑒},

(13)

where $𝐋^{𝑒}$ is a permutation-like matrix that extracts the relevant values of $𝐩$ and the relevant functions of $𝐓 (𝐱)$ and orders them such that they correspond to the local values $𝐩^{𝑒}$ and the local basis functions $𝐓^{𝑒} (𝐮)$ .

Example 3: Element localization matrix

Example 3 (Element localization matrix) For a continuous linear element $𝑒$ all three corners correspond to a value of the global vector $𝐩$ . For example, the second element in the continuous plot has local corner values given by $𝐩^{2} = {[\begin{matrix} 𝑝_{3} & 𝑝_{1} & 𝑝_{4} \end{matrix}]}^{⊤}$ . This element would have $𝐋^{𝑒}$ given as

𝐋^{2} = [\begin{matrix} 0 & 0 & 1 & 0 & \dots & 0 \\ 1 & 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & 1 & \dots & 0 \end{matrix}],

(17)

so that $𝐩^{2} = 𝐋^{2} 𝐩$ . Note that $𝐋^{𝑒}$ is only an artifact of the mathematical description. Any reasonable implementation should use indexing instead of multiplication with $𝐋^{𝑒}$ .

4. Computing the stiffness and mass matrices

Using the element description of both the geometry and the pressure function the mass matrix can be computed as follows

\begin{aligned} 𝐌 & = \int_{Ω} {𝐓 (𝐱)}^{⊤} 𝐓 (𝐱) d 𝐱 \\ \approx \sum_{𝑒 = 1}^{𝑁} \int_{Ω^{𝑒}} 𝐋_{𝑒}^{⊤} 𝐋_{𝑒} {𝐓 (𝐱)}^{⊤} 𝐓 (𝐱) 𝐋_{𝑒}^{⊤} 𝐋_{𝑒} d 𝐱 \\ = \sum_{𝑒 = 1}^{𝑁} 𝐋_{𝑒}^{⊤} (\int_{{ℒ︀}^{𝑒}} 𝐓^{𝑒} {(𝐮)}^{⊤} 𝐓^{𝑒} (𝐮) jacobian (Ω^{𝑒}, 𝐮) d 𝐮) 𝐋_{𝑒} \\ \approx \sum_{𝑒 = 1}^{𝑁} 𝐋_{𝑒}^{⊤} (\underset{𝐌_{𝑒}}{\underset{⏟}{\sum_{𝑖 = 1}^{𝑄} 𝑤_{𝑖} 𝐓^{𝑒} {(𝐮_{𝑖})}^{⊤} 𝐓^{𝑒} (𝐮_{𝑖}) jacobian (Ω^{𝑒}, 𝐮_{𝑖})}}) 𝐋_{𝑒}, \end{aligned}

(18)

where the first approximation is due to the geometry being approximated by the elements, and the second is due to the quadrature with the $𝑄$ points $𝐮_{𝑖}$ and weights $𝑤_{𝑖}$ . Furthermore, $jacobian (Ω^{𝑒}, 𝐮)$ is a function that describes the distortion from the local coordinates ${ℒ︀}^{𝑒}$ onto the global element $Ω^{𝑒}$ . This distortion represents length, area, or volume depending on the dimensionality of the problem. Note that if the elements are the same size then $𝐌_{𝑒}$ is constant with respect to each element (this is the case in the code below).

The stiffness matrix is computed similarly but with a difference in the fact that the gradient must have its coordinates changed from global to local. Using the chain rule it follows that

\nabla_{𝐮} 𝑝 (𝐱 (𝐮)) = 𝐉 (𝐮) \nabla_{𝐱} 𝑝 (𝐱 (𝐮)),

(19)

where $𝐉 (𝐮) = {(𝜕 𝐱 / 𝜕 𝐮)}^{⊤}$ , i.e. $𝐽_{𝑖 𝑗} = (𝜕 𝑥_{𝑗}) / (𝜕 𝑢_{𝑖})$ , so that $jacobian (Ω^{𝑒}, 𝐮) = | det 𝐉 (𝐮) |$ . From this we can isolate the global gradient as

\nabla_{𝐱} 𝑝 (𝐱 (𝐮)) = {𝐉 (𝐮)}^{- 1} \nabla_{𝐮} 𝑝 (𝐱 (𝐮)) .

(20)

Putting things together the stiffness matrix can therefore be computed as

\begin{aligned} 𝐊 & = \int_{Ω} {(\nabla 𝐓 (𝐱))}^{⊤} \nabla 𝐓 (𝐱) d 𝐱 \\ \approx \sum_{𝑒 = 1}^{𝑁} \int_{Ω^{𝑒}} 𝐋_{𝑒}^{⊤} 𝐋_{𝑒} {(\nabla_{𝐱} 𝐓 (𝐱))}^{⊤} \nabla_{𝐱} 𝐓 (𝐱) 𝐋_{𝑒}^{⊤} 𝐋_{𝑒} d 𝐱 \\ = \sum_{𝑒 = 1}^{𝑁} 𝐋_{𝑒}^{⊤} (\int_{{ℒ︀}^{𝑒}} {(\nabla_{𝐱} 𝐓^{𝑒} (𝐮))}^{⊤} \nabla_{𝐱} 𝐓^{𝑒} (𝐮) jacobian (Ω^{𝑒}, 𝐮) d 𝐮) 𝐋_{𝑒} \\ = \sum_{𝑒 = 1}^{𝑁} 𝐋_{𝑒}^{⊤} (\int_{{ℒ︀}^{𝑒}} {(\nabla_{𝐮} 𝐓^{𝑒} (𝐮))}^{⊤} {({𝐉 (𝐮)}^{⊤})}^{- 1} {𝐉 (𝐮)}^{- 1} \nabla_{𝐮} 𝐓^{𝑒} (𝐮) jacobian (Ω^{𝑒}, 𝐮) d 𝐮) 𝐋_{𝑒} \\ \approx \sum_{𝑒 = 1}^{𝑁} 𝐋_{𝑒}^{⊤} (\underset{𝐊_{𝑒}}{\underset{⏟}{\sum_{𝑖 = 1}^{𝑄} 𝑤_{𝑖} {(\nabla_{𝐮} 𝐓^{𝑒} (𝐮_{𝑖}))}^{⊤} {({𝐉 (𝐮_{𝑖})}^{⊤})}^{- 1} {𝐉 (𝐮_{𝑖})}^{- 1} \nabla_{𝐮} 𝐓^{𝑒} (𝐮_{𝑖}) jacobian (Ω^{𝑒}, 𝐮_{𝑖})}}) 𝐋_{𝑒} . \end{aligned}

(21)

In the following code we implement three different ways of assembling the mass and stiffness matrices. Note that the first two implementations are only for educational purposes as they are both based on dense matrices.

Example 4: Matrix assembly implementations

Example 4 (Matrix assembly implementations)

julia


1
### Importing relevant packages
2
using FastGaussQuadrature
3
using BenchmarkTools
4
using LinearAlgebra
5
using SparseArrays
6
using ForwardDiff
7
using Random
8
Random.seed!(42)
9
## Defining geometry.
10
# Impedance tube of diameter 10 cm and length 1 m (used later)
11
D = 0.1             # 100 mm diameter
12
L = 10*D            # Length of the cavity
13
ne = 600            # Number of quadratic elements
14
nnt = 2*ne+1        # Total number of nodes
15
h = L/ne            # Length of the elements
16
x = Vector(0:h/2:L) # Coordinates table
17
18
## Computing the element matrices
19
# Local basis functions (gradient by ForwardDiff: slow but easy)
20
Tᵉ(u)  = [u .* (u .- 1)/2; 1 .- u .^2; u .* (u .+ 1)/2]'
21
∇Tᵉ(u) = hcat(ForwardDiff.derivative.(Tᵉ,u)...)
22
# All elements are equal, so the Jacobian is the same for each.
23
# It maps [-1,1] onto [x_i,x_{i+1}], i.e. length 2 onto length h.
24
jacobian(u) = h/2
25
# In 1D the Jacobian function and matrix are equal (not so in 2D).
26
J(u) = h/2
27
# The element matrices, the same for all elements of equal size.
28
Q = 3  # Number of Gauss points
29
u,w = gausslegendre(Q)
30
Me = sum(i -> w[i]*Tᵉ(u[i])'*Tᵉ(u[i])*jacobian(u[i]),1:Q)
31
Ke = sum(i -> w[i]*∇Tᵉ(u[i])'*J(u[i])^(-1)*J(u[i])^(-1)*∇Tᵉ(u[i])*
32
              jacobian(u[i]),1:Q)
33
34
## Assembly 1: with the localization matrices (never do this!)
35
function assembly1(Me, Ke, nnt, ne)
36
    K = zeros(nnt, nnt)  # Dense matrix! Not ideal!
37
    M = zeros(nnt, nnt)  # Dense matrix! Not ideal!
38
    for ie = 1:ne
39
        Le = zeros(3, nnt)
40
        Le[:, ie*2-1:ie*2+1] = Diagonal(ones(3))
41
        K += Le'*Ke*Le
42
        M += Le'*Me*Le
43
    end
44
    return K, M
45
end
46
@btime K, M = assembly1(Me, Ke, nnt, ne)
47
48
## Assembly 2: with indexing instead of localization matrices
49
function assembly2(Me, Ke, nnt, ne)
50
    K = zeros(nnt, nnt)  # Dense matrix! Not ideal!
51
    M = zeros(nnt, nnt)  # Dense matrix! Not ideal!
52
    for ie = 1:ne
53
        K[ie*2-1:ie*2+1, ie*2-1:ie*2+1] += Ke
54
        M[ie*2-1:ie*2+1, ie*2-1:ie*2+1] += Me
55
    end
56
    return K, M
57
end
58
@btime K, M = assembly2(Me, Ke, nnt, ne)
59
60
## Assembly 3: sparse, using the compact support of the elements
61
function assembly3(Me, Ke, nnt, ne)
62
    I = zeros(Int64, 4nnt - 3)
63
    J = zeros(Int64, 4nnt - 3)
64
    Kd = zeros(ComplexF64, length(I))
65
    Md = zeros(ComplexF64, length(I))
66
    for ie=1:ne
67
        for i = 1:3
68
            r = (8*(ie-1)+1 + 3*(i-1)):(8*(ie-1) + 3*i)
69
            I[r]  .= ie*2-1:ie*2-1+2
70
            J[r]  .= (ie-1)*2 + i
71
            Kd[r] += Ke[:, i]
72
            Md[r] += Me[:, i]
73
        end
74
    end
75
    K = sparse(I, J, Kd)
76
    M = sparse(I, J, Md)
77
    return K, M
78
end
79
@btime K, M = assembly3(Me, Ke, nnt, ne)

5. Examples

Cite this post

bibtex


@misc{paltorp2026femintroduction,
  author = {Paltorp, Mikkel},
  title = {{The Finite Element Method}},
  year = {2026},
  month = may,
  url = {https://mipals.github.io/pdf/2026/05-11-fem-introduction/fem-introduction.pdf},
  note = {Blog post}
}

Bibliography

[1] J.-F. Allard and N. Atalla, Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials, 2nd ed. Wiley, 2009. doi: 10.1002/9780470747339.