New preconditioners for Laplace and Helmholtz integral equations on open curves: Analytical framework and Numerical results

TL;DR

This paper introduces two novel preconditioners for first-kind integral equations in 2D Helmholtz scattering problems, improving the conditioning of discretized systems through analytical and numerical analysis.

Contribution

The paper develops and analyzes new preconditioners based on square-roots of local operators for open curves, extending existing analytical preconditioners to Lipschitz scatterers.

Findings

Preconditioners significantly improve system conditioning.

Numerical examples demonstrate enhanced computational efficiency.

The analytical framework supports the design of effective preconditioners.

Abstract

The Helmholtz wave scattering problem by screens in 2D can be recast into first-kind integral equations which lead to ill-conditioned linear systems after discretization. We introduce two new preconditioners, in the form of square-roots of local operators respectively for the corresponding problems with Dirichlet and Neumann conditions on the arc. They generalize the so-called "analytical" preconditioners available for Lipschitz scatterers. We introduce a functional setting adapted to the singularity of the problem and enabling the analysis of those preconditioners. The efficiency of the method is demonstrated on several numerical examples.

Figures

Tables

Table 1. Table 1 : Computing time and number of GMRES iterations for the numerical resolution of the weighted Laplace single-layer integral equation on the segment respectively with the square root preconditioner and without preconditioner.

	with Prec.		without Prec.
$N$	$n_{it}$	t(s)	$n_{it}$	t(s)
500	8		79
2000	8		128	0.3
8000	7	$0.5$	218	11.5
32000	8	$2.2$	347	89

Table 2. Table 2 : Computing time and number of GMRES iterations for the numerical resolution of the weighted Laplace hypersingular integral equation on the segment respectively with the square root preconditioner and without preconditioner.

	with Prec.		without Prec.
$N$	$n_{it}$	t(s)	$n_{it}$	t(s)
500	5		333	0.3
2000	5		$>500$	2
8000	7	1.1	$>500$	60
32000	6	4	$>500$	725

Table 3. Table 3 : Computing time and number of GMRES iterations for the numerical resolution of the Helmholtz weighted single-layer integral equation on the segment respectively with the square root preconditioner and without preconditioner.

	with Prec.		without Prec.
$k\left|\Gamma \right|$	$n_{it}$	t(s)	$n_{it}$	t(s)
50$\pi$	8		88	0.2
200$\pi$	10	0.2	123	1.4
400$\pi$	13	4	145	40
800$\pi$	16	15	155	110
1600$\pi$	20	70	199	642

Table 4. Table 4 : Computing time and number of GMRES iterations for the numerical resolution of the Helmholtz weighted hypersingular integral equation on the segment respectively with the square root preconditioner and without preconditioner.

	with Prec.		without Prec.
$k\left|\Gamma \right|$	$n_{it}$	t(s)	$n_{it}$	t(s)
50$\pi$	10		$>500$	3.1
200$\pi$	13	0.3	$>500$	12
400$\pi$	14	7	$>500$	260
800$\pi$	18	34	-	$>15$min
1600$\pi$	25	270	-	$>15$min

Table 5. Table 5 : Computing time and number of GMRES iterations for the numerical resolution of the Helmholtz weighted single-layer integral equation on the spiral-shaped arc respectively with the square-root preconditioner and without preconditioner

	With prec.		Without prec.
$k\left|\Gamma \right|$	$n_{it}$	t(s)	$n_{it}$	t(s)
$50\pi$	19	0.1	93	0.2
$200\pi$	24	0.55	136	1.77
$400\pi$	27	4.7	160	30
$800\pi$	30	16.5	190	92
$1600\pi$	32	71	217	456

Table 6. Table 6 : Computing time and number of GMRES iterations for the numerical resolution of the Helmholtz weighted hypersingular integral equation on the spiral-shaped arc respectively with the square-root preconditioner and without preconditioner

	With prec.		Without prec.
$k\left|\Gamma \right|$	$n_{it}$	t(s)	$n_{it}$	t(s)
$50\pi$	22	0.15	$>500$	3
$200\pi$	31	0.7	$>500$	15.4
$400\pi$	34	13	$>500$	209
$800\pi$	35	46	-	$>15$min
$1600\pi$	42	283	-	$>15$min

Table 7. Table 7 : Number of GMRES iterations for the numerical resolution of the Helmholtz weighted single-layer integral equation with the square-root preconditioner, from left to right, one the flat segment, and on V-shaped arcs with increasing singularities ( θ = 3 ​ π 4 , π 2 , π 6 𝜃 3 𝜋 4 𝜋 2 𝜋 6 \theta=\frac{3\pi}{4},\frac{\pi}{2},\frac{\pi}{6} ). In all cases, the parameeter k ​ | Γ | 𝑘 Γ k\left\lvert\Gamma\right\rvert is fixed to 50 50 50 and the mesh is progressively refined. In this test, we use N p = 60 subscript 𝑁 𝑝 60 N_{p}=60 Padé approximants and u D = u i ​ n ​ c | Γ subscript 𝑢 𝐷 subscript 𝑢 conditional 𝑖 𝑛 𝑐 Γ u_{D}=u_{inc|\Gamma} where u i ​ n ​ c ​ ( x , y ) = e i ​ k ​ y subscript 𝑢 𝑖 𝑛 𝑐 𝑥 𝑦 superscript 𝑒 𝑖 𝑘 𝑦 u_{inc}(x,y)=e^{iky} .

$N/\left(k\left|\Gamma \right|\right)$	flat seg.	$\theta =\frac{3\pi }{4}$	$\theta =\frac{\pi }{2}$	$\theta =\frac{\pi }{6}$
2.5	8	9	10	17
5	7	8	9	17
7.5	7	8	10	17
10	7	8	10	17
12.5	7	8	9	17
15	7	8	10	17

Table 8. Table 8 : Computing time and number of GMRES iterations for the numerical resolution of the Helmholtz weighted single-layer integral equation on the V-shaped arc, with θ = π 2 𝜃 𝜋 2 \theta=\frac{\pi}{2} respectively with the square root preconditioner and without preconditioner.

	With prec.		Without prec.
$k\left|\Gamma \right|$	$n_{it}$	$t\left(s\right)$	$n_{it}$	$t\left(s\right)$
50$\pi$	9		97	$0.25$
200$\pi$	10	0.3	157	3.1
400$\pi$	11	3.1	190	41
800$\pi$	14	8	231	138
1600$\pi$	18	48	-	$>15$min

Table 9. Table 9 : Number of GMRES iterations and H − 1 / 2 superscript 𝐻 1 2 H^{-1/2} error for the resolution of the singe- layer integral equation S k ​ λ = u D subscript 𝑆 𝑘 𝜆 subscript 𝑢 𝐷 S_{k}\lambda=u_{D} on the flat segment, with k = 10 ​ π 𝑘 10 𝜋 k=10\pi , N = 80 𝑁 80 N=80 . In the first five lines, the equation is solved using a standard Galerkin method on a graded mesh with parameter β 𝛽 \beta ranging from 1 (uniform mesh) to 5. In all those 5 cases, we precondition the linear system using a discrete version of P k ′ subscript superscript 𝑃 ′ 𝑘 P^{\prime}_{k} eq. ( 19 ). The last line refers to the weighted Galerkin method described in this work, with the new square root preconditioner. Even for this small problem, the Galerkin matrices associated to graded meshes become very close to singular, and the preconditioners are no longer accurately evaluated using the Padé approximation as in section 5. The square root is instead computed directly using an eigenvalue decomposition of the matrices. The H − 1 / 2 superscript 𝐻 1 2 H^{-1/2} error is obtained by computing a reference solution on a mesh which is refined until the error estimates stabilize.

Numer. method	$n_{it}$	relative $H^{-1/2}$ error
Unif. mesh ($\beta =1$)	10	0.088
Graded $\beta =2$	12	0.020
Graded $\beta =3$	13	0.0066
Graded $\beta =4$	17	0.0036
Graded $\beta =5$	21	0.0030
Weighted Galerk.	7	2.2e-5

