The Boundary Element Method

12 May 2026 · Download PDF · Download code

In the previous post the Helmholtz equation was solved using the Finite Element Method (FEM), where the whole domain $\mathrm{\Omega }$ is divided into elements. The Boundary Element Method (BEM) instead first rewrites the Helmholtz equation as an integral equation that lives only on the boundary $\mathrm{\Gamma }=𝜕\mathrm{\Omega }$ [1]. As a result, only the boundary has to be divided into elements, which is particularly attractive for exterior problems, such as the sound radiated from a vibrating object into free space, where the domain $\mathrm{\Omega }$ in theory is infinite.

The building blocks are the same as in the previous post: the unknown functions are approximated as $𝑝\left(𝐱\right)\approx 𝐓\left(𝐱\right)𝐩$, the elements are parameterized by local coordinates $𝐮\in {\mathrm{ℒ︀}}^{𝑒}$, the element localization matrices ${𝐋}^{𝑒}$ pick out the local values ${𝐩}^{𝑒}={𝐋}^{𝑒}𝐩$, and integrals over the elements are approximated by a quadrature rule. The only difference in these building blocks is that the elements are now surface elements ${\mathrm{\Gamma }}^{𝑒}$ instead of volume elements ${\mathrm{\Omega }}^{𝑒}$. What changes is the equation being discretized and the approach used to discretize it.

1. The Kirchhoff-Helmholtz integral equation

For time-harmonic problems the pressure can be written as $𝑝\left(𝐱,𝑡\right)=\text{Re}\left\{𝑝\left(𝐱\right){𝑒}^{𝛼𝑖𝜔𝑡}\right\}$, where $𝛼=\pm 1$ is the sign of the chosen time dependency. Inserting this into the wave equation results in the Helmholtz equation

which looks independent of $𝛼$. This can cause confusion when the time dependency is not stated explicitly, as the relation between the normal derivative of the pressure and the normal velocity ${𝑣}_{𝐧}$ does depend on $𝛼$

where the normal $𝐧$ points into the domain $\mathrm{\Omega }$, i.e. out of a vibrating object, and $𝑠\left(𝛼\right)$ is introduced for simplicity, similarly to [2]. Note that this is opposite to the previous post, where the boundary term of the integration by parts used the outward normal of $\mathrm{\Omega }$.

In three dimensions the Green’s function, which will come in handy soon enough, is given by

It is the solution of the Helmholtz equation with a point source

where $𝛿$ is the Dirac delta function, while the sign $𝛼$ makes the waves of $𝐺$ travel outwards from the source. The Green’s function takes two inputs, $𝐱$ and $𝐲$, commonly referred to as the target location ($𝐱$) and the source location ($𝐲$), or simply the target(s) and source(s).

The derivation of the integral equation starts in the same way as the FEM: we multiply the Helmholtz equation with a test function $𝑞$ and integrate over the domain

The difference is that we now integrate by parts twice, moving both derivatives from $𝑝$ onto $𝑞$. Writing the normal derivative of $𝑝$ in terms of the normal velocity, it follows that

where the signs of the boundary terms are due to $𝐧$ pointing into $\mathrm{\Omega }$. For exterior problems the boundary also includes a sphere at infinity, whose integral vanishes when $𝑝$ satisfies the Sommerfeld radiation condition below. Choosing the test function as the Green’s function, $𝑞\left(𝐲\right)=𝐺\left(𝐱,𝐲\right)$, the volume integral collapses due to (4), and we are left with the Kirchhoff-Helmholtz integral equation

where $𝜁\left(𝐱\right)$ is the so-called integral free term depending on the geometry at $𝐱$. For $𝐱$ in the interior of $\mathrm{\Omega }$ it is equal to one, in which case (6) gives the pressure at $𝐱$ from the values of $𝑝$ and ${𝑣}_{𝐧}$ on the boundary. For $𝐱$ on a smooth part of the boundary it is equal to $1/2$, while at corners and edges it takes values between 0 and 1. Setting $𝐱\in \mathrm{\Gamma }$ makes (6) an equation involving only functions on the boundary: a boundary integral equation.

