<EMAIL_ADDRESS>(E. Aldahham<EMAIL_ADDRESS>(L. Banjai)

# A modified convolution quadrature combined with the method of fundamental
solutions and Galerkin BEM for acoustic scattering

Ebraheem Aldahham  1  and Lehel Banjai  1  11affiliationmark:  Maxwell
Institute for Mathematical Sciences, School of Mathematical and Computer
Sciences, Heriot-Watt University, Edinburgh EH14 4AS, UK

###### Abstract

We describe a numerical method for the solution of acoustic exterior
scattering problems based on the time-domain boundary integral representation
of the solution. As the spatial discretization of the resulting time-domain
boundary integral equation we use either the method of fundamental solutions
(MFS) or the Galerkin boundary element method (BEM). In time we apply either a
standard convolution quadrature (CQ) based on an A-stable linear multistep
method or a modified CQ scheme. It is well-known that the standard low-order
CQ schemes for hyperbolic problems suffer from strong dissipation and
dispersion properties. The modified scheme is designed to avoid these
properties. We give a careful description of the modified scheme and its
implementation with differences due to different spatial discretizations
highlighted. Numerous numerical experiments illustrate the effectiveness of
the modified scheme and dramatic improvement with errors up to two orders of
magnitude smaller in comparison with the standard scheme.

###### keywords:

acoustic wave scattering, convolution quadrature, modified convolution
quadrature, method of fundamental solutions, boundary integral equation.

45E10, 65M80, 65L60, 65T50

##  1 Introduction

In this paper we investigate a class of numerical methods for the scattering
problem: Find $u(t)\in H^{1}(\Omega^{+})$ such that

$\displaystyle\partial_{t}^{2}u-\Delta u$ $\displaystyle=0$
$\displaystyle\text{for }(t,x)\in[0,T]\times\Omega^{+}$ (1.1) $\displaystyle
u(0)=\partial_{t}u(0)$ $\displaystyle=0$ $\displaystyle\text{for
}x\in\Omega^{+}$ $\displaystyle u$ $\displaystyle=g$ $\displaystyle\text{for
}(t,x)\in[0,T]\times\Gamma,$

where $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, is an open, bounded Lipschitz
domain with exterior $\Omega^{+}=\mathbb{R}^{d}\setminus\overline{\Omega}$ and
boundary $\Gamma=\partial\Omega$. If $g$ is the trace of an incident wave
$-u^{\text{inc}}$ on $\Gamma$ with
$\operatorname{supp}u^{\text{inc}}(0)\subset\Omega^{+}$ then $u$ is the
scattered wave and $u^{\text{tot}}=u+u^{\text{inc}}$ is the total wave
scattered by the obstacle under the sound-soft boundary condition
$u^{\text{tot}}|_{\Gamma}=0$.

There are a number of ways to tackle problem (1.1). A popular option is to
consider the problem on a bounded domain containing $\Gamma$ facilitated by an
introduction of a transparent boundary condition. The exact transparent
boundary condition is non-local in time and space and can be computed fast on
special domains such as balls [1, 24, 26, 30]. Alternatively, an approximate
local boundary condition can be used, such as local absorbing boundary
condition [20, 27], methods based on the pole condition [32, 33] and perfectly
matched layers [14]. All of these methods apply to convex domains of a special
shape, e.g., circle or rectangle in 2D. This can result in an unnecessarily
expensive method if $\Gamma$ is such that a very large convex domain is needed
to encompass it, e.g., the elongated shape of an airplane. In such situations
methods based on time-domain boundary integral equations are of advantage [17,
35, 10]. This particularly holds when the accurate computation of the far
field potential is needed as the boundary integral potentials have the exact
far-field behaviour encoded in their kernels.

Thus in this work we represent the solution as a single layer boundary
integral potential

$u(t,x)=S(\partial_{t})\varphi(t,x):=\int_{0}^{t}\int_{\Gamma}k(t-\tau,|x-y|)\varphi(\tau,y)d\Gamma_{y}d\tau,$
(1.2)

where $\varphi\colon[0,T]\times\Gamma\to\mathbb{R}$ is an unknown density and
$k$ is the fundamental solution which depends on the spatial dimension

$k(t,r)=\left\\{\begin{array}[]{cc}\displaystyle\frac{H(t-r)}{2\pi\sqrt{t^{2}-r^{2}}}&d=2\vspace{.1cm}\\\
\displaystyle\frac{\delta(t-r)}{4\pi r}&d=3\end{array}\right.,$ (1.3)

where $H(\cdot)$ is the Heaviside function and $\delta(\cdot)$ the Dirac delta
distribution.

Taking the trace onto $\Gamma$ of (1.2) we obtain the boundary integral
equation for the unknown density $\varphi$: Find $\varphi$ such that

$V(\partial_{t})\varphi(t,x):=\int_{0}^{t}\int_{\Gamma}k(t-\tau,|x-y|)\varphi(\tau,y)d\Gamma_{y}d\tau=g(t,x)\quad(t,x)\in[0,T]\times\Gamma.$
(1.4)

Once the density $\varphi$ is obtained, the solution $u$ can be recovered from
(1.2). All the numerical methods in this paper will be based on this boundary
formulation of the scattering problem.

Note that we have implicitly defined two operators above. The single layer
potential $S(\partial_{t})\cdot$ and single layer boundary integral operator
$V(\partial_{t})\cdot$· The motivation for the notation and the mapping
properties are described in [29, 10] and will be briefly explained in the next
section.

A popular method for the discretization of (1.4) is the space-time Galerkin
method originating in the work of Bamberger and Ha Duong [2, 3].
Alternatively, any spatial discretization can be combined with convolution
quadrature (CQ) discretization in time as introduced by Lubich in [29].
Combined with Galerkin discretization in space, the analysis of CQ
discretization for standard problems of acoustics is nearly complete; see
books [35, 10]. Recently, a combination with the method of fundamental
solutions (MFS) was investigated in [28]. While CQ discretization has a number
of advantages – powerful toolbox of analysis techniques and fast methods,
relatively simple implementation, a vast literature on various applications,
broadness of applicability – a recognised weakness of the approach is its
dissipative and dispersive properties [15, 10]. While high-order methods based
on Runge-Kutta can tackle this deficiency very effectively [6, 8], it is still
of interest to develop accurate (and cheaper) low-order variants of CQ.
Despite the effectiveness of high-order methods based on Runge-Kutta time-
stepping, low order methods are of interest due to their ease of
implementation and ease of coupling with other solvers in bounded,
inhomogeneous regions; see works on FEM-BEM coupling in the time-domain [7,
10]. This has led to works on modified CQ schemes, see e.g., the works of
Davies and Duncan [18, 19], Weile [36], and [12] by the second author. Here we
make use of the method introduced by Weile. Our description will be closer to
the description given in [4] by the second author. The novelty of the current
work is the description of the FFT based implementation of the scheme, the
combination with MFS and numerous numerical experiments. The numerical results
will show up to two orders of magnitude improvement in error when compared
with the standard, low-order, schemes. We should say that the modification we
describe in this paper is applicable to high-er order schemes and are likely
to improve accuracy even for high order Runge-Kutta schemes. The investigation
of higher order schemes is however beyond the scope of the current work. We
next give a brief content of the paper.

For the spatial discretization of (1.2) and (1.4) we will consider two
approaches: the Galerkin boundary element method (BEM) and the method of
fundamental solutions (MFS). For the time discretization we will also consider
two approaches: standard convolution quadrature (CQ) based on A-stable linear
multistep methods and a class of modified CQ methods specially designed for
wave propagation problems. We begin by describing standard CQ together with
the modification mentioned above. Next, we describe spatial discretization
approaches. Finally, we explain how the discretizations combine together and
describe efficient algorithms for the implementation. The paper is finished by
a section with extensive numerical experiments.

##  2 One-sided convolutions and CQ

Convolution quadrature is a numerical method designed for the discretization
of convolution

$K(\partial_{t})g:=\int_{0}^{t}k(t-\tau)g(\tau)d\tau,$

where $k$ is a given kernel function with $K$ its Laplace transform

$K(s)=\\{\mathscr{L}k\\}(s):=\int_{0}^{\infty}e^{-st}k(t)dt\qquad\operatorname{Re}s>0.$

If $k$ is not an integrable function but potentially a distribution with
Laplace transform $K$ analytic for $\operatorname{Re}s>0$ and bounded as

$|K(s)|\leq C(\sigma)|s|^{\mu}\qquad\forall\operatorname{Re}s\geq\sigma>0,$

for some $\mu\in\mathbb{R}$ and a positive function $C$.

we define the convolution $K(\partial_{t})g$ via the inverse Laplace transform

$K(\partial_{t})g(t):=\mathscr{L}^{-1}\\{KG\\}(t)=\frac{1}{2\pi\mathrm{i}}\int_{\sigma-\mathrm{i}\infty}^{\sigma+\mathrm{i}\infty}e^{st}K(s)G(s)ds,$

where $G=\\{\mathscr{L}g\\}$. As long as $G$ decays faster than
$-\mu-1-\varepsilon$ for a fixed $\varepsilon>0$, $K(\partial_{t})g(t)$ is a
continuous function that can be continuously extended to $t<0$ by zero.
Alternatively, if $g\in C^{m}(\mathbb{R}_{\geq 0})$ with
$g(0)=g^{\prime}(0)=\dots=g^{(m-1)}(0)=0$ we have

$\left|\mathscr{L}\\{g\\}(s)\right|\leq|s|^{-m}\int_{0}^{\infty}e^{-\sigma
t}|g^{(m)}(t)|dt\qquad\operatorname{Re}s\geq\sigma>0,$

and hence $K(\partial_{t})g$ is continuous and $K(\partial_{t})g(0)=0$ if
$m>\mu+1$.

Given a time-step ${\Delta t}>0$ and a generating function of an A-stable
linear multistep method, for instance

$\delta(\zeta)=(1-\zeta)+\frac{1}{2}(1-\zeta)^{2}\;(\text{BDF2})\qquad\delta(\zeta)=2\frac{1-\zeta}{1+\zeta}\;(\text{trapezoidal
rule}),$

convolution quadrature (CQ) of the above convolution is defined by

$K(\partial_{t}^{\Delta
t})g(t):=\frac{1}{2\pi\mathrm{i}}\int_{\sigma-\mathrm{i}\infty}^{\sigma+\mathrm{i}\infty}e^{st}K(s^{\Delta
t})G(s)ds,$ (2.1)

where

$s^{\Delta t}=\frac{\delta(e^{-s{\Delta t}})}{{\Delta t}}.$

The two linear multistep methods are second order and A-stable, which in terms
of the generating functions is stated as

$\delta(e^{-z})=z+\mathcal{O}(z^{3})$ (2.2)

and

$\operatorname{Re}\delta(e^{-z})>0\qquad\text{for }\operatorname{Re}z>0.$
(2.3)

For sufficiently smooth $g$ it is shown in [29, 21, 10] that

$\left|K(\partial_{t}^{\Delta
t})g(t)-K(\partial_{t})g(t)\right|=\mathcal{O}({\Delta t}^{2})$

for the above linear multistep methods.

Let us note that $K(\partial_{t}^{\Delta t})g(t)$ is a discrete convolution,
namely with $t_{n}=n{\Delta t}$

$K(\partial_{t}^{\Delta t})g(t_{n})=\sum_{j=0}^{n}\omega_{n-j}g(t_{j})$

where $\omega_{j}$ are the convolution weights defined by the generating
function

$K\left(\frac{\delta(\zeta)}{{\Delta
t}}\right)=\sum_{j=0}^{\infty}\omega_{j}\zeta^{j}.$

A crucial property for the stability of CQ is the composition rule which
implies that solving the discrete equation

$K(\partial_{t}^{\Delta t})\varphi(t_{n})=g(t_{n})\qquad n=0,\dots,N$

is equivalent to evaluating the discrete convolution

$\varphi_{n}=K^{-1}(\partial_{t}^{\Delta t})g(t_{n})\qquad n=0,\dots,N.$

##  3 A modified CQ

A fundamental property of solutions of the wave equation is the fact that they
travel at finite speed. This is nicely visible in the definition of the
fundamental solutions in (1.3) which are zero for $t<r$, i.e., are zero before
the wave has arrived from the source to the receiver. This property is also
visible in the Laplace domain since a simple calculation shows that if $u$ is
polynomially bounded and locally integrable then for $r\geq 0$

$u(t)\equiv 0\text{ for
}t<r\;\implies\;\left|e^{sr}\mathscr{L}\\{u\\}(s)\right|\leq\int_{0}^{\infty}e^{-\sigma
t}|u(t+r)|dt\text{ for }\operatorname{Re}s\geq\sigma>0.$ (3.1)

By the inverse Laplace transform and the Cauchy integral formula the reverse
direction holds as well: For $r\geq 0$

$|K(s)e^{sr}|\leq C|s|^{\mu}\text{ for }\operatorname{Re}s>0\text{ and
}\mu<-1\implies k(t)\equiv 0\text{ for }t\leq r.$ (3.2)

In particular if we know that the kernel $K(s)$ is such that

$\left|e^{st_{m}}K(s)\right|\leq
C(\sigma)|s|^{\tilde{\mu}}\qquad\operatorname{Re}s\geq\sigma>0,$

for some $\tilde{\mu}\in\mathbb{R}$ and a positive function $C$, then
$K(\partial_{t})g(t)\equiv 0$ for $t\leq t_{m}$ for sufficiently smooth $g$.
This motivates a modified CQ, where the quadrature is only applied after
$t>t_{m}$:

$K(\partial_{t}^{\Delta
t};t_{m})g(t):=\left\\{\begin{array}[]{cc}\displaystyle\frac{1}{2\pi\mathrm{i}}\int_{\sigma-\mathrm{i}\infty}^{\sigma+\mathrm{i}\infty}e^{s(t-t_{m})}e^{s^{\Delta
t}t}K(s^{\Delta t})G(s)ds,&t>t_{m}\\\ 0&t\leq t_{m}.\end{array}\right.$

This is again a discrete convolution

$K(\partial_{t}^{\Delta
t};t_{m})g(t_{n})=\sum_{j=0}^{n}\omega_{n-j;m}g(t_{j})$ (3.3)

where the generating function for the weights is given by

$\zeta^{m}e^{t_{m}\frac{\delta(\zeta)}{{\Delta
t}}}K\left(\frac{\delta(\zeta)}{{\Delta
t}}\right)=\zeta^{m}e^{m\delta(\zeta)}K\left(\frac{\delta(\zeta)}{{\Delta
t}}\right)=\sum_{j=0}^{\infty}\omega_{j;m}\zeta^{j}.$ (3.4)

Note that $\omega_{j;m}=0$ for $j<m$ and hence $K(\partial_{t}^{\Delta
t};t_{m})g(t)=0$ for $t<t_{m}$. This modified CQ avoids approximating values
that we know are zero and furthermore, as we will see later, it has the
potential to give much more accurate solutions in certain cases.

##  4 Spatial discretization of integral operators

Time-domain boundary integral operators fit seamlessly into the setting of
Section 2. Namely, the Laplace domain single layer potential is given by

$S(s)\varphi(x):=\int_{\Gamma}K(s,|x-y|)\varphi(y)d\Gamma_{y}\qquad
x\in\mathbb{R}^{d}\setminus\Gamma,$

where $K(s,r)$ is the Laplace transform of the free space Green’s function for
wave equation in (1.3)

$K(s,r)=\left\\{\begin{array}[]{cc}\displaystyle\frac{1}{2\pi}K_{0}(sr)&d=2\vspace{.1cm}\\\
\displaystyle\frac{e^{-sr}}{4\pi r}&d=3\end{array}\right.,$

with $K_{0}(\cdot)$ the modified Bessel function of the second kind. From the
definition of $K$ for $d=3$ and from the asymptotic behaviour of $K_{0}$ [31]
we have that for $r>0$ there exists a constant $C(r)>0$ such that

$|e^{sr}K(s,r)|\leq C(r)\qquad\text{for }\operatorname{Re}s>0.$ (4.1)

Recall that this corresponds to the time-domain property $k(t,r)=0$ for $t<r$;
see (1.3).

The single layer potential is continuous across the boundary, hence in the
frequency domain we are interested in the single-layer boundary integral
operator $V(s)\colon H^{-1/2}(\Gamma)\to H^{1/2}(\Gamma)$

$V(s)\varphi(x):=\int_{\Gamma}K(s,|x-y|)\varphi(y)d\Gamma_{y}\qquad
x\in\Gamma.$

By now, [2, 35, 10], it is well-known that $V(s)$ and $V^{-1}(s)$ are bounded
for $\operatorname{Re}s>0$ as

$\begin{split}\|V(s)\|_{H^{-1/2}(\Gamma)\to H^{1/2}(\Gamma)}&\leq
C\max(1,|s|^{-2})\frac{|s|}{\operatorname{Re}s}\\\
\|V^{-1}(s)\|_{H^{1/2}(\Gamma)\to H^{-1/2}(\Gamma)}&\leq
C\max(1,|s|^{-1})\frac{|s|^{2}}{\operatorname{Re}s}.\end{split}$

We now describe two standard approaches to discretizing these operators.

###  4.1 Galerkin boundary element methods (BEM)

Let $\Gamma$ be subdivided into $M$ non-overlapping panels $\Gamma_{i}$,
$i=1,\dots,M$. We denote by $X_{h}\subset H^{-1/2}(\Gamma)$ the finite-
dimensional space of piecewise constant functions. A basis of the space is
given by $\varphi_{1},\dots,\varphi_{M}$ defined as

$\varphi_{i}(x)=\left\\{\begin{array}[]{cc}1&x\in\Gamma_{i}\\\
0&\text{otherwise}.\end{array}\right.$

Then we define the matrix representation of the Galerkin discretization of the
single-layer operator by

$\begin{split}\left(\mathbf{V}_{\text{G}}(s)\right)_{ij}&:=\int_{\Gamma}\int_{\Gamma}K(s,|x-y|)\varphi_{i}(y)\varphi_{j}(x)d\Gamma_{y}d\Gamma_{x}\\\
&=\int_{\Gamma_{i}}\int_{\Gamma_{j}}K(s,|x-y|)d\Gamma_{y}d\Gamma_{x}.\end{split}$
(4.2)

An important property of this matrix, implied by (4.1), is that

$\left|e^{sr_{ij}}\left(\mathbf{V}_{\text{G}}(s)\right)_{ij}\right|\leq const$

for $\operatorname{Re}s\geq\sigma\geq 0$, $r_{ij}>0$, and

$r_{ij}\leq\operatorname{dist}(\Gamma_{i},\Gamma_{j}).$

Again, we can interpret this in the time-domain as the fact that information
travelling from one panel, e.g., $\Gamma_{i}$, needs time greater than
$r_{ij}$ to reach the second panel, $\Gamma_{j}$.

With this spatial discretization in place, the semi-discretization of (1.4)
reads: Find $\boldsymbol{\varphi}(t)\in\mathbb{R}^{M}$ such that

$\mathbf{V}_{\text{G}}(\partial_{t})\boldsymbol{\varphi}(t)=\mathbf{g}(t),$
(4.3)

where the right-hand side is the projection $\mathbf{g}(t)\in\mathbb{R}^{M}$
of the boundary data defined by

$\left(\mathbf{g}(t)\right)_{i}=\int_{\Gamma}g(t,x)\varphi_{i}(x)d\Gamma_{x}.$

###  4.2 Method of Fundamental Solutions

As we have seen in the previous sections, the Laplace domain solution
$U(s)=\mathscr{L}\\{u\\}(s)$ can be represented by the single layer potential:

$U(s,x)=\int_{\Gamma}K(s,|x-y|)\Phi(s,y)d\Gamma_{y},$

where we have now explicitly given the dependence of the density $\Phi$ on the
Laplace domain parameter $s$. The Galerkin discretization described above,
while accurate, stable and well-understood is expensive and non-trivial to
implement due to the singular, double integrals. A simple alternative is the
Method of Fundamental Solutions (MFS) [13, 22]. In the MFS, we approximate
$U(s)$ by a sum of sources located at $y_{j}$, $j=1,\dots,K$

$U(s,x)\approx\sum_{j=1}^{K}\Phi_{j}(s)K(s,|x-y_{j}|).$ (4.4)

The unknown coefficients are determined by solving the least squares problem
derived from the boundary condition: Find $\Phi_{j}(s)$ such that

$\sum_{j=1}^{K}\Phi_{j}(s)K(s,|x_{i}-y_{j}|)=G(s,x_{i}),\quad i=1,\dots,M,$
(4.5)

where $x_{i}\in\Gamma$ are the collocation points and $G=\mathscr{L}\\{g\\}$.
Note that the source points $y_{j}$ are located in the exterior of the
computational domain and certainly away from the boundary, hence no
singularities occur in the system. A careful numerical analysis of this method
for connected planar domains with analytic boundary has been given in [13].

Again, we can define a system matrix
$\mathbf{V}_{\text{MFS}}(s)\in\mathbb{C}^{M\times K}$

$\left(\mathbf{V}_{\text{MFS}}(s)\right)_{ij}=K(s,|x_{i}-y_{j}|),$ (4.6)

which in contrast with the Galerkin method is rectangular and trivial to
compute. Importantly again we have the property

$\left|e^{sr_{ij}}\left(\mathbf{V}_{\text{MFS}}(s)\right)_{ij}\right|\leq
const,$

where now $r_{ij}=|x_{i}-y_{j}|$.

The semi-discrete system again reads: Find
$\boldsymbol{\varphi}(t)\in\mathbb{R}^{K}$ such that

$\mathbf{V}_{\text{MFS}}(\partial_{t})\boldsymbol{\varphi}(t)=\mathbf{g}(t),$
(4.7)

where the right-hand side is now simply

$\left(\mathbf{g}(t)\right)_{i}=g(t,x_{i}),\qquad i=1,\dots,M.$

##  5 Accuracy of the two time-discretizations

The approximation due to the standard CQ, according to (2.1) consists of
replacing in the Laplace domain the kernel $K(s,r)$ by its approximation
$K(s^{{\Delta t}},r)$ with $s^{{\Delta t}}=\delta(e^{-s{\Delta t}})/{\Delta
t}$. In particular, in 3D this means that we are replacing $e^{-sr}$ by
$e^{-s^{{\Delta t}}r}$. The error of this approximation can be bounded as

$\left|e^{-sr}-e^{-s^{{\Delta t}}r}\right|\leq Cr|s-s^{{\Delta t}}|\leq
C|rs||s{\Delta t}|^{2},$

for some constant $C>0$ that is allowed to change from one step to another and
$|s{\Delta t}|<c$ for some small enough constant $c>0$. Since $s=i\omega,$ if
we denote by $\omega_{\text{max}}$ the largest frequency present in the
system, the estimate above implies that to get a fixed accuracy we expect that
we need $r{\Delta t}^{2}\omega_{\text{max}}^{3}$ to be small enough, i.e., the
time-step ${\Delta t}$ needs to be chosen proportional to
$\omega_{\text{max}}^{-3/2}$ resulting in the number of degrees of freedom per
wavelength not being fixed.

Instead, in the modified method the approximation is of the form
$e^{-st_{m}}e^{-s^{{\Delta t}}(r-t_{m})}$, where $t_{m}$ is the largest time-
step smaller than $r_{ij}\leq r$. The approximation is now

$\left|e^{-sr}-e^{-st_{m}}e^{-s^{{\Delta
t}}(r-t_{m})}\right|\leq\left|e^{-s(r-t_{m})}-e^{-s^{{\Delta
t}}(r-t_{m})}\right|\leq C|(r-t_{m})s||s{\Delta t}|^{2}.$

In both the case of the Galerkin and MFS discretizations, $t_{m}$ can be
chosen so that $r-t_{m}=\mathcal{O}(h)+\mathcal{O}({\Delta t})$, where $h$ is
the spatial mesh-width. If, as is usually the case, we have $h\propto{\Delta
t}$, we obtain that the above error remains at a fixed accuracy if ${\Delta
t}^{3}\omega_{\text{max}}^{3}$ is small enough, i.e., a fixed number of
degrees of freedom is sufficient.

Similar arguments hold for the 2D case, recalling the asymptotic behaviour of
$K_{0}$ for large arguments [31, 10.25.3]:

$K_{0}(z)\sim\sqrt{\frac{\pi}{2z}}e^{-z}\qquad\text{as }z\to\infty,\;|\arg
z|<\pi.$

##  6 Fully discrete system and algorithmic implementation

In the following, $\mathbf{V}(s)$ will denote either the Galerkin system
matrix $\mathbf{V}_{\text{G}}$ or the corresponding matrix
$\mathbf{V}_{\text{MFS}}$ obtained via MFS. Hence, the semi-discrete system
reads: Find $\boldsymbol{\varphi}(t)$ such that

$\mathbf{V}(\partial_{t})\boldsymbol{\varphi}(t)=\mathbf{g}(t).$

The fully discrete system can directly be obtained by standard convolution
quadrature discretization in time: Find $\boldsymbol{\varphi}^{{\Delta
t}}_{j}$, $j=0,\dots,N$, such that

$\mathbf{V}(\partial_{t}^{\Delta t})\boldsymbol{\varphi}^{{\Delta
t}}(t_{n})=\mathbf{g}(t_{n}),\qquad n=0,\dots,N,$ (6.1)

where $t_{j}=j{\Delta t}$, ${\Delta t}>0$ is the time-step and $t_{N}=T$ the
final time.

To describe the modified scheme, recall that we are given distances $r_{ij}$
such that

$\left|e^{sr_{ij}}\left(\mathbf{V}(s)\right)_{ij}\right|\leq const.$

Next we define

$m_{ij}=\left\lfloor\frac{r_{ij}}{{\Delta t}}\right\rfloor$

so that $0\leq t_{m_{ij}}\leq r_{ij}$ and is the largest time-step satisfying
this inequality.

The modified scheme then reads: Find
$\widetilde{\boldsymbol{\varphi}}^{{\Delta t}}_{\ell}$, $\ell=0,\dots,N$, such
that

$\mathbf{V}(\widetilde{\partial}^{\Delta
t}_{t})\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}(t_{n})=\mathbf{g}(t_{n}),\qquad n=0,\dots,N,$ (6.2)

where

$\mathbf{V}(\widetilde{\partial}^{\Delta
t}_{t})\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}(t_{n}):=\sum_{\ell=0}^{n}\widetilde{\omega}_{n-\ell}\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}_{\ell}$ (6.3)

and the modified weights are given by the following generating function, c.f.
(3.4),

$\displaystyle\widetilde{\mathbf{V}}(\zeta)$
$\displaystyle=\sum_{n=0}^{\infty}\widetilde{\omega}_{n}\zeta^{n}.$ (6.4a)
$\displaystyle\left(\widetilde{\mathbf{V}}(\zeta)\right)_{ij}$
$\displaystyle:=\zeta^{m_{ij}}e^{m_{ij}\delta(\zeta)}\left(\mathbf{V}\left(\frac{\delta(\zeta)}{{\Delta
t}}\right)\right)_{ij}.$ (6.4b)

If $\tilde{\omega}_{0}$ is invertible in the least squares sense, then the
solution is given by the iteration: Find
$\widetilde{\boldsymbol{\varphi}}^{{\Delta t}}_{n}$ such that

$\widetilde{\omega}_{0}\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}_{n}=\mathbf{g}(t_{n})-\sum_{\ell=0}^{n-1}\widetilde{\omega}_{n-\ell}\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}_{\ell}$ (6.5)

in the least squares sense.

Note that the weights in (6.3) are defined by the generating function (6.4a)
as in the standard CQ. Hence, all the FFT based algorithms designed for CQ are
still available. For completeness, we describe the most basic of these
algorithms. More details and other algorithms can be found in [10, 11].

First of all, notice that (6.4a) is a Taylor expansion of the function
$\widetilde{\mathbf{V}}(\zeta)$ analytic for $|\zeta|<1$. Hence, we can
represent the coefficients by the Cauchy integral formula

$\begin{split}\widetilde{\omega}_{j}&=\frac{1}{2\pi\mathrm{i}}\oint_{\mathcal{O}}\widetilde{\mathbf{V}}\left(\zeta\right)\zeta^{-j-1}d\zeta\\\
&=\lambda^{-j}\int_{0}^{1}\widetilde{\mathbf{V}}\left(\lambda
e^{2\pi\mathrm{i}\theta}\right)e^{-2\pi\mathrm{i}j\theta}d\theta,\end{split}$

where we have chosen the disk of radius $0<\lambda<1$ as the contour. Applying
the composite trapezoidal rule to this (periodic) integral we obtain an
approximation of the weights

$\widetilde{\omega}_{j}\approx\frac{\lambda^{-j}}{N+1}\sum_{\ell=0}^{N}\widetilde{\mathbf{V}}\left(\lambda\zeta_{N+1}^{-\ell}\right)\zeta_{N+1}^{\ell
j},$ (6.6)

where $\zeta_{N+1}=e^{\frac{2\pi\mathrm{i}}{N+1}}$. This approximation is
valid for $j=0,\dots,N$ and the error is of $\mathcal{O}(\lambda^{N+1})$. In
finite precision arithmetic, due to the multiplication with the factor
$\lambda^{-j}$, $\lambda$ cannot be chosen too small. In practice
$\lambda=\varepsilon^{1/2(N+1)}$, where $\varepsilon$ is close to machine
precision, is a good choice giving as error $\sqrt{\varepsilon}$; for details
see [29] and for how to improve on this accuracy see [10].

Importantly, in (6.6) we recognise the inverse discrete Fourier transform,
which can be computed in $\mathcal{O}(N\log N)$ time for all $j=0,\dots,N$ by
the fast Fourier transform (FFT). In a similar way, further FFT based
algorithms as described in [10] are available for the modified scheme. Note,
however, that data sparse techniques such as FMM [16, 23] and
$\mathcal{H}$-matrices [25] cannot be directly applied to the computation of
the spatial operator $\mathbf{V}(s)$ due to the discontinuous changes in the
spatial kernel resulting from (6.4b). An alternative would be to use the exact
shifts for groups of panels and their interaction, rather than for each pair
in the Galerkin matrix. In this way, the kernel would still be smooth allowing
for data sparse techniques and the new method would still significantly
improve the accuracy. Detailed research would be needed to determine the
optimal way of choosing these groups of panels. This is similar to the
approach taken in [5], where the near field was treated in the time-domain and
the far field was treated in the Fourier domain with data sparse techniques.
Another possibility would be to smooth out the shifts similar to what was done
in [34] in the case of the space-time Galerkin method.

In the numerical experiments of this paper, we will use a parallel FFT
algorithm introduced in [9], whose details we describe next. First of all we
note that the approximation (6.6) of the weights $\widetilde{\omega}_{j}$ is
valid also for $j=-N,\dots,-1$ if we define $\widetilde{\omega}_{k}=0$ for
$k<0$. Furthermore, extending the sum in (6.3) to $N$ we have that the system
to be solved is

$\sum_{\ell=0}^{N}\widetilde{\omega}_{n-\ell}\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}_{\ell}=\mathbf{g}_{n}$

for $n=0,\dots,N$. Substituting now the approximation (6.6), we have the
following (approximate) system to be solved

$\sum_{\ell=0}^{N}\left[\frac{\lambda^{-(n-\ell)}}{N+1}\sum_{k=0}^{N}\widetilde{\mathbf{V}}\left(\lambda\zeta_{N+1}^{-k}\right)\zeta_{N+1}^{k(n-\ell)}\right]\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}_{\ell}\approx\mathbf{g}_{n}.$

Rearranging gives

$\frac{\lambda^{-n}}{N+1}\sum_{k=0}^{N}\widetilde{\mathbf{V}}\left(\lambda\zeta_{N+1}^{-k}\right)\left[\sum_{\ell=0}^{N}\lambda^{\ell}\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}_{\ell}\zeta_{N+1}^{-\ell k}\right]\zeta_{N+1}^{kn}\approx\mathbf{g}_{n}.$

Multiplying both sides by $\lambda^{n}$ and applying the discrete Fourier
transform we obtain $N+1$ decoupled linear systems to be solved:

$\widetilde{\mathbf{V}}\left(\lambda\zeta_{N+1}^{-k}\right)\left[\sum_{\ell=0}^{N}\lambda^{\ell}\widetilde{\boldsymbol{\varphi}}^{{\Delta
t}}_{\ell}\zeta_{N+1}^{-\ell
k}\right]\approx\sum_{n=0}^{N}\lambda^{n}\mathbf{g}_{n}\zeta_{N+1}^{-kn},$
(6.7)

$k=0,\dots,N$. Thus, once the (scaled) discrete Fourier transform of the data
$\mathbf{g}_{n}$ is computed, $N+1$ linear systems need to be solved in
parallel. The solution $\varphi_{\ell}^{{\Delta t}}$ is then recovered by
another application of the inverse discrete Fourier transform. Note that using
symmetry, roughly half of the problems in (6.7) are conjugates of the other
half and hence need not to be solved.

######  Remark 6.1.

Note that if ${\Delta t}<\min_{i,j}r_{i,j}$ then $m_{ij}>0$ for all $i,j$ and
from the definition of the weights $\tilde{\omega}_{n}$ in (6.4) it follows
that

$\left(\tilde{\omega}_{n}\right)_{ij}=\left(\widetilde{\mathbf{V}}(0)\right)_{ij}=0$

for all $i,j$. In this case we are not able to solve (6.2) by using the
iteration (6.5). However, as we will see in the numerical experiments, the
formula (6.7) is feasible even in this case and gives good results.

##  7 Numerical experiments

###  7.1 Transient wave scattering and MFS

We first consider exterior scattering by either the unit disk
$\mathcal{D}=\\{x\ \in\mathbb{R}^{2}\;:\;|x|<1\\}$ or the union of two open
ellipses $\mathcal{E}=\mathcal{E}_{1}\cup\mathcal{E}_{2}$, where the ellipses
are defined using the identification of $\mathbb{R}^{2}$ with the complex
plane via

$\mathcal{E}_{1}:=\left\\{f(z)-2\;:\;|z|<1\right\\}\qquad\mathcal{E}_{2}:=\left\\{f(z)+2\;:\;|z|<1\right\\},$

with $f$ the conformal map

$f(z)=\frac{1}{2}e^{\mathrm{i}\frac{\pi}{2}}(z+\frac{1}{5}z^{-1}).$

As the Dirichlet data in both cases we use

$g(t,x)=\sin(\omega(t-x\cdot\alpha))g(t-4-x\cdot\alpha)$

with

$g(t)=e^{-(t/0.7)^{2}},\qquad\alpha=(0,-1)^{T}$

and $\omega$ a parameter that we will choose as either $\omega=1$ or
$\omega=5$.

For the Galerkin boundary element method, we let $X_{h}$ be the space of
piecewise constant boundary element functions. Whereas for the method of
fundamental solutions in the case of the unit disk we let the collocation
points be

$x^{\mathcal{D}}_{i}=e^{2\pi\mathrm{i}i/M},\qquad
i=0,\dots,M-1\qquad(\mathcal{D})$

and the source points

$y_{j}^{\mathcal{D}}=Re^{2\pi\mathrm{i}j/K},\qquad
j=0,\dots,K-1\qquad(\mathcal{D})$

for some $R<1$. For all exterior scattering problems we use $R=0.9$. At the
end of the section, we will also show results for the interior problem for the
unit disk where we choose $R=1.1$.

For each of the two ellipses we make use of the conformal map $f$ and define
the collocation points to be

$x^{\mathcal{E}_{1}}_{i}=f(x^{\mathcal{D}}_{i})-2,\qquad
i=0,\dots,M-1\qquad(\mathcal{E}_{1})$

and source points

$y^{\mathcal{E}_{1}}_{j}=f(y^{\mathcal{D}}_{j})-2,\qquad
j=0,\dots,M-1\qquad(\mathcal{E}_{1})$

and analogously for $\mathcal{E}_{2}$. For a plot of the source and
collocation points in each of the two cases see Figure 1.

Figure 1: We show the position of the collocation points (solid lines),
position of the source points (dashed lines) and the points $X_{j}$ where we
evaluate the error. On the left the setting for the exterior scattering by the
unit disk is shown, whereas on the right by two ellipses.

This explains the spatial discretization of the problems. In time, we either
use the standard CQ based either on BDF2 or the trapezoidal rule or its
modified version. The fully discrete system is solved using the formula (6.7),
where each linear system is solved as a least squares problem.

Finally, to measure the error we choose some test points $X_{\ell}$,
$\ell=1,\dots,L$ placed in the exterior domain. As we do not have an exact
solution at hand, we will use a finer mesh to obtain an accurate approximation
denoted by $u_{\text{ex}}$. The error measure is the maximum error over all
test points and all time-steps

$\text{error}=\max_{\ell,n}|u^{h}(t_{n},X_{\ell})-u_{\text{ex}}(t_{n},X_{\ell})|,$
(7.1)

where $u^{h}$ denotes the numerical solution with a coarser mesh in time and
space.

Figure 2: Convergence of the MFS with BDF2 based methods on the left and MFS
with trapezoidal rule based methods on the right for the case of exterior
scattering by the unit disk.

In Figure 2, in the case of the exterior scattering by the unit disk, we plot
the convergence of the error for the standard CQ based on BDF2 and trapezoidal
rule and its modified counterparts. In these calculations we used $M=2000$
collocation points and $K=1000$ source points for the MFS to compute $u^{h}$.
To compute the accurate solution $u^{\text{ex}}$ we used $M^{\text{ex}}=3000$
collocation points and $K^{\text{ex}}=1500$ source points and in time the
modified CQ based on BDF2 with $N^{\text{ex}}=2^{12}$ time-steps. We see a
significant improvement when using the modified method.

Figure 3: Convergence of the MFS with BDF2 based methods on the left and MFS
with trapezoidal rule based methods on the right for the case of exterior
scattering by two ellipses.

For the exterior scattering by two ellipses we plot the convergence of the
error for the standard CQ based on BDF2 and trapezoidal rule and its modified
counterparts in Figure 3. In these calculations we used $M=4000$ collocation
points and $K=2000$ source points for the MFS to compute $u^{h}$. To compute
the accurate solution $u^{\text{ex}}$ we used $M^{\text{ex}}=6000$ collocation
points and $K^{\text{ex}}=3000$ source points and in time the modified CQ
based on BDF2 with $N=2^{12}$ time-steps. We see a significant improvement
when using the modified method. Note, however, that for the modified scheme
there is a transient region where the error is very large; see Figure 3. We
believe, though we have no proof, that the reason for this is that the new
scheme does not damp high frequencies whereas CQ does. This makes CQ stable
even when it is inaccurate, whereas the modified scheme produces controlled
results only once it starts converging.

(a) At t=2.5

(b) At t=3.75

(c) At t=5

(d) At t=6.25

(e) At t=7.5

(f) At t=8.75

Figure 4: Snapshots of the total field for the scattering by two ellipses.

Again, we see that a much larger time-step is sufficient to obtain good
accuracy when using the modified scheme compared with the standard CQ.
However, here we also see a curious behaviour for very large time-steps, where
the standard method, while inaccurate, also remains reasonably bounded,
whereas the modified scheme produces very large errors. We have seen this
effect often in numerical experiments whenever the time-step is too large.
This is not a real deficiency of the method as with such a large time-step the
results would anyway be inaccurate but is a property of the modified scheme
that one should be aware of. In Figure 4 we show snapshots of the total field
for the scattering by two ellipses with the frequency $\omega=5$.

###  7.2 Transient wave scattering and the Galerkin method in space

In this section we perform experiments with the Galerkin BEM discretization in
space as described in Section 4.1, i.e., we use the space of piecewise
constant boundary elements. The time discretization is as in the previous
section and the implementation is based on the formula (6.7). We let the
domain $\Omega$ be the non-convex domain seen in Figure 5. Its boundary is
composed of four semi-circles defined in the complex plane by

$\displaystyle\Gamma_{1}$
$\displaystyle:=\left\\{e^{\mathrm{i}\theta}\;:\;-\frac{\pi}{2}\leq\theta<\frac{\pi}{2}\right\\}$
(7.2) $\displaystyle\Gamma_{2}$
$\displaystyle:=\left\\{\frac{1}{4}e^{\mathrm{i}\theta}+\frac{3}{4}\mathrm{i}\;:\;\frac{\pi}{2}\leq\theta<\frac{3\pi}{2}\right\\}$
$\displaystyle\Gamma_{3}$
$\displaystyle:=\left\\{\frac{1}{2}e^{-\mathrm{i}\theta}\;:\;-\frac{\pi}{2}\leq\theta<\frac{\pi}{2}\right\\}$
$\displaystyle\Gamma_{4}$
$\displaystyle:=\left\\{\frac{1}{4}e^{\mathrm{i}\theta}-\frac{3}{4}\mathrm{i}\;:\;\frac{\pi}{2}\leq\theta<\frac{3\pi}{2}\right\\}.$

Figure 5: We show the non-convex domain whose boundary is described in (7.2)
and the points $X_{j}$ where the error will be computed.

As the incident wave we use

$u^{\text{inc}}(t,x)=-\sin(\omega(t-x\cdot\alpha))g(t-4-x\cdot\alpha)$ (7.3)

with

$g(t)=e^{-(t/0.7)^{2}},\qquad\alpha=\frac{1}{\sqrt{2}}(1,1)^{T}$

and $\omega$ a parameter that we will choose as either $\omega=1$ or
$\omega=5$. The Dirichlet data is then given by $-u^{\text{inc}}$ and the
final time is set to $T=10$. Snapshots of the total field with $\omega=5$ are
shown in Figure 8.

Figure 6: Convergence of the error for Galerkin BEM discretization in space
and CQ and modified CQ discretization in time. The scatterer is the non-convex
domain in Figure 5 and the incident wave is the plane-wave (7.3). An
instability is seen in the modified scheme based on the trapezoidal scheme.
This as instability is removed when finer discretization is used in space; see
Figure 7. Figure 7: Convergence of the error for Galerkin BEM discretization
in space and CQ and modified CQ discretization in time. The scatterer is the
non-convex domain in Figure 5 and the incident wave is the plane-wave (7.3)
with $M=2000$ and $M_{e}=3000$. The finer discretization in space removes the
instability seen in Figure 6.

The exact solution is not known hence we use a more accurate approximation
instead and denote this by $u_{\text{ex}}$. As the error measure we use the
maximum error in space at points $X_{\ell}$, see Figure 7.2, and all time-
steps; see (7.1). For the computation of $u^{\text{ex}}$ a uniform boundary
element mesh is used with $M^{\text{ex}}=1500$ degrees of freedom and
$N^{\text{ex}}=2^{12}$. For the computation of the approximate solution
$u^{h}$ we use again a uniform boundary element mesh with $M=1000$ degrees of
freedom, whereas the number of time-steps is increased from $N=2^{3}$ to
$N=2^{10}$. Again we perform experiments with both BDF2 and trapezoidal based
schemes. The convergence results are shown in Figure 6 and again show
considerably better performance of the modified scheme except that for the
smallest time-step, the modified scheme based on the trapezoidal scheme
becomes unstable. This is not an unknown problem even for standard convolution
quadrature. Its stability is only assured if the quadrature is sufficiently
accurate. This problem is easily remedied either by improving the quadrature
rule or by using a finer mesh in space; see Figure 7 where we increased the
spatial degrees of freedom to $M=2000$ and $M^{\text{ex}}=3000$. It is by now
well-understood that the CQ based on the trapezoidal is more difficult to
discretize in space. The reason behind this is that the frequencies

$s_{\ell}=\frac{\delta(\lambda\zeta_{N+1}^{-\ell})}{{\Delta t}}$

are of size $\mathcal{O}({\Delta t}^{-1})$ for the BDF2 scheme and of size
$\mathcal{O}({\Delta t}^{-2})$ for the trapezoidal scheme due to the
singularity of $\delta(\zeta)$ at $\zeta=-1$ in the case of the latter; for
details see [21, 10]. As the numerical results for the modified BDF2 scheme
are comparable to those of the modified trapezoidal scheme, in practice one
may prefer to restrict using the BDF2 scheme.

(a) At t=2.5

(b) At t=3.75

(c) At t=5

(d) At t=6.25

(e) At t=7.5

(f) At t=8.75

Figure 8: Snapshots of the total field for the scattering by the non-convex
domain.

###  7.3 Interior scattering (MFS)

To conclude the numerical experiments, we consider an interior scattering
problem. The domain $\Omega$ is the unit disk and the incident wave is given
by its initial data

$u^{\text{inc}}(0)=e^{-\frac{1}{2}a^{2}|x-x_{0}|^{2}}\qquad\partial_{t}u^{\text{inc}}(0)=0,$

with

$a=10\qquad x_{0}=(\tfrac{1}{4},0).$

Note that the choice of parameters is such that $|u^{\text{inc}}(x,0)|\lesssim
10^{-13}$ for $x\in\Gamma=\partial\Omega$, i.e., to almost machine precision
the data is initially zero at the boundary. Due to the radial form of the
initial data, the incident wave can be computed for any $x$, and $t$ using the
Hankel transform. Namely,

$u^{\text{inc}}(x,t)=\int_{0}^{\infty}F(k)J_{0}(kr)k\cos(kt)dt\qquad
r=|x-x_{0}|,$

where $J_{0}$ is a Bessel function of the 1st kind and

$F(k)=\frac{1}{a^{2}}e^{-\frac{k^{2}}{2a^{2}}}$

is the Hankel function of the function $e^{-\frac{1}{2}a^{2}r^{2}}$. Thus, the
Dirichlet data is given by

$g(x,t)=-u^{\text{inc}}(x,t)\qquad x\in\Gamma=\partial\Omega,\;t\geq 0.$

We discretize in space using MFS, where use $M=2000$ collocation points and
$K=1000$ source points. The setting is shown in Figure 9. To compute the
accurate solution $u^{\text{ex}}$ we used $M^{\text{ex}}=3000$ collocation
points and $K^{\text{ex}}=1500$ source points and in time the modified CQ
based on BDF2 with $N^{\text{ex}}=2^{11}$ time-steps. The convergence of the
modified and standard methods is shown in Figure 10. Again, we see a
significant improvement when using the modified method. Similar results were
obtained using the Galerkin method. To save on space we do not show the
results here.

Figure 9: We show the setting for the interior scattering by the unit disk,
where the position of the collocation points (solid lines), position of the
source points (dashed lines) and the points $X_{j}$ where we evaluate the
error. Figure 10: Convergence of the MFS with BDF2 based methods on the left
and MFS with trapezoidal rule based methods on the right for the case of
interior scattering by the unit disk.

## References

  * [1] B. Alpert, L. Greengard, and T. Hagstrom.  Nonreflecting boundary conditions for the time-dependent wave equation.  J. Comput. Phys., 180(1):270–296, 2002.
  * [2] A. Bamberger and T. H. Duong.  Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d’une onde acoustique. I.  Math. Methods Appl. Sci., 8(3):405–435, 1986.
  * [3] A. Bamberger and T. H. Duong.  Formulation variationnelle pour le calcul de la diffraction d’une onde acoustique par une surface rigide.  Math. Methods Appl. Sci., 8(4):598–608, 1986.
  * [4] L. Banjai.  Dissipation free low order convolution quadrature for TDBIE.  In 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA), pages 1210–1213, 2015.
  * [5] L. Banjai and M. Kachanovska.  Fast convolution quadrature for the wave equation in three dimensions.  J. Comput. Phys., 279:103–126, 2014.
  * [6] L. Banjai, C. Lubich, and J. M. Melenk.  Runge-Kutta convolution quadrature for operators arising in wave propagation.  Numer. Math., 119(1):1–20, 2011.
  * [7] L. Banjai, C. Lubich, and F.-J. Sayas.  Stable numerical coupling of exterior and interior problems for the wave equation.  Numer. Math., 129(4):611–646, 2015.
  * [8] L. Banjai, M. Messner, and M. Schanz.  Runge-Kutta convolution quadrature for the boundary element method.  Comput. Methods Appl. Mech. Engrg., 245/246:90–101, 2012.
  * [9] L. Banjai and S. Sauter.  Rapid solution of the wave equation in unbounded domains.  SIAM J. Numer. Anal., 47(1):227–249, 2008/09.
  * [10] L. Banjai and F.-J. Sayas.  Integral equation methods for evolutionary PDE.  Springer Series in Computational Mathematics. Springer, To be published in 2022.
  * [11] L. Banjai and M. Schanz.  Wave propagation problems treated with convolution quadrature and BEM.  In Fast boundary element methods in engineering and industrial applications, volume 63 of Lect. Notes Appl. Comput. Mech., pages 145–184. Springer, Heidelberg, 2012.
  * [12] L. Banjai and Y. Zhang.  A family of efficient numerical solvers of time domain boundary integral equations.  In Proceedings of Forum Acusticum, Aalborg, 2011.
  * [13] A. H. Barnett and T. Betcke.  Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains.  J. Comput. Phys., 227(14):7003–7026, 2008.
  * [14] J.-P. Berenger.  A perfectly matched layer for the absorption of electromagnetic waves.  J. Comput. Phys., 114(2):185–200, 1994.
  * [15] Q. Chen, P. Monk, X. Wang, and D. Weile.  Analysis of convolution quadrature applied to the time-domain electric field integral equation.  Commun. Comput. Phys., 11(2):383–399, 2012.
  * [16] R. Coifman, V. Rokhlin, and S. Wandzura.  The fast multipole method for the wave equation: a pedestrian prescription.  IEEE Antennas and Propagation Magazine, 35(3):7–12, 1993.
  * [17] M. Costabel and F.-J. Sayas.  Time-dependent problems with boundary integral equation method.  In E. e. a. Stein, editor, Encyclopedia of computational mechanics Second Edition. Part 2, pages 1–24. John Wiley & Sons, Ltd., 2017\.
  * [18] P. J. Davies and D. B. Duncan.  Convolution-in-time approximations of time domain boundary integral equations.  SIAM J. Sci. Comput., 35(1):B43–B61, 2013.
  * [19] P. J. Davies and D. B. Duncan.  Convolution spline approximations for time domain boundary integral equations.  J. Integral Equations Appl., 26(3):369–410, 2014.
  * [20] B. Engquist and A. Majda.  Absorbing boundary conditions for the numerical simulation of waves.  Math. Comp., 31(139):629–651, 1977.
  * [21] H. Eruslu and F.-J. Sayas.  Polynomially bounded error estimates for trapezoidal rule convolution quadrature.  Comput. Math. Appl., 79(6):1634–1643, 2020.
  * [22] G. Fairweather and A. Karageorghis.  The method of fundamental solutions for elliptic boundary value problems.  Adv. Comput. Math., 9(1-2):69–95, 1998.
  * [23] L. Greengard.  The rapid evaluation of potential fields in particle systems.  ACM Distinguished Dissertations. MIT Press, Cambridge, MA, 1988.
  * [24] M. J. Grote and J. B. Keller.  Nonreflecting boundary conditions for time-dependent scattering.  J. Comput. Phys., 127(1):52–65, 1996.
  * [25] W. Hackbusch.  Hierarchical matrices: algorithms and analysis, volume 49 of Springer Series in Computational Mathematics.  Springer, Heidelberg, 2015.
  * [26] T. Hagstrom.  Radiation boundary conditions for the numerical simulation of waves.  In Acta numerica, 1999, volume 8 of Acta Numer., pages 47–106. Cambridge Univ. Press, Cambridge, 1999.
  * [27] T. Hagstrom, A. Mar-Or, and D. Givoli.  High-order local absorbing conditions for the wave equation: extensions and improvements.  J. Comput. Phys., 227(6):3322–3357, 2008.
  * [28] I. Labarca and R. Hiptmair.  Acoustic scattering problems with convolution quadrature and the method of fundamental solutions.  Commun. Comput. Phys., 30(4):985–1008, 2021.
  * [29] C. Lubich.  On the multistep time discretization of linear initial-boundary value problems and their boundary integral equations.  Numer. Math., 67(3):365–389, 1994.
  * [30] C. Lubich and A. Schädle.  Fast convolution for nonreflecting boundary conditions.  SIAM J. Sci. Comput., 24(1):161–182, 2002.
  * [31] NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, Release 1.1.1 of 2021-03-15.  F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * [32] D. Ruprecht, A. Schädle, and F. Schmidt.  Transparent boundary conditions based on the pole condition for time-dependent, two-dimensional problems.  Numer. Methods Partial Differential Equations, 29(4):1367–1390, 2013.
  * [33] D. Ruprecht, A. Schädle, F. Schmidt, and L. Zschiedrich.  Transparent boundary conditions for time-dependent problems.  SIAM J. Sci. Comput., 30(5):2358–2385, 2008.
  * [34] S. Sauter and A. Veit.  A Galerkin method for retarded boundary integral equations with smooth and compactly supported temporal basis functions.  Numer. Math., 123(1):145–176, 2013.
  * [35] F.-J. Sayas.  Retarded potentials and time domain boundary integral equations, volume 50 of Springer Series in Computational Mathematics.  Springer, [Cham], 2016.  A road map.
  * [36] D. S. Weile.  A hybrid Runge-Kutta convolution quadrature-temporal Galerkin approach to the solution of the time domain integral equations of electromagnetics.  In 2014 International Conference on Electromagnetics in Advanced Applications (ICEAA), pages 391–394, 2014.