AMATH 843: Integral Equation Methods for PDEs

These notes synthesize material from R. Kress, Linear Integral Equations (3rd ed., Springer); K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind (Cambridge); S. Sauter and C. Schwab, Boundary Element Methods (Springer); the foundational papers of L. Greengard and V. Rokhlin on the Fast Multipole Method; and NYU Courant MATH-GA 2011 lecture materials on integral equations.

---

Chapter 1: Potential Theory and Green’s Functions

1.1 Fundamental Solutions

The theory of integral equations for partial differential equations rests on a simple but powerful idea: rather than solving a PDE directly in a volume, one can reformulate the problem in terms of integral operators defined on the boundary of the domain. This reformulation reduces the dimension of the problem by one and often yields equations with superior analytical and numerical properties. The starting point for the entire enterprise is the concept of a fundamental solution.

A fundamental solution for a linear partial differential operator \(L\) is a distribution \(\Phi(\mathbf{x}, \mathbf{y})\) satisfying \(L_{\mathbf{x}} \Phi(\mathbf{x}, \mathbf{y}) = \delta(\mathbf{x} - \mathbf{y})\) in the sense of distributions, where \(\delta\) denotes the Dirac delta function. Physically, \(\Phi\) represents the response at the point \(\mathbf{x}\) due to a unit point source located at \(\mathbf{y}\). From the fundamental solution, one constructs layer potentials and representation formulas that form the backbone of boundary integral equation methods.

The logarithmic singularity in two dimensions and the Coulomb-type singularity in three dimensions reflect the different geometries of the respective spaces. The sign convention varies across the literature; we follow Kress, where the two-dimensional fundamental solution carries a negative sign so that the Laplacian (not its negative) equals the delta function.

The choice of the Hankel function \(H_0^{(1)}\) (rather than \(H_0^{(2)}\)) enforces the Sommerfeld radiation condition, which selects outgoing waves at infinity. This is the physically relevant choice for scattering problems. As \(k \to 0\), the Helmholtz fundamental solution degenerates (up to constants) to the Laplace fundamental solution, a fact that provides a useful consistency check and connects the static and time-harmonic theories.

1.2 Green’s Representation Formula

The bridge between volume-based PDE theory and boundary integral equations is the Green’s representation formula, which expresses a solution of a PDE inside a domain entirely in terms of boundary data and the fundamental solution. The derivation relies on Green’s second identity, a classical tool from vector calculus.

By choosing \(v(\mathbf{x}) = \Phi(\mathbf{x}, \mathbf{y})\) in Green’s second identity and carefully excising a small ball around the singularity at \(\mathbf{y}\), one obtains the celebrated Green’s representation formula. This formula is the Rosetta Stone of the boundary integral equation method: it translates between volume PDE information and boundary integral data.

This representation expresses any harmonic function inside \(\Omega\) in terms of its boundary values \(u|_\Gamma\) and its normal derivative \(\frac{\partial u}{\partial \nu}\big|_\Gamma\). For a well-posed boundary value problem, one of these quantities is given as data and the other is unknown; the task of the boundary integral equation is to determine the unknown boundary quantity.

1.3 Single-Layer and Double-Layer Potentials

The two integrals appearing in Green’s representation formula motivate the definition of two fundamental objects: the single-layer and double-layer potentials. These potentials, defined for arbitrary densities on the boundary, are the building blocks from which all boundary integral formulations are constructed.

Physically, the single-layer potential represents the field generated by a continuous distribution of point sources of strength \(\varphi\) on the boundary, while the double-layer potential represents the field generated by a distribution of dipoles. The single-layer potential produces a continuous field across \(\Gamma\) but with a discontinuous normal derivative, while the double-layer potential itself is discontinuous across \(\Gamma\) but has a continuous normal derivative. These complementary regularity properties are precisely captured by the jump relations.

1.4 Jump Relations Across Boundaries

The behavior of layer potentials as the evaluation point \(\mathbf{y}\) approaches the boundary \(\Gamma\) is one of the central results in potential theory. The jump relations govern this limiting behavior and are essential for deriving boundary integral equations from representation formulas.

The \(\pm \frac{1}{2}\) terms arise from the solid angle subtended by the boundary at a point on itself. For smooth boundaries, this solid angle is exactly \(2\pi\) (in 2D) or \(2\pi\) steradians (in 3D, corresponding to a half-sphere), leading to the factor of \(\frac{1}{2}\). At corners or edges of non-smooth boundaries, the solid angle differs from \(2\pi\) and the coefficient \(\frac{1}{2}\) is replaced by a geometry-dependent value \(\frac{\omega(\mathbf{x}_0)}{2\pi}\), where \(\omega(\mathbf{x}_0)\) is the interior solid angle at the boundary point \(\mathbf{x}_0\). For a two-dimensional domain with a corner of interior angle \(\alpha\), the coefficient becomes \(\frac{\alpha}{2\pi}\).