The time dependency still enters (6) through $𝐺$ and $𝑠\left(𝛼\right)$. Changing $𝛼$ conjugates the Green’s function and flips the sign of $𝑠\left(𝛼\right)$, so a sign error in either of them results in a solution of the wrong problem. Finally, the derivation assumes that $𝑝$ satisfies the Sommerfeld radiation condition

where $𝑑$ is the dimension, meaning that no waves come in from infinity. In the FEM this condition has to be imposed by truncating the domain with an artificial boundary, whereas in the BEM it is built into the Green’s function, so that the solution of (6) satisfies it automatically.

2. Collocation

As in the FEM, we multiply (6) with a test function $𝜑$, this time integrating over the boundary,

Choosing $𝜑\left(𝐱\right)\approx {𝐚}^{\top }{𝐓\left(𝐱\right)}^{\top }$ would result in a Galerkin BEM. The integrals of the Galerkin BEM are, however, double surface integrals with a singular integrand, which makes them expensive to compute. Instead, the most common approach is the so-called collocation approach, where the test function is a sum of Dirac delta functions placed at the collocation points ${𝐭}_1,\dots ,{𝐭}_{𝑛}\in \mathrm{\Gamma }$

where $𝐚\in {ℂ}^{𝑛}$ is again a vector of arbitrary coefficients. The collocation points are in most cases chosen as the interpolation nodes of the basis functions, so that the cardinal property $𝐓\left({𝐭}_{𝑗}\right)𝐩={𝑝}_{𝑗}$ can be used. Inserting $𝜑$ into (8), the sifting property of the Dirac delta function evaluates the outer integral at the collocation points, so that (8) becomes ${𝐚}^{\top }𝐫=0$ with

In its current form the unknowns are still functions, which a computer cannot be tasked to find. Instead, we insert the approximations $𝑝\left(𝐲\right)\approx 𝐓\left(𝐲\right)𝐩$ and ${𝑣}_{𝐧}\left(𝐲\right)\approx 𝐓\left(𝐲\right){𝐯}_{𝐧}$ and handle the three vectors one at a time. For the first vector, the cardinal property of the basis functions gives

where $𝜻={\left[\begin{array}{ccc}𝜁\left({𝐭}_1\right)& \dots & 𝜁\left({𝐭}_{𝑛}\right)\end{array}\right]}^{\top }$. For the two integrals, the coefficients $𝐩$ and ${𝐯}_{𝐧}$ do not depend on $𝐲$ and can be moved outside of the integrals

where each entry of the vectors is a row vector, as $𝐓\left(𝐲\right)$ is, so that $𝐅$ and $𝐆$ are $𝑛\times 𝑛$ matrices. Collecting the terms, ${𝐚}^{\top }𝐫=0$ becomes

As this has to hold for all $𝐚\in {ℂ}^{𝑛}$, the term inside the parentheses must vanish, and the collocation BEM therefore reads

where we simply defined $𝐇=\text{diag}\left(𝜻\right)-𝐅$ for the sake of convenience. In a typical problem ${𝐯}_{𝐧}$ is known (e.g. from the vibration of a structure) and (14) is solved for $𝐩$.

3. Computing the matrices

While (14) is easily stated, it does not mention how to actually compute $𝐅$ and $𝐆$. The procedure is the same as for the mass and stiffness matrices in the previous post: the boundary integral is split into a sum over the elements, each element integral is moved to local coordinates, and the result is approximated by a quadrature rule. On element $𝑒$ the pressure is interpolated as $𝑝\left({𝐲}^{𝑒}\left(𝐮\right)\right)\approx {𝐓}^{𝑒}\left(𝐮\right){𝐋}^{𝑒}𝐩$, so that the $𝑘$th row of $𝐅$ becomes

