Lecture 10. Integral equations: discretization

Previous lecture

• Introduction into the integral equations: the basic one, the concept of non-locality
• Exterior and interior problems

Todays lecture

• How to discretize: method of moments (MoM) and the Galerkin method
• Computation of the singular integrals and exponential sums
• Fast methods: the idea of the Barnes-Hut method

Books

Books on integral equations

• Integral equation methods in scattering theory, D. Colton R. Kress
• Sauter, Schwab, Boundary element methods (2011)

Method of moments (collocation method)

Recall the basic integral equation from electrostatics:

$$\int_{\Omega} \frac{q(y) dy}{\Vert x - y \Vert} = f(x).$$

1. First we mesh the domain in simplices.
2. Then, we introduce the discrete solution space (piecewise-linear, piecewise-constant functions).
3. Then, we approximate the solution as a linear combination of the basis functions $$ q(y) = \sum_i q_i \phi_i(y).$$

MoM (continued)

Then, we get the system of equations $$ f(x) = \sum_i q_i \int_{\Pi_i} \frac{\phi_i}{\Vert x - y \Vert } dy. $$ The method of moments is obtained by the selection of collocation points (aka staggerred grids), typically in the middle of triangles/rectangles

Piecewise-constant basis functions

For a piecewise constant basis functions, we have

$$\int_{\Pi} \frac{1}{\Vert x - y \Vert} dy, $$

where this is a two-dimensional integral.

Computing the simplest integral

We can actually compute the anti-derivative of the integral in question, but taking successive integrals

$$ \frac{1}{\sqrt{x^2 + y^2}} $$\begin{equation} \begin{split} F(x, y) = \int \int \frac{1}{\sqrt{x^2 + y^2}} dx dy = \\ -y + y \log(x + \sqrt{x^2 + y^2}) + x \log(y + \sqrt{x^2 + y^2}), \end{split} \end{equation}

and the final integral is just $$ I(x_0, y_0, x_1, y_1) = F(x_0, y_0) + F(x_1, y_1) - F(x_0, y_1) - F(x_1, y_0), $$

(generalization of the Newton-Leibniz formula)

Nystrom method

A Nystrom method is a popular choice: replace the integral by quadrature

probably by the simplest one: the rectangle formula,

which in 1D has the form

$$I(f) = \sum_{k=0}^n w_k f(x_k), $$

where $x_k$ is the uniform grid on the interval, $w_0 = w_n = h/2, \quad w_k = h$ for all other $k$.

For our equation, it looks like

$$A_{ij} = \frac{1} {\Vert x_i - y_j \Vert} h^2,$$

the set $x_i$ can be called receivers and the set $y_j$ can be called sources ;

A typical way is to shift the grid by half-step, i.e. in 1D it will be

$$ A_{ij} = \frac{1}{|i-j| + \frac{1}{2}}. $$ By the way, do you remember, what is the rank of this matrix?

Galerkin method

1. Collocation method gives non-symmetric matrices, even if the original operator was symmetric;
2. Nystrom method adds logarithmic factor into the accuracy estimate
3. Galerkin method is the method of choice in many cases

Galerkin method

A Galerkin method is defined for an arbitrary linear (and even non-linear) equations

$$ A u = f, $$

you select the test space and trial space .

For simplicity, let them be equal,

i.e. $$ u = \sum_i u_i \phi_i, $$

and

$$A_{ij} = (A \phi_i, \phi_j) $$

Integrals in the Galerkin method

We have to compute 4D integrals (here it is written for piecewise-constant functions)

$$ A_{ij} = \int_{\Omega_i} \int_{\Omega_j}\frac{1}{\Vert x - y \Vert} dx dy $$

Analytic expression are cumbersome!

Semi-numerical method

You can be a little bit smarter: approximate $1/r$ by something that can be well-integrated.

By a sum of separable functions!

$$\frac{1}{\Vert x - y\Vert } \approx \sum_{\alpha=1}^r u_{\alpha}(x) u_{\alpha}(y)$$

We have used the symmetry.

But how to compute these functions $u_{\alpha}$? (why it is useful?)

Exponential sums

There is an identity:

$$\frac{1}{\sqrt{x}} = \frac{1}{\sqrt{\pi}} \int^{\infty}_0 \frac{e^{-px}}{\sqrt{p}} dp.$$

It can be very useful in separating variables.

A quadrature

Suppose we have a quadrature for the integral:

$$ \frac{1}{\sqrt{x}} = \frac{1}{\sqrt{\pi}} \int^{\infty}_0 \frac{e^{-px}}{\sqrt{p}} dp. \approx \sum_{k} w_k e^{-p_k x}.$$

Then put

$$x := x^2 + y^2$$