These jump relations were first established rigorously by Carl Neumann in the 1870s for smooth boundaries, and their extension to Lipschitz domains was completed in the 1980s by Verchota, Dahlberg, and Kenig. They remain the cornerstone of the boundary integral equation method.

The jump relations immediately yield boundary integral equations. For instance, if \(u\) is harmonic in \(\Omega\) and we represent it as a double-layer potential \(u = \mathcal{D}\psi\), then taking the interior limit and using the jump relation gives

\[ u(\mathbf{y}_0) = -\frac{1}{2}\psi(\mathbf{y}_0) + \int_\Gamma \frac{\partial \Phi(\mathbf{x}, \mathbf{y}_0)}{\partial \nu(\mathbf{x})} \psi(\mathbf{x}) \, dS(\mathbf{x}), \quad \mathbf{y}_0 \in \Gamma. \]

If \(u|_\Gamma = f\) is prescribed (Dirichlet problem), this becomes a second-kind Fredholm integral equation for the unknown density \(\psi\). This observation, due to Carl Neumann, launched the entire field.

Chapter 2: Fredholm Theory

2.1 Compact Operators on Banach Spaces

The boundary integral equations that arise from potential theory are, in favorable circumstances, Fredholm integral equations of the second kind. The analysis of such equations relies on the spectral theory of compact operators, a beautiful and complete theory developed by Ivar Fredholm (1903), Frigyes Riesz (1918), and Juliusz Schauder (1930). This chapter develops the abstract functional-analytic framework before we return to its application to boundary integral equations.

Compact operators can be thought of as “almost finite-dimensional” operators: they are limits (in the operator norm) of finite-rank operators when \(X\) is a Hilbert space, and they enjoy spectral properties reminiscent of finite-dimensional linear algebra. The prototypical example in our setting is an integral operator with a continuous or weakly singular kernel.

More generally, integral operators with weakly singular kernels — kernels that blow up at the diagonal but are integrable — are also compact. This is the case for the single-layer boundary operator on smooth boundaries.

The compactness of weakly singular integral operators follows from a refinement of the Arzela-Ascoli argument: although the kernel is unbounded, the integrability of the singularity ensures that the operator still maps bounded sets to equicontinuous sets. This result, due to various authors including Kress and Atkinson, underpins the applicability of Fredholm theory to boundary integral operators.

This result is important because it allows us to pass between an integral equation and its adjoint (or transpose) equation, a maneuver that is frequently necessary in indirect boundary integral formulations.

2.2 Fredholm Integral Equations

The integral equations that arise from boundary value problems fall into two main classes. Understanding their distinct characters is essential for choosing appropriate solution strategies.

First-kind equations are inherently ill-posed in the sense of Hadamard: small perturbations in \(f\) can lead to large changes in \(\varphi\), because the compact operator \(T\) has eigenvalues accumulating at zero. Second-kind equations, by contrast, enjoy the full power of the Fredholm alternative and are generally well-posed. This is why the indirect and direct formulations of boundary integral equations strive to produce second-kind equations whenever possible.

In operator notation, a second-kind equation reads \((I - \lambda T)\varphi = f\), where \(T\) is the compact integral operator and \(I\) is the identity. The structure \(I - \lambda T\) (identity minus a compact perturbation) is the key that unlocks the Fredholm theory.

2.3 The Fredholm Alternative

The Fredholm alternative is the central result of the theory. It provides a complete dichotomy: either the equation \((I - \lambda T)\varphi = f\) has a unique solution for every right-hand side, or the homogeneous equation has nontrivial solutions and solvability of the inhomogeneous equation is governed by a finite-dimensional compatibility condition. There is no middle ground.

The Fredholm alternative reduces the question of solvability of an integral equation to a finite-dimensional problem: one must check whether the homogeneous equation has nontrivial solutions and, if so, verify a finite number of orthogonality conditions. For the boundary integral equations of potential theory, uniqueness of the underlying PDE boundary value problem often (but not always) guarantees that the first alternative holds.

2.4 Spectral Theory of Compact Operators

The spectral theory of compact operators on Banach (and especially Hilbert) spaces provides detailed information about the eigenvalues and eigenfunctions that govern the behavior of integral equations.

For self-adjoint compact operators on Hilbert spaces, the spectral theorem takes an especially elegant form that is directly relevant to the analysis of symmetric kernels.

The rate of decay of eigenvalues is intimately connected to the smoothness of the kernel: smoother kernels have faster-decaying eigenvalues, and this in turn controls the conditioning of the integral equation and the convergence rates of numerical methods.

2.5 Riesz-Schauder Theory and the Resolvent

The resolvent operator \(R(\lambda; T) = (I - \lambda T)^{-1}\), when it exists, is a bounded operator that depends analytically on \(\lambda\) in the resolvent set. The Riesz-Schauder theory provides a complete description of the singularities of this analytic operator-valued function.

This theorem is essential in the analysis of integral equations depending on a parameter (such as the wavenumber \(k\) in the Helmholtz equation). It guarantees that failure of unique solvability can only occur at isolated parameter values, and the behavior near these exceptional values is completely controlled by finite-dimensional linear algebra.