where the first approximation is due to the geometry being approximated by the elements, and the second is due to the quadrature. The $𝑘$th row of $𝐆$ is computed in the same way, with $𝐺\left({𝐭}_{𝑘},𝐲\right)$ in place of its normal derivative. The jacobian describes the local deformation of the map from local to global coordinates. For the curve elements of 2D problems this is a change of length, while for the surface elements of 3D problems it is a change of area

The new difficulty compared to the FEM is that the Green’s function is singular at $𝐲={𝐭}_{𝑘}$, so a standard quadrature rule is inaccurate on the elements that contain, or lie close to, the collocation point. This is why (15) highlights that the number of quadrature points $𝑄\left({𝐭}_{𝑘},𝑒\right)$ can depend on both the collocation point and the element. For the elements containing ${𝐭}_{𝑘}$, special techniques are needed, such as splitting the element at ${𝐭}_{𝑘}$ and using a transformation whose jacobian cancels the $1/𝑟$ singularity.

Example 1: Pulsating sphere

Example 1 (Pulsating sphere)   A sphere of radius $𝑎$ whose surface moves uniformly with the normal velocity ${𝑣}_0$ radiates the pressure $𝑝\left(𝑟\right)=𝐴exp\left(-𝛼𝑖𝑘𝑟\right)/𝑟$. Matching the radial velocity at $𝑟=𝑎$ gives the analytical surface pressure

$$𝑝\left(𝑎\right)={𝜌}_0{𝑐}_0{𝑣}_0\frac{𝛼𝑖𝑘𝑎}{1+𝛼𝑖𝑘𝑎}.$$

(17)

To compute it with the BEM we use BoundaryIntegralEquations.jl, which assembles $𝐅$ and $𝐆$ as in (15), including the special quadrature for the elements close to the collocation points. The sphere is meshed with continuous quadratic triangles, both for the geometry (${𝐍}^{𝑒}$) and for the pressure and normal velocity (${𝐓}^{𝑒}$), and the collocation points are the interpolation nodes. The package uses the time dependency $exp\left(-𝑖𝜔𝑡\right)$, i.e. $𝛼=-1$, and the normals of the mesh point into the sphere. The latter flips the sign of the matrix returned as $𝐅$, and the package returns $1-𝜁$ instead of $𝜁$.

julia

Copy codeCopyCopied!

1
using BoundaryIntegralEquations, LinearAlgebra
2
## Physical parameters. The sphere pulsates with the normal velocity v₀.
3
a = 1.0             # Radius of the sphere
4
ρ = 1.21            # Density of air
5
c = 343.0           # Speed of sound
6
v₀ = 1.0            # Normal velocity of the surface
7
s = im*ρ*c          # s(α) for the time dependency exp(-iωt), i.e. α = -1
8
## Mesh: 464 quadratic triangles with 930 nodes
9
pkg_dir = pkgdir(BoundaryIntegralEquations)
10
mesh_file = joinpath(pkg_dir, "examples", "meshes", "sphere_1m_coarse")
11
mesh = load3dTriangularComsolMesh(mesh_file)
12
vₙ = fill(v₀, size(mesh.sources, 2))
13
## Solving (diag(ζ) - F)p = -s k G vₙ for each wavenumber.
14
# The normals of the mesh point into the sphere, so the F of the package
15
# is -F, and the integral free term of the exterior problem is ζ = 1 - C.
16
ka = range(0.05, 5, length=100)
17
p_bem = map(ka) do kai
18
    k = kai/a
19
    F, G, C = assemble_parallel!(mesh, k, mesh.sources; progress=false)
20
    H = Diagonal(1 .- C) + F
21
    p = H \ (-s*k*G*vₙ)
22
    sum(p)/length(p)    # The surface pressure is constant; average it.
23
end
24
## The analytical surface pressure of the pulsating sphere.
25
p_analytical = @. ρ*c*v₀*(-im*ka)/(1 - im*ka)

Figure 1: The surface pressure of the pulsating sphere normalized by ${𝜌}_0{𝑐}_0{𝑣}_0$, analytical and computed with the BEM using 464 quadratic elements.