and we get the sum of Gaussians on the right.

Sum of Gaussians

$$\frac{1}{r} \approx \sum_{k} w_k e^{-p_k (x^2 + y^2)},$$

and the integrals reduce to the 1D integrals

$$ \int_a^b e^{-p_k x^2} dx, $$

and this can be expressed via the Error functions (and there is a standard numpy erfc function)

But who you can approximate?

We reduced the problem to the problem of creating a good quadrature rule for the integral

$$\int^{\infty}_0 \frac{e^{-px}}{\sqrt{p}} dp$$

You have to do the $p = e^t$ (why do we do so?) and we the new integral

$$ \int^{\infty}_{-\infty} e^{-e^{t} x + t/2} dt $$

and this integral is then approximate by trapezoidal rule $t_k = a_t + k h_t, \quad h_t = (b_t - a_t)/K$

The integral is doubly exponentially decaying, and also its Fourier transform is decaying exponentially fast.

Thus, trapezoidal rule has exponential convergence in this setting

A short demo...

A sidenote about fitting by sum of exponentials

• Approximation of a given function by a sum of exponentials is a ill-posed problem
• Original method by Prony (1795)
• Recent work by G. Beylkin, et. al, Approximation by exponential sums revisited

Efficient quadratures (papers)

• C. Schwab, S. Sauter Efficient automatic quadrature in 3-d Galerkin BEM
• A. Polimeridis, Traianos V Yioultsis On the direct evaluation of weakly singular integrals in Galerkin mixed potential integral equation formulations

Green functions and fundamental solutions

Terminology:

• A fundamental solution is a "point charge" in free space: no boundary conditions
• Green function can be for the boundary conditions/layered media/etc. (i.e. "inverse" of the differential operator)

Fundamental solution: Helmholtz

The Helmholtz equation has the form

$$\Delta u + k^2 u = f,$$

and it describes the propagation of scalar waves (acoustic wave).

$k$ is called wavenumber.

The fundamental solution (the one that gives $\delta$-function after the application of the operator) has the form

$$K(r) = \frac{\exp(i k r)}{r}.$$

So it oscillates.

Laplace in 2D

What is the fundamental solution to the Laplace equation in 2D?

$$G(x, y) = \log \Vert x - y \Vert.$$

It plays a role in the theory of conformal mappings

Helmholtz in 2D is solved by the Hankel function

Stokes problem

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 p = 0 \end{split} \end{equation}

Maxwell equations

$$ \begin{split} \nabla \times H = \frac{\partial D}{\partial t} + J, \\\quad \nabla \times E = -\frac{\partial B}{\partial t}, \quad \nabla \cdot D = \rho, \quad \nabla \cdot B = 0. \\ D = \varepsilon E, \quad H = \mu B. \end{split} $$

How to make it fast

The memory is $\mathcal{O}(N^2)$, the complexity is $\mathcal{O}(N^3)$.

Need to make it $\mathcal{O}(N \log N)$ both in memory and complexity of matrix-by-vector product.

Matrix structure

The integral operator matrix:

1. Matrix is dense
2. Matrix is not of low-rank
3. Some blocks are of low rank!

A demo

Off-diagonal blocks correspond to "far" interaction

$$y_i = \sum_{j} a_{ij} q_j,$$

and the sources are separated from the receivers.

This gives a hint: separated geometrically -- low-rank property.

Geometric separability and low rank

For $\frac{1}{r}$ we can have the exponential sums to separate variables, if $r > r_0$.

What happens, if $r_0$ goes to zero?

The separability property vanishes!

One-dimensional case

$$A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22}\end{bmatrix},$$

where $A_{12}$ and $A_{21}$ correspond to the "far" interaction, and are of low-rank.

Remember, that a low-rank matrix can be factored as

$$A = U V^{\top},$$

where $U$ is $n \times r$, $V$ is $m \times r$,

Memory estimate

Thus we have

$$Memory(n) = 2 Memory(n/2) + 2 n/2 r$$

and that gives

$Memory(n) = 2 n/2 * r + 4 n/4 * r + \ldots = \mathcal{O}( n \log n)$

Similar estimates hold for the matrix-by-vector product procedure (why do we care about that?)

Two-dimensional case

In two dimension, the splitting is different.

Plotting the singular values of the off-diagonal block shows that they do not decay that fast.

Why? Because they correspond to the non-separated sets of sources and receivers.

The general scheme

• Construct trees for the sources and receivers
• Check the sets for geometric separability.
• If they are not separated, call the recursion
• Otherwise, do low-rank approximation.

This is the basis for the Barnes-Hut method (original work in astronomy)

and the preliminary approach to the Fast Multipole Method , proposed by Greengard and Rokhlin