# Nyström methods for high-order CQ solutions of the wave equation in two
dimensions

Peter Petropoulos, Catalin Turc, Erli Wind-andersen

###### Abstract

We investigate high-order Convolution Quadratures methods for the solution of
the wave equation in unbounded domains in two dimensions that rely on Nyström
discretizations for the solution of the ensemble of associated Laplace domain
modified Helmholtz problems. We consider two classes of CQ discretizations,
one based on linear multistep methods and the other based on Runge-Kutta
methods, in conjunction with Nyström discretizations based on Alpert and QBX
quadratures of Boundary Integral Equation (BIE) formulations of the Laplace
domain Helmholtz problems with complex wavenumbers. We present a variety of
accuracy tests that showcase the high-order in time convergence (up to and
including fifth order) that the Nyström CQ discretizations are capable of
delivering for a variety of two dimensional scatterers and types of boundary
conditions.
Keywords: wave equation, convolution quadratures, Nyström discretizations.

AMS subject classifications: 65N38, 35J05, 65T40,65F08

##  1 Introduction

Numerical solutions of the wave equation based on retarded potentials boundary
integral equations [22, 36, 7, 8, 4] enjoy certain advantages over volumetric
solvers such as FDTD and FEM, especially when the wave equation is posed in
unbounded domains. In particular, Convolution Quadratures (CQ) have emerged in
the last two decades as a powerful machinery for the numerical solution of
time domain boundary integral equations. Indeed, the CQ provide a means to
extend boundary integral equation (BIE) based frequency domain solvers to time
domain BIE solvers via Laplace transforms. We present in this paper a CQ based
methodology to extend high-order Nyström BIE frequency domain solvers to high-
order solvers of the time domain equation in the exterior of two dimensional
obstacles that exhibit geometric singularities (e.g. corners, arc tips).

Whenever applicable, BIE solvers are attractive choices for the solution of
wave equations in the frequency domain (e.g. Helmholtz equation), owing to the
dimensional reduction they provide and the explicit enforcement of radiation
conditions when dealing with unbounded regions. In addition, high-order
Nyström BIE solvers for Helmholtz equations in geometrically complex two
dimensional geometries (including non-smooth geometries) have been proposed in
the literature [19, 16, 17, 20, 25, 12, 15, 23, 26] and can be more efficient
than BEM solvers. Most of these BIE solvers were designed to handle the case
of Helmholtz equations with _real_ wavenumbers, which is the most relevant in
physical applications. CQ methods for the solution of the two dimensional wave
equation, on the other hand, require the solution of an ensemble of _modified_
Helmholtz boundary value problems, that is Helmholtz equations with _complex_
wavenumbers, and the size of the ensemble grows with the order in time that CQ
achieve. We focus in this paper on the second order in time BDF2 CQ [5] and
the third and fifth order in time Runge Kutta CQ [6], whose algorithmic
details are presented in those references. BIE solvers can be incorporated
into the CQ machinery to solve wave equations provided their underpinning
quadratures can handle integral kernels that feature complex wavenumbers. As
it turns out, the wavenumbers associated with the CQ modified Helmholtz
equations may have very large (positive) imaginary parts which become larger
as the order of the CQ time integrator increases.

While some of the high-order Helmholtz BIE quadratures proposed in the
literature are oblivious to the nature of the wavenumbers, others struggle to
resolve the additional challenges that complex wavenumbers with large
imaginary parts pose. For instance, and as illustrated in this paper, the
popular Kussmaul Martensen singularity splitting high-order quadratures [29,
32], while very effective in the case of real wavenumbers, cannot be extended
in a straightforward manner to handle the complex wavenumbers countered in the
CQ context. A method that does not suffer from this issue is the delta-bem
method presented in [19, 16, 17, 20] which produces third order quadratures
schemes for the discretization of the four boundary integral operators needed
for BIE solutions of Helmholtz boundary value problems in two dimensions,
irrespective of whether the wavenumbers are real or complex. Owing to this
desirable feature, the delta-bem quadratures have been incorporated into high-
order CQ solvers of the wave equation [16, 17, 20], and the code has been
freely available for almost a decade [33] serving as an important
computational platform in the CQ community.

We focus in this paper on two different types of high-order Nyström
quadratures that can also handle seamlessly discretizations of Helmholtz BIOs
with complex wavenumbers: (1) one that relies on Alpert quadratures [1] for
weakly singular kernels, and (2) the other which is the increasingly popular
Quadratures by Expansion (QBX) introduced in [26]. Unlike the delta-bem
methods, both aforementioned quadratures (i) do not rely on global
interpolation and therefore are amenable to panel discretizations and (ii) can
be extended to three dimension configurations (albeit the Alpert quadrature
can only handle axisymmetric scatterers [23]). Most of the delta-bem and
Alpert quadrature literature is concerned with smooth boundaries, and
therefore we focus mostly in this paper on two dimensional Lipschitz
boundaries. High-order Nyström discretizations of the Helmholtz BIOs
associated with Lipschitz curves have been proposed [12, 25, 18], and they
resolve the singularities of the BIO functional densities through the use of
graded meshes that allow for the definition of weighted densities that are
more regular. We follow this program and we present in this paper a simple
extension of Alpert quadratures that produces high-order Nyström
discretizations of Helmholtz BIOs corresponding to Lipschitz boundaries via
sigmoid graded meshes [15]. In addition, we present a new version of QBX
methods based on Chebyshev meshes and Féjer quadratures that leads to high-
order discretizations of Helmholtz BIOs in non smooth (including open)
boundaries. An attractive feature of QBX methods based on Chebyshev meshes is
their natural treatment of the corner and/or open end singularities of the
Dirichlet and Neumann traces of Helmholtz solutions on non-smooth boundaries.

The high-order discretizations of Helmholtz BIOs with complex wavenumbers lead
to high-order BIE solvers for Helmholtz problems, which can then be
incorporated in the (blackbox) CQ machinery to deliver high-order in time
solutions of the wave equation. While the choice of the BIE formulations can
affect the accuracy of CQ schemes [10], we focus on the simplest possible
formulations, that is the single layer and respectively the double layer
formulations in the case of Dirichlet and respectively Neumann boundary
conditions. Besides the fact that both of these formulations are well posed
for the CQ Laplace domain problems, we prefer them because they can be equally
applied to both closed and open boundaries. In the case of Neumann boundary
conditions the formulation we chose entails the Helmholtz hypersingular
operator whose discretization is challenging for Nyström discretizations. We
investigate an alternative formulation of the hypersingular operator which is
he Current and Charge (CC) formulation [35, 24]. The CC formulation uses two
functional densities instead of one in order to bypass the need to effect
numerical differentiation, and in addition it involves Cauchy Principal Value
BIOs. Nevertheless, we show that the QBX Nyström discretization can still
deliver high-order discretizations of the Current and Charge (CC) formulations
on closed and open Lipschitz boundaries.

The paper is organized as follows: in Section 2 we present the details of the
BDF2 and RK CQ applied to the wave equation in two dimensions; in Section 3 we
present Nyström discretizations of the ensemble of Laplace domain modified
Helmholtz problems related to BDF2 and RK CQ methods for the wave equation;
finally in Section 4 we present numerical results that illustrate the high
order in time (up and including fifth order) achieved by CQ solutions of the
wave equation in the exterior of geometrically complex two dimensional
obstacles.

##  2 Convolution Quadratures for the solution of the Wave Equation

We are interested in solving the wave equation in a two-dimensional infinite
domain $\Omega^{+}:=\mathbb{R}^{2}\setminus\overline{\Omega}$, where $\Omega$
is a bounded domain in $\mathbb{R}^{2}$:

$\begin{cases}c^{-2}\frac{\partial^{2}u}{\partial t^{2}}&=\Delta u\quad{\rm
in}\ \Omega^{+}\times(0,\infty)\\\ u(x,0)&=\frac{\partial u(x,0)}{\partial
t}=0\quad{\rm in}\ \Omega^{+}\\\ u(x,t)&=g(x,t)\quad{\rm on}\
\Gamma\times(0,\infty),\end{cases}$ (1)

and we denoted $\Gamma:=\partial\Omega$. Using the fundamental solution of the
wave equation in two dimensions

$k_{2D}(x,t):=\frac{H(t-c^{-1}|x|)}{2\pi\sqrt{t^{2}-c^{-2}|x|^{2}}}$

we seek for a solution $u(x,t)$ of the wave equation in the form of a Time
Domain Single-Layer (TDSL) potential:

$\displaystyle u(x,t)$ $\displaystyle=$
$\displaystyle(S*\lambda)(x,t)=\int_{\Gamma}\int_{0}^{t}k_{2D}(t-\tau,x-y)\varphi(y,\tau)ds(y)d\tau$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{\Gamma}\int_{0}^{t-c^{-1}|x-y|}\frac{\varphi(y,\tau)}{\sqrt{(t-\tau)^{2}-c^{-2}|x-y|^{2}}}ds(y)d\tau$

in terms of the unknown density function $\varphi$ defined on
$\Gamma\times(0,\infty)$. Imposing the Dirichlet boundary conditions leads to
the following Time Domain Boundary Integral Equation (TDBIE)

$\int_{\Gamma}\int_{0}^{t}k_{2D}(t-\tau,x-y)\varphi(y,\tau)ds(y)d\tau=g(x,t)\quad{\rm
on}\ \Gamma\times(0,\infty).$ (2)

We assume in what follows that the wave speed is $c=1$ and we solve the TDBIE
(2) using the Convolution Quadrature (CQ) methodology which we describe in
what follows. Our presentation is concerned exclusively with the algorithmic
aspects of CQ and is based on [7, 10, 33].

###  2.1 Convolution Quadratures for TDBIE

CQ are hybrid time stepping methods that combine Laplace transforms and
A-stable ODE solvers to produce approximations of convolution integrals

$(f*g)(t)\coloneqq\int_{0}^{t}f(\tau)g(t-\tau)d\tau,$ (3)

for values of $t$ on an equidistant mesh $0=t_{0}\leq\dots\leq t_{N}=N\Delta
t:=T$.  For a given function $f$ we define its Laplace transform

$F(s)\coloneqq\int_{0}^{\infty}e^{-st}f(t)dt,$

which we assume to exist for $s\in\mathbb{C},\ \Re{s}>0$ and recast the
convolution $f*g$ in the form

$(f*g)(t)=\frac{1}{2\pi
i}\int_{\sigma-\infty}^{\sigma+\infty}F(s)\underbrace{\left\\{\int_{0}^{t}e^{s\tau}g(t-\tau)d\tau\right\\}}_{:=y(t;s)}ds.$

In a nutshell, multistep CQ methods apply A-stable multistep solvers (e.g. the
second order in time multistep method BDF2 which is considered in this
section) to the ODE $\frac{dy}{dt}=sy+g,\ y(0)=0$ satisfied by the expression
$y(t;s)$ defined in the curly brackets above and ultimately deliver a discrete
convolution

$(f*g)(t_{n})\approx\sum_{m=0}^{n}\omega_{m}(\Delta t)g_{n-m}\quad\text{ where
}g_{n}=g(t_{n})$

that approximates the continuous convolution $(f*g)(t_{n}),0\leq n\leq N$. The
convolution weights $\omega_{m}(\Delta t)$ in the formula above are defined as

$F\left(\frac{P(\zeta)}{\Delta t}\right)=\sum_{m=0}^{\infty}\omega_{m}(\Delta
t)\zeta^{m}\quad\text{ where }\omega_{m}(\Delta
t)=\frac{1}{m!}\frac{d^{m}}{d\zeta^{m}}\left(F\left(\frac{P(\zeta)}{\Delta
t}\right)\right)\Big{|}_{\zeta=0},$ (4)

where where $P(\zeta)=\frac{1}{2}(\zeta^{2}-4\zeta+3)$ is the characteristic
polynomial of BDF2. The convolution in time present in the Kirchhoff integral
on the left hand side of equation (2) is recast through Laplace transforms in
the form

$\int_{0}^{t}\int_{\Gamma}k_{2D}(x-y,t-\tau)\varphi(y,\tau)ds(y)d\tau=\frac{1}{2\pi
i}\int_{\sigma-i\infty}^{\sigma+i\infty}\int_{0}^{t}e^{\tau
s}\left(V(s)\varphi(\cdot,t-\tau)\right)(x)d\tau ds,$

where $V(s)$ denotes the single layer BIO associated with a modified Helmholtz
equation

$\left(V(s)\varphi\right)(x)=\int_{\Gamma}K(x-y,s)\varphi(y)ds(y),\quad
x\in\Gamma.$ (5)

In equation (5), the kernel $K(x,s)$ denotes the Laplace transform of the wave
equation Green function $k_{2D}(x,t)$ with respect to the variable $t$, that
is

$K(x,s)=\frac{i}{4}H_{0}^{(1)}(is|x|).$

We note that the single layer potential

$w(x;s):=\int_{\Gamma}K(x-y,s)\varphi(y)ds(y),\quad
x\in\mathbb{R}^{2}\setminus\Gamma$

is a radiative solution of the modified Helmholtz equation

$\Delta_{x}w(x;s)-s^{2}w(x;s)=0\quad x\in\mathbb{R}^{2}\setminus\Omega.$ (6)

The CQ machinery produces discrete approximations of the continuous
convolutions on the left hand side of the TDBIE

$\frac{1}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}\int_{0}^{t}e^{\tau
s}\left(V(s)\varphi(\cdot,t-\tau)\right)(x)d\tau ds=g(x,t)\quad t>0\quad
x\in\Gamma.$ (7)

Specifically, denoting by

$\varphi^{\Delta t}_{j}(x)=\varphi^{\Delta t}(x,t_{j}),\quad 0\leq j\leq N$

the discrete CQ approximations lead us to the following linear system of BIEs

$\sum_{j=0}^{n}\left(\omega_{n-j}^{\Delta
t}\varphi_{j}\right)(x)=g_{n}(x),\quad n=1,\cdots,N\quad x\in\Gamma,$ (8)

where according to equations (4) the convolution weights are defined as Taylor
series operator coefficients

$V\left(\frac{P(\zeta)}{\Delta t}\right)=\sum_{n=0}^{\infty}\omega_{n}^{\Delta
t}(V)\zeta^{n},\quad|\zeta|\leq 1.$ (9)

The operator convolution weights $\omega_{j}^{\Delta t}(V)$ are computed from
the Cauchy formula using the trapezoidal rule

$\omega_{j}^{\Delta t}(V)=\frac{1}{2\pi
i}\int_{|\zeta|=\lambda}\frac{V(\frac{P(\zeta)}{\Delta
t})}{\zeta^{j+1}}d\zeta\approx\frac{\lambda^{-j}}{N+1}\sum_{\ell=0}^{N}V(s_{\ell})\zeta_{N+1}^{\ell
j},\quad s_{\ell}=\frac{P(\lambda\zeta^{-\ell}_{N+1})}{\Delta t}.$ (10)

The lower triangular linear system of frequency domain BIEs (8) can be fully
decoupled using the Discrete Fourier Transform (DTF). Indeed, straightforward
DFT manipulations lead us to the following $N+1$ decoupled BIEs

$\left(V(s_{\ell})\hat{\varphi_{\ell}}\right)(x)=\hat{g_{\ell}}(x)\quad
x\in\Gamma\quad\ell=0,\ldots,N,$ (11)

where

$\displaystyle\hat{g_{\ell}}(x)$
$\displaystyle:=\sum_{n=0}^{N}\lambda^{n}g_{n}(x)\zeta^{-\ell n}_{N+1}$
$\displaystyle\hat{\varphi_{l}}(x)$
$\displaystyle:=\sum_{n=0}^{N}\lambda^{n}\varphi_{n}(x)\zeta^{-ln}_{N+1}.$

Once the densities $\hat{\varphi_{\ell}}(x)$ are computed by solving the
sequence of frequency domain BIE for all wavenumbers $s_{\ell}$, we retrieve
the quantities $\varphi_{\ell}(x)$ which are approximations of
$\varphi(x,\ell\Delta t)$ using the inverse DFT:

${\varphi_{\ell}}(x)=\sum_{j=0}^{N}\lambda^{n}\hat{\varphi}_{n}(x)\zeta^{\ell
j}_{N+1},\quad\ell=0,\ldots,N.$ (12)

Due to the Dalquist barrier, the maximum achievable order of convergence of
linear multistep CQ methods is second order in time. Higher order in time CQ
methods can be obtained using Runge-Kutta methods which we will review next.

###  2.2 Runge-Kutta Convolution Quadrature Methods

Runge-Kutta based CQ methods can achieve orders of convergence in time higher
than three when applied to the wave equation [5, 6]. Our exposition of these
methods follows closely the one in [10]. Applying an $m$ stage Runge-Kutta
(RK) scheme with the associated Butcher Tableau consisting of the ensemble of
matrices and vectors

$A=[a_{ij}]_{1\leq i,j\leq m}\quad b=[b_{j}]_{1\leq j\leq m}\quad
c=[c_{j}]_{1\leq j\leq m},$

to the first order system reformulation of the wave equation, one arrives via
the $\zeta$-transforms at the following ensemble of Laplace domain modified
Helmholtz equations

$\begin{cases}\Delta_{x}W_{j}(x;\zeta)-\left(\frac{\gamma_{j}(\zeta)}{\Delta
t}\right)^{2}W_{j}(x;\zeta)=0\quad&x\in\Omega^{+}\\\
W_{j}(x;\zeta)=\sum_{\ell=1}^{m}\mathds{P}^{-1}_{j\ell}(\zeta)G_{\ell}(x;\zeta)\quad
x\in\Gamma\end{cases}$ (13)

where $G_{\ell}(x;\zeta)$ is the following $\zeta$-transform

$G_{\ell}(x;\zeta):=\sum_{n\geq 0}g(x;t_{n}+c_{\ell}\Delta t)\zeta^{n}\quad
x\in\Gamma.$ (14)

Once the solutions $W_{j}(\cdot;\zeta)$ of the modified Helmholtz problems
above are obtained, the $\zeta$-transform of the solution of the wave equation

$U_{d}(x;\zeta):=\sum_{n\geq 0}u_{d}(x;t_{n})\zeta^{n}$

is in turn immediately retrieved from the formula

$U_{d}(x;\zeta)=\zeta\sum_{j=1}^{m}\mathds{P}_{mj}(\zeta)W_{j}(x;\zeta).$ (15)

In equations (13) the wavenumbers $\gamma_{j}(\zeta)$ are related to the
diagonalization of the matrix valued function $\Delta(\zeta)$ defined by

$\Delta(\zeta)=\left(A+\frac{\zeta}{1-\zeta}\mathds{1}b^{T}\right)^{-1}$

where $\mathds{1}=[1,\ldots,1]^{T}$. Specifically, $\mathds{P}(\zeta)$ is the
matrix consisting of the eigenvectors of $\Delta(\zeta)$ and
$\mathds{D}(\zeta)=diag(\gamma_{1}(\zeta),\ldots,\gamma_{m}(\zeta))$ is the
diagonal matrix with the corresponding eigenvalues of $\Delta(\zeta)$ so that
$\Delta(\zeta)=\mathds{P}(\zeta)\mathds{D}(\zeta)\mathds{P}(\zeta)^{-1}$.
Finally, reverting from the Laplace domain to the physical domain is performed
in the same manner as in the case of linear multistep methods.

We will use in our numerical solutions two-stage and three-stage Implicit RK
Radau IIa CQ methods. The two-stage RK Radau IIa method gives rise to third
order in time CQ and will be denoted by the acronym RK3. The three-stage RK
Radau IIa method gives rise to fifth order in time CQ and will be denoted by
the acronym RK5. We note that the higher order in time convergence that can be
achieved by RK CQ methods entail commensurately more solutions of frequency
domain modified Helmholtz equations.

##  3 Nyström discretizations for the solution of modified Helmholz BIE

As we have already remarked, a key component of CQ methods for the solution of
TDBIE is solving ensembles of modified Helmholtz equations (13) in unbounded
domains. We will use Nyström discretizations for the solution of the BIE
formulations of those modified Helmholtz problems that involve complex
wavenumbers. Whereas in the case of BIE associated with Helmholtz equations
with real wavenumbers several Nyström discretization strategies exist in the
literature (see for example [23] for an excellent review of some of those
strategies), the case of Helmholtz equations with complex wavenumbers has
received less attention. We discuss in what follows the challenges that are
encountered in extending Nyström discretizations to the case of BIOs with
complex wavenumbers and we present several methods that can deal seamlessly
with ensemble evaluation of Helmholtz BIOs for a wide range of complex
wavenumbers, and for Lipschitz and smooth boundaries (both closed and open).

###  3.1 Dirichlet Boundary Conditions

We note that the BDF2 CQ solution of TDBIE (2) presented in Section 2.1
requires in turn the solution of Helmholtz problems of the type (6) via
associated single layer operators (9). Furthermore, each of the Helmholtz
problems (13) related to the RK CQ methods presented in Section 2.2 can be
again solved via single layer potentials and boundary integral operators.
Therefore, we discuss in what follows the single layer potential approach for
the solution of Helmholtz problems in unbounded domains with Dirichlet
boundary conditions. Specifically, we seek to solve the following Helmholtz
problem with wavenumber $k\in\mathbb{C}$

$\begin{cases}\Delta u+k^{2}u=0\quad&x\in\Omega^{+}\\\
u(x)=g(x)\quad&x\in\Gamma\\\
u(x)=\sum_{m=-\infty}^{\infty}u_{m}H_{m}^{(1)}(k|x|)e^{im\theta},\quad&|x|>r_{0},\quad
x=|x|(\cos{\theta},\sin{\theta}).\end{cases}$ (16)

We note that the last condition at infinity imposed in equations (16) is
equivalent to the Sommerfeld radiation condition
$\lim_{|x|\to\infty}\sqrt{|}x|(\frac{\partial}{\partial|x|}-ik)u(|x|,\theta)=0$
in the case when $k>0$.

For wavenumbers $k\in\mathbb{C}$ such that $\Im{k}>0$—which turn out to be all
of the wavenumbers involved in the Helmholtz problems (6) and respectively
(13), we seek the solution $u$ of the Helmholtz problem (16) in the form of
the single layer potential

$u(x)=SL(k)[\varphi](x):=\int_{\Gamma}G_{k}(x-y)\varphi(y)ds(y),x\in\mathbb{R}^{2}\setminus\Gamma,\quad
G_{k}(z)=\frac{i}{4}H_{0}^{(1)}(k|z|),\ z\neq 0,$

and thus the boundary density $\varphi$ solves the BIE of the first kind

$V(k)[\varphi](x):=\int_{\Gamma}G_{k}(x-y)\varphi(y)ds(y)=g(x),\quad
x\in\Gamma.$ (17)

We note that the BIE (17) is uniquely solvable for all wavenumbers $k$ such
that $\Im{k}>0$, and all types of boundaries $\Gamma$, be the boundaries
closed (Lipschitz) or open. Given its versatility, the single layer BIE
formulation (17) seems to be the preferred method of solution in the CQ
literature [7, 5, 33]. Furthermore, coercivity estimates on the operators
$V(k):H^{-1/2}(\Gamma)\to H^{1/2}(\Gamma)$ are available when $\Im{k}>0$ [22]
and they greatly facilitate the error analysis [7, 33] of CQ solution of TDBIE
(2). Nevertheless, in the case when $\Gamma$ is closed, a variety of other BIE
formulations are available for the solution of the Helmholtz problems (16),
such as the double layer or the combined field approach, and their influence
on the accuracy of the CQ schemes in connection with the location of the
complex poles of the BIE formulation was discussed in [10].

The main challenge faced by Nyström discretizations of the single layer
formulations of the type (17) is the resolution of the logarithmic singularity
of the integral kernel. To this end, several singular quadratures have been
proposed in the literature. For instance, the Kussmaul-Martensen quadrature
[29, 32] assumes a $2\pi$ periodic parametrization
$x:[0,2\pi]\to\mathbb{R}^{2}$ of the closed curve $\Gamma$, global
trigonometric interpolation of the periodic functional density $\varphi$, and
a singular splitting of the kernels in the form

$\frac{i}{4}H_{0}^{(1)}(k|z|)=-\frac{1}{2\pi}J_{0}(k|z|)\log{|z|}+C_{k}(z),$
(18)

together with trapezoidal rules for treatment of regular and periodic
integrands (we note that the function $C_{k}(z)$ in equation (21) is smooth)
as well as closed form formulas for the evaluation of the integrals of
products of the kernel
$\log\left(4\sin^{2}\left(\frac{t-\tau}{2}\right)\right)$ and the Fourier
harmonics (see [15] for a detailed description). Unfortunately, the Kussmaul-
Martensen quadratures cannot be applied directly in the case when the
wavenumber $k$ is complex because the Bessel function $J_{0}(k|z|)$ in the
logarithmic splitting (21) grows exponentially as $|z|\to\infty$ in the case
when the wavenumber $k$ is complex with positive imaginary part. This setback
can be remedied in various ways, either using cutoff functions to mollify the
exponential growth of the Bessel functions [11] or by replacing the Bessel
function $J_{0}(k|z|)$ in equation (21) with a function that has the same
asymptotic behavior as $|z|\to 0$ but decreases exponentially as
$|z|\to\infty$ [31]. The former approach amounts to localizing the logarithmic
splitting only in a neighborhood of the singularity, that is

$\frac{i}{4}H_{0}^{(1)}(k|z|)=\chi(k|z|^{4})\\{-\frac{1}{2\pi}J_{0}(k|z|)\log{|z|}+C_{k}(z)\\}+(1-\chi(k|z|^{4}))\\{\frac{i}{4}H_{0}^{(1)}(k|z|)\\},$
(19)

where $\chi\in C_{0}^{\infty}(\mathbb{R})$ such as $\chi(t)\equiv 1,\ |t|\leq
1/2$ and $\chi(t)\equiv 0,\ |t|\geq 1$. On the other hand, the latter approach
[31] starts with the Taylor expansion

$e^{|k|z}J_{0}(k|z|)=\sum_{\ell=0}^{L}c_{\ell}|kz|^{\ell}+\mathcal{O}(|kz|^{L+1}),\quad
z\to 0^{+}$

where

$c_{\ell}=\sum_{m=0}^{[\ell/2]}\frac{e^{(2\arg{k}+\pi)mi}}{(m!)^{2}(\ell-2m)!}$

We define then

$\widetilde{J}_{0}(k|z|)=e^{-|k|z}\sum_{\ell=0}^{L}c_{\ell}|kz|^{\ell}$ (20)

and the modified splitting

$\frac{i}{4}H_{0}^{(1)}(k|z|)=-\frac{1}{2\pi}\widetilde{J}_{0}(k|z|)\log{|z|}+\widetilde{C}_{k}(z),$
(21)

where $\widetilde{J}_{0}(k|z|)$ decays exponentially as $|z|\to\infty$ when
$\Im{k}>0$.

_Alpert quadratures._ In contrast to the Kusmaul-Martensen quadratures that
rely on the logarithmic splitting techniques described above, the Alpert
quadratures [1] can handle the same kernel singularities in a manner that is
agnostic to the argument of the Hankel functions. We note that the use of
Alpert quadratures in the CQ context was already advocated in [30] where it
was shown that these quadratures can handle in a seamless manner the
discretization of the BIOs involved in the transmission analogues of the
ensemble of Laplace domain Helmholtz problems (6) and respectively (13). We
describe briefly in what follows the details of Alpert quadratures, following
the exposition in [23]. In essence, these quadratures are designed to handle
logarithmic singularities in the following manner

$\displaystyle\int_{0}^{T}\rho(x_{i},x^{\prime})\mu(x^{\prime})dx^{\prime}$
$\displaystyle\approx
h\sum_{p=0}^{N-2a}\rho(x_{i},x_{i}+ah+ph)\mu(x_{i}+ah+ph)$
$\displaystyle+h\sum_{p=1}^{m}w_{p}\rho(x_{i},x_{i}+\chi_{p}h)\mu(x_{i}+\chi_{p}h)$
$\displaystyle+h\sum_{p=1}^{m}w_{p}\rho(x_{i},x_{i}+T-\chi_{p}h)\mu(x_{i}+T-\chi_{p}h),$

where $\\{x_{i}:=ih\\}$ is an equispaced mesh on the interval $[0,T]$, the
kernel $\rho(x_{i},x^{\prime})$ has a logarithmic singularity as
$x^{\prime}\to x_{i}$, $\mu$ is assumed to be a regular enough density
function, and the weights $w_{p}$ and the nodes $\chi_{p}$ are selected so
that the ensuing quadratures achieves a prescribed rate of convergence
provided that the integer quantities $a$ and $m$ are chosen appropriately [1].
The endpoint correction nodes $\chi_{p}$ are typically not integers, and as
such the Alpert quadrature requires evaluation of the density $\mu$ outside of
the equispaced mesh $\\{x_{i}\\}$. This is achieved by resorting to Lagrange
interpolation of order $m+3$ with $m+4$ equispaced nodes around $x_{i}$ and
$x_{i}+T$

$\sigma(x)=\sum_{q=0}^{m+3}L_{q}^{(x_{i})}(x)\mu(x_{i}+qh)\quad\quad\text{where}\quad\quad
L_{q}^{(x_{i})}(x)=\prod_{r=0}\frac{x-(x_{i}+rh)}{(x_{i}+qh)-(x_{i}+rh)}.$

Consequently, the Alpert quadrature takes on the explicit form

$\displaystyle\int_{0}^{T}\rho(x_{i},x^{\prime})\mu(x^{\prime})dx^{\prime}$
$\displaystyle\approx
h\sum_{p=0}^{N-2a}\rho(x_{i},x_{i}+ah+ph)\mu(x_{i}+ah+ph)$
$\displaystyle+h\sum_{q=0}^{m+3}\left(\sum_{p=1}^{m}w_{p}\rho(x_{i},x_{i}+\chi_{p}h)L_{q}^{(x_{i})}(x_{i}+\chi_{p})\right)\mu(x_{i}+qh)$
$\displaystyle+h\sum_{q=0}^{m+3}\left(\sum_{p=1}^{m}w_{p}\rho(x_{i},x_{i}+T-\chi_{p})L_{q}^{(x_{i}+T)}(x_{i}+T-\chi_{p}h)\right)\mu(x_{i}+T-qh),$

which can be readily re-expressed in matrix-vector product form [23]. The
Alpert quadratures described above can be applied directly to obtain high-
order Nyström discretizations of the parametrized versions of the BIE (17)
(because the kernels of the latter exhibit a logarithmic singularity) for any
value of the wavenumber $k$ such that $\Im{k}\geq 0$. Clearly, the high-order
of Alpert quadratures Nyström discretizations of the BIE (17) also depends on
the regularity of the density function $\varphi$. Extensive comparisons of
Alpert quadrature Nyström discretizations with other high-order Nyström
discretizations were performed in [23] in the case when the boundary $\Gamma$
is a smooth closed curve. We present in what follows a simple extension of
Alpert quadrature Nyström discretizations of the single layer BIE (17) to the
case when the boundaries $\Gamma$ are either (a) Lipschitz and closed or (b)
open arcs. To this end we rely on sigmoid transforms [27] and weighted
formulations [18] as we explain next. We mention that this procedure is
readily extended to the case of the double layer (and combined field) BIE
formulations of the Helmholtz equation (16).

We assume that the curve $\Gamma$ is piecewise smooth has corners at
$x_{1},x_{2},\ldots,x_{P}$; in the case when $\Gamma$ is an open arc, we
denote by $x_{1}$ and $x_{P}$ its end points. We assume that the boundary
curve $\Gamma$ has a $2\pi$ periodic parametrization so that each of the
curved segments $[x_{j},x_{j+1}]$ is parametrized by
$\gamma(t)=(\gamma_{1}(w_{j}(t)),\gamma_{2}(w_{j}(t)))$ with
$t\in[T_{j},T_{j+1}]$ (so that $x_{j}=\gamma(T_{j})$, and
$\gamma_{j}:\mathbb{R}\to\mathbb{R}$ are $2\pi$ periodic with
$(\gamma_{1}^{\prime}(t))^{2}+(\gamma_{2}^{\prime}(t))^{2}>0$ for all $t$)
where $0=T_{0}<T_{1}<T_{2}<\ldots<T_{P}<T_{P+1}=2\pi$ in the case when
$\Gamma$ is closed and $0=T_{1}<T_{2}<\ldots<T_{P}=2\pi$ when $\Gamma$ is an
open arc. The maps $w_{j}:[T_{j},T_{j+1}]\to[T_{j},T_{j+1}],\ 0\leq j\leq P$
are the sigmoid transforms introduced by Kress [27]

$\displaystyle w_{j}(s)$ $\displaystyle=$
$\displaystyle\frac{T_{j+1}[v_{j}(s)]^{\sigma}+T_{j}[1-v_{j}(s)]^{\sigma}}{[v_{j}(s)]^{\sigma}+[1-v_{j}(s)]^{\sigma}},\
T_{j}\leq s\leq T_{j+1},\ 0\leq j\leq P$ (22) $\displaystyle v_{j}(s)$
$\displaystyle=$
$\displaystyle\left(\frac{1}{\sigma}-\frac{1}{2}\right)\left(\frac{T_{j}+T_{j+1}-2s}{T_{j+1}-T_{j}}\right)^{3}+\frac{1}{\sigma}\
\frac{2s-T_{j}-T_{j+1}}{T_{j+1}-T_{j}}+\frac{1}{2}$

where $\sigma>2$. The functions $w_{j}$ are smooth and increasing bijective
functions on each of the intervals $[T_{j},T_{j+1}]$ for $0\leq j\leq P$, with
$w_{j}^{(k)}(T_{j})=w_{j}^{(k)}(T_{j+1})=0$ for $1\leq k\leq\sigma-1$ and all
$0\leq j\leq P$. With the aid of the graded meshes just introduced, we define
_weighted_ densities $\varphi^{w}$ in the following manner

$\varphi^{w}(t):=|\gamma^{\prime}(t)|\varphi(\gamma(t)),\quad 0\leq t\leq
2\pi$

as well as _weighted_ parametrized single layer formulations

$\frac{i}{4}\int_{0}^{2\pi}H_{0}^{(1)}(k|\gamma(t)-\gamma(\tau)|)\varphi^{w}(\tau)d\tau=g(\gamma(t)),\quad
t\in[0,2\pi].$ (23)

On account of the vanishing property of the sigmoid transforms at the
breakpoints $T_{j}$, the weighted density $\varphi^{w}$ are $2\pi$ periodic
functions whose regularity can be control by the value of $p$ in the graded
mesh used. As such, we apply the Alpert quadrature described above to produce
high-order Nyström discretizations of the weighted single layer formulations
(23). Finally, in order to avoid the evaluation of the density function at
corner points, we employ the equispaced mesh
$t_{j}=(j-1/2)\left(\frac{2\pi}{N}\right),1\leq j\leq N$ in the Alpert Nyström
method. We present in Table 1 the errors in the near field achieved by the
Alpert discretization of the weighted formulation (23) for real wavenumbers
with parameters $a=2$ and $m=3$ for various Lipschitz and open arc $\Gamma$
configurations. Specifically, we considered the cases when $\Gamma$ is (i) the
teardrop geometry $x(t)=\left(2\sin{\frac{t}{2}},-\beta\sin{t}\right),\
\beta=\tan{\frac{\alpha\pi}{2}},\ \alpha=1/2,\ 0\leq t\leq 2\pi$, (ii) the
boomerang geometry
$x(t)=\left(-\frac{2}{3}\sin{\frac{3t}{2}},-\sin{t}\right),\ \ 0\leq t\leq
2\pi$ [15]—both domains have diameters equal to 2, and (iii) the strip
connecting $(-1,0)$ to $(1,0)$, and we report errors in the near field
evaluated at $512$ equispaced points placed on circles located at distance $1$
from the scatterers. We focus on the case when the wavenumber is $k=8$ under
normal plane wave incidence, because in the case of real wavenumbers extremely
accurate reference solutions can be produced via alternative spectral Nyström
solvers (e.g. those described in [15] in the case of closed curves and
respectively in [13] in the case of open arcs). We report the estimated orders
of convergence achieved by the Alpert methods under the acronym “e.o.c.”. We
mention that similar convergence rates were observed in the case when the
wavenumber $k$ is complex in the formulation (23).

$N$ | Teardrop | Boomerang | Strip | V-shaped
---|---|---|---|---
| $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c.
64 | 4.9 $\times$ $10^{-4}$ |  | 1.9 $\times$ $10^{-3}$ |  | 4.1 $\times$ $10^{-6}$ |  | 8.5 $\times$ $10^{-5}$ |
128 | 1.2 $\times$ $10^{-5}$ | 4.63 | 6.1 $\times$ $10^{-5}$ | 4.95 | 1.9 $\times$ $10^{-7}$ | 4.37 | 3.2 $\times$ $10^{-6}$ | 4.71
256 | 7.3 $\times$ $10^{-7}$ | 4.75 | 2.3 $\times$ $10^{-6}$ | 4.72 | 1.2 $\times$ $10^{-8}$ | 4.02 | 1.3 $\times$ $10^{-7}$ | 4.57
512 | 2.8 $\times$ $10^{-8}$ | 4.66 | 9.0 $\times$ $10^{-8}$ | 4.67 | 7.9 $\times$ $10^{-10}$ | 3.93 | 6.7 $\times$ $10^{-9}$ | 4.35
Table 1: Errors in the near field and estimated orders of convergence
corresponding to the Alpert discretization of the weighted formulation (23)
for real wavenumber $k=8$ and plane wave normal incidence. We used the
following Alpert quadrature parameters: $\sigma=4$ in the sigmoid transform,
and respectively $a=2,m=3$.

_QBX quadratures._ The quadrature by expansion (QBX) method is another
powerful method to discretize the boundary layer potentials and operators that
arise in Helmholtz problems [26]. Just like the Alpert quadratures, QBX can be
applied without fuss to discretizations of layer potentials whose associated
Green functions feature either real or complex wavenumbers. We review in what
follows the details of the QBX method in the case of evaluations of single
layer potentials and BIOs. The QBX method is based on the continuity of the
single layer potential across the boundary $\Gamma$

$V(k)[\varphi](x)=\lim_{x^{+}\to x\in\Gamma}SL(k)[\varphi](x^{+}),$

and it relies on smooth extensions of the single layer potentials
$SL(k)[\varphi](x^{+}),\ x^{+}\in\Omega^{+}$ to $\overline{\Omega^{+}}$ via
certain series expansions. Specifically, the method uses expansion centers,
that is for $x\in\Gamma$ we define $x^{\pm}=x\pm\varepsilon n(x)$, where the
vector $n(x)$ is the exterior unit normal at $x\in\Gamma$ and the parameter
$\varepsilon=\varepsilon(x)>0$ is small enough, and relies on the addition
theorem

$H_{0}^{(1)}(k|x-x^{\prime}|)=\sum_{\ell=-\infty}^{\infty}H_{\ell}^{(1)}(k|x^{\prime}-x^{+}|)e^{i\ell\theta^{\prime}}J_{\ell}(k|x-x^{+}|)e^{-i\ell\theta},\
x^{\prime}\in\Gamma$ (24)

where $\theta$ and $\theta^{\prime}$ are the angular coordinates of $x$ and
respectively $x^{\prime}$ in the polar coordinate system centered at $c=x^{+}$
to produce the following series representation for the single layer potential

$V(k)[\varphi](x)=\sum_{\ell=-\infty}^{\infty}\alpha_{\ell}[\varphi]J_{\ell}(k|x-x^{+}|)e^{-i\ell\theta},\quad\alpha_{\ell}[\varphi]:=\frac{i}{4}\int_{\Gamma}H_{\ell}^{(1)}(k|x^{\prime}-x^{+}|)e^{i\ell\theta^{\prime}}\varphi(x^{\prime})ds(x^{\prime}).$
(25)

In practice a truncation parameter $p$ is chosen in the series representation
(25) leading to the following QBX approximation

$V(k)[\varphi](x)\approx\sum_{\ell=-p}^{p}\alpha_{\ell}[\varphi]J_{\ell}(k|x-x^{+}|)e^{-i\ell\theta}$
(26)

where the coefficients $\alpha_{\ell}$ have to be evaluated numerically for
all $-p\leq\ell\leq p$. We note that the integrands in the definition (25) of
the coefficients $\alpha_{\ell}$ do not contain kernel singularities as
$x^{+}\notin\Gamma$. Owing to the fact that the evaluation of the expansion
coefficients $\alpha_{\ell}$ requires integration on the whole boundary
$\Gamma$, this version of QBX is referred to as global; local versions have
also been designed [34], but we will not consider them in this paper.

We will use in what follows panel Chebyshev meshes and Clenshaw-Curtis
quadratures to evaluate numerically the coefficients $\alpha_{\ell}$ in
equation (25). Indeed, we assume that the piecewise smooth curve $\Gamma$ can
be represented as the union of (disjoint) smooth panels $\Gamma_{m}$ in the
form $\Gamma=\bigcup_{m=1}^{M}\Gamma_{m}$ so that the corner points
$x_{1},x_{2},\ldots,x_{P}$ of $\Gamma$ are end points of these panels. We thus
have

$\alpha_{\ell}[\varphi]=\sum_{m}\alpha_{\ell,m}[\varphi],\qquad\alpha_{\ell,m}[\varphi]:=\frac{i}{4}\int_{\Gamma_{m}}H_{\ell}^{(1)}(k|x^{\prime}-x^{+}|)e^{i\ell\theta^{\prime}}\varphi(x^{\prime})ds(x^{\prime}).$

We apply Clenshaw-Curtis quadratures for the evaluation of each of the
integrals in the definition of $\alpha_{\ell,m}$, that is

$\alpha_{\ell,m}[\varphi]\approx\frac{i}{4}\sum_{j=1}^{N_{m}}\omega_{j}H_{\ell}^{(1)}(k|\gamma_{m}(t_{j})-x^{+}|)e^{i\ell\theta^{\prime}_{j}}\widetilde{\varphi}(\gamma_{m}(t_{j}))$
(27)

where we assume that $\Gamma_{m}$ is parametrized in the form
$\Gamma_{m}=\\{\gamma_{m}(t):t\in[-1,1]\\}$ and
$\gamma_{m}:[-1,1]\to\mathbb{R}^{2}$ is smooth,
$\widetilde{\varphi}(\gamma_{m}(t_{j}))=\varphi(\gamma_{m}(t_{j}))|\gamma^{\prime}_{m}(t_{j})|$,
the quadrature points $t_{j}$ are the Chebyshev zero points

$t_{j}:=\cos(\vartheta_{j}),\quad\vartheta_{j}:=\frac{(2j-1)\pi}{2N_{m}},\quad
j=1,\ldots,N_{m}$

and the Fejér quadrature weights $\omega_{j}$ are given by

$\omega_{j}:=\frac{2}{N_{m}}\left(1-2\sum_{q=1}^{[N_{m}/2]}\frac{1}{4q^{2}-1}\cos(2q\vartheta_{j})\right),\quad
j=1,\ldots,N_{m}.$

The steps of the QBX discretization algorithm of the single layer BIE (17) are
summarized in what follows: (1) choose a grid on $\Gamma$ such that
$x_{j,m}=\gamma_{m}(t_{j}),\ 1\leq m\leq M,\ 1\leq j\leq n_{m}$ and
$t_{j}:=\cos(\vartheta_{j}),\quad\vartheta_{j}:=\frac{(2j-1)\pi}{2n_{m}},\quad
j=1,\ldots,n_{m}$; (2) for each $x_{j}$ on $\Gamma$ define the expansion
centers $x_{j}^{\pm}:=x_{j}\pm\varepsilon(x_{j})n(x_{j})$; (3) use the QBX
approximation (26) for a given truncation parameter $p$; and (4) evaluate the
ensuing coefficients $\alpha_{\ell,m}[\varphi]$ according to the Clenshaw-
Curtis quadrature (27) using oversampling, that is choose $N_{m}=\beta n_{m},\
\beta\geq 1$ in equation (27) and use Chebyshev interpolation to access the
values of the density $\varphi^{w}$ on the fine Chebyshev grid with $N_{m}$
nodes on each $\Gamma_{m}$. Finally, we choose
$\varepsilon(x_{j})=\min(|x_{j}-x_{j-1}|,|x_{j}-x_{j+1}|)$ in the definition
of the centers $x_{j}^{\pm}$.

It is also possible to consider weighted versions of the single layer BIO
$V(k)$ by defining weighted functional densities that are more regular in the
neighborhood of the geometric singularities of the curve $\Gamma$.
Specifically we define the weighted densities

$\varphi_{m}^{w}(t):=\varphi(\gamma_{m}(t))|\gamma^{\prime}_{m}(t)|\sqrt{1-t^{2}},\quad-1\leq
t\leq 1,\quad 1\leq m\leq M$ (28)

whose square root weights are chosen to resolve the worst case single layer
density singularities that arise in the neighborhood of open endpoints of
$\Gamma$ (in the case when $\Gamma$ is an open arc the solutions of the single
layer BIE (17) behave asymptotically as $\varphi\sim\mathcal{O}(1/\sqrt{d})$
in the neighborhood of the open ends of the arc, where $d$ is the distance to
the open end of the arc). Incidentally, the incorporation of weighted
densities (28) in the definition of single layer BIO leads to a simpler Fejér
quadrature rule for the evaluation of the expansion coefficients
$\alpha_{\ell,m}$ in the form

$\alpha_{\ell,m}[\varphi^{w}]\approx\frac{i\pi}{4N_{m}}\sum_{j=1}^{N_{m}}H_{\ell}^{(1)}(k|\gamma_{m}(t_{j})-x^{+}|)e^{i\ell\theta^{\prime}_{j}}\varphi^{w}_{m}(t_{j}).$
(29)

In addition, the condition numbers of the QBX Nyström matrices corresponding
to the weighted expansion coefficients $\alpha_{\ell,m}[\varphi^{w}]$
according to equations (29) are about one order of magnitude smaller than
those corresponding to the expansion coefficients $\alpha_{\ell,m}[\varphi]$
according to equations (27). Also, should iterative solvers such as GMRES be
used for the solution of the QBX Nyström linear systems, the iterations counts
are more favorable in the case of the single layer formulation that
incorporates the weighted unknown $\varphi^{w}$ defined in equation (28). For
these reasons, the use of weighted single layer formulations is preferred.

We report in Table 2 numerical results concerning QBX quadrature Nyström
discretizations of the single layer formulation of the Helmholtz equation (16)
using the weighted unknown (28), a number of Chebyshev panels that coincides
with the number of corners in the case of closed scatterers and respectively
the number of corners plus one in the case of open arcs (that is one global
Chebyshev mesh in the case of the teardrop, boomerang, and strip, and two
Chebyshev meshes in the case of the V-shaped scatterer), as well as the
Clenshaw-Curtis quadratures (29) to evaluate the QBX expansion coefficients
$\alpha_{\ell,m}[\varphi^{w}]$ in the expansion (25). We focused on the case
of real wavenumbers, so that highly accurate reference solutions can be
produced by alternative high-order discretizations [15, 13]; qualitatively
similar results are obtained in the cases when the wavenumber is complex. In
these numerical experiments we fixed the QBX expansion order parameter $p=8$
and the upsampling factor $\beta=6$, and we present results as the
discretization Chebyshev mesh is increasingly refined. We note that in the
case of open arcs the convergence of the Nyström QBX discretizations with
respect to the size of the mesh appears to be second order, which can be
accounted for by the fact that the endpoints of the open arc are the nearest
singularities of the smooth extension onto $\Gamma$ of single layer potentials
evaluated at nearby exterior points.

$n_{1}$ | Teardrop | Boomerang | Strip | V-shaped
---|---|---|---|---
| $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c.
32 | 1.3 $\times$ $10^{-2}$ |  | 4.4 $\times$ $10^{-2}$ |  | 2.6 $\times$ $10^{-4}$ |  | 4.7 $\times$ $10^{-4}$ |
64 | 6.5 $\times$ $10^{-6}$ | 10.9 | 1.0 $\times$ $10^{-5}$ | 12.01 | 5.9 $\times$ $10^{-5}$ | 2.18 | 1.1 $\times$ $10^{-4}$ | 1.99
128 | 1.0 $\times$ $10^{-6}$ | 2.70 | 4.9 $\times$ $10^{-8}$ | 7.73 | 1.3 $\times$ $10^{-5}$ | 2.10 | 2.9 $\times$ $10^{-5}$ | 1.99
256 | 1.6 $\times$ $10^{-7}$ | 2.64 | 2.2 $\times$ $10^{-10}$ | 7.80 | 3.3 $\times$ $10^{-6}$ | 2.07 | 7.4 $\times$ $10^{-6}$ | 1.99
Table 2: Errors in the near field and estimated orders of convergence obtained
from the QBX Nyström quadrature discretization of the single layer formulation
for real wavenumber $k=8$ and plane wave normal incidence. We used the
weighted unknown (28) and the Clenshaw-Curtis quadratures (29) in the
expansion (25), as well as the QBX parameters $p=8$ and $\beta=6$.

We continue to investigate in Table 3 the role played by nearby singularities
in QBX field expansions. Specifically, we present near field accuracy errors
achieved by Nyström QBX discretizations of the single layer formulation of the
Helmholtz equation with real wavenumber $k=8$ in the case of teardrop
scatterers with very acute apertures whose parametrizations are given by
$x(t)=\left(2\sin{\frac{t}{2}},-\beta\sin{t}\right),\
\beta=\tan{\frac{\alpha\pi}{2}},\ \alpha=1/32,\ 0\leq t\leq 2\pi$. We used in
our numerical experiments both one panel (a global Chebyshev mesh—second
column) as well as 8 Chebyshev panels dyadically refined near the corner
(third column), we fixed the same expansion parameter $p=6$ and the up
sampling factor $\beta=4$ for each panel, and we refined the size $n_{m},\
1\leq m\leq M$ of the Chebyshev meshes. We observe that the QBX
discretizations give rise to second order convergent Nyström solvers (with
respect to the size $n_{m}$ of Chebyshev meshes), and we note that Alpert
solvers (first column in Table 3) do not exhibit a deterioration of the order
of convergence upon Chebyshev mesh refinement for very acute corners.
Nevertheless, the accuracy levels of the Alpert and QBX Nyström
discretizations appear to be quite similar for the same mesh size. On the
other hand, following [26], we present in the fourth column of Table 3 the
convergence properties of QBX Nyström discretizations from a different
perspective, that is we considered 64 Chebyshev panels dyadically refined near
the corner so that the smallest size of the panels is of the order $10^{-9}$,
the size of the Chebyshev meshes on each panel was fixed at $n_{m}=8$, and we
progressively increased the values of the parameters $(p,\beta)$. In this
setting the QBX Nyström discretizations exhibit high order convergence with
respect to the parameters $(p,\beta)$. We applied in the experiments presented
in Table 4 the same dyadic panel refinement QBX strategy in the case of
scattering off strips and high order convergence with respect to the
parameters $(p,\beta)$ is observed again. However, since the single layer
formulation is a BIE of the first kind, the condition numbers and numbers of
GMRES iterations required for convergence of the QBX Nyström matrices tend to
grow with the number $M$ of panels used for the same global size of the
discretization. One possibility to mitigate this issue is to resort to
recursively compressed inverse preconditioning (RCIP) [25], which will be
pursued elsewhere. For the scatterers considered in this paper and the levels
of accuracy of CQ methods for the solution of the time domain scattering (e.g.
$10^{-8}$), the use of geometrically large panels in connection to QBX appears
to be adequate as either the direct or iterative solvers of the linear systems
involving QBX Nyström matrices are well behaved.

Alpert | QBX $M=1$ | QBX $M=8$ | QBX $M=64,n_{m}=8$
---|---|---|---
$N$ | $\varepsilon_{\infty}$ | e.o.c. | $n_{1}$ | $\varepsilon_{\infty}$ | e.o.c. | $n_{m}$ | $\varepsilon_{\infty}$ | e.o.c. | $(p,\beta)$ | $\varepsilon_{\infty}$
64 | 9.8 $\times$ $10^{-3}$ |  | 64 | 2.0 $\times$ $10^{-4}$ |  | 8 | 4.3 $\times$ $10^{-4}$ |  | (1,1) | 5.8 $\times$ $10^{-2}$
128 | 3.5 $\times$ $10^{-4}$ | 4.81 | 128 | 4.9 $\times$ $10^{-5}$ | 2.02 | 16 | 1.1 $\times$ $10^{-4}$ | 2.04 | (2,2) | 7.7 $\times$ $10^{-3}$
256 | 1.0 $\times$ $10^{-5}$ | 5.01 | 256 | 1.1 $\times$ $10^{-5}$ | 2.02 | 32 | 2.6 $\times$ $10^{-5}$ | 2.00 | (4,2) | 3.9 $\times$ $10^{-4}$
512 | 5.7 $\times$ $10^{-7}$ | 4.25 | 512 | 2.9 $\times$ $10^{-6}$ | 2.01 | 64 | 6.4 $\times$ $10^{-6}$ | 2.01 | (6,4) | 2.5 $\times$ $10^{-5}$
Table 3: Errors in the near field and estimated orders of convergence obtained from the Alpert and QBX Nyström quadrature discretizations of the single layer formulation for real wavenumber $k=8$ and plane wave normal incidence impinging on a teardrop with a very acute aperture. We used the parameters $\sigma=4,\ a=2,\ m=3$ for the Alpert quadrature, the single layer formulation with the weighted unknown (28), $M=1$ and $M=8$ panels with QBX parameters $p=6$ and $\beta=4$ and increasingly refined Chebyshev meshes, as well as $M=64$ dyadically refined panels around the corner with Chebyshev meshes of a fixed size $n_{m}=8$ and various values of the $(p,\beta)$ QBX parameters. QBX $M=16,n_{m}=16$ | QBX $M=64,n_{m}=8$
---|---
$(p,\beta)$ | $\varepsilon_{\infty}$ | $(p,\beta)$ | $\varepsilon_{\infty}$
(1,1) | 1.4 $\times$ $10^{-2}$ | 1,1) | 4.9 $\times$ $10^{-2}$
(2,2) | 7.5 $\times$ $10^{-4}$ | (2,2) | 4.9 $\times$ $10^{-3}$
(6,4) | 1.1 $\times$ $10^{-5}$ | (6,4) | 1.3 $\times$ $10^{-4}$
(8,4) | 6.5 $\times$ $10^{-7}$ | (8,4) | 5.3 $\times$ $10^{-6}$
Table 4: Errors in the near field obtained from the QBX Nyström quadrature
discretization of the single layer formulation for real wavenumber $k=8$ and
plane wave normal incidence impinging on the strip scatterer. We used the
formulation involving the weighted unknown (28), 16 and respectively 64
dyadically refined panels around the end points of the strip, Chebyshev grids
of size $n_{m}=16$ and respectively $n_{m}=8$ on each panel, Clenshaw-Curtis
quadratures (29) in the expansion (25), and various QBX parameters
$(p,\beta)$.

###  3.2 Neumann Boundary Conditions

The CQ methodology can be also applied to the case when the wave equation (1)
has Neumann boundary conditions. In this case we would have to solve an
ensemble of Helmholtz problems with Neumann boundary conditions

$\begin{cases}\Delta u(x)+k^{2}u(x)=0\quad x\in\Omega\\\ \frac{\partial
u(x)}{\partial n}=f(x)\quad x\in\Gamma\\\
u(x)=\sum_{m=-\infty}^{\infty}u_{m}H_{m}^{(1)}(k|x|)e^{im\theta},\quad&|x|>r_{0},\quad
x=|x|(\cos{\theta},\sin{\theta}).\end{cases}$ (30)

for a family of complex wave numbers $k$ such that $\Im{k}\geq 0$.

_Single layer formulation._ We begin with the case when $\Gamma$ is a closed
boundary and we seek the solution of the Helmholtz problem (30) in the form of
the single layer potential

$u(x)=\int_{\Gamma}G_{k}(x-y)\varphi(y)ds(y),\ x\in\Omega^{+},$

and we obtain the following BIE of the second kind

$-\frac{\varphi}{2}+K^{\top}(k)[\varphi]=f\quad{\rm on}\ \Gamma,\quad
K^{\top}(k)[\varphi](x):=\int_{\Gamma}\frac{\partial G_{k}(x-y)}{\partial
n(x)}\varphi(y)ds(y),\ x\in\Gamma.$ (31)

We note that the BIE formulation (31) is uniquely solvable whenever
$\Im{k}>0$. In order to make the BIE formulation (31) amenable to Nyström
discretizations based on Alpert quadratures, we use the parametrization
$\gamma$ of the boundary curve $\Gamma$ that incorporates the sigmoid
transforms and we employ the following weighted formulation

$-\frac{1}{2}\varphi^{w}(t)-\frac{ik}{4}\int_{0}^{2\pi}H_{1}^{(1)}(k|\gamma(t)-\gamma(\tau)|)\frac{(\gamma(t)-\gamma(\tau))\cdot(\gamma^{\prime}(t))^{\perp}}{|\gamma(t)-\gamma(\tau)|}\varphi^{w}(\tau)d\tau=f(\gamma(t))|\gamma^{\prime}(t)|$
(32)

where
$(\gamma^{\prime}(t))^{\perp}=(\gamma^{\prime}_{2}(t),-\gamma^{\prime}_{1}(t))$
and again here $\varphi^{w}(t):=|\gamma^{\prime}(t)|\varphi(\gamma(t))$. On
the other hand, the application of QBX quadratures to the evaluation of the
BIO $K^{\top}(k)[\varphi]$ relies on the classical jump formulas for the
gradients of single layer potentials

$\nabla_{x}SL(k)[\varphi](x)=\lim_{x^{\pm}\to
x}\nabla_{x^{\pm}}SL(k)[\varphi](x^{\pm})\pm\frac{1}{2}n(x)\varphi(x),\
x\in\Gamma,\ x^{\pm}=x\pm\varepsilon n(x),$

which amounts to differentiating term by term the expressions

$SL(k)[\varphi](x^{\pm})\approx\sum_{\ell=-p}^{p}\alpha_{\ell}[\varphi]J_{\ell}(k|x^{\pm}-x|)e^{-i\ell\theta^{\pm}}$

in the polar coordinate systems around the expansion centers [26].
Furthermore, per the prescriptions in [26], we average the derivatives of the
expansions above with respect to both $x^{+}$ and $x^{-}$ in order to evaluate
the BIO expressions in equation (31). Specifically, the QBX quadrature takes
on the form

$\displaystyle K^{\top}(k)[\varphi](x)$ $\displaystyle\approx$
$\displaystyle-\frac{1}{2}\sum_{\ell=-p}^{p}k\
\alpha_{\ell}[\varphi]\left(-J_{\ell+1}(k|x^{+}-x|+\frac{\ell}{k|x^{+}-x|}J_{\ell}(k|x^{+}-x|)\right)e^{-i\ell\theta^{+}}$
$\displaystyle+$ $\displaystyle\frac{1}{2}\sum_{\ell=-p}^{p}k\
\alpha_{\ell}[\varphi]\left(-J_{\ell+1}(k|x^{-}-x|+\frac{\ell}{k|x^{-}-x|}J_{\ell}(k|x^{-}-x|)\right)e^{-i\ell\theta^{-}}$

where the coefficients $\alpha_{\ell}[\varphi]$ are in turn computed via the
formulas (27). It is straightforward to modify the formulas above to
accommodate for a weighted formulation whose unknown is defined as in formula
(28). We present in Tables 5 and 6 numerical results concerning Alpert and QBX
discretizations of the second kind formulation (31) using weighted unknowns in
the case when the wavenumber $k=8$ using reference solutions produced by the
high-order kernel-splitting Nyström discretization in [2]. It is also possible
to apply the dyadic refinement around corners strategy in connection to QBX
discretizations in the case of the second kind formulation (31), and we
obtained qualitatively similar results to those in Table 3. On account of the
BIE formulation used being of the second kind, both the condition numbers of
the Nyström QBX matrices and their spectral properties in connection with
GMRES iteration counts are virtually independent of the number of panels $M$
(and thus their geometrical size).

$N$ | Teardrop | Boomerang
---|---|---
| $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c.
64 | 2.4 $\times$ $10^{-4}$ |  | 6.2 $\times$ $10^{-3}$ |
128 | 5.7 $\times$ $10^{-6}$ | 5.40 | 1.4 $\times$ $10^{-4}$ | 5.40
256 | 1.9 $\times$ $10^{-7}$ | 4.87 | 4.7 $\times$ $10^{-6}$ | 4.91
512 | 7.7 $\times$ $10^{-9}$ | 4.67 | 1.8 $\times$ $10^{-7}$ | 4.71
Table 5: Errors in the near field and estimated orders of convergence corresponding to the Alpert discretization of the weighted formulation (32) for real wavenumber $k=8$ and plane wave normal incidence. We used the following Alpert quadrature parameters: $\sigma=4$ in the sigmoid transform, and respectively $a=2,m=3$. $n_{1}$ | Teardrop $M=1$ | Boomerang $M=1$
---|---|---
| $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c.
32 | 2.5 $\times$ $10^{-2}$ |  | 3.3 $\times$ $10^{-1}$ |
64 | 3.0 $\times$ $10^{-3}$ | 3.07 | 2.6 $\times$ $10^{-2}$ | 3.65
128 | 3.6 $\times$ $10^{-4}$ | 3.07 | 1.6 $\times$ $10^{-3}$ | 4.03
256 | 3.1 $\times$ $10^{-5}$ | 3.54 | 8.9 $\times$ $10^{-5}$ | 4.16
Table 6: Errors in the near field and estimated orders of convergence obtained
with the QBX Nyström quadrature discretization of the weighted single layer
formulation for real wavenumber $k=8$ and plane wave normal incidence. We used
a weighted unknown, one global Chebyshev grid, the QBX expansion (3.2) and the
Clenshaw-Curtis quadratures (29) in the expansion (25), as well as the QBX
parameters $p=4$ and $\beta=4$.

_Double layer formulation._ We note that the BIE formulation (31) is
restricted to closed boundaries $\Gamma$ only. An alternative formulation that
is applicable to both closed and open boundaries can be derived if we seek the
solution of (30) in the form of a double layer potential

$u(x)=\int_{\Gamma}\frac{\partial G_{k}(x-y)}{\partial n(y)}\varphi(y)ds(y),\
x\in\Omega^{+}.$

Then the enforcement of Neumann boundary conditions results in the following
BIE

$N(k)[\varphi](x)=f(x),\ x\in\Gamma$ (34)

where the operator $N(k):H^{1/2}(\Gamma)\to H^{-1/2}(\Gamma)$ is a BIO which
has to be understood in the sense of Hadamard finite parts when $\Gamma$ is a
closed curve

$N(k)[\varphi](x)=f.p.\int_{\Gamma}\frac{\partial^{2}G_{k}(x-y)}{\partial
n(x)\partial n(y)}\varphi(y)ds(y),\ x\in\Gamma.$ (35)

We can recast the operator $N(k)$ in the following form

$\displaystyle N(k)[\varphi](x)$
$\displaystyle=\int_{\Gamma}\frac{\partial^{2}(G_{k}-G_{0})(x-y)}{\partial
n(x)\partial n(y)}\varphi(y)ds(y)+N(0)[\varphi](x)$ $\displaystyle
N(0)[\varphi](x)$
$\displaystyle:=f.p.\int_{\Gamma}\frac{\partial^{2}G_{0}(x-y)}{\partial
n(x)\partial n(y)}\varphi(y)ds(y)$ (36)

where $G_{0}$ is the Green’s function of the Laplace equation:

$G_{0}(x)=-\frac{1}{2\pi}\log{|x|}.$

The integral kernel of the operator $N(k)-N(0)$ has a logarithmic singularity
[2], and as such, the Nyström discretization of the operator $N(k)-N(0)$ is
amenable to Alpert quadratures. On the other hand, the Hadamard finite parts
integral operator term in the definition of the hypersingular BIO $N(0)$ can
be manipulated using integration by parts techniques starting from the
identity

$\frac{\partial^{2}G_{0}(x-y)}{\partial n(x)\partial
n(y)}=-\frac{\partial^{2}G_{0}(x-y)}{\partial t(x)\partial t(y)},$

where $t(x)$ is the unit tangent on $\Gamma$ at $x$. In the case when $\Gamma$
is at least $C^{3}$, the discretization of the Hadamard finite parts integral
operator $N(0)$ can be performed using the fact that the kernel

$N_{0}(t,s):=\frac{\partial^{2}}{\partial t\ \partial
s}\left(\log(|\gamma(t)-\gamma(s)|^{2})-\log\left(4\sin^{2}\left(\frac{t-s}{2}\right)\right)\right)$

is a bounded function in the variables $(t,s)$, assuming a
$\gamma:[0,2\pi]\to\mathbb{R}^{2}$ parametrization of the closed curve
$\Gamma$ [28]. Indeed, denoting

$\widetilde{N}(t,s)=\frac{(\gamma(t)-\gamma(s))\cdot\gamma^{\prime}(t)\
(\gamma(t)-\gamma(s))\cdot\gamma^{\prime}(s)}{|\gamma(t)-\gamma(s)|^{2}}$

we have

$N_{0}(t,s)=-\frac{2}{|\gamma(t)-\gamma(s)|^{2}}(\gamma^{\prime}(t)\cdot\gamma^{\prime}(s)-2\widetilde{N}(t,s))-\frac{1}{2\sin^{2}\left(\frac{t-s}{2}\right)}$

with

$N_{0}(t,t)=-\frac{(\gamma^{\prime}(t)\cdot\gamma^{\prime\prime}(t))^{2}}{(\gamma^{\prime}(t)\cdot\gamma^{\prime}(t))^{2}}+\frac{1}{2}\frac{\gamma^{\prime\prime}(t)\cdot\gamma^{\prime\prime}(t)}{\gamma^{\prime}(t)\cdot\gamma^{\prime}(t)}+\frac{1}{3}\frac{\gamma^{\prime}(t)\cdot\gamma^{\prime\prime\prime}(t)}{\gamma^{\prime}(t)\cdot\gamma^{\prime}(t)}-\frac{1}{6}.$

Finally, the Hilbert transform operators

$H[g](t):=\frac{1}{2\pi}\int_{0}^{2\pi}\cot{\frac{t-s}{2}}{g}^{\prime}({x}(s))ds$

can be discretized using $2\pi$ trigonometric interpolation [28]. In the case
when $\Gamma$ is an open smooth arc, the numerical scheme above has to be
slightly modified in order to accommodate for the different mapping property
of the BIO $N(k)$, that is $N(k):\widetilde{H}^{1/2}(\Gamma)\to
H^{-1/2}(\Gamma)$. Thus, the density function $\varphi$ has to vanish at the
end points of the arc $\Gamma$, which has to be explicitly accounted for in
the numerical scheme. Following the procedure in [14], we assume a
parametrization of $\Gamma$ of the form $\gamma:[-1,1]\to\mathbb{R}^{2}$,
followed by the cosine change of variables and an odd periodic extension of
the density and the parametrized integral equations on the interval
$[0,2\pi]$. Finally, the ensuing discretized integral equations are projected
back onto the grid points on the interval $(0,\pi)$. We present in Table 7
numerical results concerning Alpert quadratures Nyström discretizations of the
hypersingular BIE equation (34) for smooth closed (specifically the kite
geometry in [15] and the starfish geometry [23]) and open boundaries $\Gamma$
in the case of real wavenumbers; the errors were computed using reference
solutions produced via the high-order kernel-splitting Nyström discretization
presented in [2]. Similar orders of convergence are observed in the case when
the wavenumbers are complex.

In the case when $\Gamma$ is Lipschitz, the discretization above used in
conjunction with sigmoid transforms (22) does not result in high-order
convergence. The lack of high-order convergence is attributable to the fact
that the expressions $N_{0}(t,s)$ become unbounded for $s=t$. A similar
situation occurs when Maue’s type integration by parts formulas are employed
to recast the operators $N(k)$ in terms of Cauchy Principal Value integral
operators. We mention that in the case when $\Gamma$ is closed, the evaluation
of hypersingular operators can be entirely bypassed through the use of
Calderón identities when dealing with the BIO $N(k)$ [3]. However, the
extension of the Calderón identities to the case when $\Gamma$ is an open arc
requires that the singularities of density functions near the end points be
explicitly accounted for [13]. We present in what follows a reformulation of
the hypersingular BIE (34) whose discretization does not require numerical
differentiation, and which is amenable to QBX discretizations that are
oblivious to whether $\Gamma$ is open or closed.

$N$ | Kite | Starfish | Strip
---|---|---|---
| $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c.
64 | 9.8 $\times$ $10^{-4}$ |  | 3.6 $\times$ $10^{-4}$ |  | 1.5 $\times$ $10^{-5}$ |
128 | 2.7 $\times$ $10^{-6}$ | 5.07 | 1.3 $\times$ $10^{-5}$ | 4.71 | 5.0 $\times$ $10^{-7}$ | 4.94
256 | 1.1 $\times$ $10^{-6}$ | 4.70 | 5.4 $\times$ $10^{-7}$ | 4.66 | 1.9 $\times$ $10^{-8}$ | 4.71
512 | 7.7 $\times$ $10^{-9}$ | 4.63 | 2.2 $\times$ $10^{-8}$ | 4.62 | 7.6 $\times$ $10^{-10}$ | 4.66
Table 7: Errors in the near field and estimated orders of convergence obtained
when the Alpert discretization is applied to the weighted formulation (34) in
the case of real wavenumber $k=8$ and normal plane wave incidence. We used the
following Alpert quadrature parameteres: $\sigma=4$ in the sigmoid
transform,and respectively $a=2,m=3$.

_Current and Charge formulations and QBX discretizations._ We start with the
Maue’s integration by parts recasting of the BIO $N(k)$ in the form

$N(k)[\varphi](x)=k^{2}\int_{\Gamma}G_{k}(x-y)\hat{n}(x)\cdot\hat{n}(y)\varphi(y)ds(y)+\partial_{s}^{x}\int_{\Gamma}G_{k}(x-y)\partial_{s}\varphi(y)ds(y),\
x\in\Gamma$ (37)

where $\partial_{s}$ denote tangential differentiation, where $\varphi\in
H^{1/2}(\Gamma)$ when $\Gamma$ is closed and
$\varphi\in\widetilde{H}^{1/2}(\Gamma)$ when $\Gamma$ is open (we not that the
integration by parts is still justified in the latter case since $\varphi$
vanishes at the endpoints of $\Gamma$). We shall use in what follows the
$\hat{\cdot}$ notation to denote vector quantities. We employ in what follows
the Current and Charge formulation of (34) starting from equation (37) which
is a variant of the formulations pioneered in [35] for 3D Maxwell equations
and subsequently studied in the excellent contributions [9, 21]. We mention
that extended BIE formulations that use charges as unknowns were recently
analyzed in the context of transmission problems with complex wavenumbers
[24]. First, we introduce the _charge_ unknown $\rho$ which is related to the
_current_ unknown $\varphi$ in the following manner

$\partial_{s}\varphi=ik\rho.$ (38)

Replacing the tangential derivative $\partial_{s}\varphi$ by the new unknown
$\rho$ in the Maue formula (37) and enforcing the constraint (3.2) in an
integral form, we arrive at the system of BIE

$\displaystyle-
ik\int_{\Gamma}G_{k}(x-y)\hat{n}(x)\cdot\hat{n}(y)\varphi(y)ds(y)+\partial_{s}\int_{\Gamma}G_{k}(x-y)\rho(y)ds(y)$
$\displaystyle=$ $\displaystyle\frac{1}{ik}f(x)$
$\displaystyle\int_{\Gamma}G_{k}(x-y)\partial_{s}\varphi(y)ds(y)-ik\int_{\Gamma}G_{k}(x-y)\rho(y)ds(y)$
$\displaystyle=$ $\displaystyle 0.$

Performing an integration by parts on the integral term above that contains
the quantity $\partial_{s}\varphi$ (this is justified even in the case when
$\Gamma$ is an open arc since $\varphi$ vanishes at the open ends) we obtain
the following Current and Charge (CC) BIE formulation of the Helmholtz problem
(30) whose unknowns are the densities $(\varphi,\rho)\in H^{1/2}(\Gamma)\times
H^{-1/2}(\Gamma)$ in the case when $\Gamma$ is closed and respectively
$(\varphi,\rho)\in\widetilde{H}^{1/2}(\Gamma)\times\widetilde{H}^{-1/2}(\Gamma)$
when $\Gamma$ is an open arc

$\displaystyle-
ik\int_{\Gamma}G_{k}(x-y)\hat{n}(x)\cdot\hat{n}(y)\varphi(y)ds(y)+\partial_{s}\int_{\Gamma}G_{k}(x-y)\rho(y)ds(y)$
$\displaystyle=$ $\displaystyle\frac{1}{ik}f(x)$
$\displaystyle\nabla\cdot\int_{\Gamma}G_{k}(x-y)\hat{t}(y)\varphi(y)ds(y)-ik\int_{\Gamma}G_{k}(x-y)\rho(y)ds(y)$
$\displaystyle=$ $\displaystyle 0.$ (39)

We note that the CC formulation (3.2) above features BIOs whose singularities
are no longer integrable (i.e they exhibit Cauchy p.v. operators), and as
such, they cannot be discretized directly using Alpert quadratures. On the
other hand, the QBX methods can handle relatively seamlessly these new
operators. In the QBX approach, the derivatives of the layer potentials that
feature in the CC formulation (3.2) are performed by differentiating term by
term the expansions around the centers, via the formulas

$\partial_{s}SL(k)[\rho](x)\approx-\sum_{\ell=-p}^{p}\ell\alpha_{\ell}[\rho]\frac{J_{\ell}(k|x^{+}-x|)}{|x^{+}-x|}e^{-i\ell\theta^{+}}$

and respectively

$\displaystyle\nabla^{x}\cdot SL(k)[\hat{t}\varphi](x)$ $\displaystyle\approx$
$\displaystyle-\sum_{\ell=-p}^{p}k\hat{\alpha}_{\ell}[\varphi]\cdot\hat{n}(x)\left(-J_{\ell+1}(k|x^{+}-x|+\frac{\ell}{k|x^{+}-x|}J_{\ell}(k|x^{+}-x|)\right)e^{-i\ell\theta^{+}}$
$\displaystyle-$ $\displaystyle
i\sum_{\ell=-p}^{p}\ell\hat{\alpha}_{\ell}[\varphi]\cdot\hat{t}(x)\frac{J_{\ell}(k|x^{+}-x|)}{|x^{+}-x|}e^{-i\ell\theta^{+}}$

where

$\hat{\alpha}_{\ell}[\varphi]:=\frac{i}{4}\int_{\Gamma}H_{\ell}^{(1)}(k|x^{\prime}-x^{+}|)e^{i\ell\theta^{\prime}}\hat{t}(x^{\prime})\varphi(x^{\prime})ds(x^{\prime}).$

Clearly, we can express the QBX expansion coefficients
$\alpha_{\ell}[\rho]=\sum_{m}\alpha_{\ell,m}[\rho]$ and
$\hat{\alpha}_{\ell}[\varphi]=\sum_{m}\hat{\alpha}_{\ell,m}[\varphi]$ by
accounting for the panel contributions. Just like in the Dirichlet case, the
incorporation of weighted unknowns in the CC formulation is beneficial with
regards to both the condition numbers and the spectral properties vis a vis
iterative solvers of the QBX Nyström matrices. Indeed, instead of using the
more singular functional density $\rho$, we favor using its weighted version
of defined on each panel $\Gamma_{m}$ as

$\rho_{m}^{w}(t):=\rho(\gamma_{m}(t))|\gamma_{m}^{\prime}(t)|\sqrt{1-t^{2}},\quad-1\leq
t\leq 1,\ 1\leq m\leq M,$ (40)

together with the Crenshaw-Curtis quadratures (29) for the evaluation of the
expansion coefficients $\alpha_{\ell,m}[\rho^{w}]$. Given that the functional
density $\varphi$ is more regular than $\rho$, we do not choose to use a
weighted version of it in the case when $\Gamma$ is closed, and thus we use
the Crenshaw-Curtis quadratures (27) to evaluate the coefficients
$\hat{\alpha}_{\ell,m}[\varphi]$.

In the case when $\Gamma$ is an open arc, we can incorporate the explicit
asymptotic behavior of the more regular density $\varphi$ in the neighborhood
of the end points into our numerical scheme. We start with the simpler case
when $\Gamma$ is a smooth open arc with the usual
$\gamma:[-1,1]\to\mathbb{R}^{2}$, and we employ regularized unknowns that take
into account the singular behavior of the functional densities
$(\varphi,\rho)$ in the BIE (3.2) as proposed in [13]. Specifically, we use
the asymptotic behavior $\varphi\sim\mathcal{O}(\sqrt{d})$ and
$\rho\sim\mathcal{O}(1/\sqrt{d})$ where $d$ denotes the distance to the open
ends, to define weighted densities that are globally regular

$\varphi(\gamma(t))|\gamma^{\prime}(t)|=\varphi^{w}(t)\sqrt{1-t^{2}},\
\rho^{w}(t):=\rho(\gamma(t))|\gamma^{\prime}(t)|\sqrt{1-t^{2}},\quad-1\leq
t\leq 1.$ (41)

However, the weighted approach (41) has to be slightly modified if we
represent the open arc $\Gamma$ via panels, that is
$\Gamma=\bigcup_{m=1}^{M}\Gamma_{m}$. In this case we make a distinction
between the treatment of densities defined on panels that contain the end
points of $\Gamma$—we assume without loss of generality that those are
$\Gamma_{1}$ and respectively $\Gamma_{M}$, and densities defined on panels
that do not contain the open ends of $\Gamma$. The latter case is
straightforward, as we assume $[-1,1]$ parametrizations of those panels
equipped with Chebyshev meshes, and we use the weighted version
$\rho_{m}^{w},2\leq m\leq M-1$ defined in equation (40), and the unweighted
density $\varphi$. In the former case, since the densities defined on the
panels $\Gamma_{1}$ and $\Gamma_{M}$ exhibit the asymptotic behavior
$\varphi\sim\mathcal{O}(\sqrt{d})$ and $\rho\sim\mathcal{O}(1/\sqrt{d})$ where
$d$ denotes the distance to the open ends of $\Gamma$ (and implicitly of
$\Gamma_{1}$ and $\Gamma_{M}$), we consider in turn $[0,1]$ parametrizations
of the open end panels $\Gamma_{1}$ and $\Gamma_{M}$ such that the open ends
correspond in the parameter space to $t=0$. We then use on $\Gamma_{1}$ and
$\Gamma_{M}$ Chebyshev quadrature points $v_{n}$ defined as

$v_{n}:=\frac{1}{2}(1+\cos(\vartheta_{n})),\quad\vartheta_{n}:=\frac{(2n-1)\pi}{2N},\quad
n=1,\ldots,N,$

weighted densities $\varphi_{m}^{w}$ and $\rho_{m}^{w}$ for $m\in\\{1,M\\}$
that incorporate explicitly the asymptotic behavior of the densities
$(\varphi,\rho)$ around the open end points, and the following modified Fejér
quadratures to evaluate the QBX expansion coefficients

$\int_{0}^{1}\frac{f(v)}{\sqrt{v}}dv\approx\sum_{n=1}^{N}{}^{{}^{\prime}}\omega_{n}^{(1)}f(v_{n}),\
\omega_{n}^{(1)}:=\frac{2}{N}\sum_{q=1}^{N}\frac{4}{1-4(q-1)^{2}}\cos((q-1)\vartheta_{n}),\quad
n=1,\ldots,N.$ (42)

and respectively

$\int_{0}^{1}f(v)\sqrt{v}\
dv\approx\sum_{n=1}^{N}{}^{{}^{\prime}}\omega_{n}^{(2)}f(v_{n}),\
\omega_{n}^{(2)}:=\frac{2}{N}\sum_{q=1}^{N}\left(\frac{1}{1-4(q-1)^{2}}+\frac{3}{9-4(q-1)^{2}}\right)\cos((q-1)\vartheta_{n}).$
(43)

where the prime notation denotes that the first term is halved. We note that
the panel approach presented above can also handle the case when the open arc
$\Gamma$ is piecewise smooth.

We present in Table 8 numerical results concerning QBX discretizations of the
CC formulations (3.2) in the case when the wavenumber $k$ is real; reference
solutions were obtained via (a) the high-order kernel-splitting Nyström
methods in [2, 13] in the case of closed Lipschitz curves and smooth arcs and
(b) refined QBX discretizations of the CC formulations with the modified Fejér
quadratures (42) and (43). Similar orders of convergence were observed in the
case when the wavenumbers are complex. We used one panel $M=1$ for the first
three scatterers and two panels $M=2$ for the V-shaped arc, a fixed QBX
expansion parameter $p=4$ and up sampling parameter $\beta=4$ and we present
errors in the near field together with estimated orders of convergence as the
size of the Chebyshev meshes is refined. We have observed in practice that
using larger numbers of panels (of smaller size) leads to QBX matrices with
larger condition numbers that also require larger GMRES iterations counts for
the solution of ensuing linear systems. These findings can be attributed to
the fact that the CC formulation is a formulation of the first kind; the
condition numbers and GMRES iteration counts do not vary with the size of the
panels on the case of the QBX Nyström matrices corresponding to the second
kind formulation (31). Nevertheless, the CC formulations (3.2) are
advantageous because they are universally applicable to all types of
boundaries $\Gamma$, including piecewise smooth closed boundaries as well as
open arcs. In addition, we mention that analytical preconditioners can be
employed in order to speed up the iterative convergence of the discrete linear
systems corresponding to Nyström discretizations of the first kind CC
formulations. Indeed, assuming that the Nyström CC system is expressed in the
block matrix form

$\begin{bmatrix}\mathcal{C}\mathcal{C}_{11}^{n}&\mathcal{C}\mathcal{C}_{12}^{n}\\\
\mathcal{C}\mathcal{C}_{21}^{n}&\mathcal{C}\mathcal{C}_{22}^{n}\end{bmatrix}\begin{bmatrix}\varphi_{n}\\\
\rho_{n}\end{bmatrix}=\frac{1}{ik}\begin{bmatrix}f_{n}\\\ 0\end{bmatrix}$

we use the following preconditioned CC formulation

$\begin{bmatrix}\mathcal{C}\mathcal{C}_{22}^{n}&-\mathcal{C}\mathcal{C}_{12}^{n}\\\
-\mathcal{C}\mathcal{C}_{21}^{n}&\mathcal{C}\mathcal{C}_{11}^{n}\end{bmatrix}\begin{bmatrix}\mathcal{C}\mathcal{C}_{11}^{n}&\mathcal{C}\mathcal{C}_{12}^{n}\\\
\mathcal{C}\mathcal{C}_{21}^{n}&\mathcal{C}\mathcal{C}_{22}^{n}\end{bmatrix}\begin{bmatrix}\varphi_{n}\\\
\rho_{n}\end{bmatrix}=\frac{1}{ik}\begin{bmatrix}\mathcal{C}\mathcal{C}_{22}^{n}&-\mathcal{C}\mathcal{C}_{12}^{n}\\\
-\mathcal{C}\mathcal{C}_{21}^{n}&\mathcal{C}\mathcal{C}_{11}^{n}\end{bmatrix}\begin{bmatrix}f_{n}\\\
0\end{bmatrix}.$ (44)

The preconditioned CC formulation (44) can be justified heuristically on
account of the pseudodifferential calculus, whose details we leave out for
now. We have observed in practice that this preconditioned formulations
appears to be of the second kind.

$n_{1}$ | Teardrop $M=1$ | Boomerang $M=1$ | Strip $M=1$ | V-shaped $M=2$
---|---|---|---|---
| $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c. | $\varepsilon_{\infty}$ | e.o.c.
32 | 2.1 $\times$ $10^{-1}$ |  | 8.4 $\times$ $10^{-2}$ |  | 1.2 $\times$ $10^{-2}$ |  | 8.5 $\times$ $10^{-3}$ |
64 | 2.6 $\times$ $10^{-2}$ | 3.03 | 7.6 $\times$ $10^{-3}$ | 3.46 | 3.4 $\times$ $10^{-3}$ | 1.89 | 2.2 $\times$ $10^{-3}$ | 1.95
128 | 2.5 $\times$ $10^{-3}$ | 3.32 | 5.1 $\times$ $10^{-4}$ | 3.88 | 9.1 $\times$ $10^{-4}$ | 1.89 | 5.4 $\times$ $10^{-4}$ | 2.01
256 | 2.4 $\times$ $10^{-4}$ | 3.42 | 6.1 $\times$ $10^{-5}$ | 3.08 | 2.4 $\times$ $10^{-4}$ | 1.89 | 1.3 $\times$ $10^{-4}$ | 1.99
Table 8: Errors in the near field and estimated orders of convergence obtained
when the QBX Nyström discretization is applied to the CC formulation (3.2)
with real wavenumber $k=8$ and normal plane wave incidence. We used the QBX
parameters $p=4$ and $\beta=4$.

Having described various discretization strategies of BIE formulations of
Helmholtz equations that are oblivious to the nature of the wavenumber (be it
real or complex), we present in the final section numerical results concerning
CQ discretizations of the wave equation.

##  4 CQ Numerical Experiments

We start with an example in Table 9 that showcases the advantages that can be
garnered from using Nyström discretizations based on Alpert and QBX
quadratures when dealing with layer potentials featuring complex wavenumbers.
Specifically, we considered the ensemble of wavenumbers
$s_{\ell}=\frac{P(\lambda\zeta^{-\ell}_{N+1})}{\Delta t}$ where
$P(\zeta)=\frac{1}{2}(\zeta^{2}-4\zeta+3)$ for $0\leq\ell\leq N$ with
$\lambda=\max\\{(\Delta t)^{3/N},{\rm{eps}}^{1/2N}\\}$, $N=1024$ and $\Delta
t=\frac{2}{N}$ associated with the CQ BDF2 method (here $\rm{eps}$ is the
machine precision). In the case when $\Gamma$ is the unit circle, the single
layer operators are diagonalizable, that is

$V(s_{\ell})e^{imx}=\frac{i\pi}{2}J_{|m|}(s_{\ell})H_{m}^{(1)}(s_{\ell})e^{imx},\quad
m\in\mathbb{Z}.$ (45)

We used the relation (45) in the case $m=3$ to test the accuracy of various
Nyström discretizations of the single layer BIO $V(s_{\ell}),\ 0\leq\ell\leq
N$. Specifically, we considered equispaced discretizations $x^{M}$ of size $M$
of the unit circle and various Nyström discretizations $V^{M}(s_{\ell})$ end
we report the errors on the boundary

$\varepsilon_{M}=\max_{0\leq\ell\leq N}\max_{x\in
x^{M}}\left|V^{M}(s_{\ell})e^{imx}-\frac{i\pi}{2}J_{|m|}(s_{\ell})H_{m}^{(1)}(s_{\ell})e^{imx}\right|$

as the size $M$ of the discretization is increased. We investigated the
behavior of two classes of Nyström discretization: (i) the Alpert and QBX
quadratures that are applied directly to the kernel of the single layer BIO;
and (ii) the kernel splitting techniques that involve cutoff localizations
(19) and modifications of the Bessel function according to formula (20) (with
a truncation parameter $L=4$ in the latter). We considered discretizations
that can deliver accuracies at the level of ${\rm eps}^{\frac{1}{2}}$ which is
the level achievable by CQ quadratures [7]. As it can be seen from the results
in Table 9, the Nyström discretizations based on Alpert and QBX quadratures
respectively achieve uniformly small errors which are orders of magnitude more
accurate than those produced by Nyström discretizations based on kernel
splitting techniques.

$M$ | Alpert $a=6,\ m=10$ | QBX $M=1,\ p=12,\ \beta=4$ | Cutoff splitting (19) |  $\widetilde{J}_{0}$ (20)
---|---|---|---|---
256 | 7.4 $\times$ $10^{-6}$ | 7.1 $\times$ $10^{-7}$ | 3.4 $\times$ $10^{-3}$ | 3.6 $\times$ $10^{-3}$
512 | 2.7 $\times$ $10^{-7}$ | 1.1 $\times$ $10^{-7}$ | 7.9 $\times$ $10^{-4}$ | 8.4 $\times$ $10^{-4}$
1024 |  |  | 1.7 $\times$ $10^{-4}$ | 8.6 $\times$ $10^{-5}$
2048 |  |  | 3.1 $\times$ $10^{-6}$ | 1.9 $\times$ $10^{-6}$
Table 9: Boundary errors in the eigenvalue test for circular geometries and a
CQ BDF ensemble of complex wavenumbers.

For all the time domain scattering experiments we consider the incident plane
wave

$g(x,t)=\cos(5t-x\cdot(1,0))\exp(-1.5(5t-x\cdot(1,0)-5)^{2}),$

and we present results concerning the convergence of CQ numerical schemes that
evaluate the scattered fields at $T=2$ at $512$ points placed on circles
located at distance $1$ from the scatterers. The incident field considered is
smooth enough (in the causality sense, that is the field and its higher order
derivatives vanish on the scatterers for negative times all the way to time
$t=0$) so that the predicted theoretical rates of convergence rates of CQ
methods [5, 6] are validated in our numerical experiments. The smoothness of
the incident fields also manifests itself in the sparse frequency content of
the boundary data in the Laplace domain [7], e.g. most boundary data defined
in equations (12) and respectively (14) are small in the infinity norm. As a
result, we follow the common thresholding practice in CQ and we solve only
those Laplace domain Helmholtz equations whose boundary data $\hat{g}_{\ell}$
in equations (12) is larger than $10^{-10}$ when measured in the infinity
norm. In all the time domain CQ experiments based on QBX quadratures for the
discretization of the Laplace domain modified Helmholtz problems we considered
representations of the scatterers in terms of geometrically large Chebyshev
panels, whose number depended only on the number of corners or open ends
present on the scatterers’ boundaries. We show results in Figure 1–3
convergence plots of the BDF2 and RK3/RK5 CQ methods in the case of Dirichlet
boundary conditions for Lipschitz closed scatterers as well as open scatterers
(strips). The Laplace domain problems that enter the CQ methods were solved
using the weighted single layer formulation whose associated BIOs were
discretized using both Alpert quadratures and QBX (based on the weighted
unknown defined in equation (28) and the Fejér quadratures (29)) of
appropriately high order. As it can be seen from the results presented in
Figure 1–3 the expected second, third, and respectively fifth orders of
convergence were observed for the BDF2, RK3 and respectively RK5 CQ methods.
We mention that qualitatively similar results can be observed for smooth
scatterers; we simply chose to report time domain results for non-smooth
scatterers in this paper as the behavior of CQ solvers for such problems was
not discussed in great detail in the delta-bem literature [33].

Figure 1: Orders of convergence of the near field CQ solutions of the wave
equation in the exterior of the teardrop domain at time $T=2$ with Dirichlet
boundary conditions. We used Alpert quadratures with the same parameters
$\sigma=4$, $a=2$ and $m=3$ for the solution of the ensemble of weighted
single layer BIE formulations (23) in the Laplace domain using $N=256$
discretization points on the boundary. The reference solution was produced
using the CQ-RK5 discretization with $2048$ time steps and $N=512$ boundary
discretization points.

Figure 2: Orders of convergence of the near field CQ solutions of the wave
equation in the exterior of the boomerang domain at time $T=2$ with Dirichlet
boundary conditions. We used Alpert quadratures with the same parameters
$\sigma=4$, $a=2$ and $m=3$ for the solution of the ensemble of weighted
single layer BIE formulations (23) in the Laplace domain using $N=256$
discretization points on the boundary. The reference solution was produced
using the CQ-RK5 discretization with $2048$ time steps and $N=512$ boundary
discretization points.

Figure 3: Orders of convergence of the near field CQ solutions of the wave
equation in the exterior of the strip domain at time $T=2$ with Dirichlet
boundary conditions. We used (i) Alpert quadratures with the same parameters
$\sigma=4$, $a=2$ and $m=3$ and (ii) Chebyshev based QBX discretizations (29)
with parameters $M=1$ (that is one global Chebyshev mesh) $p=8$ and $\beta=4$
for the solution of the ensemble of BIE formulations (23) in the Laplace
domain using $N=256$ discretization points on the boundary for both
quadratures. The reference solution was produced using the CQ-RK5
discretization with $2048$ time steps and $N=512$ boundary discretization
points.

We continue in Figure 4 with an illustration of the high orders of convergence
achieved by various CQ methods considered in this paper—up to and including
fifth order, in the case of a time domain scattering problem off a closed non-
smooth obstacle with Neumann boundary conditions. We considered BIE solvers
based on both the weighted single layer and the double layer formulations of
the CQ Laplace domain problems with Neumann boundary conditions, leading to
BIE of the second kind (32) in the former case and to the hypersingular BIE
alternative CC formulation (3.2) in the latter case. The ensemble of CQ second
kind formulations (32) was discretized using both Alpert quadratures and QBX
methods of high enough orders, whereas the CC formulations (32) were
discretized using QBX. We continue in Figure 5 with an illustration of the
high-order convergence achieved by CQ solvers in the case of a time domain
scattering problem off a strip with Neumann boundary conditions. In this case
only the hypersingular BIE formulation (34) is available for the solution of
Laplace domain problems, and we used both Alpert quadratures and QBX for its
discretization—the former in conjunction with cosine changes of variables that
resolve the end point singularities and the latter in the context of the CC
reformulation (3.2) of (34). Finally, we present in Figure 6 and Figure 7 time
domain scattering experiments involving non-smooth open arcs with Neumann
boundary conditions: one (top panel) and two V-shaped (bottom panel) obstacles
in the former figure and a W-shaped obstacle in the latter figure. We used the
modified Fejér quadratures (42) and (43) in the QBX discretizations for the
solutions of the frequency domain CQ modified Helmholtz problems. Since the CC
formulation (3.2) is of the first kind, we investigated the efficacy of the
simple preconditioner (44) for the Laplace domain problems corresponding to
the BDF2 frequencies. As it can be seen from the plots in Figure 7 and Figure
8 of GMRES iteration numbers required to reach $10^{-7}$ residuals for the
unpreconditioned and preconditioned CC formulations, the preconditioner (44)
is highly effective.

Figure 4: Orders of convergence of the near field CQ solutions of the wave
equation in the exterior of the teardrop domain at time $T=2$ with Neumann
boundary conditions. We used (i) Alpert quadratures with the same parameters
$\sigma=4$, $a=2$ and $m=3$ for the solution of the ensemble of second kind
BIE formulations (32) in the Laplace domain using $N=256$ discretization
points on the boundary and (ii) Chebyshev based QBX discretizations with
parameters $M=1,\ n_{1}=256$, $p=8$, and $\beta=4$ for the solution of the
ensemble of CC BIE formulations (3.2). The reference solution was produced
using the CQ-RK5 discretization with $2048$ time steps and $N=512$ boundary
discretization points.

Figure 5: Orders of convergence of the near field CQ solutions of the wave
equation in the exterior of the strip at time $T=2$ with Neumann boundary
conditions. We used (i) Alpert quadratures with the same parameters
$\sigma=4$, $a=2$ and $m=3$ for the solution of the ensemble of the
hypersingular BIE formulations (34)in the Laplace domain using $N=256$
discretization points on the boundary and(ii) Chebyshev based QBX
discretizations with parameters $M=1,\ n_{1}=256$, $p=8$ and $\beta=4$ for the
solution of the ensemble of CC BIE formulations (3.2). The reference solution
was produced using the CQ-RK5 discretization with $2048$ time steps and
$N=512$ boundary discretization points.

Figure 6: Orders of convergence of the near field CQ solutions of the wave
equation in the exterior of one V-shaped scatterer (top panel) and two
V-shaped scatterers arranged in a chevron pattern (bottom panel) at time $T=2$
with Neumann boundary conditions. We used the Chebyshev based QBX
discretizations with parameters $M=2,\ n_{1}=n_{2}=128$, $p=8$ and $\beta=4$
that incorporated the modified Fejér quadratures (42) and (43) on each
V-shaped scatterer for the solution of the ensemble of CC BIE formulations
(3.2). The reference solution was produced using the CQ-RK5 discretization
with $2048$ time steps and $N=512$ boundary discretization points per
scatterer.

Figure 7: Left panel: Orders of convergence of the near field CQ solutions of
the wave equation in the exterior of a W-shaped scatterer at time $T=2$ with
Neumann boundary conditions. We used the Chebyshev based QBX discretizations
with parameters $M=4,\ n_{j}=128,1\leq j\leq 4$, $p=8$ and $\beta=4$ that
incorporated the modified Fejér quadratures (42) and (43) for the solution of
the ensemble of CC BIE formulations (3.2). The reference solution was produced
using the CQ-RK5 discretization with $2048$ time steps and $N=1024$ Chebyshev
points on the W-shaped boundary. Right panel: GMRES iterations required by the
CC formulation and its preconditioned version to reach residuals of $10^{-7}$
for the Laplace domain modified Helmholtz problems corresponding to BDF2
frequencies for $T=2$ and $1024$ time steps; we considered only those
frequencies whose Laplace domain BDF2 boundary data has a norm larger than
$10^{-10}$.

Figure 8: Numbers of GMRES iterations to reach $10^{-7}$ residuals of the
Nyström CC and preconditioned CC (44) systems for the V-shaped scatterer and
the whole ensemble of Laplace domain frequencies associated with the CQ BDF2
formulation with $N=1024$ time steps and final time $T=2$. We considered only
those frequencies whose Laplace domain BDF2 boundary data has a norm larger
than $10^{-10}$

Finally, we present in Figures 9 and 10 time domain multiple scattering
experiments involving a slotted cylinder, that is a cavity like structure
comprised of a circular scatterer enclosed by an open arc of small aperture,
with both Dirichlet and Neumann boundary conditions over long periods of time
(e,g, final time is $T=16$). The numerical simulations presented in these
Figures were produced with CQ-RK3 methods using $2048$ time steps (that is
$\Delta t=7.8\times 10^{-3}$) and $512$ discretization points per boundary for
the discretization of the BIE formulation of the modified Helmholtz equations
in the Laplace domain, parameters which result in accuracies at the level of
$10^{-4}$ in the near field.

Figure 9: Time sequences starting left to right, top to bottom, of the total
fields scattered by a multiple scattering configuration consisting of a
circular scatterer enclosed by an open arc, with Dirichlet boundary conditions
on both scatterers.

Figure 10: Time sequences starting left to right, top to bottom, of the total
fields scattered by a multiple scattering configuration consisting of a
circular scatterer enclosed by an open arc, with mixed Dirichlet and Neumann
boundary conditions on the scatterers.

##  5 Conclusions

We have shown that CQ methods used in conjunction with Nyström discretizations
for the BIE solution of their associated Laplace domain Helmholtz problems can
deliver high-order numerical solutions of two dimensional wave equations in
unbounded domains with non smooth open/closed boundaries. In particular, the
CQ-RK5 methods require the largest time steps to reach desired accuracy
levels, and more importantly the gains in accuracy over the CQ-BDF2 and CQ-RK3
are not outweighed by the increased number of Laplace domain problems that
ought to be solved in this CQ version. We focused on two classes of high-order
Nyström discretizations, one based on Alpert quadratures, and one on QBX
methods. Our choice was motivated by the fact that both of these quadratures
are rather straightforward to implement, are applicable seamlessly to real as
well as complex wavenumbers, are compatible with fast algorithms, and can be
extended to three dimensional settings. The extension of the CQ Nyström
methodology presented in this paper to the solution of three dimensional wave
equations and the acceleration of the solvers via fast algorithms are
currently underway.

## Acknowledgments

Catalin Turc gratefully acknowledges support from NSF through contract
DMS-1908602.

## References

  * [1] Bradley K Alpert.  Hybrid gauss-trapezoidal quadrature rules.  SIAM Journal on Scientific Computing, 20(5):1551–1584, 1999.
  * [2] A. Anand, J. S. Ovall, and C. Turc.  Well-conditioned boundary integral equations for two-dimensional sound-hard scattering problems in domains with corners.  J. Integral Equations Appl., 24(3):321–358, 2012.
  * [3] Akash Anand, Jeff Ovall, Catalin Turc, et al.  Well-conditioned boundary integral equations for two-dimensional sound-hard scattering problems in domains with corners.  Journal of Integral Equations and Applications, 24(3):321–358, 2012\.
  * [4] Thomas G Anderson, Oscar P Bruno, and Mark Lyon.  High-order, dispersionless fast-hybrid wave equation solver. part i: O (1) sampling cost via incident-field windowing and recentering.  SIAM Journal on Scientific Computing, 42(2):A1348–A1379, 2020.
  * [5] Lehel Banjai.  Multistep and multistage convolution quadrature for the wave equation: algorithms and experiments.  SIAM Journal on Scientific Computing, 32(5):2964–2994, 2010.
  * [6] Lehel Banjai, Christian Lubich, and Jens Markus Melenk.  Runge–kutta convolution quadrature for operators arising in wave propagation.  Numerische Mathematik, 119(1):1–20, 2011.
  * [7] Lehel Banjai and Stefan Sauter.  Rapid solution of the wave equation in unbounded domains.  SIAM Journal on Numerical Analysis, 47(1):227–249, 2009.
  * [8] Alex Barnett, Leslie Greengard, and Thomas Hagstrom.  High-order discretization of a stable time-domain integral equation for 3d acoustic scattering.  Journal of Computational Physics, 402:109047, 2020.
  * [9] Abderrahmane Bendali, Francis Collino, Mbarek Fares, and Bassam Steif.  Extension to nonconforming meshes of the combined current and charge integral equation.  IEEE transactions on antennas and propagation, 60(10):4732–4744, 2012.
  * [10] Timo Betcke, Nicolas Salles, and Wojciech Smigaj.  Overresolving in the laplace domain for convolution quadrature methods.  SIAM Journal on Scientific Computing, 39(1):A188–A213, 2017.
  * [11] Y. Boubendir and C. Turc.  Wave-number estimates for regularized combined field boundary integral operators in acoustic scattering problems with neumann boundary conditions.  IMA Journal of Numerical Analysis, 33(4):1176–1225, 2013.
  * [12] James Bremer.  On the nyström discretization of integral equations on planar curves with corners.  Applied and Computational Harmonic Analysis, 32(1):45–64, 2012\.
  * [13] Oscar P Bruno and Stéphane K Lintner.  Second-kind integral solvers for te and tm problems of diffraction by open arcs.  Radio Science, 47(06):1–13, 2012.
  * [14] Roman Chapko, Rainer Kress, and Lars Mönch.  On the numerical solution of a hypersingular integral equation for elastic scattering from a planar crack.  IMA journal of numerical analysis, 20(4):601–619, 2000.
  * [15] D. Colton and R. Kress.  Inverse acoustic and electromagnetic scattering theory, volume 93 of Applied Mathematical Sciences.  Springer-Verlag, Berlin, second edition, 1998.
  * [16] Victor Dominguez, Sijiang L Lu, and Francisco-Javier Sayas.  A fully discrete calderón calculus for two dimensional time harmonic waves.  arXiv preprint arXiv:1210.7017, 2012.
  * [17] Víctor Domínguez, Sijiang L Lu, and Francisco-Javier Sayas.  A nyström flavored calderón calculus of order three for two dimensional waves, time-harmonic and transient.  Computers & Mathematics with Applications, 67(1):217–236, 2014\.
  * [18] Victor Dominguez, Mark Lyon, and Catalin Turc.  Well-posed boundary integral equation formulations and nystr$\backslash$” om discretizations for the solution of helmholtz transmission problems in two-dimensional lipschitz domains.  Journal of Integral Equations and Applications, 28(3):395–440, 2016\.
  * [19] Víctor Domínguez, M-L Rapún, and F-J Sayas.  Dirac delta methods for helmholtz transmission problems.  Advances in Computational Mathematics, 28(2):119–139, 2008.
  * [20] Víctor Domínguez, Tonatiuh Sánchez-Vizuet, and Francisco-Javier Sayas.  A fully discrete calderón calculus for the two-dimensional elastic wave equation.  Computers & Mathematics with Applications, 69(7):620–635, 2015\.
  * [21] Charles L Epstein and Leslie Greengard.  Debye sources and the numerical solution of the time harmonic maxwell equations.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 63(4):413–463, 2010.
  * [22] A Bamberger et T. Ha Duong and JC Nedelec.  Formulation variationnelle espace-temps pour le calcul par potentiel retardé de la diffraction d’une onde acoustique (i).  Mathematical methods in the applied sciences, 8(1):405–435, 1986\.
  * [23] Sijia Hao, Alex H Barnett, Per-Gunnar Martinsson, and P Young.  High-order accurate methods for nyström discretization of integral equations on smooth curves in the plane.  Advances in Computational Mathematics, 40(1):245–272, 2014.
  * [24] Johan Helsing and Anders Karlsson.  An extended charge-current formulation of the electromagnetic transmission problem.  SIAM Journal on Applied Mathematics, 80(2):951–976, 2020.
  * [25] Johan Helsing and Rikard Ojala.  Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning.  Journal of Computational Physics, 227(20):8820–8840, 2008.
  * [26] Andreas Klöckner, Alexander Barnett, Leslie Greengard, and Michael O’Neil.  Quadrature by expansion: A new method for the evaluation of layer potentials.  Journal of Computational Physics, 252:332–349, 2013.
  * [27] R. Kress.  A Nyström method for boundary integral equations in domains with corners.  Numer. Math., 58(2):145–161, 1990.
  * [28] R. Kress.  On the numerical solution of a hypersingular integral equation in scattering theory.  J. Comput. Appl. Math., 61(3):345–360, 1995.
  * [29] R. Kussmaul.  Ein numerisches Verfahren zur Lösung des Neumannschen Neumannschen Aussenraumproblems für die Helmholtzsche Schwingungsgleichung.  Computing (Arch. Elektron. Rechnen), 4:246–273, 1969.
  * [30] Ignacio Labarca, Luiz M Faria, and Carlos Pérez-Arancibia.  Convolution quadrature methods for time-domain scattering from unbounded penetrable interfaces.  Proceedings of the Royal Society A, 475(2227):20190029, 2019.
  * [31] Wangtao Lu and Ya Yan Lu.  Efficient high order waveguide mode solvers based on boundary integral equations.  Journal of Computational Physics, 272:507–525, 2014.
  * [32] E. Martensen.  Über eine Methode zum räumlichen Neumannschen Problem mit einer Anwendung für torusartige Berandungen.  Acta Math., 109:75–135, 1963.
  * [33] Francisco-Javier Sayas.  Retarded potentials and time domain boundary integral equations.  2013\.
  * [34] Michael Siegel and Anna-Karin Tornberg.  A local target specific quadrature by expansion method for evaluation of layer potentials in 3d.  Journal of Computational Physics, 364:365–392, 2018.
  * [35] Matti Taskinen and Pasi Yla-Oijala.  Current and charge integral equation formulation.  IEEE Transactions on Antennas and Propagation, 54(1):58–67, 2006\.
  * [36] Ali E Yilmaz, Jian-Ming Jin, and Eric Michielssen.  Time domain adaptive integral method for surface integral equations.  IEEE Transactions on Antennas and Propagation, 52(10):2692–2708, 2004.