Derivation and analysis of computational methods for fractional Laplacian equations with absorbing layers

Abstract

This paper is devoted to the derivation and analysis of accurate and efficient perfectly matched layers (PMLs) or efficient absorbing layers for solving fractional Laplacian equations within initial boundary value problems (IBVP). Two main approaches are derived: we first propose a Fourier-based pseudospectral method, and then present a real space method based on an efficient computation of the fractional Laplacian with PML. Some numerical experiments and analytical results are proposed along the paper to illustrate the presented methods.

Appendices

Appendix A: Padé approximant-based PML

In the following, we detail the procedure for deriving PMLs by using Padé’s approximants with α = p/2^k, $k \in {\mathbb {N}}^{*}$ and $p\in {\mathbb {N}}^{*}$.

Case α = 1/2^k, $k \in {\mathbb {N}}^{*}$.:

The idea developed for α = 1/2 can easily be extended to coefficients of the form 1/2^k. We can iteratively repeat the process described above, by simply using:

$$ \begin{array}{@{}rcl@{}} (-\triangle_{\text{PML}})^{1/2^{k}} = \sqrt{(-\triangle_{\text{PML}})^{1/2^{k-1}}} . \end{array} $$

$$ \begin{array}{@{}rcl@{}} \text{Op}\Big(\sqrt{\sigma\big((-\triangle_{\text{PML}})^{1/2^{k}}\big)}\Big) \approx \text{Op}\Big(\sum\limits_{k=0}^{M}a_{k}^{(M)} - \sum\limits_{k=1}^{M}\frac{a_{k}^{(M)}d_{k}^{(M)}}{\sigma\big((-\triangle_{\text{PML}})^{1/2^{k-1}}\big) + d_{k}^{(M)}} \Big) . \end{array} $$

This leads to long calculations, which however have to be done once for all for any given α = 1/2^k for $k \in {\mathbb {N}}^{*}$.

Case $\alpha \in {\mathbb {N}}^{*}/2^{k}$, $k\in {\mathbb {N}}^{*}$.:

We extend the above ideas to rational numbers α in the form p/2^k, for $p\in {\mathbb {N}}^{*}$. In fact, thanks to the above discussion, we simply need to detail the case α = p/2, for $p \in {\mathbb {N}}^{*}$. Although the expressions look quite complex, in practice, simplifications and approximations are possible:

(61)

In the above system, we assume that designates some smooth real- or purely complex-valued functions and that is a finite set of strictly positive numbers in ${\mathbb {N}}^{*}/2$. We then consider the corresponding IBVP:

(62)

We can formally rewrite the symbol of (−△_PML)^α [1] as:

$$ \begin{array}{@{}rcl@{}} \sigma\big((-\triangle_{\text{PML}})^{p/2}\big) = \sigma\big(\sqrt{\!-\triangle_{\text{PML}}}^{p}\big) = \sigma\big(\sqrt{\!-\triangle_{\text{PML}}}\big) \# \sigma\big(\sqrt{\!-\triangle_{\text{PML}}}\big) {\cdots} \# \sigma\big(\sqrt{-\triangle_{\text{PML}}}\big) , \end{array} $$

where we recall that:

Proposition 11

[1] For two pseudodifferential operators A and B with $C^{\infty }$-coefficients, α = (α₁,α₂) with |α| = α₁ + α₂ and α! = α₁!α₂!, the symbol to the composed operator AB is given by:

$$ \begin{array}{@{}rcl@{}} \sigma(AB) = \sigma(A) \#\sigma(B) \sim \sum\limits_{|\alpha|=0}^{\infty}\frac{(-\texttt{i})^{|\alpha|}}{\alpha!}\partial^{\alpha_{1}}_{x}\partial^{\alpha_{2}}_{y}\sigma(A)\partial^{\alpha_{1}}_{\xi_{x}}\partial^{\alpha_{2}}_{\xi_{y}}\sigma(B) . \end{array} $$

From a practical point of view, the computation of these symbols and the approximation of the corresponding operators can be complex. Instead, we can proceed as follows:

• If $p \in 2{\mathbb {N}}^{*}$, and denoting $q=p/2\in {\mathbb {N}}^{*}$, then the corresponding differential operator simply reads:

  $$ \begin{array}{@{}rcl@{}} (-\triangle_{\text{PML}})^{p/2} & = & \text{Op}\Big(\Big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} +\texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} + \texttt{i}\frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} \Big)^{q}\Big) \end{array} $$

  and can easily be analytically computed and numerically approximated, as a standard differential operator.