4. Acceleration methods

A considerable downside of the BEM is that $𝐇,𝐆\in {ℂ}^{𝑛\times 𝑛}$ are dense, unlike the sparse mass and stiffness matrices of the FEM. The memory used to store them therefore scales as $\mathrm{𝒪︀}\left({𝑛}^2\right)$, and a direct solve of (14) is $\mathrm{𝒪︀}\left({𝑛}^3\right)$, rendering the direct application of the method unusable for large $𝑛$. This is somewhat offset by the fact that only the boundary is discretized, but for large problems it is a serious limitation. However, various methods have been developed that reduce the memory and the cost of a matrix-vector product to $\mathrm{𝒪︀}\left(𝑛log𝑛\right)$ or even $\mathrm{𝒪︀}\left(𝑛\right)$ by exploiting the rank structure of the problem. As these methods only define the matrix-vector products, (14) must then be solved by an iterative method, the most common being the generalized minimal residual method (GMRES) [3], [4], [5].

The aim of this section is to show how the matrix-vector products with $𝐆$ and $𝐅$ can be computed in $\mathrm{𝒪︀}\left(𝑛\right)$ (or $\mathrm{𝒪︀}\left(𝑛log𝑛\right)$) time using the fast multipole method (FMM) or hierarchical matrices ($\mathrm{\mathscr{H}︀}$-matrices). To do so, we start by assuming that all elements are integrated using the same $𝑄$ quadrature points.¹¹ This results in inaccurate integrals for the elements close to the collocation point, which are fixed by a near-field correction below. Multiplying the $𝑘$th row of $𝐆$ with a vector $𝐳$ it follows that

where the subscript $𝑗$ refers to an ordering of the quadrature points of all elements, and $𝑒\left(𝑗\right)$ returns the element that quadrature point $𝑗$ lies on. The matrix $𝐂\in {ℝ}^{𝑁𝑄\times 𝑛}$ is extremely sparse and can be thought of as a map from the values $𝐳$ at the collocation points to the coefficients $𝐜=𝐂𝐳$ at the quadrature points ${𝐲}_{𝑗}$. The same approach applies to $𝐅$, where the normal derivative is written as ${𝐧\left({𝐲}_{𝑗}\right)}^{\top }{\nabla }_{𝐲}𝐺\left({𝐭}_{𝑘},{𝐲}_{𝑗}\right)$

Applying (18) or (19) directly for all $𝑛$ rows results in products that scale as $\mathrm{𝒪︀}\left(𝑁𝑄𝑛\right)$, which could hurt more than we gain, as $𝑁𝑄\gtrsim 𝑛$. The point is that the dense part of the product, a sum of Green’s functions over the points ${𝐲}_{𝑗}$, is exactly what the FMM and $\mathrm{\mathscr{H}︀}$-matrices accelerate. Doing so blindly, however, keeps the inaccurate integration of the elements close to the collocation points. To fix this, a sparse near-field correction $𝐒$ is added, which subtracts the inaccurate integration using only $𝑄$ quadrature points and adds the accurate one. Applying (18) or (19) for all rows at the same time can therefore be thought of as

where $𝐁$ is the dense part approximated by either the FMM or an $\mathrm{\mathscr{H}︀}$-matrix (see (22) and (23) below), $𝐂$ is the coefficient map, and $𝐒$ is the near-field correction. As $𝐒$ and $𝐂$ are highly sparse, both their assembly and their products scale as $\mathrm{𝒪︀}\left(𝑛𝜏\right)$ with $𝜏\ll 𝑛$. Using an approximation of $𝐁$ whose product scales linearly in time and memory therefore results in a representation of $𝐀$ that scales similarly.

To compute $𝐒$ efficiently, its sparsity pattern must be known in advance. For collocation this could, for example, be the elements neighboring each interpolation node, as shown in Figure 2. For discontinuous elements, whose nodes lie inside the elements, the support of each basis function is a single element (panel (b)). While this works well in many cases, it fails when elements that are not neighbors are geometrically close, as is the case for e.g. thin structures.