Table 10. Table 10 : Computing time and number of GMRES iterations for the Helmholtz weighted single-layer integral equation on the flat segment, respectively with the Calderón preconditioner and with our new square-root preconditioners.

	Calderón Prec.		Square root Prec.
$k\left|\Gamma \right|$	$n_{it}$	t(s)	$n_{it}$	t(s)
50$\pi$	15		8
200$\pi$	15	0.45	10	0.35
400$\pi$	15	11	13	5
800$\pi$	15	42	16	18

Table 11. Table 11 : Computing time and number of GMRES iterations for the Helmholtz weighted hypersingular integral equation on the flat segment, respectively with the Calderón preconditioner and with our new square-root preconditioners.

	Calderón Prec.		Square root Prec.
$k\left|\Gamma \right|$	$n_{it}$	t(s)	$n_{it}$	t(s)
50$\pi$	15		10
200$\pi$	16	0.3	13	0.3
400$\pi$	17	18	15	7
800$\pi$	17	68	18	34

Equations

$$G_k:=S_k{\tilde{\Lambda }}_k+\frac{I_d}{2}+D_k$$

$${\tilde{\Lambda }}_k:=\sqrt{-{\Delta }_{\Gamma }-k^2{I}_d},$$

$$S_k\lambda =u_D{N}_k\mu =u_N,$$

$$P_k=\frac{1}{\omega }\sqrt{-(\omega {\partial }_{\tau }{)}^2-k^2{\omega }^2},Q_k=\omega {\left[-({\partial }_{\tau }\omega {)}^2-k^2{\omega }^2\right]}^{-\frac{1}{2}},$$

$$\Gamma =\left\{\gamma (t)|t\in [-1,1]\right\}$$