• If $p \in 2{\mathbb {N}}+1$, with p = 2q + 1 and $q\in {\mathbb {N}}$, then we rewrite:

  $$ \begin{array}{@{}rcl@{}} \left. \begin{array}{lcl} \sigma\big((-\triangle_{\text{PML}})^{p/2}\big) & =&\displaystyle \sigma\big((-\triangle_{\text{PML}})^{q}\sqrt{-\triangle_{\text{PML}}} \big)\\ & = & \sigma\big((-\triangle_{\text{PML}})^{q}\big) \# \sigma\big(\sqrt{-\triangle_{\text{PML}}}\big)\\ & = & \displaystyle \sum\limits_{|\upbeta|=0}^{\infty}\frac{(-\texttt{i})^{|\upbeta|}}{\upbeta!}\partial^{{\upbeta}_{1}}_{x}\partial^{{\upbeta}_{2}}_{y} \sigma\big((-\triangle_{\text{PML}})^{q}\big)\partial^{{\upbeta}_{1}}_{\xi_{x}}\partial^{{\upbeta}_{2}}_{\xi_{y}}\sigma\big(\sqrt{-\triangle_{\text{PML}}}\big) , \end{array} \right. \end{array} $$

  where $\upbeta =({\upbeta }_{1},{\upbeta }_{2}) \in {\mathbb {N}}^{2}$ denotes a 2-index. Regarding $\sigma \big (\sqrt {-\triangle _{\text {PML}}}\big )$, we use Padé’s approximants, so that:

  $$ \begin{array}{@{}rcl@{}} \left. \begin{array}{lcl} \sigma\big((-\triangle_{\text{PML}})^{q+1/2}\big) \!& \approx &\! \displaystyle \sum\limits_{|\upbeta|=0}^{\infty}\frac{(-\texttt{i})^{|\upbeta|}}{\upbeta!}\partial^{{\upbeta}_{1}}_{x}\partial^{{\upbeta}_{2}}_{y} \Big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} -\texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} - \texttt{i}\frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} \Big)^{q} \\ & & \!\times \partial^{{\upbeta}_{1}}_{\xi_{x}}\partial^{{\upbeta}_{2}}_{\xi_{y}}\Big({\sum}_{k=0}^{M}a_{k}^{(M)} - {\sum}_{k=1}^{M}a_{k}^{(M)}d_{k}^{(M)}\big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} \\ & & \!+ \texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \texttt{i} \frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} + d_{k}^{(M)}\big)^{-1}\Big) . \end{array} \right. \end{array} $$

  From a practical point of view, we define:

  $$ \begin{array}{@{}rcl@{}} \sigma\big((-\triangle_{\text{PML}})_{m}^{q+1/2}\big) : = \sum\limits_{|\upbeta|=0}^{m}\lambda_{\upbeta}^{(m)}(x,y,\xi_{x},\xi_{y}) , \end{array} $$

  where

  $$ \begin{array}{@{}rcl@{}} \left. \begin{array}{lcl} \lambda_{\upbeta}^{(m)} & := & \partial^{{\upbeta}_{1}}_{x}\partial^{{\upbeta}_{2}}_{y} \Big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} -\texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} - \texttt{i}\frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} \Big)^{q} \\ & & \times \partial^{{\upbeta}_{1}}_{\xi_{x}}\partial^{{\upbeta}_{2}}_{\xi_{y}}\Big({\sum}_{k=0}^{M}a_{k}^{(M)} - {\sum}_{k=1}^{M}a_{k}^{(M)}d_{k}^{(M)}\big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} \\ & & + \texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \texttt{i} \frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} + d_{k}^{(M)}\big)^{-1}\Big) . \end{array} \right. \end{array} $$

  Tedious computations allow for an explicit expression of $\{\lambda _{\upbeta }^{(m)}\}_{\upbeta }$. We get:

  $$ \begin{array}{@{}rcl@{}} \left. \begin{array}{lcl} \lambda_{(0,0)}^{(m)} & = & \Big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} -\texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} - \texttt{i}\frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} \Big)^{q}\\ & & \displaystyle\times \Big(\sum\limits_{k=0}^{M}a_{k}^{(M)} - \sum\limits_{k=1}^{M}a_{k}^{(M)}d_{k}^{(M)}\big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} + \texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} \\&&+ \texttt{i} \frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} + d_{k}^{(M)}\big)^{-1}\Big) . \end{array} \right. \end{array} $$

  Next, we obtain:

  $$ \begin{array}{@{}rcl@{}} \left. \begin{array}{lcl} \lambda_{(1,0)}^{(m)} \!& = &\! -q\Big(\partial_{x}\Big(\frac{1}{{S_{x}^{2}}}\Big)|\xi_{x}|^{2} - \texttt{i}\partial_{x}\Big(\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\Big)\xi_{x}\Big)\Big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} - {\tt i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} - {\tt i}\frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} \Big)^{q-1}\\ & & \!\times \frac{{\sum}_{k=1}^{M}a_{k}^{(M)}d_{k}^{(M)}\big(\frac{1}{{S_{x}^{2}}}\partial_{\xi_{x}}|\xi_{x}|^{2} + \texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\big)}{{\sum}_{k=0}^{M}a_{k}^{(M)} - {\sum}_{k=1}^{M}a_{k}^{(M)}d_{k}^{(M)}\big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} + \texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \texttt{i} \frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} + d_{k}^{(M)}} , \end{array} \right. \end{array} $$

  and

  $$ \begin{array}{@{}rcl@{}} \left. \begin{array}{lcl} \lambda_{(0,1)}^{(m)} \!& = &\! -q\Big(\partial_{y}\Big(\frac{1}{{S_{y}^{2}}}\Big)|\xi_{y}|^{2} -\texttt{i}\partial_{y}\Big(\frac{S_{y}^{\prime}}{{S^{3}_{y}}}\Big)\xi_{x}\Big)\Big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} - {\tt i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} - {\tt i}\frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} \Big)^{q-1}\\ & & \!\times \frac{{\sum}_{k=1}^{M}a_{k}^{(M)}d_{k}^{(M)}\big(\frac{1}{{S_{x}^{2}}}\partial_{\xi_{y}}|\xi_{y}|^{2} + \texttt{i}\frac{S_{y}^{\prime}}{{S^{3}_{y}}}\big)}{{\sum}_{k=0}^{M}a_{k}^{(M)} - {\sum}_{k=1}^{M}a_{k}^{(M)}d_{k}^{(M)}\big(\frac{1}{{S_{x}^{2}}}|\xi_{x}|^{2} + \frac{1}{{S_{y}^{2}}}|\xi_{y}|^{2} + \texttt{i}\frac{S_{x}^{\prime}}{{S^{3}_{x}}}\xi_{x} + \texttt{i} \frac{S^{\prime}_{y}}{{S^{3}_{y}}}\xi_{y} \!+ d_{k}^{(M)}} . \end{array} \right. \end{array} $$

  For |β| = 1, we have to construct $\lambda ^{(m)}_{(1,1)}$, $\lambda ^{(m)}_{(2,0)}$, and $\lambda ^{(m)}_{(0,2)}$.

Appendix B: Cauchy integral approximation

In this Appendix, we discuss the approximation of A^α by using the Cauchy integral representation, for $\alpha \in {\mathbb {R}}$ and $A \in {\mathbb {R}}^{N\times N}$. We recall that:

$$ \begin{array}{@{}rcl@{}} A^{\alpha} & = & (2\pi \texttt{i})^{-1}A{\int}_{{\Gamma}_{A}}z^{\alpha-1}(zI-A)^{-1}dz , \end{array} $$

(63)

where Γ_A is a closed contour in the complex plane enclosing the spectrum of matrix A, where the latter is assumed to have its spectrum in ${\mathbb {C}} \backslash {\mathbb {R}}_{-}$. This approach can be quite inefficient if the spectrum of the matrix A has a large radius. This leads to a straightforward approximation of the Cauchy integral based on a quadrature rule:

$$ \begin{array}{@{}rcl@{}} A_{h}^{\alpha} & = & (2\pi \texttt{i})^{-1}A\sum\limits_{j}{\Delta} z_{j} \theta_{j}z_{j}^{\alpha-1}(z_{j}I-A)^{-1} , \end{array} $$

where {𝜃_j}_j are interpolation weights and $\{z_{j}\}_{j}\in {\Gamma }_{A} \subset {\mathbb {C}}$ are the interpolation nodes on Γ_A. There are many ways to reduce the computational complexity [6, 7, 25, 31]. Among others, we propose in [6] the following possible approach based on the use of a traditional preconditioner M for the linear system. Typically, M ≈ A^− 1, and MA has a spectrum clustering at the point (1,0) in the complex plane. Thus,

$$ \begin{array}{@{}rcl@{}} (MA)^{\alpha} & = & (2\pi \texttt{i})^{-1}MA{\int}_{{\Gamma}_{MA}}z^{\alpha-1}(zI-MA)^{-1}dz , \end{array} $$

(64)

where ℓ(Γ_M) ≪ ℓ(Γ), ℓ denoting the length of a curve in the complex plan. In particular computing (64) is cheaper than (63). However, the connection between (MA)^α and A^α is not necessarily simple.

Proposition 12

Assume that A is symmetric and M is a preconditioner commuting with A. Then, we have [6]:

$$ \begin{array}{@{}rcl@{}} A^{\alpha} = M^{-\alpha}(MA)^{\alpha}. \end{array} $$

In other words, the polynomial preconditioning allows for an efficient computation of matrix powers.

Practically, the proposed preconditioning allows for a reduction of the length of the contour enclosing the spectrum of the preconditioned matrix MA, as long as M^−α can be efficiently computed. We refer to [6] for additional details.