Figure 2: The support of the basis function of a node (red point) is shown in red. In some cases the near-field correction is extended to also include the blue elements.

4.1. The fast multipole method

The FMM can be used to accelerate sums of the Green’s function ($𝐆$) as well as of its normal derivative ($𝐅$) [6], [7], [8]. Throughout the years, many good resources explaining the intricacies of the FMM have been written. For a gentle introduction I recommend the short course by Greengard and the pedestrian introduction by Ying [9], while the books by Gumerov and Duraiswami [10] and Liu [11] give a more detailed look with an emphasis on the BEM. Ying also has a wonderful introductory talk. Given these resources, the details of the FMM are left out here, and we focus on the connection to the BEM. For this, it is enough to think of the FMM as a way of speeding up sums of the form

where the $𝑗$th term is left out when ${𝐭}_{𝑘}={𝐲}_{𝑗}$. Looking back at (18) and (19), this is exactly what is computed for each row of $𝐆$ and $𝐅$, with the coefficients $𝐜=𝐂𝐳$. For a fixed precision, applying the FMM to all $𝑛$ collocation points scales as $\mathrm{𝒪︀}\left(𝑛+𝑁𝑄\right)$ at low frequencies, and as $\mathrm{𝒪︀}\left(\left(𝑛+𝑁𝑄\right)log\left(𝑛+𝑁𝑄\right)\right)$ for the wideband methods needed at high frequencies. Many FMM libraries exist, such as ExaFMM, ScalFMM and FMM3D, the last of which has a Julia interface. For 2D problems, the companion library FMM2D can be used from Julia through FMM2D.jl.²² I maintain this Julia package.

4.2. Hierarchical matrices

The $\mathrm{\mathscr{H}︀}$-matrix approach differs from the FMM in that it does not accelerate the sums directly [12], [13]. Instead, it approximates the matrix that corresponds to computing the sums for all $𝑛$ collocation points at the same time. For $𝐆$ this is the matrix

where the approximation is based on the fact that subblocks of the matrix corresponding to points far from each other are well approximated by low-rank matrices. For $𝐅$ the matrix is

How the hierarchical low-rank approximation is computed is outside the scope of this post, but the interested reader can look at [14], [15], [16], of which [14] is my preferred introduction. Furthermore, a wonderful lecture series by Alex Barnett, Adrianna Gillman and Per-Gunnar Martinsson from the 2014 CBMS-NSF Conference: Fast Direct Solvers for Elliptic PDEs can be found online. There are many implementations of hierarchical matrices, such as HLIBpro, HLib and H2Lib, and in Julia KernelMatrices.jl [17], HssMatrices.jl and HMatrices.jl implement various formats with different properties. A longer list can be found here. An often mentioned strength of $\mathrm{\mathscr{H}︀}$-matrices is that their LU factorization can be computed efficiently, which is useful when iterative solvers converge slowly, e.g. near resonances, although the ranks, and thereby the cost, grow with the frequency. However, an LU factorization cannot be applied directly to the form of (20).

4.3. The interpolated factored Green function

A more recent approach is the interpolated factored Green function (IFGF) method [18]. An advantage of the IFGF is that it does not need the split into low and high frequencies required by the FMM, while at the same time being easily parallelized, as it does not rely on the fast Fourier transform [19]. However, as opposed to the purely algebraic $\mathrm{\mathscr{H}︀}$-matrices, the IFGF is based on interpolation. It therefore requires that the matrix is generated by a function that can be evaluated at any pair of points $\left(𝐱,𝐲\right)$. Unfortunately, this is not the case for $𝐅$, as its entries depend on the normal, which is computed from the local coordinates $𝐮$ of the elements. A workaround is to approximate the full gradient of the Green’s function with the IFGF, as it can be evaluated at any $\left(𝐱,𝐲\right)$, and then compute the normal derivative from it. A simple implementation of the method can be found in IFGF.jl.

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