Chapter 3: Boundary Integral Equations

3.1 Direct and Indirect Methods

With the functional-analytic machinery of Fredholm theory in hand, we now return to the boundary integral equations that arise from potential theory. There are two fundamentally different approaches to formulating boundary integral equations for a given PDE boundary value problem: the direct method and the indirect method. Each has its own advantages, and the choice between them depends on the problem at hand.

The direct method starts from Green’s representation formula (Theorem 1.5) and uses the prescribed boundary data directly. The resulting integral equation involves the physical boundary data (values of \(u\) and \(\frac{\partial u}{\partial \nu}\)), and the unknown is the missing piece of boundary data. This approach is natural and yields equations with clear physical interpretations.

The indirect method, by contrast, represents the solution as a single-layer or double-layer potential with an unknown density function that typically has no direct physical meaning. The boundary condition then yields an integral equation for this fictitious density. The indirect method often produces simpler integral equations at the cost of a less transparent physical interpretation.

These four operators — \(S\), \(D\), \(D'\), and \(N\) — are collectively known as the Calderon projector operators or boundary integral operators. They satisfy important algebraic relations (the Calderon identities) that connect the interior and exterior problems and provide consistency checks for numerical implementations.

The Calderon identities have important implications for the conditioning of boundary integral equations and for the construction of preconditioners. For instance, the relation \(D^2 = \frac{1}{4}I + SN\) suggests that \(S\) can serve as a preconditioner for \(N\) (and vice versa), an idea that underlies the operator preconditioning approach of Steinbach and Wendland.

3.2 Interior Dirichlet Problem for the Laplace Equation

We begin with the most classical boundary value problem: the interior Dirichlet problem for Laplace’s equation.

The existence and uniqueness of solutions to this problem are guaranteed by classical potential theory (the maximum principle gives uniqueness; Perron’s method gives existence for general domains).

Direct method. Green’s representation formula, combined with the jump relation for the double-layer potential, yields the boundary integral equation

\[ \frac{1}{2} u(\mathbf{x}) + \int_\Gamma \frac{\partial \Phi(\mathbf{y}, \mathbf{x})}{\partial \nu(\mathbf{y})} u(\mathbf{y}) \, dS(\mathbf{y}) = \int_\Gamma \Phi(\mathbf{y}, \mathbf{x}) \frac{\partial u}{\partial \nu}(\mathbf{y}) \, dS(\mathbf{y}), \quad \mathbf{x} \in \Gamma. \]

Here \(u|_\Gamma = f\) is known, and the unknown is the Neumann data \(g = \frac{\partial u}{\partial \nu}\big|_\Gamma\). Rearranging, this is a first-kind equation \(Sg = \frac{1}{2}f + Df\) for \(g\).

Indirect method (double-layer ansatz). Represent the solution as a double-layer potential:

\[ u(\mathbf{y}) = \int_\Gamma \frac{\partial \Phi(\mathbf{x}, \mathbf{y})}{\partial \nu(\mathbf{x})} \psi(\mathbf{x}) \, dS(\mathbf{x}), \quad \mathbf{y} \in \Omega. \]

Taking the interior limit and applying the jump relation gives the second-kind Fredholm equation

\[ -\frac{1}{2}\psi(\mathbf{x}) + (D\psi)(\mathbf{x}) = f(\mathbf{x}), \quad \mathbf{x} \in \Gamma. \]

3.3 Interior Neumann Problem

For the Neumann problem, the indirect method uses a single-layer ansatz \(u = \mathcal{S}\varphi\). Taking the interior normal derivative and using the jump relation gives

\[ \frac{1}{2}\varphi(\mathbf{x}) + (D'\varphi)(\mathbf{x}) = g(\mathbf{x}), \quad \mathbf{x} \in \Gamma. \]

This is again a second-kind Fredholm equation, but now the homogeneous equation \(\frac{1}{2}\varphi_0 + D'\varphi_0 = 0\) has the one-dimensional solution space spanned by the constant function (reflecting the non-uniqueness of the Neumann problem up to constants). By the Fredholm alternative, solvability requires \(g \perp \ker(\frac{1}{2}I + D) = \text{span}\{1\}\), which is precisely the compatibility condition \(\int_\Gamma g \, dS = 0\).

3.4 Exterior Problems and Uniqueness

Exterior boundary value problems require special care because the unbounded domain introduces complications not present in the interior case. For the Laplace equation, one must impose a decay condition at infinity to ensure uniqueness; for the Helmholtz equation, one imposes the Sommerfeld radiation condition.

The indirect double-layer formulation for the exterior Dirichlet problem leads to the equation

\[ \frac{1}{2}\psi(\mathbf{x}) + (D\psi)(\mathbf{x}) = f(\mathbf{x}), \quad \mathbf{x} \in \Gamma, \]

where the sign in front of the \(\frac{1}{2}\) has changed relative to the interior problem due to the reversal of the limiting direction. However, this equation has a nontrivial kernel: the constant function satisfies the homogeneous equation in two dimensions. This reflects the fact that the exterior Laplace Dirichlet problem in 2D requires careful handling of the logarithmic capacity.

The idea of the combined-layer approach is to break the symmetry that causes the nullspace of the double-layer operator to be nontrivial. This idea, due to Brakhage and Werner (1965) and independently to Leis (1965) and Panich (1965), is one of the most important practical devices in the boundary integral equation method.

3.5 Integral Equations for the Helmholtz Equation

The Helmholtz equation \(\Delta u + k^2 u = 0\) arises in the study of time-harmonic wave propagation at frequency \(\omega = ck\), where \(c\) is the wave speed. The boundary integral formulations for the Helmholtz equation are closely analogous to those for the Laplace equation, but with the Helmholtz fundamental solution \(\Phi_k\) replacing \(\Phi\).

For the exterior Dirichlet problem (sound-soft scattering), the indirect double-layer formulation yields

\[ \frac{1}{2}\psi(\mathbf{x}) + (D_k\psi)(\mathbf{x}) = f(\mathbf{x}), \quad \mathbf{x} \in \Gamma, \]

where \(D_k\) is the double-layer operator formed with the Helmholtz kernel. This equation suffers from the notorious problem of irregular frequencies (also called spurious resonances or fictitious eigenvalues): at wavenumbers \(k\) corresponding to Dirichlet eigenvalues of the Laplacian in the interior domain \(\Omega\), the integral equation fails to be uniquely solvable even though the exterior scattering problem itself is well-posed.

Chapter 4: Numerical Methods — Nystrom and Collocation

4.1 The Nystrom Method

We turn now to the numerical solution of the boundary integral equations developed in Chapters 1-3. The Nystrom method, named after E.J. Nystrom who introduced it in 1930, is one of the simplest and most elegant approaches: it replaces the integral in the integral equation by a quadrature rule. The method is particularly effective for equations with smooth or mildly singular kernels on smooth geometries.

The beauty of the Nystrom method is that it simultaneously provides both the discrete approximation at the nodes and, through interpolation, an approximation everywhere on \(\Gamma\). Once the values \(\varphi_N(\mathbf{x}_j)\) are computed, the Nystrom interpolation formula

\[ \varphi_N(\mathbf{x}) = f(\mathbf{x}) + \sum_{j=1}^N w_j K(\mathbf{x}, \mathbf{x}_j) \varphi_N(\mathbf{x}_j) \]

gives the approximate solution at any point \(\mathbf{x} \in \Gamma\).

The key message of this theorem is that the convergence rate of the Nystrom method is controlled by the convergence rate of the underlying quadrature rule applied to the specific integrand \(K(\mathbf{x}, \cdot)\varphi(\cdot)\). For smooth kernels and smooth solutions on smooth boundaries, using a high-order quadrature rule yields spectral (superalgebraic) convergence.

4.2 Quadrature Rules for Smooth and Singular Kernels

The choice of quadrature rule is the heart of the Nystrom method. For smooth periodic integrands (as arise on smooth closed curves in 2D), the trapezoidal rule is optimal in a precise sense.

For singular kernels, one must modify the quadrature to handle the singularity. The most important case is the logarithmic singularity that arises from the single-layer operator in two dimensions. After parametrizing the boundary, the kernel takes the form \(K(s,t) = A(s,t) \log|s - t| + B(s,t)\), where \(A\) and \(B\) are smooth.

4.3 Collocation Methods

Collocation methods provide an alternative to the Nystrom approach, particularly suited to integral equations on surfaces in three dimensions where the geometry is typically represented by a mesh of panels (triangles or quadrilaterals).

The resulting linear system is \((I_N - K_N)\boldsymbol{\alpha} = \mathbf{f}\), where \((I_N)_{ij} = \phi_j(\mathbf{x}_i)\) and \((K_N)_{ij} = \int_\Gamma K(\mathbf{x}_i, \mathbf{y}) \phi_j(\mathbf{y}) \, dS(\mathbf{y})\). The matrix entries \((K_N)_{ij}\) require evaluating integrals of the kernel against the basis functions, which must be done carefully when the kernel is singular.

The most common choice of approximation space is the space of piecewise constant or piecewise linear functions on a triangulation of \(\Gamma\). For piecewise constant basis functions, the collocation method is equivalent to the Nystrom method with a particular quadrature rule (the midpoint rule on each panel).

This superconvergence phenomenon, first discovered by Sloan and Chandler in the 1980s, is one of the pleasant surprises of the collocation method and partially explains its excellent performance in practice.

4.4 Galerkin Boundary Element Method

The Galerkin method is the most mathematically principled approach to discretizing boundary integral equations. It provides a variational framework that yields optimal error estimates and is the method of choice when rigorous error analysis is required.

The resulting linear system is \((\mathbf{M} - \mathbf{K})\boldsymbol{\alpha} = \mathbf{f}\), where the mass matrix \(\mathbf{M}_{ij} = \langle \phi_j, \phi_i \rangle\) and the stiffness matrix \(\mathbf{K}_{ij} = \langle T\phi_j, \phi_i \rangle = \int_\Gamma \int_\Gamma K(\mathbf{x}, \mathbf{y}) \phi_j(\mathbf{y}) \phi_i(\mathbf{x}) \, dS(\mathbf{y}) \, dS(\mathbf{x})\). The Galerkin matrix entries involve double integrals over the boundary, making them more expensive to compute than collocation or Nystrom entries but yielding superior error estimates.

The constant \(C\) in the error estimate depends on \(\|(I-T)^{-1}\|\), and can be large when the integral equation is poorly conditioned. This occurs, for instance, near irregular frequencies in the Helmholtz case or for domains with corners or edges where the solution has reduced regularity.

Chapter 5: Singular and Hypersingular Integrals

5.1 Cauchy Principal Value Integrals

The boundary integral operators introduced in Chapter 3 involve kernels with various types of singularities. The single-layer kernel is weakly singular (integrable), the double-layer kernel has a Cauchy-type singularity, and the hypersingular kernel is more singular still. The rigorous definition and practical computation of these singular integrals is one of the main technical challenges of the boundary element method.

The double-layer kernel \(\frac{\partial \Phi(\mathbf{x}, \mathbf{y})}{\partial \nu(\mathbf{y})}\) in three dimensions has a singularity of order \(|\mathbf{x} - \mathbf{y}|^{-2}\) times the surface measure, which makes the integral conditionally convergent. The Cauchy principal value is the correct interpretation: the contributions from opposite sides of the singularity cancel due to the antisymmetry of the kernel, and the limit exists for Holder-continuous densities.

5.2 Hadamard Finite-Part Integrals

The hypersingular operator \(N\) involves the kernel \(\frac{\partial^2 \Phi(\mathbf{x}, \mathbf{y})}{\partial \nu(\mathbf{x}) \partial \nu(\mathbf{y})}\), which has a singularity of order \(|\mathbf{x} - \mathbf{y}|^{-3}\) in \(\mathbb{R}^3\) — too singular even for the Cauchy principal value to exist. The correct interpretation is the Hadamard finite-part integral, introduced by Jacques Hadamard in his 1923 treatise on Cauchy’s problem.

The Hadamard finite part is not an integral in the Lebesgue sense; it is a distribution acting on smooth densities. Nevertheless, it defines a bounded operator between appropriate Sobolev spaces. The hypersingular operator \(N: H^1(\Gamma) \to L^2(\Gamma)\) is bounded, and its Galerkin discretization is well-defined because the test and trial functions in \(H^1(\Gamma)\) provide enough regularity to make the double integral convergent.

5.3 Regularization of Singular Kernels

A powerful approach to handling singular integrals in practice is regularization: rewriting the singular integral in a form that is amenable to standard quadrature. Several regularization techniques exist, each with its own advantages.

Maue’s identity transfers the hypersingularity from the kernel to the density through integration by parts on the surface. The resulting integrals involve only the weakly singular kernel \(\Phi_k\) and are amenable to standard quadrature techniques. This regularization is the basis for many practical BEM implementations. The price paid is that the density \(\varphi\) must now have surface gradients, which means the trial space must be at least \(H^1(\Gamma)\) (e.g., piecewise linear rather than piecewise constant). In the Galerkin framework, this is natural: both trial and test functions are taken from \(H^1(\Gamma)\), and the double surface integral converges in the classical sense after Maue’s regularization.

An alternative and widely used regularization in two dimensions exploits the fact that the logarithmic kernel can be decomposed into a smooth part and a standard singular part.

5.4 Numerical Quadrature for Singular Integrals

The practical computation of weakly singular, singular, and hypersingular integrals requires specialized quadrature methods. The choice of quadrature is one of the most critical implementation decisions in any boundary element code, as it directly determines both the accuracy and the efficiency of the method. We survey the main approaches, beginning with classical techniques and progressing to modern high-order methods.

For weakly singular integrals in 2D with a logarithmic singularity, the standard approach is the log-weighted Gaussian quadrature: one uses Gaussian rules for the weight function \(w(t) = -\log t\) on \([0,1]\), which integrate \((-\log t) p(t)\) exactly for polynomials \(p\) of degree up to \(2n - 1\). These rules were tabulated by Stroud and Secrest and are now computed to arbitrary precision using the Golub-Welsch algorithm.

5.5 Panel-Based Quadrature

For three-dimensional problems on general geometries, the boundary \(\Gamma\) is discretized into a mesh of panels (typically flat or curved triangles). Near-field interactions — where the source and target panels are close together or coincide — require careful treatment of the kernel singularity.

The Duffy transform is one of the most widely used tools in practical BEM implementations. For edge-adjacent and vertex-adjacent panel pairs, analogous coordinate transformations reduce the singular integrals to smooth integrals on the unit hypercube, though the number of required subdivisions increases with the type of adjacency. A self-interaction on a single triangle requires one Duffy transform (splitting into three sub-triangles, each with the singularity at a vertex), while an edge-adjacent pair requires four sub-integrals and a vertex-adjacent pair requires two.

Chapter 6: Fast Algorithms

6.1 The Need for Fast Algorithms

The boundary integral equation method reduces a PDE in \(d\) dimensions to an integral equation on a \((d-1)\)-dimensional boundary, but the resulting dense linear system of \(N\) equations has an \(O(N^2)\) storage cost and an \(O(N^3)\) cost for direct solution by Gaussian elimination (or \(O(N^2)\) per iteration of an iterative solver). For large-scale problems with \(N = 10^5\) to \(10^8\) unknowns, these costs are prohibitive. Indeed, storing a dense \(10^6 \times 10^6\) matrix of double-precision numbers requires approximately 8 terabytes of memory, far exceeding the capacity of any single workstation. Fast algorithms exploit the structure of the kernel to reduce the cost of matrix-vector products to \(O(N)\) or \(O(N \log N)\), making iterative solution feasible even for very large problems.

The key structural property is that the kernel \(K(\mathbf{x}, \mathbf{y})\) is smooth (and hence numerically low-rank) when \(\mathbf{x}\) and \(\mathbf{y}\) are well-separated. This means that the interaction between a cluster of sources and a well-separated cluster of targets can be accurately approximated by a low-rank factorization, and the cost of evaluating this interaction is proportional to the rank rather than the product of cluster sizes.

6.2 The Fast Multipole Method

The Fast Multipole Method (FMM), introduced by Greengard and Rokhlin in their landmark 1987 paper, is a systematic algorithm for exploiting the low-rank structure of far-field interactions. It was selected by the journal Computing in Science and Engineering as one of the “Top 10 Algorithms of the 20th Century.”

The FMM organizes the computation using a hierarchical tree decomposition of space (a quadtree in 2D, an octree in 3D). At each level of the tree, it computes multipole expansions for each box, translates them to local (Taylor) expansions at well-separated target boxes, and passes information between levels. The key operations are:

(a) Multipole-to-multipole (M2M) translation: shifting a multipole expansion from a child box center to its parent.

(b) Multipole-to-local (M2L) translation: converting a multipole expansion at a source box to a local expansion at a well-separated target box.

(c) Local-to-local (L2L) translation: shifting a local expansion from a parent box to its child.

The \(O(N)\) complexity is achieved because each source and target interacts with only a bounded number of boxes at each level of the tree, and the total number of translations is proportional to the number of boxes, which is \(O(N)\) for a quasi-uniform distribution.

The extension of the FMM to the Helmholtz equation requires different expansion types depending on the distance between source and target. At low frequencies (\(k \cdot \text{box size} \ll 1\)), the multipole expansion is similar to the Laplace case but with complex coefficients. At high frequencies, the multipole expansion must be replaced by a plane wave expansion or a diagonal translation operator, leading to the “wideband FMM” that handles all frequency regimes.

6.3 Hierarchical Matrices

An alternative framework for exploiting the low-rank structure of integral equation matrices is provided by hierarchical matrices (\(\mathcal{H}\)-matrices), developed primarily by Wolfgang Hackbusch and collaborators beginning in the late 1990s.

The \(\mathcal{H}\)-matrix framework provides not only fast matrix-vector products (in \(O(N \log N)\) or \(O(N)\) time) but also approximate matrix arithmetic: addition, multiplication, inversion, and LU factorization can all be performed in nearly linear complexity. This makes \(\mathcal{H}\)-matrices a powerful tool for constructing fast direct solvers, which can be far more robust than iterative methods for ill-conditioned problems.

6.4 The Proxy Point Method

The proxy point method is an elegant technique for compressing far-field interactions without explicitly computing analytic expansions. It is a key ingredient in modern fast direct solvers and kernel-independent FMM variants.

The proxy point method works because any potential generated by sources inside \(B\) that is observed at well-separated targets can be equivalently generated by an equivalent source distribution on the proxy surface. This is a consequence of the classical equivalence principle: from the perspective of well-separated observers, the detailed source distribution inside \(B\) is indistinguishable from a suitably chosen distribution on the proxy surface. The number of proxy points needed depends on the expansion order (accuracy) but not on the number of original sources, which is precisely what makes the method useful for compression. This makes the proxy point method particularly useful for kernel-independent algorithms that cannot exploit analytic expansions, including algorithms for kernels arising from oscillatory (Helmholtz), tensorial (Stokes, elasticity), or layered-medium Green’s functions.

6.5 O(N) and O(N log N) Solvers

Combining the fast matrix-vector product (via FMM or \(\mathcal{H}\)-matrices) with an iterative linear algebra method yields an \(O(N \log N)\) or \(O(N)\) solver for boundary integral equations. The typical workflow is:

(1) Discretize the integral equation to obtain a dense \(N \times N\) linear system \(A\mathbf{x} = \mathbf{b}\).

(2) Apply GMRES (or another Krylov subspace method) to solve the system iteratively, using the FMM to compute each matrix-vector product in \(O(N)\) time.

(3) The total cost is \(O(N \cdot n_{\text{iter}})\), where \(n_{\text{iter}}\) is the number of GMRES iterations.

This is one of the great advantages of second-kind integral equations over first-kind equations: the condition number of the discrete system is bounded as \(N \to \infty\), so iterative methods converge in \(O(1)\) iterations. Combined with an \(O(N)\) fast matrix-vector product, the total solution cost is \(O(N)\) — optimal complexity.

Chapter 7: Applications

The preceding chapters have developed the theoretical and computational foundations of integral equation methods. We now turn to applications, surveying the major areas where boundary integral equations have proven most effective. In each application, the key ingredients are: (i) a fundamental solution for the governing PDE, (ii) a boundary integral formulation that is well-posed and well-conditioned, and (iii) efficient numerical methods (quadrature, fast algorithms, iterative solvers) to handle practical problem sizes.

7.1 Acoustic Scattering

The prototypical application of boundary integral equations for the Helmholtz equation is acoustic scattering: the determination of how a sound wave interacts with an obstacle. The total field \(u = u^i + u^s\) consists of an incident field \(u^i\) (a known solution of the Helmholtz equation, typically a plane wave) and a scattered field \(u^s\) that satisfies the Helmholtz equation in the exterior domain, the boundary condition on the scatterer, and the Sommerfeld radiation condition at infinity.

Using the Burton-Miller formulation (Theorem 3.9), one obtains a well-posed second-kind integral equation for all wavenumbers. After discretization by a Galerkin or Nystrom method and acceleration by the FMM, the resulting system can be solved for problems involving millions of unknowns. State-of-the-art BIE solvers for acoustic scattering can handle problems with \(N > 10^8\) unknowns, corresponding to scattering objects hundreds of wavelengths in diameter. The far-field pattern (or scattering amplitude) is then obtained by evaluating the representation formula in the asymptotic regime.

The far-field pattern contains all the information about the scattering process that is observable at large distances. It is the primary quantity of interest in radar cross-section computation, sonar target classification, and non-destructive testing. The far-field pattern also determines the total scattering cross-section via the optical theorem:

\[ \sigma_{\text{total}} = \frac{4\pi}{k} \operatorname{Im}\left[ u_\infty(\hat{\mathbf{d}}) \right], \]

where \(\hat{\mathbf{d}}\) is the direction of the incident wave. This remarkable identity, which connects the total power scattered in all directions to the imaginary part of the forward-scattering amplitude, provides both a physical check on numerical solutions and a mathematically elegant consequence of energy conservation.

The computational cost of solving the acoustic scattering problem by BIE methods is governed by two factors: the number of unknowns \(N\) (which must be proportional to \(k^{d-1}\) to resolve the oscillatory surface field, where \(d\) is the spatial dimension) and the number of GMRES iterations (which grows with \(k\) for the standard second-kind formulation). Research on high-frequency methods seeks to reduce both costs. Hybrid asymptotic-numerical methods partition the surface into illuminated, shadow, and transition regions, using asymptotic approximations in regions where the field varies predictably and BIE discretization only where needed.

7.2 Stokes Flow

Boundary integral equations are particularly natural for Stokes flow (creeping or low-Reynolds-number flow), where the governing equations are linear and the free-space Green’s function is known explicitly. Applications include microorganism locomotion, microfluidics, and sedimentation of particles.

The Stokeslet has two key properties that distinguish it from the Laplace fundamental solution: it is a tensor (reflecting the vector nature of the velocity field), and it decays more slowly at infinity (\(O(1/r)\) versus \(O(1/r)\) for Laplace in 3D, but the Stokeslet has an additional algebraic structure that couples different components of the force and velocity). The associated pressure field \(p_j\) decays as \(O(1/r^2)\), faster than the velocity.

The boundary integral formulation for Stokes flow is entirely analogous to the scalar Laplace case, with the scalar kernel replaced by the Stokeslet tensor. The single-layer and double-layer Stokes potentials are defined as

\[ (\mathcal{S}\mathbf{f})_i(\mathbf{x}) = \int_\Gamma G_{ij}(\mathbf{x}, \mathbf{y}) f_j(\mathbf{y}) \, dS(\mathbf{y}), \quad (\mathcal{D}\boldsymbol{\mu})_i(\mathbf{x}) = \int_\Gamma T_{ijk}(\mathbf{x}, \mathbf{y}) \mu_j(\mathbf{y}) \nu_k(\mathbf{y}) \, dS(\mathbf{y}), \]

where the stress tensor kernel \(T_{ijk}\) (or “stresslet”) is

\[ T_{ijk}(\mathbf{x}, \mathbf{y}) = -\frac{3}{4\pi} \frac{r_i r_j r_k}{|\mathbf{r}|^5}. \]

One of the most striking applications of Stokes BIE methods is to the simulation of suspensions of many particles. When \(M\) rigid particles are suspended in a Stokes flow, the boundary integral formulation involves surfaces of all \(M\) particles simultaneously. With FMM acceleration, the cost per time step is \(O(M \cdot N_p)\), where \(N_p\) is the number of unknowns per particle, enabling simulations with thousands to millions of particles. This capability has transformed the study of blood flow (red blood cells as deformable capsules), sedimentation, and microfluidic sorting.

7.3 Electromagnetic Scattering

The boundary integral formulation for Maxwell’s equations is considerably more complex than for scalar equations, owing to the vector nature of the fields and the need for surface vector calculus. Nevertheless, boundary integral methods are among the most effective tools for computational electromagnetics, particularly for scattering by perfectly conducting or dielectric bodies.

The free-space dyadic Green’s function for Maxwell’s equations is constructed from the scalar Helmholtz Green’s function:

\[ \overline{\overline{G}}(\mathbf{x}, \mathbf{y}) = \left(I + \frac{1}{k^2}\nabla\nabla\right) \Phi_k(\mathbf{x}, \mathbf{y}). \]

This leads to the electric field integral equation (EFIE) and the magnetic field integral equation (MFIE) for the surface current \(\mathbf{J} = \hat{\nu} \times \mathbf{H}\) on a perfect conductor.

The EFIE is a first-kind equation (and hence ill-conditioned), while the MFIE is second-kind but only valid for closed surfaces. For open surfaces (such as thin plates or slot antennas), only the EFIE applies. The combined field integral equation (CFIE), a weighted average of the EFIE and MFIE, is both uniquely solvable at all frequencies and applicable to closed surfaces. The standard CFIE takes the form

\[ \alpha \cdot \text{EFIE} + (1 - \alpha) \cdot \text{MFIE}, \]

with the conventional choice \(\alpha = 0.5\). The CFIE is the electromagnetic analog of the Burton-Miller formulation for the Helmholtz equation (Theorem 3.9), and it eliminates the interior resonance problem that plagues both the EFIE and MFIE individually. The conditioning of the CFIE grows only as \(O(k)\) with frequency, compared to \(O(k^2)\) for the EFIE, making it far more amenable to iterative solution.

7.4 Elasticity

Linear elasticity is another classical domain of boundary integral equations, with a long history going back to the work of Betti (1872) and Somigliana (1886). The boundary integral equation method extends naturally to the equations of linear elastostatics and elastodynamics. The fundamental solution for the Navier-Cauchy equations of linear elasticity (the Kelvin solution) plays the same role that the Laplace fundamental solution plays for potential theory.

The structure of the elasticity boundary integral equation is identical in spirit to the Stokes case (indeed, the Stokes equations are the incompressible limit \(\nu \to 1/2\) of the elasticity equations). The single-layer potential involves the Kelvin solution \(U_{ij}\), and the double-layer kernel involves the associated traction tensor.

The Somigliana identity expresses the displacement at interior points in terms of boundary displacements and tractions. For the interior Dirichlet problem (prescribed displacement), the unknown traction is determined by a boundary integral equation obtained by taking the limit \(\mathbf{y} \to \Gamma\) and applying the jump relations for the elasticity layer potentials.

7.5 Comparison with FEM and FDM

To conclude these notes, we briefly compare the boundary integral equation / boundary element method (BIE/BEM) with the two dominant volume discretization approaches: the finite element method (FEM) and the finite difference method (FDM). This comparison is not about declaring a winner — each method has its natural domain of applicability — but about understanding which tools are best suited to which problems.

The principal advantage of the BIE/BEM approach is dimensional reduction: only the boundary \(\Gamma\) is discretized, reducing a 3D problem to a 2D one. This is particularly beneficial for exterior problems in unbounded domains, where volume methods require artificial truncation of the domain and absorbing boundary conditions. BIE methods automatically enforce the correct behavior at infinity (radiation condition, decay conditions) through the choice of fundamental solution.

A second important advantage of BIE methods is the high accuracy achievable on the boundary itself. Since the integral equation is solved directly on \(\Gamma\), boundary quantities (such as surface stresses, heat fluxes, or surface currents) are computed as primary unknowns rather than being recovered by differentiation of a volume solution. This avoids the loss of accuracy inherent in numerical differentiation and makes BIE methods the preferred choice for problems where boundary quantities are the primary interest, such as fracture mechanics (stress intensity factors) and electromagnetic compatibility (surface currents and charges).

However, the BIE/BEM approach has significant disadvantages as well. The system matrix is dense (\(O(N^2)\) storage for \(N\) unknowns, compared to \(O(N)\) for sparse FEM matrices), although fast algorithms reduce this to \(O(N)\) or \(O(N \log N)\). The method is restricted to linear, homogeneous-medium problems for which a fundamental solution is known (or can be efficiently computed). Non-linear problems, variable coefficients, and multi-physics coupling are difficult to handle. Finally, the implementation of BEM codes is considerably more complex than FEM codes, owing to the singular integrals, specialized quadrature, and fast algorithms required.