The main problem with integral equation methods (we talked about boundary integral equations, but there are also volume integral equations is that we have to deal with non-local interactions.
After discretization, we have to compute the sum
$$v_i = \sum_{j=1}^N A_{ij} q_j,$$and all $A_{ij}$ are not small (i.e., they can not be approximated by sparse matrices).
If we have two-dimensional PDE, the integral equation will be defined on a curve.
For a second-order discretization, 1000 elements are needed, i.e. the matrix will be
$10^{3} \times 10^{3}$, and storage is just 8 megabytes, i.e. fast methods are not very needed.
For 3D problems and surface integral equations, the storage problem becomes much more complicated:
For $h = 10^{-3}$ we have $10^6$ unknowns and the matrix is $10^6 \times 10^6$, more than 8 Terabytes.
One of the ways to deal with the storage problem is to solve the system using matrix-by-vector product.
Apply iterative method, and compute
$$v_i = \sum_{j=1}^N A_{ij} y_j $$on the fly.
The complexity, however, stays $\mathcal{O}(N^2)$ as well.
Moreover, the evaluation of matrix elements can be expensive, especially if Galerkin/collocation methods are used.
Consider the integral equation on a square:
$$\int_{\Omega} \frac{q(y)}{\Vert x - y \Vert} dy = f(x), $$Let us discretize using shifts of the same basis functions,
$$q(y) = \sum_{ij} \phi_{ij} q_{ij}, $$where $$\phi_{ij}(y_1, y_2) = \phi( y_1 - h i, y_2 - h j).$$
What would be the structure of the matrix?
Recall from numerical linear algebra, that BTTB matrix can be multiplied by a vector in $\mathcal{O}(N \log N)$ operations.
The idea is a generalization of the idea for one level Toeplitz matrix: any $N \times N$ matrix can be embedded into a $2N \times 2N$ circulant matrix.
For two-level model problem on a square, the matrix is not Toeplitz, but block Toeplitz with Toeplitz blocks.
It can be embedded into block circulant with circulant blocks, which can be diagonalized in $\mathcal{O}(N)$ operations using two-dimensional FFT:
$$F_2 = F \otimes F.$$For a translation-invariant case with a square uniform grid we have a matrix with BTTB structure, which can be:
But if the domain is not a square (i.e. a sphere or a complicated domain) using Toeplitz structure becomes very difficult.
What to do? The idea of precorrected FFT comes to mind, where 2D structure is embedded into a regular 3D grid.
Consider a 3D parallelepiped splitted into $k \times l \times m$ array of small cubes in such a way that each cube contains only a small number of panels.
The key idea is that for all evaluation points, distant for the particular cell, it can be represented as a weighted sum of point charges located on a uniform grid.
In a cell, we want to match the potential, created by $m$ panels, by a uniform grid of $p \times p \times p$ point charges. Then, we just select a test surface of radius $r_c$, select test points $z_i$, $i=1, \ldots, N_{test}$ and require that
$$\psi_{uniform}(z_i) \approx \psi_{panels}(z_i), \quad i=1,\ldots,N_{test}$$on the test points.
We have $p^3$ unknowns and $N_{test}$ equations, and it is solved as a least squares.
To get the value at evaluation point, it is sufficient to use the transpose matrix to the matrix that interpolates the charges from a non-uniform grid to a uniform grid.
The projection/interpolation have similar accuracy.
For collocation, the result can be written as
$$\psi_{c} = V^{\top} H W q, $$and for Galerkin method:
$$\psi_{G} = W^{\top} H W q,$$i.e. pFFT preserves the symmetry of the Galerkin method.
The main diffuculty is that the projection-FFT-interpolation step does not accurately represent the nearby interactions.
Let $(k, l)$ be a pair of nearby cells.
Then,
$$\psi_{c (k, l)} = V^{\top}(k) H(k, l) W(l) q_l,$$adding and subtracting the wrong contribution we have
$$\psi_{(k, l)} = \psi_{c (k, l)} - V^{\top}(k) H(k, l) W(l) q_l + P_{k, l} q_l.$$The simplicity of pFFT is that it can be applied very easily to other integral equation kernels with a similar efficiency, and few additional programming.
A disadvantage is that for simple sufraces many of the panel charges will be zero, i.e. the representation can be redundant.
Finally, let us present some other integral equation kernels that are often used.
Stokes equation define the laminar flow in low Reynolds number regime:
\begin{equation} \begin{split} - \nabla p + \mu \Delta v = -F \delta(r) \\ \nabla \cdot v = 0 \end{split} \end{equation}Here $p$ is pressure, $v$ is velocity.
The solution is so-called Stokeslet, and it is a matrix:
$$\mathbb{J}(r) = \frac{1}{8 \pi \mu} \left( \frac{I}{r} + \frac{r_i r_j}{r^3} \right), \quad u = F J , \quad p = \frac{(F, \mathbf{r})}{r^3}.$$If we consider time-harmonic case, i.e. the dependence from time is $e^{i w t}$, we get time-harmonic Maxwell equations.
However, for Maxwell equations the most natural setting is scattering from inhomogenious media, and volume integral equations (or volume integro-differential) equations are used instead.
We will illustrate how the BEM++ package works, which is quite general package for solving boundary integral equations (unfortunately, installation can be painful).
In [2]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/custom.css", "r").read()
return HTML(styles)
css_styling()
Out[2]: