Lecture 2: Integral equations/quadratures

Previous lecture

• Intro
• IE for external problems

Todays lecture

• IE discretization
• Nystrom, collocation, Galerkin method

Simplest integral equation

The simplest integral equation we talked about reads

$$\int_{\partial \Omega} \frac{q(y)}{\Vert x - y \Vert} dy = f(x), \quad x \in \partial \Omega,$$

where $\Omega$ is a certain domain in 3D.

The general form of a first-kind integral equation is

$$ \int_{\partial \Omega} q(y) G(x, y) dy = f(x), \quad x \in \partial \Omega, $$

where $G(x, y)$ is a Green function of the PDE.

Typical Green functions

The Green function of a PDE operator $A$ satisfies $$A G(x, y) = \delta(x-y), $$ i.e. it plays the role of the "inverse" of the operator without boundary conditions.

It is also called free-space Green function. The Green function can also incorporate the boundary conditions.

Green functions with boundary conditions and reflection method

Suppose we want to find a function $G(x, y)$ such that

$$\Delta_x G(x, y) = \delta (x - y), \; G(x, y) \mid_{x \in \partial \Omega} = 0. $$

If we found such a function, solution of

$$\Delta u = f,$$

is given by the explicit integral

$$u(x) = \int_{\Omega} G(x, y) f(y) dy.$$

Reflection method

In a few cases, you can write down such a Green function, typically based on the method of reflections.

For a Dirictlet bouundary conditions, we can use not only charge located at $y$, but also many other charges that lie outside $\Omega$ and give the required b.c.

• A plane - two charges are ok.
• A cube - many reflections, infinite series

Discretization

Now to discretization of the equation $$\int_{\partial \Omega} \frac{q(y)}{\Vert x - y \Vert} dy = f(x), \quad x \in \partial \Omega.$$

Steps of discretization:

1. Select a finite-dimensional subspace, spanned by basis function $\phi_i, i = 1, \dots, N$. Typically, $\phi_i$ are local basis functions defined on a certain mesh $\mathcal{T}$ over the surface.
2. Look for the solution in the form $$q(y) \approx \sum_{i=1}^N q_i \phi_i(y)$$

Note, that the basis functions can be even piecewise-constant, compared to PDEs, where generally they have to be at least piecewise-linear for Finite Element Method (they can be piecewise-constant for Discontinious Galerkin, but that is more an exception)

Discretization (2)

$$\int_{\partial \Omega} \frac{q(y)}{\Vert x - y \Vert} dy = f(x), \quad x \in \partial \Omega,$$

Suppose that $\mathrm{supp}(\phi_i) = T_i$ (i.e., local support). Then,

$$ \sum_{i=1}^N q_i \int_{T_i} \frac{\phi_i(y)}{\Vert x - y \Vert} dy = f(x).$$

This has to hold for all $x$, i.e. it is infinite amount of equations. There are two ways of dealing with this situation:

1. Galerkin method
2. Collocation method.

Elements $T_i$ are elements of mesh $\mathcal{T}$. Usually they are triangles, but other configuration cam be used.

Galerkin method

The Galerkin method takes the equation $$ \sum_{i=1}^N q_i \int_{T_i} \frac{\phi_i(y)}{\Vert x - y \Vert} dy = f(x),$$ and multiplies it by $\phi_j(x)$ and integrates over $T_j$. Then we get

$$ \sum_{i=1}^N q_i \int_{T_j} \int_{T_i} \frac{\phi_j(x) \phi_i(y)}{\Vert x - y \Vert} dx dy = \int_{T_j} f(x) \phi_j(x) dx. $$

This is a square system of linear equations of the form $$ A q = f, $$ where $$ A = [a_{ji}], \; a_{ji} = \int_{T_j} \int_{T_i} \frac{\phi_j(x) \phi_i(y)}{\Vert x - y \Vert} dx dy, \quad f_j = \int_{T_j} f(x) \phi_j(x) dx $$
Everything is fine, the matrix is symmetric (it can be also shown that it is positive definite, let us do that), but the evaluation of the element requires 4D integral to be evaluated. That is main computational problem, thus engineers often use a much more simple collocation method.

Collocation method

Collocation method consists of taking the equation $$ \sum_{i=1}^N q_i \int_{T_i} \frac{\phi_i(y)}{\Vert x - y \Vert} dy = f(x),$$ at some collocation points $x_j$.