$$G_k(x)=\{\begin{array}{ll}-\frac{1}{2\pi }\ln |x|& \text{if }k=0,\\ \frac{i}{4}{H}_0^{(1)}(k|x|)& \text{if }k>0,\end{array}$$

$$\forall x\in \Gamma ,S_k\varphi (x)={\int }_{\Gamma }{G}_k(x-y)\varphi (y)d\sigma (y)$$

$$⟨N_k\varphi ,\psi ⟩:={\iint }_{\Gamma \times \Gamma }{G}_k(x-y)\left({\varphi }^{\prime }(x){\psi }^{\prime }(y)-k^2\varphi (x)\psi (y)\right)d\sigma (x)d\sigma (y),$$

$$S_k:{\tilde{H}}^{-\frac{1}{2}}(\Gamma )\to H^{1/2}(\Gamma )$$

$$N_k:{\tilde{H}}^{\frac{1}{2}}(\Gamma )\to H^{-\frac{1}{2}}(\Gamma )$$

$$S_k\lambda =u_D,N_k\mu =u_N.$$

$$\lambda =\frac{\alpha }{\omega },\mu =\omega \beta$$

$$\omega (x)=\sqrt{d(x,\partial \Gamma )}.$$

$$S_0\lambda =u_D\text{ and }{N}_0\mu =u_N,$$

$$\forall x\in [-1,1]T_n(x)=\cos (n\arccos (x)),U_n(x)=\frac{\sin ((n+1)\arccos (x))}{\sqrt{1-x^2}},$$

$$(1-x^2){\partial }_{xx}{T}_n-x{\partial }_x{T}_n+n^2{T}_n=0$$

$$(1-x^2){\partial }_{xx}{U}_n-3x{\partial }_x{U}_n+n(n+2)U_n=0$$

$$(\omega {\partial }_x{)}^2{T}_n$$

$$({\partial }_x\omega {)}^2{U}_n$$

$$L_{\frac{1}{\omega }}^2:=\left\{u\in L_{\text{loc}}^1(-1,1)|{\int }_{-1}^1\frac{f^2(x)}{\sqrt{1-x^2}}dx<+\infty \right\}$$

$$L_{\omega }^2:=\left\{u\in L_{\text{loc}}^1(-1,1)|{\int }_{-1}^1{f}^2(x)\sqrt{1-x^2}dx<+\infty \right\},$$

$$\mathcal{C}f(\theta ):=f(\cos \theta ),\mathcal{S}g(\theta ):=\sin (\theta )g(\cos (\theta ))$$

$$\mathcal{C}{T}_n(\theta )=\cos (n\theta )\text{and}\mathcal{S}{U}_n(\theta )=\sin ((n+1)\theta )$$

$$u(x)=\sum _{n=0}^{+\infty }{\hat{u}}_n{T}_n(x)$$

$${\left(u,v\right)}_{\frac{1}{\omega }}:=\frac{1}{\pi }{\int }_{-1}^1\frac{u(x)\overline{v(x)}}{\omega (x)}dx.$$

$$v(x)=\sum _{n=0}^{+\infty }{\check{v}}_n{U}_n(x)$$

$${\left(u,v\right)}_{\omega }={\int }_{-1}^{1}u(x)v(x)\omega (x)dx.$$

$$T^s=\left\{u\in L_{\frac{1}{\omega }}^2|\sum _{n=0}^{+\infty }(1+n^2{)}^s{\left|{\hat{u}}_n\right|}^2<+\infty \right\}.$$

$${\left(u,v\right)}_{T^s}={\hat{u}}_0\overline{{\hat{v}}_0}+\frac{1}{2}\sum _{n=1}^{+\infty }(1+n^2{)}^s{\hat{u}}_n\overline{{\hat{v}}_n},$$

$$U^s=\left\{u\in L_{\omega }^2|\sum _{n=0}^{+\infty }(1+n^2{)}^s{\left|{\check{u}}_n\right|}^2\right\}$$

$${\left(u,v\right)}_{U^s}=\frac{1}{2}\sum _{n=0}^{+\infty }(1+(n+1{)}^2{)}^s{\check{u}}_n\overline{{\check{v}}_n}.$$

Full text

Introduction

For the numerical resolution of wave scattering problems, a by-now well-established approach is the boundary elements method, which involves the discretization of integral equations leading to linear systems of smaller size when compared to finite element methods, but the corresponding matrices are fully populated. They are solved efficiently by combining iterative methods such as GMRES [38], with acceleration methods (e.g. the fast multipole method, see [19] and references therein). The number of iterations in GMRES depends on the condition number of the linear system. The question of preconditioning, by which we mean minimizing this number of iterations, either by formulating well-conditioned integral equations or by finding preconditioners for the linear systems, has attracted a lot of attention since more than two decades.

In the case of an obstacle with a $C^{\infty}$ smooth boundary $\Gamma$, general methods to build efficient preconditioners are known [1, 2, 3, 11, 23, 41], which in most cases can also be applied in theory and/or in practice if the scatterer is only assumed to have Lipschitz regularity. Among them, one efficient approach is the Generalized Combined Source Integral Equation (GCSIE) method [2] initiated by Levadoux in [28]. It involves the inversion of the discretized GCSIE operator

[TABLE]

where $S_{k}$ and $D_{k}$ stand for the classical single- and double-layer potentials on $\Gamma$, $I_{d}$ is the identity operator and $\tilde{\Lambda}_{k}$ is an approximation of the exterior Dirichlet-to-Neumann (DtN) map $\Lambda_{k}$ for the Helmholtz equation. In the ideal case ${\tilde{\Lambda}_{k}=\Lambda_{k}}$, $G_{k}$ becomes the identity and, consequently, when $\tilde{\Lambda}_{k}$ is a compact perturbation of $\Lambda_{k}$, $G_{k}$ is a compact perturbation of the identity. It is well-known that the Galerkin discretization of such operators leads to well-conditioned linear systems for which the GMRES method converges supra-linearly [8, 22, 41]. To build suitable candidates for $\tilde{\Lambda}_{k}$, a generic tool is the theory of pseudo-differential operators (see e.g. [25]). This has been succesfully applied by Antoine and Darbas in [3] who proposed the choice

[TABLE]

based on low-order expansions of the symbol of $\Lambda_{k}$. The authors also introduce an efficient scheme to approximate the square root operator in (1) relying on Padé approximants. Their numerical results clearly demonstrate that the preconditioning performances are independent of the discretization parameters and robust in $k$.

When the scattering object is a “screen”, that is a curve in 2D and a surface in $3D$ with edges (therefore not a Lipschitz domain), the robustness of the precondtioners in the discretization parameters and in $k$ is lost in practice for all of the aforementioned methods. In this work, we propose a generalized version of the square-root preconditioners of Antoine and Darbas on screens in 2D that overcomes this limitation. To fix the notation, let us write the classical first-kind integral equations corresponding to the Dirichlet and Neumann problems as

[TABLE]

where $S_{k}$ and $N_{k}$ are the classical single- and hypersingular layer potentials on the screen $\Gamma$ and $u_{D}$ and $u_{N}$ are smooth right-hand sides. We define the operators

[TABLE]

where $\partial_{\tau}$ is the tangential derivative on $\Gamma$ and $\omega(x)$ is defined for $x\in\Gamma$ as the square root of the distance from $x$ to the edges of $\Gamma$. We present some numerical evidence showing that the preconditioning performances of $P_{k}$ (resp. $Q_{k}$) for $S_{k}$ (resp. $N_{k}$) are independent of the discretization parameters, and only very lightly sensitive to the wavenumber $k$. The careful analysis of $P_{k}$ and $Q_{k}$, and the proof of the fact that they provide parametrices of $S_{k}$ and $N_{k}$ respectively relies on two new classes of pseudo-differential operators on open curves, and is treated in full details in [5].

To our knowledge, the recent literature mainly contains two other approaches for the problem of preconditioning for Helmholtz screen scattering problems, and the present work is connected to both of them. The first one, see e.g. [23, 24, 26, 35, 44], builds on recently found closed-form expressions for the inverses of the Laplace layer potentials ($k=0$) on flat screens in 2D and 3D (recovered by a unified approach in [18]). Using compact perturbations arguments, those operators are expected to provide efficient preconditioners for the Helmholtz ($k>0$) layer potentials on arbitrary smooth screens when the wavenumber $k$ is reasonably small. Here, we show that when $\Gamma$ is a flat segment, $P_{0}$ and $Q_{0}$ coincide with those exact inverses (up to minor modifications). Our numerical results show that using $P_{k}$ and $Q_{k}$ instead of $P_{0}$ and $Q_{0}$ for the Helmholtz problem improves greatly the preconditioning performances.

The second one, developed by Bruno and Lintner [10] is a generalization of the Calderón preconditioners available for Lipschitz scatterers. In their work, some weighted versions of the Helmholtz layer potentials, $S_{k,\omega}$ and $N_{k,\omega}$, are considered and it is shown that $S_{k,\omega}N_{k,\omega}$ is a second-kind operator, therefore, $S_{k,\omega}$ gives a good preconditioner for $N_{k,\omega}$ and reciprocally. The weighted layer potentials $S_{k,\omega}$ and $N_{k,\omega}$ also appear in the present study as the fundamental objects on which the arguments are developed. Furthermore, we compare numerically our preconditioning method to that of Bruno and Lintner. In our implementation, we find that their method seems slightly more robust in $k$. However, our preconditioners are significantly faster to evaluate due to their quasi-local form.

The outline of the paper is as follows. We use the first section to fix notation and recall some important results concerning integral equations on screens. In the second section, we introduce the preconditioners for the Laplace problem on a flat screen. The formulas are generalized to the Helmholtz in the third section. In section 4, we describe the weighted Galerkin setup that we use to discretize the integral equations. In section 5, we describe how we build the preconditioners operators introduced in sections 2 and 3. Finally, in section 6, we show the performance of our preconditioners in a variety of cases and compare them to other ideas of the literature.

1 First kind integral equations

In the following, we consider a smooth, non-intersecting open arc $\Gamma$, that is, a set of the form

[TABLE]

where $\gamma:[-1,1]\to\mathbb{R}^{2}$ is an injective and $C^{\infty}$ function. The arc $\Gamma$ models a $1$-dimensional screen in a 2D setting. Let $k\geq 0$ be the wavenumber and $G_{k}$ the Green kernel defined by

[TABLE]

where $H_{0}^{(1)}$ is the Hankel function of first kind and [math] order. The classical single-layer potential, denoted by $S_{k}$, is defined for $\phi$ is the space $C^{\infty}_{c}(\Gamma)$ of compactly supported functions smooth functions on $\Gamma$ by

[TABLE]

where $d\sigma$ is the uniform measure on $\Gamma$. Moreover, the Helmholtz hypersingular layer potential $N_{k}:C^{\infty}_{c}(\Gamma)\to\mathcal{D}^{\prime}(\Gamma)$ is defined by

[TABLE]

where $\phi^{\prime}$ and $\psi^{\prime}$ denote the arclength derivatives of $\phi$ and $\psi$ and $\mathcal{D}^{\prime}(\Gamma)$ is the set of distributions on $\Gamma$.

Let us now recall the definitions of the spaces $H^{s}(\Gamma)$ and $\tilde{H}^{s}(\Gamma)$ for $s\in\mathbb{R}$, following [30, chap. 3]. For this purpose, we consider a closed curve $\tilde{\Gamma}$ such that $\Gamma\subset\tilde{\Gamma}$. A distribution $u$ on $\Gamma$ is said to be in $H^{s}(\Gamma)$ if there exists an $U\in H^{s}(\tilde{\Gamma})$ such that $u=U_{|\Gamma}$. Furthermore, $u\in\tilde{H}^{s}(\Gamma)$ if $u$ is in ${H}^{s}(\tilde{\Gamma})$ and $\textup{supp}\,u\subset\Gamma$.

Proposition 1** (see [42, Thm 1.8] and [43, Thm 1.4]).**

The operator $S_{k}$ has a continuous extension

[TABLE]

which is bicontinuous when $k\neq 0$. When $k=0$, it is bicontinuous if and only if the logarithmic capacity of $\Gamma$ (see e.g. [16]) is not equal to $1$. Similarly, the operator $N_{k}$ can be extended continuously as

[TABLE]

and this extension is bicontinuous.

It will be assumed in the following that the logarithmic capacity of $\Gamma$ is not $1$. This can always be satisfied by properly rescaling the problem. In this work, we are concerned with the resolution of the integral equations

[TABLE]

where $u_{D}\in H^{1/2}(\Gamma)$ and $u_{N}\in H^{-1/2}(\Gamma)$. They are related to the Helmholtz problem in $\mathbb{R}^{2}\setminus\Gamma$ with prescribed Dirichlet boundary data $u_{D}$ (“sound-soft”) or Neumann boundary data $u_{N}$ (“sound-hard”) respectively.

Due to the singularity of the manifold $\Gamma$, the solutions $\lambda$ and $\mu$ of the previous equations have edge singularities even if $u_{D}$ and $u_{N}$ are arbitrarily smooth. More precisely, Costabel et al. have shown

Proposition 2** ([14, Cor. A.5.1]).**

Assume that $u_{D}$ and $u_{N}$ are in $C^{\infty}(\Gamma)$. Then the solutions $\lambda$ and $\mu$ of eq. (2) can be expressed as

[TABLE]

where $\alpha$ and $\beta$ are in $C^{\infty}(\Gamma)$ and where

[TABLE]

In addition to the singular nature of the problem, the Galerkin discretization of first-kind integral equations is known to produce ill-conditioned linear systems, see [39, Sec. 4.5]. The aim of this paper is to introduce a formalism that enables to resolve the singularity and provide efficient preconditioners for the linear systems.

2 Laplace equation on a flat segment

We start with the particular case where $\Gamma$ is the flat segment $\Gamma=[-1,1]\times\{0\}$ and furthermore $k=0$. The associated integral equations,

[TABLE]

with logarithmic kernels, are the object of a considerable number of papers in the 1990’s, for instance [4, 17, 27, 32, 40, 45].

More recently, Jerez-Hanckes and Nédélec [26] have exhibited exact inverses of those operators through explicit variational forms (Prop 3.1 and 3.3) and identities involving tangential square-roots of differential operators with boundary conditions at the end points (Prop 3.10 and subsequent remark). In Gimperlein et al. [18], those results are recovered with a different method, and it is shown that $S_{0}$ and $N_{0}$ are related fractional powers of the Laplace-Beltrami operator on the line $\left\{(x,0)\mathrel{}\middle|\mathrel{}x\in\mathbb{R}\right\}$, by considering the natural extensions of functions by [math] outside $\Gamma$.

Here we give a new expression of the inverses of $S_{0}$ and $N_{0}$ (Theorem 1 and Theorem 2), in the form of square roots of suitably weighted tangential differential operators. Those formulas do not involve boundary conditions at the end points nor extensions by [math] outside $\Gamma$. The proofs of the identities are very simple, nevertheless, it seems, to the best of our knowledge, that they have remained unnoticed so far. Their generalization to the case $k>0$ (section 3) is the main object of this work.

2.1 Analytical setting

We introduce the Chebyshev polynomials of first and second kinds [29], respectively given by

[TABLE]

and we denote by $\omega$ the operator $\omega:u(x)\mapsto\omega(x)u(x)$ with $\omega(x)=\sqrt{1-x^{2}}$. We also denote by $\partial_{x}$ the derivation operator. The Chebyshev polynomials satisfy the ordinary differential equations

[TABLE]

and

[TABLE]

which can be rewritten in divergence form as

[TABLE]

where we emphasize that $\partial_{x}\omega f$ should be understood as the composition of the operators $\partial_{x}$ and $\omega$ applied to $f$, that is $\partial_{x}\omega f(x)=\frac{d}{dx}\left(\omega(x)f(x)\right)$. Both $T_{n}$ and $U_{n}$ are polynomials of degree $n$, and form complete orthogonal families respectively of the Hilbert spaces

[TABLE]

and

[TABLE]

see e.g. [29] and in particular Thm. 5.2. This fact can also be verified directly by remarking that for $(f,g)\in L^{2}_{\frac{1}{\omega}}\times L^{2}_{\omega}$, the functions

[TABLE]

are respectively even and odd functions in $L^{2}(-\pi,\pi)$. This can be verified simply using the change of variables $x=\cos\theta$. The completeness of the polynomials $T_{n}$ and $U_{n}$ in $L^{2}_{\frac{1}{\omega}}$ and $L^{2}_{\omega}$ are deduced from the completeness of

[TABLE]

in the space of $2\pi$-periodic $L^{2}$ even or odd functions respectively.

As a consequence, any $u\in L^{2}_{\frac{1}{\omega}}$ can be decomposed through the first kind Chebyshev series

[TABLE]

where the Fourier-Chebyshev coefficients $\hat{u}_{n}$ are given by $\hat{u}_{n}\mathrel{\mathop{:}}=\dfrac{\left(u_{n},T_{n}\right)_{\frac{1}{\omega}}}{\left(T_{n},T_{n}\right)_{\frac{1}{\omega}}}\,.$ with the inner product

[TABLE]

Similarly, any function $v\in L^{2}_{\omega}$ can be decomposed along the $(U_{n})_{n\geq 0}$ as

[TABLE]

where the coefficients $\check{v}_{n}$ are given by $\check{v}_{n}\mathrel{\mathop{:}}=\dfrac{\left(v,U_{n}\right)_{\omega}}{\left(U_{n},U_{n}\right)_{\omega}}$ with the inner product

[TABLE]

Those properties can be used to define Sobolev-like spaces.

Definition 1**.**

For all $s\geq 0$, we define

[TABLE]

Endowed with the scalar product

[TABLE]

$T^{s}$* is a Hilbert space for all $s\geq 0$. Similarly, we set*

[TABLE]

which is a Hilbert space for the scalar product

[TABLE]

One can extend the definitions of $T^{s}$ ad $U^{s}$ for $s\in\mathbb{R}$, in which case they form interpolating scales of Hilbert space, see [5, Section 1] for details. Letting $H^{s}_{e}\oplus H^{s}_{o}$ be the partition of the Sobolev space of $2\pi$-periodic functions $H^{s}_{\textup{per}}$ into even and odd functions, it is easy to check that the pullbacks

[TABLE]

are isomorphisms. It follows that the inclusions

[TABLE]

are compact. Denoting by $T^{\infty}=\cap_{s\geq 0}T^{s}$ and similarly for $U^{\infty}$, it is shown in [5, Lem. 7] that

Lemma 1**.**

[TABLE]

For $s=\pm\frac{1}{2}$, those spaces have been analyzed (with different notation) e.g. in [26] and verify

[TABLE]

[TABLE]

2.2 Laplace single-layer equation

We start with the single-layer integral equation $S_{0}\lambda=g$, with $\lambda\in\tilde{H}^{-1/2}(\Gamma)$, that is

[TABLE]

The following result is the fundamental tool for studying this equation:

Lemma 2** (See e.g. [29, Thm 9.2]).**

For all $n\in\mathbb{N}$, we have

[TABLE]

where

[TABLE]

Using the decomposition of $g$ on the basis $(T_{n})_{n}$, we see at once that the solution $\lambda$ to equation (8) admits the expansion

[TABLE]

As a corollary, we obtain by an alternative proof the result of Costabel et al. cited in Proposition 2 in this particular case:

Corollary 1**.**

If the data $g$ is in $C^{\infty}([-1,1])$, the solution $\lambda$ to the equation

[TABLE]

is of the form

[TABLE]

with $\alpha\in C^{\infty}([-1,1])$.

Proof.

Let $\alpha(x)=\sqrt{1-x^{2}}\,\lambda(x)$ where $\lambda$ is the solution of $S_{0}\lambda=g$. By 1, if $g\in C^{\infty}([-1,1])$, then $g\in T^{\infty}$, and by equation (9),

[TABLE]

from which we deduce that $\alpha$ also belongs to $T^{\infty}=C^{\infty}([-1,1])$. ∎

Following [10], we introduce the weighted single layer operator as the operator that appears in 2.

Definition 2**.**

Let $S_{0,\omega}$ be the weighted single layer operator defined by

[TABLE]

Lemma 2 can be restated as saying that the Chebyshev polynomials $T_{n}$ are eigenfunctions of $S_{0,\omega}$. To obtain the solution of the single-layer( integral equation (8), we thus solve the weighted single-layer integral equation

[TABLE]

and let $\lambda=\frac{\alpha}{\omega}$, which indeed belongs to $\tilde{H}^{-1/2}(\Gamma)$ by eq. (6).

Comparing the eigenvalues of $S_{0,\omega}$ and $-(\omega\partial_{x})^{2}$, we directly obtain the following result:

Theorem 1**.**

There holds

[TABLE]

or equivalently,

[TABLE]

where $\pi_{0}$ is the $L^{2}_{\frac{1}{\omega}}$ orthogonal projection on $T_{0}$, that is

[TABLE]

The operator $-(\omega\partial_{x})^{2}$ being self-adjoint (with domain $T^{1}$), its square-root can be defined via the spectral theorem, see for example [21, Def. 10.5]. In later occurences, the square-root will also be applied to self-adjoint operators whose spectrum may contain negative values. In this case, we take the standard definition of the square root, with $\sqrt{-r}=ir$ for all $r>0$.

Proof.

We simply notice that for all $n\in\mathbb{N}$, $-(\omega\partial_{x})^{2}T_{n}=n^{2}T_{n}$, thus

[TABLE]

Therefore,

[TABLE]

∎

Since the operator $2\sqrt{-(\omega\partial_{x})^{2}}+\frac{2}{\ln(2)}\pi_{0}$ is the inverse of $S_{0,\omega}$, it can be used as an efficient preconditioner for the weighted integral equation (10).

2.3 Laplace hypersingular equation

We now turn our attention to the Laplace hypersingular equation

[TABLE]

Similarly to the previous section and again following [10], we consider the weighted operator $N_{0,\omega}\mathrel{\mathop{:}}=N_{0}\omega$. We can get the solution to equation (11) by solving the weighted hypersingular integral equation

[TABLE]

and letting $\mu=\omega\beta$. We now show that $N_{0,\omega}$ can be analyzed using this time the spaces $U^{s}$.

Lemma 3**.**

For any $\beta$, $\beta^{\prime}\in U^{\infty}$, one has

[TABLE]

Proof.

We use the well-known integration by parts formula

[TABLE]

valid when $u$ and $v$ are regular enough and vanish at the endpoints of the segment. For a smooth $\beta$, we thus have

[TABLE]

which implies the claimed identity. ∎

Lemma 4**.**

For all $n\in\mathbb{N}$, there holds

[TABLE]

Proof.

From the identity $\partial_{x}T_{n+1}=(n+1)U_{n}$ and eq. (4), we obtain

[TABLE]

Therefore, by 3

[TABLE]

∎

The identity $-(\partial_{x}\omega)^{2}U_{n}=(n+1)^{2}U_{n}$ now leads to the following result:

Theorem 2**.**

[TABLE]

or equivalently,

[TABLE]

Here again the inverse square root is defined by functional calculus, which is possible since $-(\partial_{x}\omega)^{2}$ is self-adjoint (with domain $U^{1}$).

3 Helmholtz equation

In this section, we aim at generalizing the preceding analysis to the case of Helmholtz equation ($k>0$) on $\mathbb{R}^{2}\setminus\Gamma$ with $\Gamma=[-1,1]\times\{0\}$, based on the explicit formulas presented in the previous section. In an analogous fashion as above, let $S_{k,\omega}\mathrel{\mathop{:}}=S_{k}\frac{1}{\omega}$ and $N_{k,\omega}\mathrel{\mathop{:}}=N_{k}\omega$. The following commutation holds:

Theorem 3**.**

[TABLE]

Proof.

Let us compute

[TABLE]

where we use the notation $\omega_{y}$ and $\partial_{y}$ to emphasize the dependence in the variable $y$. Since $(\omega\partial_{x})^{2}$ is symmetric with respect to the bilinear form $\left(\cdot,\cdot\right)_{\frac{1}{\omega}}$, we have

[TABLE]

Thus,

[TABLE]

where $D_{k}(x,y)\mathrel{\mathop{:}}=\left[(\omega_{y}\partial_{y})^{2}-(\omega_{x}\partial_{x})^{2}\right]\left[G_{k}(x-y)\right]$. A simple computation leads to

[TABLE]

Since $G_{k}$ is a solution of the Helmholtz equation, we have for all $(x\neq y)\in\mathbb{R}$

[TABLE]

thus

[TABLE]

A careful analysis shows that no Dirac mass appears in the previous formula. Note that $y^{2}-x^{2}=\omega_{x}^{2}-\omega_{y}^{2}$ so the first term vanishes and we find

[TABLE]

as claimed. ∎

There also holds the following identity:

[TABLE]

The proof can be found in [7, Chap. 2, Thm 2.2].

Those commutations imply that the operators $S_{k,\omega}$ and $N_{k,\omega}$ share the same eigenvectors as, respectively, $\left[-(\omega\partial_{x})^{2}-k^{2}\omega^{2}\right]$ and $\left[-(\partial_{x}\omega)^{2}-k^{2}\omega^{2}\right]$. The eigenfunctions of the operator $\left[-(\omega\partial_{x})^{2}-k^{2}\omega^{2}\right]$ thus provide us with a diagonal basis for $S_{k,\omega}$. They are the solutions of another Sturm-Liouville problem

[TABLE]

Once we set $x=\cos\theta$, $\tilde{y}(\theta)=y(x)$, $q=\frac{k^{2}}{4}$, $a=\lambda+2q$, $\tilde{y}$ is a solution of the standard Mathieu equation

[TABLE]

There exists a discrete set of values $a_{2n}(q)$ for which this equation possesses even and $2\pi$ periodic solutions, which are known as the Mathieu cosine functions, and usually denoted by $\textup{ce}_{n}$. Here, we use the notation $\textup{ce}^{k}_{n}$ to emphasize the dependency in the parameter $k=\sqrt{2q}$ of those functions. The normalization is taken as

[TABLE]

The Mathieu cosine functions are $L^{2}$ orthogonal:

[TABLE]

so that any even $2\pi$ periodic function in $L^{2}(-\pi,\pi)$ can be expanded along the functions $\textup{ce}_{n}$, with the coefficients obtained by orthonormal projection. Setting

[TABLE]

in analogy to the zero-frequency case, we have

[TABLE]

For large $n$, using the general results from the theory of Hill’s equations (see e.g. [33, eqs. (21), (28) and (29)]), we have the following asymptotic formula for $\lambda_{n,k}$:

[TABLE]

The first commutation established in Theorem 3 implies that the Mathieu cosine functions are also the eigenfunctions of the single-layer operator. (An equivalent statement is given in [9, Thm 4.2], if we allow the degenerate case $\mu=0$.)

A similar analysis can be applied to the hypersingular operator. The eigenfunctions of $\left[-(\partial_{x}\omega)^{2}-k^{2}\omega^{2}\right]$ are given by

[TABLE]

where $\textup{se}_{n}^{k}$ are the so-called Mathieu sine functions, which also satisfy the Mathieu differential equation (13), but with the condition that they are $2\pi$ periodic and odd functions.

One could furthermore expect that the operators

[TABLE]

provide compact perturbations of the inverses of $S_{k,\omega}$ and $N_{k,\omega}$ respectively, and the knowledge of the eigenvalues of $S_{k,\omega}$ and $N_{k,\omega}$ would allow for a simple proof of this fact. Those eigenvalues are not known but the result can be obtained by a more involved analysis. To state precisely the result, let us introduce the following terminology:

Definition 3**.**

A linear operator $A:T^{-\infty}\to T^{-\infty}$ (resp. $U^{-\infty}\to U^{-\infty}$) is of order $\alpha$ in the scale $(T^{s})_{s}$ (resp. $(U^{s})_{s}$) if for all real $s$, $A$ maps continuously $T^{s}$ to $T^{s-\alpha}$ (resp. $U^{s}$ to $U^{s-\alpha}$).

Theorem 4** ([5, Thms. 4 and 6]).**

The operators

[TABLE]

are well defined and satisfy

[TABLE]

respectively in $T^{-\infty}$ and $U^{-\infty}$, where $I_{d}$ is the identity operator and $K_{1}$ and $K_{2}$ are of order $-4$ and $-3$ respectively in the scales $T^{s}$ and $U^{s}$.

Since $K_{1}$ and $K_{2}$ are of negative order, they are in particular compact endomorphisms of $L^{2}_{\frac{1}{\omega}}$ and $L^{2}_{\omega}$ respectively, due to the compact inclusions that hold for

[TABLE]

Therefore, the operators appearing in (14) are of second-kind in those spaces. The extent to which the dependence in $k$ in $P_{k}$ and $Q_{k}$ is optimal is reflected by the next theorem.

Theorem 5** ([5, Cor. 10 ]).**

Let $K$ be an operator of order [math] in the scale $T^{s}$ and let $\Delta_{K}$ be defined by

[TABLE]

Then $\Delta_{k}$ is of order $-2$ in the scale $T^{s}$. Moreover, $\Delta_{K}$ is of order $-4$ if and only if

[TABLE]

where $L$ is an operator of order $-2$ in the scale $T^{s}$.

An analogous result holds for the hypersingular operator, see [5, Cor. 12] Taking $K=0$, we see that $P_{0}$ and $Q_{0}$ are also compact equivalent inverses of $S_{k,\omega}$ and $N_{k,\omega}$, however up to a less regularizing remainder than $P_{k}$ and $Q_{k}$. A clear link between the order of the remainder and the numerical performance of the preconditioner remains to be elucidated. We conjecture that a smoother remainder leads to better performances and especially, more robustness with respect to the parameter $k$. This is strongly supported by our numerical results exposed in section 6.

We note that all the previous (except the commutations) carries over to the more general case of a $C^{\infty}$ non-intersecting open curve $\Gamma$. In this case, we define a weight $\omega_{\Gamma}$ on the curve by replacing $\partial_{x}$ by $\omega_{\Gamma}(r(t))=\frac{\left\lvert\Gamma\right\rvert}{2}\omega(t)$, where $\left\lvert\Gamma\right\rvert$ is the length of the curve and $r:[-1,1]\to\Gamma$ is such that for all $x$, $\left\lvert\partial_{x}r(x)\right\rvert=\frac{\left\lvert\Gamma\right\rvert}{2}$. Then, the derivation operator $\partial_{x}$ is replaced by $\partial_{\tau}$ the tangential derivative on $\Gamma$, and the weighted layer potentials by $S_{k,\omega_{\Gamma}}\mathrel{\mathop{:}}=S_{k}\frac{1}{\omega_{\Gamma}}$, and $N_{k,\omega_{\Gamma}}\mathrel{\mathop{:}}=N_{k}\omega_{\Gamma}$,

The previous theoretical analysis suggests to use $P_{k}$ and $Q_{k}$ as operator preconditioners for $S_{k,\omega}$ and $N_{k,\omega}$. The remainder of this paper is dedicated to testing this idea in practice.

4 Galerkin method

It is known that the naive piecewise polynomial Galerkin discretization of the non-weighted integral equations with a uniform mesh converges very slowly in terms of the mesh size (the error in energy norm converges in $O(\sqrt{h})$, see e.g. [34, Thm 1.2]). Several alternative discretization schemes for the integral equations on open arcs have been proposed in the literature, including mesh grading [34, 13], Galerkin method with special singular functions [42, 12], cosine change of variables [45], and high-order Nyström methods [10].

Here, we describe a simple Galerkin setting, suited to the spaces $T^{s}$ and $U^{s}$. We use standard piecewise linear functions defined on a non-uniform mesh, which is refined towards the edges as follows. Let $X:[0,\left\lvert\Gamma\right\rvert]$ be the parametrization of $\Gamma$ by the arclength, where $\left\lvert\Gamma\right\rvert$ is the length of the curve. We choose the breakpoints $(X_{i})_{1\leq i\leq N}$ as

[TABLE]

where $(s_{i})_{1\leq i\leq N}$ are such that the value

[TABLE]

be (approximately) constant in $i$.

Such a mesh turns out to be analogous to an algebraically graded mesh with a grading parameter $\beta=2$. That is to say, near an edge of the curve, the width of the $i-th$ interval is approximately $(ih)^{2}$. Notice that this modification alone, i.e. using the h-BEM method with a polynomial order $p=1$, is not sufficient to get an optimal rate of convergence. Indeed, it only leads to a convergence rate in $O(h)$ for the energy norm (cf. [34, Theorem 1.3]) instead of the expected $O(h^{5/2})$ behavior (to reach such an order of convergence would require $\beta=5$).

The key ingredient to recover optimal convergence, beside the graded mesh, is to use a weighted $L^{2}$ scalar product (with weight $\frac{1}{\omega}$ or $\omega$ depending on the considered equation), in order to assemble the operators in their natural spaces. We state here the orders of convergence that one gets with this new method, and refer the reader to [7, Chap. 3, Sec. 2.3] for the proofs. To keep the exposition simple, we also restrict our presentation to the case where $\Gamma=[-1,1]\times\{0\}$ and $k=0$.

Dirichlet problem.

For the resolution of the single-layer equation (8) we use a variational formulation of (10) to compute an approximation $\alpha_{h}$ of $\alpha$. Namely, let $V_{h}$ the Galerkin space of (discontinuous) piecewise affine functions defined on the mesh $(x_{i})_{0\leq i\leq N}$ defined above, and $\alpha_{h}$ the unique solution in $V_{h}$ to

[TABLE]

We then compute $\lambda_{h}=\frac{\alpha_{h}}{\omega}$. Using the notation $C$ to denote any constant that does not depend on the parameter $h$, we then have

Theorem 6** ([7, Thm. 3.1]).**

If the data $u_{D}$ is in $T^{s+1}$ for some $-1/2\leq s\leq 2$, then there holds:

[TABLE]

In particular, when $u_{D}$ is smooth, the solution $\alpha=\omega\lambda$ belongs to $T^{\infty}$, and we get the optimal rate of convergence of the error in $O(h^{5/2})$.

Neumann problem.

For the numerical resolution of (11), we use a variational form for equation (12) to compute an approximation $\beta_{h}$ of $\beta$, and solve it using a Galerkin method with continuous piecewise affine functions. Introducing $W_{h}$ the space of continuous piecewise affine functions on the mesh defined by the points $(x_{i})_{0\leq i\leq N}$, we denote by $\beta_{h}$ the unique solution in $W_{h}$ to the variational equation:

[TABLE]

Then, the proposed approximation for $\mu$, given by $\mu_{h}=\omega\beta_{h}$, satisfies the following error estimate.

Theorem 7** ([7, Thm. 3.2]).**

If $u_{N}\in U^{s-1}$, for some $\frac{1}{2}\leq s\leq 2$, there holds

[TABLE]

Numerical validation.

In fact, estimates in the local weighted $L^{2}$ norms can be derived from the previous, namely:

[TABLE]

for the Dirichlet problem and

[TABLE]

[TABLE]

for the Neumann problem. We verify those rates numerically. For the Dirichlet problem, we solve two test cases $S_{0,\omega}\alpha_{1}=u_{1}$ and $S_{0,\omega}\alpha_{2}=u_{2}$ having the explicit solutions $\alpha_{1}(x)=\omega(x)$ and $\alpha_{2}=\omega(x)^{3}$, for adequately chosen right hand sides (rhs) $u_{1}$ and $u_{2}$. One can check that $\alpha_{1}\in T^{s}$ for $s<\frac{3}{2}$ and $\alpha_{1}\notin T^{3/2}$, while $\alpha_{2}\in T^{2}$. The $L^{2}_{\frac{1}{\omega}}$ error is plotted in Figure 1 in each case as a function of the mesh size $h$. We find that the expected rates $O(h^{3/2})$ and $O(h^{2})$ predicted by the theory are precisely observed in practice.

Similarly, for the Neumann case, we solve a a test case $N_{0,\omega}\beta=u_{N}$ where the solution $\beta$ is explicit. We take $u_{N}=U_{2}$ the second Chebyshev polynomial of the second kind. The corresponding solution $\beta$ is proportional to $U_{2}$ and thus belongs to $U^{\infty}$. The theory therefore predicts a convergence rate of the error in the $L^{2}_{\omega}$ and $U^{1}$ norms respectively in $O(h^{2})$ and $O(h)$. This behavior is again confirmed by our numerical results, exposed in Figure 2.

5 Building the preconditioners

Let $X_{h}$ the considered finite element space ($X_{h}=V_{h}$ or $W_{h}$), and $(\phi_{i})_{i}$ the basis functions. For an operator $A$, we denote by $\left[A\right]_{p}$ the Galerkin matrix of the operator for the relevant weight $p(x)=\frac{1}{\omega(x)}$ or $\omega(x)$, defined by

[TABLE]

where $\phi_{i}$ and $\phi_{j}$ are the basis functions of the Galerkin space. When the operator $BA$ is a compact perturbation of the identity (either in $T^{s}$ or $U^{s}$) then, following [22, 41], we precondition the linear system $\left[A\right]_{p}x=b$ by the matrix $\left[I_{d}\right]^{-1}_{p}\left[B\right]_{p}\left[I_{d}\right]_{p}^{-1}$, which amounts to solve

[TABLE]

When $B$ is the inverse of a local operator $C$, then it may be more convenient to compute $\left[C\right]_{p}$, and solve instead

[TABLE]

The operators

[TABLE]

introduced in section 2 for $k=0$ (Theorem 1 and Theorem 2) and section $3$ for $k>0$ are at the base of our preconditioning strategy, as $P_{k}S_{k,\omega}$ and $Q_{k}N_{k,\omega}$ are compact perturbations of the identity (Theorem 4). To define preconditioners for the linear systems, following the above remark, we need to compute the Galerkin matrices of those operators. For $Q_{k}$, we rewrite

[TABLE]

This brings us back to computing the Galerkin matrix of the square root of a differential operator. When the frequency is [math], we use the method exposed in [20], relying on the discretization of contour integrals in the complex plane. When the frequency is non-zero, the previous method fails since the spectrum of the matrix contains negative values. We then follow Antoine and Darbas [3] by using a Padé approximation of the square root with regularization and rotation of the branch cut. Let us reproduce here some details of the method, for the reader’s convenience. Consider the classical Padé approximation

[TABLE]

where the coefficients $c_{0}$, $a_{i}$ and $b_{i}$ are given by the formulas

[TABLE]

It is preferable to use a “rotated” version of this approximation to avoid the singularity related to the branch cut for $X<-1$:

[TABLE]

This yields the new approximation

[TABLE]

where

[TABLE]

[TABLE]

This provides a good approximation of the square root in any region of the real line away from $X=-1$, as described by the next result.

Lemma 5**.**

Let $\theta\in\mathbb{R}$ and let $z$ a complex number. Let $r=\left\lvert z+1\right\rvert$. One has

[TABLE]

where

[TABLE]

As a consequence, it is not difficult to check that, if $\theta\in(-\pi,\pi)$, then when $N_{p}\to\infty$, the Padé approximants converge exponentially in the uniform error in any region of the form $\delta<|z+1|<R$ where $0<\delta<1$ and $R>1$. The previous scheme can be exploited to approximate an operator $\sqrt{X-k^{2}I_{d}}$ where $X$ is a positive self-adjoint operator. Writing

[TABLE]

we have

[TABLE]

Using 5, one can conclude that, when $N_{p}$ is large enough, this yields a good approximation in the eigenspaces that are associated to eigenvalues $\lambda$ such that

[TABLE]

In our context, the eigenvalues $\lambda\approx k^{2}$ correspond to the so-called “grazing modes”. To deal with them, Marion Darbas introduced in her thesis [15] a regularization recipe, which consists in adding to the wavenumber some damping. Namely, letting $\varepsilon>0$, the approximation is replaced by

[TABLE]

Based on those considerations, we can approximate the Galerkin matrix of the operator appearing on the left by

[TABLE]

When $X$ is a local operator, this formula involves sparse matrix products and sparse linear system resolutions, which can be performed efficiently. To get the Galerkin matrices of $P_{k}$ and $Q_{k}$, we apply this strategy with $X=-(\omega\partial_{x})^{2}+k^{2}(I_{d}-\omega^{2})$ in the case of the Dirichlet problem and $X=-(\partial_{x}\omega)^{2}+k^{2}(I_{d}-\omega^{2})$ for the Neumann problem. We use the parameters $N_{p}=15$, $\theta=\frac{\pi}{3}$ and $\varepsilon=0.05k^{1/3}$. Our numerical results do not depend crucially on those choices, see subsection 6.4. The choice of $N_{p}$ is dictated by the ratio $\frac{N_{dof}}{k}$ where $N_{dof}$ is the dimension of the Galerkin space. In our tests, this ratio will be held fixed unless stated otherwise, and thus, a fixed value of $N_{p}$ yields good results for all $k$. An informal way to explain this fact is that for our choices of operator $X$, the largest eigenvalues of the Galerkin matrix behave as $\lambda_{M}\leq C_{1}N_{dof}^{2}$. Using a fixed number of points per wavelength gives in turn $N_{dof}=C_{2}k$ and thus

[TABLE]

Choosing $N_{p}$ such that (16) holds with $R=C_{1}C_{2}$, and forgetting about the grazing modes, we thus ensure a small error in (17) independently of $k$. On the other hand, when $k$ is fixed and the mesh is refined, it is necessary to increase $N_{p}$ to maintain accuracy.

6 Numerical results

In this section, we present some numerical results concerning the efficiency of the preconditioners defined above for solving the linear systems arising from the Galerkin method detailed in section 4. The linear systems are solved with the GMRES method [38] with no restart and a tolerance of $\varepsilon=10^{-8}$. Execution times are reported only to show in which case the preconditioned linear system is solved faster than the non-preconditioned system. We report the best timing of three successive runs when the time is less than $5$ seconds. If the number of iterations is greater than 500, we stop the calculations and report the time to reach the 500th iteration. The computations are also interrupted when they last more than $15$ minutes and in this case no iteration number is reported. All the simulations are performed on a personal laptop running on an eight cores intel i7 processor with a clock rate of 2.8GHz. The method is implemented in the language Matlab R2018a. In this section, $k$ stands for the wave number, $N$ for the number of mesh points and $\left\lvert\Gamma\right\rvert$ for the length of a curve $\Gamma$. Moreover, $\omega_{\Gamma}$ is the weight defined on any curve $\Gamma$ by

[TABLE]

where $r:[-1,1]\to\Gamma$ is a parametrization of $\Gamma$ satisfying

[TABLE]

For problems such that $N>5\times 10^{3}$, we use the Efficient Bessel Decomposition [6] to compress the Galerkin matrices.

6.1 Laplace equation on the flat segment

We start by testing the performance of the exact inverses, characterized in section 2, as preconditioners. This is rather a validation stage.

Flat segment, Laplace-Dirichlet problem.

In Table 1, we report the timings and number of GMRES iterations for the Laplace weighted single-layer equation

[TABLE]

Two cases are considered, first without any preconditioner, and then with a preconditioner given by the exact inverse $2\sqrt{-(\omega\partial_{x})^{2}}+\frac{2}{\ln(2)}\pi_{0}$ (see Theorem 1 and the previous section for the detailed construction of the preconditioner). The rhs is chosen as

[TABLE]

A graph of the history of the GMRES relative residual is given in Figure 3 for a mesh with $N=2000$ node points.

Flat segment, Laplace-Neumann problem.

For the Laplace weighted hypersingular equation

[TABLE]

we also report in Table 2 the timings and number of iterations of the GMRES method first without preconditioner, and then with the preconditioner obtained from the operator $[-(\partial_{x}\omega)^{2}]^{-1/2}$ by the method described in the previous section. The rhs is chosen as

[TABLE]

A graph of the history of the GMRES relative residual is given in Figure 3 for a mesh mesh with $N=2000$ node points.

We observe in both Dirichlet and Neumann cases a low and stable number of iteration which is the expected behavior. In the Neumann problem, the presence of the preconditioner leads to huge speedups.

6.2 Helmholtz equation on the flat segment

We now turn our attention to the Helmholtz equation ($k>0$). From now on, unless stated otherwise, the number of segments in the discretization is set to $N\approx 5k\left\lvert\Gamma\right\rvert$ (5 points per wavelength).

Flat segment, Helmholtz-Dirichlet problem.

In Table 3 we report the number of GMRES iterations for the numerical resolution of the weighted single-layer integral equation

[TABLE]

on the flat segment $\Gamma=[-1,1]\times\{0\}$, when the using a preconditioner based on the opretor

[TABLE]

as compared to the case where no preconditioner is used. We take, for the Dirichlet data, the plane wave $u_{D}(x)=e^{ikx}$. We also provide, in 4(a), the history of the GMRES relative residual with and without preconditioner, for a problem with $k\left\lvert\Gamma\right\rvert=200\pi$. When the preconditioner is used, the number of iterations is approximately reduced by a factor $10$.

Flat segment, Helmholtz-Neumann problem.

We run the same numerical comparisons, this time for the Helmholtz weighted hypersingular integral equation

[TABLE]

on $\Gamma=[-1,1]\times\{0\}$ and taking the preconditioner based on the operator

[TABLE]

Results are given in Table 4 for different meshes and in 4(b) for the history of the GMRES relative residual in a case where $k\left\lvert\Gamma\right\rvert=200\pi$. The rhs is chosen as the normal derivative of a diagonal plane wave

[TABLE]

Huge differences, both in time and number of iterations are shown in favor of the preconditioned system.

6.3 Helmholtz equation on non-flat arcs

Recall that for any non-flat arc $\Gamma$, the Helmholtz weighted single- and hypersingular layer potentials are defined by

[TABLE]

Spiral-shaped arc

We first consider a spiral-shaped arc of equation

[TABLE]

[TABLE]

for $t\in[-1,1]$. The curve has a length of about $\left\lvert\Gamma\right\rvert=3.88$. We report in Tables 5 and 6 the number of iterations and computing times respectively for the Dirichlet and Neumann problems with the rhs given by

[TABLE]

where for all $(x,y)$, $u_{inc}(x,y)=e^{ikx}$. To illustrate the problem, the scattering pattern with Neumann boundary conditions is shown in Figure 5 for this geometry. The preconditioning performances are qualitatively similar to the case of the flat segment. This shows that the preconditioning strategy is also efficient in presence of non-zero curvature.

Non-smooth arc

We consider now the case where $\Gamma$ is the V-shaped arc given by the parametric equations

[TABLE]

where $\theta$ is a parameter. When $\theta=\pi$, $\Gamma$ is the flat segment. When $\theta\neq\pi$, this arc has a corner in the middle of angle $\theta$. Since it is not a smooth arc, the theory presented in this work does not apply. For example, the solution $\alpha$ to the weighted single-layer integral equation

[TABLE]

where $u_{D}$ is a smooth function, has a singularity at the corner of $\Gamma$. This is illustrated in Figure 6 where we plot the solution $\alpha(s)$ as a function of the arclength $s$ for a rhs $u_{D}=e^{ikx}$, when $\theta=\frac{\pi}{2}$. Despit this singularity, the number of GMRES iterations remains independent of the mesh size for a fixed frequency. This result is reported in Table 7, where we compare, for $k\left\lvert\Gamma\right\rvert=50$ fixed, the number of GMRES iterations for the resolution of the Helmholtz weighted single-layer integral equation respectively on the flat segment and on the V-shaped curve, for different values of $N$.

We show in Table 8 the influence of the frequency on the preconditioning performances for the Helmholtz weighted single-layer integral equation on the V-shaped arc of angle $\theta=\frac{\pi}{2}$. The results are qualitatively the same as in the case of a smooth curve. To illustrate the problem, the scattering pattern for a sound-hard V-shaped arc of angle $\frac{\pi}{2}$ (Neumann conditions) is shown in Figure 7 for $k\left\lvert\Gamma\right\rvert=50$.

6.4 Influence of the number of Padé approximants

The method is not very sensitive to the number of Padé approximantes, i.e. the parameter $N_{p}$ in Equation 17. We show this in Figure 8 in the case of the Dirichlet problem for the spiral-shaped arc with $k\left\lvert\Gamma\right\rvert=200\pi$ and $u_{D}=e^{ikx}$. The parameter $\varepsilon$ and the angle of the branch rotation $\theta$ remain fixed (see section 5).

6.5 Importance of the correction

It is crucial to include the correct dependence in $k$ in the preconditioners. We report here some numerical results in several situations where this dependence is not respected.

Laplace preconditioning

First, we precondition the Helmholtz weighted single-layer integral equation on the flat segment with the operator

[TABLE]

instead of

[TABLE]

An identity operator is added under the square root to $P^{\prime}_{0}$ to make it invertible. It is easy to check that $2P^{\prime}_{0}$ is spectrally equivalent to the inverse of $S_{0,\omega}$ with

[TABLE]

Since the Laplace preconditioner is also a compact perturbation of the inverse of $S_{k,\omega}$, the theory predicts (see e.g. [22]) that the number of iterations remains bounded when the frequency is fixed and the mesh is refined. This result is confirmed numerically in Figure 9. We see indeed in practice that no matter how the mesh is refined, the number of iterations remains constant. However, we see in the previous example that this number of iterations grows with $k$. Including the dependence in $k$ in the preconditioner reduces a lot this behavior, as illustrated in Figure 10.

No singularity correction

Second, we test the preconditioner without singularity correction, which is the method obtained when we take $\omega\equiv 1$. That is, we solve the non-weighted integral equation

[TABLE]

on the flat segment with a standard $\mathbb{P}^{1}$ Galerkin method, and build a preconditioner based on the operator

[TABLE]

This is the direct application of the method of Antoine and Darbas [3] to the context of an open curve. As stated at the beginning of section 4, if a uniform mesh is used, the Galerkin approximation converges at the rate $O(\sqrt{h})$ only. One remedy is to use a mesh graded towards the edges. A graded mesh of parameter $\beta$ is a mesh such that near the edge, the width of the $i$-th interval is approximately $(ih)^{\beta}$. The parameter $\beta=2$ corresponds to the mesh defined in our Galerkin method defined, and $\beta=5$ is the one that theoretically leads to the same rate of convergence as in our method [34]. In Table 9, we report the number of iterations in the GMRES method for the resolution of eq. (18) preconditioned by $P^{\prime}_{k}$ for $k=10\pi$ on graded meshes for different parameters $\beta$. We compare the results with our preconditioned weighted Galerkin method. In each case, we report the $H^{-\frac{1}{2}}$ error. One can see that mesh-refinement allows to decrease the error at the price of losing the performance of the preconditioner. This justifies the need for the method introduced in this paper.

6.6 Comparison with the generalized Calderón preconditioners

We finally adapt to our context the idea of Bruno and Lintner [10], namely to use $S_{k,\omega}$ and $N_{k,\omega}$ as mutual preconditioners. Notice that the way we discretize the problem is different from [10] where a spectral method is used. In our setting, using the notation of section 5, we define the preconditioners

[TABLE]

respectively for the weighted single-layer and weighted hypersingular integral equations. We report the number of iterations and computing times respectively for the Dirichlet and Neumann problems on the flat segment respectively in Table 10 and Table 11. The performance is compared to that of our new preconditioners. The rhs are respectively $u_{D}=u_{inc}$ and $u_{N}=\frac{\partial u_{inc}}{\partial n}$ where $u_{inc}$ is a plan wave of angle of incidence $\frac{\pi}{4}$. Our results confirm the efficiency of the Generalized Calderón preconditioners, which iteration counts remain very stable with respect to $k$. Despite a slightly larger increase of the iterations for our method, the resolution resolution remains faster for the tests presented here, particularly for the Dirichlet problem. This is due to the fact that our preconditioners are evaluated faster.

7 Conclusion

We have presented a new approach for the preconditioning of integral equations coming from the discretization of wave scattering problems in 2D by open arcs. The methodology is very effective and proven to be optimal for Laplace problems on straight segments. It generalizes the formulas mainly proposed in [3] for regular domains, by the simple addition of a suitable weight. We deeply believe that the methodology opens new perspectives for such problems. First, it is possible to generalize the approach in 3D for the diffraction by a disk (see [7, Chap. 4]). Second, the strategy that we used here seems very likely to be extended to the half line and hopefully to 2D sectors, giving, on the one hand a new pseudo-differential analysis more suitable than classical ones (see e.g. [31, 36, 37]) for handling Helmholtz-like problems on singular domains, and, on the other hand, a completely new preconditioning technique adapted to the treatment of BEM operators on domains with corners or wedges in 3D. Eventually, the weighted square root operators that appear in the present context might well be generalized to give suitable approximation of the exterior Dirichlet to Neumann map for the Helmholtz equation. Having efficient approximations of this map is of particular importance in many contexts, such as e.g. domain decomposition methods.