Consider piecewise-constant case, i.e. $\phi_i(y) = 1$ on $T_i$ and $0$ otherwise. Then the most natural choice for collocation points are centers of the triangles, and we again get a square linear system

$$A q = f,$$

where

$$ A = [a_{ji}], \; a_{ji} = \int_{T_i} \frac{dy}{\Vert x_j - y \Vert}, \quad f_j = f(x_j) $$

i.e. it has the meaning of uniformly charged triangle.

Note that this matrix is non-symmetric, so you can not use symmetric-oriented methods to solve linear system, but integrals are in 2D domain.

Also if number of collocation points is more than number of triangles $T_i$, then we have to solve least squares problem, because linear system becomes rectangular and underdetermined.

Collocation method as Petrov-Galerkin method

• In Galerkin method, function $f(x)$ is projected onto local basis functions $\phi_j(x)$: $$ f_j = \int_{T_j} f(x) \phi_j(x) dx $$
• In collocation method, function $f(x)$ is projected onto the delta functions $\delta(x - x_j)$, where $x_j$ are collocation points: $$ f_j = \int_{T_j} f(x) \delta (x - x_j) dx = f(x_j) $$

Nystrom method

The Nystrom method is just another simplication of the collocation method (as collocation method can be considered as a quadrature applied to the Galerkin method).

We just approximate

$$ A_{ij} = \int_{T_i} \frac{dy}{\Vert x_j - y \Vert} dy \approx S_i \frac{1}{\|x_j - y_i\|}. $$

If the source triangle $T_i$ is sufficiently far from the receiver $x_j$, indeed, this can be a very good approximation.

If $y_i$ lies inside the triangle, then the approximation can be really bad.

Locally corrected Nystrom

Typically, a locally corrected Nystrom is used: much more advanced quadratures are used in the close zone, where the singularity of the integrals appears. The concept is as follows:

• To evaluate the matrix-by-vector product $Aq$ it is split into the far part and close part.
• Close part is evaluated separately, far is evaluated using a simple quadrature

Far-close splitting is the key to fast methods!

Computing the integrals

The main art about Galerkin/collocation methods is computation of singular integrals, which can not be handled by standard Gauss quadratures in the close zone.

There are several ways

• Analytic evaluation
• Subtraction of singularity, singularity is integrated analytically, remainder by Gauss
• Gauss-Duffy quadrature
• Schwab-Sauter quadrature
• Rokhlin generalized quadratures

Model integral

We will study the integrals of the form

$$ \int^b_a \int^d_c \frac{\phi(x) \phi(y)}{\Vert x - y \Vert} dy, $$

where $\phi(x), \phi(y)$ are standard local basis function (constant or linear).

Computing the simplest integral

We can actually compute the anti-derivative of the integral in question, by 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)

Smarter than Nystrom

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.

Such representation is obviously useful, since in Galerkin formula we can replace 4D integration by 2D integration,

$$ \int \int \phi_i(x) \phi_j(y) u_{\alpha}(x) u_{\alpha}(y) dx dy = \int \phi_i(x) u_{\alpha}(x) dx \int \phi_j(y) u_{\alpha}(y) dy. $$

But how to compute these functions $u_{\alpha}$?

Exponential sums

A very useful approach is based on so-called exponential sums. A simple identity has the form

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

It can be very useful in separating variables.

A quadrature

Suppose we have a quadrature for the integral:

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

Then put

$$ r = x^2 + y^2 $$

and we get the sum of Gaussians on the right.

Sum of Gaussians

$$\frac{1}{\sqrt{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 how to find $w_k$ and $p_k$?

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 have the new integral

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

and this integral is then approximated 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

Integration over triangles

Exponential sums are well-applicable to the integration over rectangles. But the most typical choice is the integration over triangles.

Consider the standard triangle $$ T = \{ (x, y) \; | \; 0 \leq x \leq 1, \; 0 \leq y \leq x\}. $$

Consider the integral $$\int_{T} h dS,$$

which can be written as

$$ \int^{1}_{0} dx \int^x_0 h(x, y)dy. $$

Duffy transformation

Duffy transformation is a mapping from rectangle to triangle, given by formulas

$$x = x, \quad y = x u, \quad 0 \leq u \leq 1,$$

i.e. it is a mapping from $(x, u)$ to $(x, y)$.

Duffy transformation(2)

Then, we have

$$\int^{1}_0 \int^{1}_0 h(x, xu) x dx du.$$

Recall the singularity:

$$h(x, y) = \frac{g(x, y)}{\sqrt{x^2 + y^2}},$$

therefore the $x$ cancels out and we have

$$\int_0^1 \int_0^1 \frac{g(x, xu)}{\sqrt{1 + u^2}} dx du, $$

which is non-singular.

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
• Vico-Greengard-Ferrando quadrature

Vico-Greengard-Ferrando quadrature

Suppose we have to compute $$f(x) = \int G(x-y) \rho (y) dy, $$ where $G(x-y)$ is a spherically symmetric Green function of some differential operator, that is $$G(r) = G(|r|)$$

Idea: use Fourier transform $$\widehat{(f \ast g)} = {\hat f} {\hat g}$$

This suggests that we can compute $f(x)$ as $$f(x) = \Big(\frac{1}{2\pi}\Big)^3 \int e^{i s \cdot x} \widehat{\rho}(s) \widehat{G}(s) ds$$

However, there is a problem: $\widehat{G}(s)$ is typically very singular, e.g. for Helmholtz differential operator $\Delta + k^2$ $$G(s) = \frac{1}{k^2 - s^2},$$ there is a sphere of singular points. So how we can avoid it?

Vico-Greengard-Ferrando quadratures (2)

We can use the following simple observation. Suppose that $\rho(x)$ is $0$ outside of cube $[0, 1]^3$. Then, in this cube we have the following identity:

$$ \int G(x-y) \rho (y) dy \equiv \int G_L(x-y) \rho (y) dy, $$

with $$ G_L(r) = G(r) \chi_{[0, L]}(|r|), \quad L = \sqrt{3} $$

and $\chi_D(x)$ is a charactersitic function of a set $D$, and identity holds because in the cube $[0, 1]^3$ distance between any two points doesn't exceed $\sqrt{3}$.

So we can use $G_L$ instead of $G$. $G_L$ has compact support and it's Fourier transform is smooth!

E.g. for Helmholtz differential operator:

Vico-Greengard-Ferrando quadratures (3)

So we need to compute $$f(x) = \Big(\frac{1}{2\pi}\Big)^3 \int e^{i s \cdot x} \widehat{\rho}(s) \widehat{G_L}(s) ds, $$ with $G_L(s)$ being smooth. This can be performed using trapezoidal rule.

Rokhlin generalized quadratures

Finally, I will discuss a general quadrature idea. For a particular kernel of the integral operator, the quadrature points and weights in

$$\int K f \approx \sum_i w_i (K f)(x_i),$$

have to be defined only once.

Idea in 1D

I will illustrate the idea in 1D to compute

$$\int^1_0 f(x) dx \approx \sum_{j=1}^M f(x_j) w_j$$

and this quadrature has to integrate certain subset of functions $f_1, \ldots, f_n$,

i.e.

$$\sum_{j=1}^M f_i(x_j) w_j = \int^1_0 f_i(x) dx.$$

Assumptions

The quadrature formula with nodes $x_1, \ldots, x_m$ and weights $w_1, \ldots, w_m$:

• Integrates $f_i, i = 1, \ldots, n$ and their products
• $f_i$ form orthonormal basis set in $L_2[0, 1]$
• $m \gg n$

Linear system for the weights

Then we have a linear $n \times m$ system for the weights, which is underdetermined.

Specifically,

$$\sum_{j=1}^m f_i(x_j) \sqrt{w_j} v_j = b_n.$$

Since the quadrature formula integrates the products, the matrix has singular values $0$ or $1$ (prove!).

We can also set most of the solution coefficients to zero, and leave only $m$ non-zero entries.

If we denote the positions by $i_1, \ldots, i_n$, the final quadrature will read

$$\int^1_0 f(x) dx \approx \sum_{k=1}^n f(x_{i_k}) z_{k} \sqrt{w_{i_k}}.$$

Summary of the generalized quadrature idea

• Start with an orthonormal basis set
• Select very large quadrature with many points
• Select a "pruned" quadrature.

The approach can be generalized to singular integrals.

Summary of the lecture

• Basic discretization schemes (Galerkin, collocation, Nystrom)
• Some approaches to compute singular integrals

Next lecture

• Fast methods: computation of the convolution via FFT, idea of precorrected FFT.
• More integral equation kernels