Smoothed Combined Field Integral Equations for Exterior Helmholtz Problems

TL;DR

This paper introduces smooth combined field integral equations for exterior Helmholtz problems that are easier to evaluate numerically and are uniquely solvable across all frequencies, with a novel regularization of hypersingular operators.

Contribution

It proposes a new regularization method for hypersingular operators that avoids complex principal value integrals and extends the applicability of combined field integral equations.

Findings

Integral equations are smooth and involve only differentiable integrands.

Numerical evaluations near boundaries are accurate and feasible.

Numerical examples demonstrate effectiveness across different discretizations.

Abstract

This paper presents smoothed combined field integral equations for the solution of Dirichlet and Neumann exterior Helmholtz problems. The integral equations introduced in this paper are smooth in the sense that they only involve continuously differentiable integrands in both Dirichlet and Neumann cases. These integral equations coincide with the well-known combined field equations and are therefore uniquely solvable for all frequencies. In particular, a novel regularization of the hypersingular operator is obtained, which, unlike regularizations based on Maue's integration-by-parts formula, does not give rise to involved Cauchy principal value integrals. The smoothed integral operators and layer potentials, on the other hand, can be numerically evaluated at target points that are arbitrarily close to the boundary without severely compromising their accuracy. A variety of numerical examples in two spatial dimensions that consider three different Nystr\"om discretizations for smooth domains and domains with corners---one of which is based on direct application of the trapezoidal rule---demonstrates the effectiveness of the proposed integral approach. In certain aspects, this work extends to the uniquely solvable Dirichlet and Neumann combined field integral equations, the ideas presented in the recent contribution R. Soc. Open Sci. 2(140520), 2015.

Tables

Table 1. Table 1: Number of GMRES iterations needed to solve of the linear systems resulting from the TR, MK and KR-10 Nyström discretizations of the smoothed (SD-CFIE ( 15 ) and SN-CFIE ( 25 )) and non-smoothed (D-CFIE ( 9 ) and N-CFIE ( 10 )) integral equations. The upper panel displays the results for the Dirichlet integral equations, D-CFIE and SD-CFIE, and the lower panel displays the results for the Neumann integral equations, N-CFIE and SN-CFIE. The GMRES tolerance 10 − 6 superscript 10 6 10^{-6} was utilized in all the examples considered in this table.

$k$	$N$	TR SD-CFIE		MK D-CFIE		MK SD-CFIE		KR-10 D-CFIE		KR-10 SD-CFIE
	$\left(\frac{2\pi }{h}\right)$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$
1	30	8	$7.15\times {10}^{-4}$	8	$1.01\times {10}^{-5}$	8	$6.08\times {10}^{-4}$	58	$9.24\times {10}^{-1}$	58	$5.59\times {10}^{-1}$
4	120	10	$9.52\times {10}^{-5}$	10	$2.05\times {10}^{-5}$	10	$2.05\times {10}^{-5}$	83	$3.65\times {10}^{-3}$	83	$4.21\times {10}^{-3}$
16	480	12	$8.15\times {10}^{-5}$	12	$2.19\times {10}^{-5}$	12	$2.19\times {10}^{-5}$	214	$2.38\times {10}^{-3}$	215	$2.13\times {10}^{-3}$
64	1920	16	$1.05\times {10}^{-4}$	16	$1.39\times {10}^{-5}$	16	$1.39\times {10}^{-5}$	793	$1.73\times {10}^{-3}$	809	$1.67\times {10}^{-3}$
$k$	$N$	TR SN-CFIE		MK N-CFIE		MK SN-CFIE		KR-10 N-CFIE		KR-10 SD-CFIE
	$\left(\frac{2\pi }{h}\right)$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$	It.	${\epsilon }_{\infty }$
1	30	17	$3.35\times {10}^{-4}$	20	$2.16\times {10}^{-5}$	20	$7.32\times {10}^{-5}$	28	$2.47\times {10}^{-1}$	58	$2.88\times {10}^{-0}$
4	120	24	$1.80\times {10}^{-4}$	29	$5.49\times {10}^{-5}$	29	$4.74\times {10}^{-5}$	103	$3.32\times {10}^{-3}$	84	$3.41\times {10}^{-3}$
16	480	23	$1.00\times {10}^{-4}$	27	$3.04\times {10}^{-5}$	27	$3.03\times {10}^{-5}$	346	$1.45\times {10}^{-3}$	230	$9.61\times {10}^{-4}$
64	1920	12	$1.03\times {10}^{-4}$	14	$2.22\times {10}^{-5}$	14	$2.19\times {10}^{-5}$	998	$1.29\times {10}^{-3}$	719	$7.97\times {10}^{-4}$

Equations

$$\begin{array}{rl}{\int }_{\Gamma }K(\mathbf{x},\mathbf{y})\phi (\mathbf{y})\mathrm{d}s(\mathbf{y})=& \sum _{j=0}^N{\partial }_s^j\phi (\mathbf{x}){\int }_{\Gamma }\{K_1(\mathbf{x},\mathbf{y})p_j(\mathbf{y}|\mathbf{x})-K_2(\mathbf{x},\mathbf{y}){\partial }_n{p}_j(\mathbf{y}|\mathbf{x})\}\mathrm{d}s(\mathbf{y})\\ & +{\int }_{\Gamma }\{K_1(\mathbf{x},\mathbf{y}){\rho }_1(\mathbf{y}|\mathbf{x})-K_2(\mathbf{x},\mathbf{y}){\rho }_2(\mathbf{y}|\mathbf{x})\}\mathrm{d}s(\mathbf{y})\end{array}$$

$$\begin{array}[]{rclll}\Delta u_{D}+k^{2}u_{D}&=&0&\mbox{in}&\mathbb{R}^{2}\setminus\overline{\Omega},\vskip 6.0pt plus 2.0pt minus 2.0pt\\ u_{D}&=&-u^{\mathrm{inc}}&\mbox{on}&\Gamma,\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \lx@intercol\hfil\displaystyle\lim_{|\boldsymbol{r}|\to\infty}\sqrt{|\boldsymbol{r}|}\left(\frac{\partial u_{D}}{\partial|\boldsymbol{r}|}-iku_{D}\right)\hfil\lx@intercol&=&0,\end{array}$$

$$\begin{array}[]{rclll}\Delta u_{N}+k^{2}u_{N}&=&0&\mbox{in}&\mathbb{R}^{2}\setminus\overline{\Omega},\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \displaystyle\partial_{n}u_{N}&=&-\partial_{n}u^{\mathrm{inc}}&\mbox{on}&\Gamma,\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \lx@intercol\hfil\displaystyle\lim_{|\boldsymbol{r}|\to\infty}\sqrt{|\boldsymbol{r}|}\left(\frac{\partial u_{N}}{\partial|\boldsymbol{r}|}-iku_{N}\right)\hfil\lx@intercol&=&0,\end{array}$$

$$\mathcal{S}[\phi ](\mathbf{r})={\int }_{\Gamma }G(\mathbf{r},\mathbf{y})\phi (\mathbf{y})\mathrm{d}s(\mathbf{y})\text{mbox}and\mathcal{D}[\phi ](\mathbf{r})={\int }_{\Gamma }\frac{\partial G(\mathbf{r},\mathbf{y})}{\partial \mathbf{n}(\mathbf{y})}\phi (\mathbf{y})\mathrm{d}s(\mathbf{y}),\mathbf{r}\in {\mathbb{R}}^2\setminus \Gamma ,$$

$$\begin{array}{r}\mathcal{S}[\phi ]=S[\phi ],{\partial }_n\mathcal{S}[\phi ]=-\frac{\phi }{2}+K^{\prime }[\phi ],\mathcal{D}[\phi ]=\frac{\phi }{2}+K[\phi ]\text{mbox}and{\partial }_n\mathcal{D}[\phi ]=N[\phi ],\end{array}$$

$$S[\phi ](\mathbf{x})$$

$$K[\phi ](\mathbf{x})$$

$$K^{\prime }[\phi ](\mathbf{x})$$

$$N[\phi ](\mathbf{x})$$

$$N[\phi ](\mathbf{x})=k^2{\int }_{\Gamma }G(\mathbf{x},\mathbf{y})\mathbf{n}(\mathbf{x})\cdot \mathbf{n}(\mathbf{y})\phi (\mathbf{y})\mathrm{d}s(\mathbf{y})+p.v.{\int }_{\Gamma }{\partial }_{s}G(\mathbf{x},\mathbf{y}){\partial }_s\phi (\mathbf{y})\mathrm{d}s(\mathbf{y})$$

$$\Gamma =\{x(t):t\in [0,2\pi )\},$$

$$u(\mathbf{r})=(\mathcal{D}-i\eta \mathcal{S})[\phi ](\mathbf{r}),\mathbf{r}\in {\mathbb{R}}^2\setminus \overline{\Omega },$$

$$\left(\frac{I}{2}+K-i\eta S\right)[\phi ]=-u^{\mathrm{inc}}\text{mbox}on\Gamma$$

$$\left(\frac{i\eta }{2}I+N-i\eta K^{\prime }\right)[\phi ]=-{\partial }_n{u}^{\mathrm{inc}}\text{mbox}on\Gamma ,$$

$$\begin{array}{c}\Delta p_0(\mathbf{r}|{\mathbf{x}}_0)+k^2{p}_0(\mathbf{r}|{\mathbf{x}}_0)=0,\mathbf{r}\in {\mathbb{R}}^2,\\ p_0(\mathbf{r}|{\mathbf{x}}_0)=1,{\partial }_n{p}_0(\mathbf{r}|{\mathbf{x}}_0)=i\eta ,{\partial }_s{p}_0(\mathbf{r}|{\mathbf{x}}_0)=0\text{mbox}and{\partial }_s{\partial }_n{p}_0(\mathbf{r}|{\mathbf{x}}_0)=0\text{mbox}at\mathbf{r}={\mathbf{x}}_0,\end{array}$$

$$\begin{array}{c}\Delta p_1(\mathbf{r}|{\mathbf{x}}_0)+k^2{p}_1(\mathbf{r}|{\mathbf{x}}_0)=0,\mathbf{r}\in {\mathbb{R}}^2,\\ p_1(\mathbf{r}|{\mathbf{x}}_0)=0,{\partial }_n{p}_1(\mathbf{r}|{\mathbf{x}}_0)=0,{\partial }_s{p}_1(\mathbf{r}|{\mathbf{x}}_0)=1\text{mbox}and{\partial }_s{\partial }_n{p}_1(\mathbf{r}|{\mathbf{x}}_0)=i\eta \text{mbox}at\mathbf{r}={\mathbf{x}}_0,\end{array}$$

$$\mathcal{D}\left[p_0(\cdot |{\mathbf{x}}_0)\right](\mathbf{r})-\mathcal{S}\left[{\partial }_n{p}_0(\cdot |{\mathbf{x}}_0)\right](\mathbf{r})$$

$$\mathcal{D}\left[p_1(\cdot |{\mathbf{x}}_0)\right](\mathbf{r})-\mathcal{S}\left[{\partial }_n{p}_1(\cdot |{\mathbf{x}}_0)\right](\mathbf{r})$$

$$\begin{array}{rl}u(\mathbf{r})=& \mathcal{D}\left[\phi -\phi ({\mathbf{x}}_0)p_0(\cdot |{\mathbf{x}}_0)-{\partial }_s\phi ({\mathbf{x}}_0)p_1(\cdot |{\mathbf{x}}_0)\right](\mathbf{r})\\ & -\mathcal{S}\left[i\eta \phi -\phi ({\mathbf{x}}_0){\partial }_n{p}_0(\cdot |{\mathbf{x}}_0)-{\partial }_s\phi ({\mathbf{x}}_0){\partial }_n{p}_1(\cdot |{\mathbf{x}}_0)\right](\mathbf{r}),\mathbf{r}\in {\mathbb{R}}^2\setminus \overline{\Omega },{\mathbf{x}}_0\in \Gamma .\end{array}$$

$$\begin{array}{rl}u(\mathbf{x})=& \frac{\phi (\mathbf{x})-\phi ({\mathbf{x}}_0)p_0(\mathbf{x}|{\mathbf{x}}_0)-{\partial }_s\phi ({\mathbf{x}}_0)p_1(\mathbf{x}|{\mathbf{x}}_0)}{2}+K\left[\phi -\phi ({\mathbf{x}}_0)p_0(\cdot |{\mathbf{x}}_0)-{\partial }_s\phi ({\mathbf{x}}_0)p_1(\cdot |{\mathbf{x}}_0)\right](\mathbf{x})\\ & -S\left[i\eta \phi -\phi ({\mathbf{x}}_0){\partial }_n{p}_0(\cdot |{\mathbf{x}}_0)-{\partial }_s\phi ({\mathbf{x}}_0){\partial }_n{p}_1(\cdot |{\mathbf{x}}_0)\right](\mathbf{x})\text{mbox}forall\mathbf{x},{\mathbf{x}}_0\in \Gamma ,\end{array}$$

$$\left(K\circ R_D-S\circ R_S\right)[\phi ]=-u^{\mathrm{inc}}\text{mbox}on\Gamma ,$$

$$\begin{array}{rl}{R}_D[\phi |\mathbf{x}](\mathbf{y})=& \phi (\mathbf{y})-\phi (\mathbf{x})p_0(\mathbf{y}|\mathbf{x})-{\partial }_s\phi (\mathbf{x})p_1(\mathbf{y}|\mathbf{x})\text{mbox}and\\ R_S[\phi |\mathbf{x}](\mathbf{y})=& i\eta \phi (\mathbf{y})-\phi (\mathbf{x}){\partial }_n{p}_0(\mathbf{y}|\mathbf{x})-{\partial }_s\phi (\mathbf{x}){\partial }_n{p}_1(\mathbf{y}|\mathbf{x}).\end{array}$$

$$R_D[\phi |\mathbf{x}](\mathbf{y})=O(|\mathbf{x}-\mathbf{y}{|}^2)\text{mbox}and{R}_S[\phi |\mathbf{x}](\mathbf{y})=O(|\mathbf{x}-\mathbf{y}{|}^2)\text{mbox}as\mathbf{y}\to \mathbf{x},(\mathbf{y}\in \Gamma ),$$

$${\rho }_D(\tau |t)=R_D[\phi |x(t)](x(\tau ))=\varphi (\tau )-\varphi (t){\tilde{p}}_0(\tau |t)-|x^{\prime }(t){|}^{-1}{\varphi }^{\prime }(t){\tilde{p}}_1(\tau |t),$$

$$\begin{array}{rl}{\rho }_D(\tau |t)=& \varphi (t)+{\varphi }^{\prime }(t)(\tau -t)+\frac{{\varphi }^{\prime \prime }(t)}{2}(\tau -t{)}^2+o(|\tau -t{|}^2)\\ & -\varphi (t)\{1+O(|\tau -t{|}^2)\}-{\varphi }^{\prime }(t)\{\tau -t+O(|\tau -t{|}^2)\}=O(|\tau -t{|}^2)\text{mbox}as\tau \to t,\end{array}$$

$$\begin{array}{rl}{\tilde{p}}_0(t|t)=p_0(x(t)|x(t))=1,& {\tilde{p}}^{\prime }(t|t)=|x^{\prime }(t)|{\partial }_s{p}_0(x(t)|x(t))=0,\\ {\tilde{p}}_1(t|t)=p_1(x(t)|x(t))=0,& {\tilde{q}}^{\prime }(t|t)=|x^{\prime }(t)|{\partial }_s{p}_1(x(t)|x(t))=|x^{\prime }(t)|,\end{array}$$

$$\begin{array}{rl}{\partial }_n{\tilde{p}}_0(t|t)={\partial }_n{p}_0(x(t)|x(t))=i\eta ,& {\partial }_n{\tilde{p}}_0^{\prime }(t|t)=|x^{\prime }(t)|{\partial }_s{\partial }_n{p}_0(x(t)|x(t))=0,\\ {\partial }_n{\tilde{p}}_1(t|t)={\partial }_n{p}_1(x(t)|x(t))=0,& {\partial }_n{\tilde{p}}_1^{\prime }(t|t)=|x^{\prime }(t)|{\partial }_s{\partial }_n{p}_1(x(t)|x(t))=i\eta |x^{\prime }(t)|,\end{array}$$

$${\rho }_S(\tau |t)=R_S[\phi |x(t)](x(\tau ))=i\eta \varphi (\tau )-\varphi (t){\partial }_n{\tilde{p}}_0(\tau |t)-|x^{\prime }(t){|}^{-1}{\varphi }^{\prime }(t){\partial }_n{\tilde{p}}_1(\tau |t)$$

$${\rho }_S(\tau |t)=O(|\tau -t{|}^2)\text{mbox}as\tau \to t.$$

$$v(t)={\int }_0^{2\pi }\left\{L(t,\tau ){\rho }_D(\tau |t)-M(t,\tau ){\rho }_S(\tau |t)\right\}\mathrm{d}\tau ,t\in [0,2\pi ),$$

$$\begin{array}{rl}L(t,\tau )=& \frac{ik}{4}\frac{H_1^{(1)}(kR)}{R}(x(t)-x(\tau ))\cdot n(\tau )|x^{\prime }(\tau )|=L_1(t,\tau )\log |t-\tau |+L_2(t,\tau ),\\ M(t,\tau )=& \frac{i}{4}{H}_0^{(1)}(kR)|x^{\prime }(\tau )|=M_1(t,\tau )\log |t-\tau |+M_2(t,\tau ),\end{array}$$

$$L_1(t,\tau ):=-\frac{k}{2\pi }\frac{J_1(kR)}{R}(x(t)-x(\tau ))\cdot n(\tau )|x^{\prime }(\tau )|,L_2(t,\tau ):=L(t.\tau )-L_1(t,\tau )\log |t-\tau |$$

Full text

Keywords: Combined field integral equation, regularization, hypersingular operator, Helmholtz equation, Nyström discretization

1 Introduction

As is well known, boundary integral equation (BIE) methods, such as boundary element methods [5, 34] as well as Nyström methods [9, 23, 24, 26, 28], provide several advantages over methods based on volume discretization of the computational domain, such as finite difference [36] and finite element methods [20], for the solution of exterior Helmholtz problems. For example, BIE methods can easily handle unbounded domains and radiation conditions at infinity without recourse to approximate absorbing/transparent boundary conditions for truncation of the computational domain [16]. Additionally, BIE methods are based on discretization of the relevant physical boundaries, and they therefore give rise to linear systems of reduced dimensionality—which, although dense, can be efficiently solved by means of accelerated iterative linear algebra solvers [4, 9, 17, 32].

One of the main issues associated with the use of BIE methods is the numerical evaluation of the challenging singular, weakly-singular and nearly-singular integrals that are inherent to the integral operators and layer potentials upon which BIE methods are based on. In two spatial dimensions, for example, the single-layer, double-layer and adjoint double-layer operators feature weak $O(\log|\boldsymbol{x}-\boldsymbol{y}|)$ kernel singularities, while the hypersingular operator features a much more stronger $O(|\boldsymbol{x}-\boldsymbol{y}|^{-2})+O(\log|\boldsymbol{x}-\boldsymbol{y}|)$ kernel singularity as $\boldsymbol{y}\to\boldsymbol{x}$, where $\boldsymbol{x}$ and $\boldsymbol{y}$ denote points on the (assumed smooth) boundary. As is known, however, application of the standard regularization procedure, which was originally proposed by Maue [29], enables the hypersingular operator to be expressed in terms of a Cauchy principal value integral that exhibits a $O(|\boldsymbol{x}-\boldsymbol{y}|^{-1})+O(\log|\boldsymbol{x}-\boldsymbol{y}|)$ kernel singularity as $\boldsymbol{y}\to\boldsymbol{x}$. Nearly-singular integrals arise, on the other hand, when integral operators and layer potentials are evaluated at target points close to but not on the boundary of the domain. All these issues greatly hinder the use of BIE methods, as numerical evaluation of integral operators and layer potentials requires special treatment of the kernel singularities by means of specialized quadrature rules and/or semi-analytical techniques for which there is a vast literature that will not be reviewed here; cf. [2, 3, 5, 9, 8, 12, 14, 24, 23, 22, 33, 34].

This paper presents uniquely solvable BIEs for the solution of exterior Helmholtz problems with Dirichlet and Neumann boundary conditions, that arise from electromagnetic scattering by perfectly conducting obstacles in two spatial dimensions. The BIEs introduced in this contribution, dubbed smoothed combined field integral equations, coincide with the common combined field integral equations (CFIEs) [11, 31, 7, 27] in both Dirichlet and Nuemann cases. Unlike the common CFIEs, however, they are given in terms of operators expressed as integrals of continuously differentiable functions (they indeed exhibit a mild singularity of the form $O(|\boldsymbol{x}-\boldsymbol{y}|^{2}\log|\boldsymbol{x}-\boldsymbol{y}|)$ as $\boldsymbol{y}\to\boldsymbol{x}$). In particular, we hereby introduce a novel regularization of the hypersingular operator that involves neither Cauchy principal value nor weakly singular integrals. The proposed smoothing procedure is also utilized to regularize nearly-singular integrals that arise from evaluation of layer potentials at target points near the boundary, and from integral operators that result from integral formulations of problems involving two or more obstacles close to each other.

Our smoothing procedure relies on the existence of certain homogenous solutions $p_{j}$, $j=0,\ldots,N$ ($N>0$) of the Helmholtz equation, that we referred to as smoothing functions. Such functions allow a sufficiently smooth density $\varphi$ to be expressed as $\varphi(\boldsymbol{y})=\sum_{j=0}^{N}\partial_{s}^{j}\varphi(\boldsymbol{x})p_{j}(\boldsymbol{y}|\boldsymbol{x})+O(|\boldsymbol{x}-\boldsymbol{y}|^{N+1})$ and $i\eta\varphi(\boldsymbol{y})=\sum_{j=0}^{N}\partial_{s}^{j}\varphi(\boldsymbol{x})\partial_{n}p_{j}(\boldsymbol{y}|\boldsymbol{x})+O(|\boldsymbol{x}-\boldsymbol{y}|^{N+1})$ where $\boldsymbol{x}$ is a given point on the boundary $\Gamma$ and where $\eta>0$ is a constant (the symbols $\partial_{s}$ and $\partial_{n}$ denote tangential and normal derivatives on $\Gamma$). Calling $K=K_{1}-i\eta K_{2}$ the kernel of the Dirichlet or Neumann combined field operators, we thus can write

[TABLE]

where $\rho_{1}(\boldsymbol{y}|\boldsymbol{x})=\varphi(\boldsymbol{y})-\sum_{j=0}^{N}\partial_{s}^{j}\varphi(\boldsymbol{x})p_{j}(\boldsymbol{y}|\boldsymbol{x})$ and $\rho_{2}(\boldsymbol{y}|\boldsymbol{x})=i\eta\varphi(\boldsymbol{y})-\sum_{j=0}^{N}\partial_{s}^{j}\varphi(\boldsymbol{x})\partial_{n}p_{j}(\boldsymbol{y}|\boldsymbol{x})$. As it turns out, Green’s third identity provides closed-form expressions for the boundary integrals inside the sum. Therefore, the operator $\int_{\Gamma}K(\boldsymbol{x},\boldsymbol{y})\varphi(\boldsymbol{y})\,\mathrm{d}s(\boldsymbol{y})$ can be easily evaluated by integrating the smoothed mildly singular functions $K_{j}(\boldsymbol{x},\boldsymbol{y})\rho_{j}(\boldsymbol{y}|\boldsymbol{x})$, $j=1,2$—which satisfy $K_{j}(\boldsymbol{x},\boldsymbol{y})\rho_{j}(\boldsymbol{y}|\boldsymbol{x})=O(|\boldsymbol{x}-\boldsymbol{y}|^{N+1}\log|\boldsymbol{x}-\boldsymbol{y}|)$—wherever on $\Gamma$ the tangential derivatives $\partial_{s}^{j}\varphi(\boldsymbol{x})$, $j=1,\ldots,N,$ exist. In this paper we present a smoothing procedure that considers functions $p_{j}$, $j=0,1$, that are obtained explicitly as linear combinations of plane waves.

A smoothing procedure similar in nature to the one presented here was originally introduced in [21] for the solution of the Laplace equation and was later extended in [35] to the Helmholtz equation in three spatial dimensions. Both contributions consider integral equations derived from direct use of Green’s third identity. As such, the associated smoothed integral equations for the Helmholtz equation suffer from spurious resonances in both Dirichlet and Neumann cases [13]. The smoothing procedure introduced in those references, on the other hand, which provides a smoothing factor that turns weakly singular integrands (in three-dimensions) into bounded but discontinuous functions, does not suffice for the regularization of the hypersingular operator that appears in the combined field integral equation for the Neumann problem.

The structure of this paper is as follows. Section 2 presents the boundary value problems considered in this paper and reviews the definition and main properties of the layer potentials and boundary integral operators. Section 3, subsequently, introduces the smoothed CFIE formulations for both Dirichlet and Nuemann problems. Details on the construction of the smoothing functions are provided in Section 4. Finally, Section 5 presents a variety of numerical examples in two spatial dimensions that include three different Nyström discretizations for smooth domains and domains with corners.

2 Preliminaries

This paper considers exterior Helmholtz boundary value problems that arise as an incident TE- or TM-polarized electromagnetic wave impinges on the surface of an axially symmetric perfect electric conductor with cross section $\Omega\subset\mathbb{R}^{2}$ and boundary $\partial\Omega=\Gamma$. In TE-polarization the scattered field $u_{D}:\mathbb{R}^{2}\setminus\Omega\to\mathbb{C}$ is solution of the exterior Dirichlet problem

[TABLE]

where $k\in\mathbb{C}$, $\mathrm{Im}\,k\geq 0$, $\mathrm{Re}\,k>0$ denotes the wavenumber of the unbounded medium surrounding $\Omega$. In TM-polarization, on the other hand, the scattered field $u_{N}:\mathbb{R}^{2}\setminus\Omega\to\mathbb{C}$ is solution the exterior Neumann problem

[TABLE]

where the symbol $\partial_{n}$ in (2) denotes the exterior normal derivative on the boundary $\Gamma$.

As is well-known (cf. [13]) for a continuous boundary data the exterior Dirichlet (1) (resp. Neumann  (2)) admits a unique solution $u_{D}$ (resp. $u_{N}$) for all wavenumbers $k\in\mathbb{C}$, $\mathrm{Re}\,k>0$, $\mathrm{Im}\,k\geq 0$.

Given a density function $\varphi:\Gamma\to\mathbb{C}$ we define the single- and double-layer potentials as

[TABLE]

respectively, where $G(\boldsymbol{r},\boldsymbol{y})=\frac{i}{4}H_{0}^{(1)}(k|\boldsymbol{r}-\boldsymbol{y}|)$ is the free-space Green function for the Helmholtz equation. Evaluation of the layer potentials (3) and their exterior normal derivatives on $\Gamma$ yields the jump relations [13]

[TABLE]

which are expressed in terms of the single-layer ($S$), double-layer ($K$), adjoint double-layer $(K^{\prime})$ and hypersingular $(H)$ operators. These operators are given by the integral expressions

[TABLE]

for $\boldsymbol{x}\in\Gamma$. Note that the integral in the definition of the hypersingular operator (5d) must be understood as a Hadamard finite-part integral which, upon integration by parts, can expressed as

[TABLE]

in terms of a Cauchy principal value integral and the tangential derivative $\partial_{s}$ of the surface density $\varphi$ on $\Gamma$. The expression (6) is sometimes called Maue’s integration by parts formula [29] and can be interpreted as regularization of the integral operator in the sense that it involves integrands are “smoother” (see Section 1).

We assume, for the time being, that the boundary of the PEC obstacle $\Gamma=\partial\Omega$ admits a real analytic $2\pi$-periodic parametric representation

[TABLE]

where $|\mathbb{x}^{\prime}(t)|\neq 0$ for all $t\in[0,2\pi)$. (This smoothness assumption on $\Gamma$ is relaxed in Section 5.3 where numerical examples for domains with corners are considered.) Utilizing the boundary parameterization we define the Sobolev space $H^{s}(\Gamma)$, $s>0$, as the space of functions $\varphi\in L^{2}(\Gamma)$ such that $\varphi\circ\mathbb{x}\in H^{s}[0,2\pi]$, where $H^{s}[0,2\pi]=\{v\in L^{2}[0,2\pi]:\|v\|_{s}<\infty\}$ (see [25, Chapter 8] or [33, Section 5.3] for a more detailed definition of this space and its properties). In particular, the integral operators (5): $S:H^{s}(\Gamma)\to H^{s+1}(\Gamma)$, $K^{\prime},K:H^{s}(\Gamma)\to H^{s+3}(\Gamma)$ and $N:H^{s}(\Gamma)\to H^{s-1}(\Gamma)$ are continuous for all $s>0$ [30].

In order to solve the exterior boundary value problems we look for solutions given by the combined double- and single-layer potential [7, 27, 31]

[TABLE]

where the density function $\varphi$ can be determined by matching the potential $u$ or its normal derivative $\partial_{n}u$ with the appropriate boundary data on $\Gamma$.

Therefore, in the case of Dirichlet problem (1) we obtain the Dirichlet Combined Field Integral Equation (D-CFIE)

[TABLE]

for an unknown density function $\varphi$, which was obtained evaluating the potential (8) on $\Gamma$ using the jump conditions (4). Since $K$ and $S$ are compact operators on $H^{s}(\Gamma)$ we have that the combined field integral operator $K-i\eta S$ is also a compact operator on $H^{s}(\Gamma)$, and thus (9) is a Fredholm integral equation of the second second-kind. Therefore, the well-posedness of the D-CFIE follows from the fact that (9) admits at most one solution for all wavenumbers $\mathrm{Re}\,k>0$, $\mathrm{Im}\,k\geq 0$, provided $\eta>0$ [13, Theorem 3.3].

In the case of the Neumann problem (2), on the other hand, evaluation of the normal derivative of the combined potential (8) on $\Gamma$ yields the Neumann Combined Field Integral Equation (N-CFIE)

[TABLE]

for an unknown density function $\varphi$. Here, $\frac{i\eta}{2}I+N-i\eta K^{\prime}:H^{s}\to H^{s-1}$ is not a compact operator on $H^{s}(\Gamma)$. As is known, however, the N-CFIE admits at most one solution $\varphi$ for all wavenumbers $\mathrm{Re}\,k>0$, $\mathrm{Im}\,k\geq 0$, provided $\eta>0$ [13, Theorem 3.34].

3 Smoothed combined field integral equation formulations

In this section we introduce smoothed versions of the combined potential (8) and associated CFIEs (9) and (10).

3.1 Dirichlet problem

Let us assume that we are given two $C^{\infty}$-smooth functions $p_{0}(\>\cdot\>|\boldsymbol{x}_{0})$ and $p_{1}(\>\cdot\>|\boldsymbol{x}_{0})$ that satisfy

[TABLE]

and

[TABLE]

respectively, where $\boldsymbol{x}_{0}$ is given a point on the boundary $\Gamma$. Assume further that both functions $p_{0}(\>\cdot\>|\boldsymbol{x}_{0})$ and $p_{1}(\>\cdot\>|\boldsymbol{x}_{0})$ are given by certain linear combinations of plane waves. (Expressions for such functions are given in Section 4.) Therefore, Green’s third identity together with standard stationary phase arguments yield the relations

[TABLE]

Multiplying (13a) and (13b) through by $\varphi(\boldsymbol{x}_{0})$ and $\partial_{s}\varphi(\boldsymbol{x}_{0})$, respectively, and subtracting the resulting expressions from the combined potential (8), we obtain

[TABLE]

Evaluating the potential (14) on $\Gamma$, using the jump conditions (4), we get

[TABLE]

where the operators $K$ and $S$ are defined in (5b) and (5a), respectively. Therefore, selecting $\boldsymbol{x}_{0}=\boldsymbol{x}$ in the relation above and using the identities $p_{0}(\boldsymbol{x}|\boldsymbol{x})=1$ and $p_{1}(\boldsymbol{x}|\boldsymbol{x})=0$, we obtain the Smoothed Dirichlet Combined Field Integral Equation (SD-CFIE)

[TABLE]

for the unknown density function $\varphi$. Here, the operators $R_{D}$ and $R_{S}$ are explicitly defined in terms of the smoothing functions $p_{0}$ and $p_{1}$ by

[TABLE]

The following lemma establishes the essential property of $R_{D}$ and $R_{S}$:

Lemma 3.1**.**

Given $\boldsymbol{x}\in\Gamma$, the operators $R_{D}[\>\cdot\>|\boldsymbol{x}]$ $R_{S}[\>\cdot\>|\boldsymbol{x}]:H^{s}(\Gamma)\to H^{s}(\Gamma)$ introduced in (16), are well defined for all $s>3/2$. Furthermore, they satisfy

[TABLE]

for all $s>5/2$.

Proof.

Clearly $R_{D}$ and $R_{S}$ admit the representations $R_{D}=I-p_{0}(\,\cdot\,|\boldsymbol{x})\delta_{\boldsymbol{x}}-p_{1}(\,\cdot\,|\boldsymbol{x})\delta^{\prime}_{\boldsymbol{x}}$ and $R_{S}=i\eta I-\partial_{n}p_{0}(\,\cdot\,|\boldsymbol{x})\delta_{\boldsymbol{x}}-\partial_{n}p_{1}(\,\cdot\,|\boldsymbol{x})\delta^{\prime}_{\boldsymbol{x}}$ in terms of the Dirac’s distribution $\delta_{\boldsymbol{x}}$ (supported at $\boldsymbol{x}\in\Gamma$) and its derivative, both of which belong to $H^{-s}(\Gamma)$ for all $s>3/2$. Since $p_{0}(\>\cdot\>|\boldsymbol{x}),\partial_{n}p_{0}(\>\cdot\>|\boldsymbol{x})\in H^{s}(\Gamma)$ for all $s\in\mathbb{R}$, on the other hand, it readily follows that $R^{(1)}_{D}[\varphi|\boldsymbol{x}],R^{(1)}_{S}[\varphi|\boldsymbol{x}]\in H^{s}(\Gamma)$, $s>3/2$, for any given point $\boldsymbol{x}\in\Gamma$.

We now prove the asymptotic identities in (17). Let the parameter values $t,\tau\in[0,2\pi)$ be such that $\boldsymbol{x}=\mathbb{x}(t)$ and $\boldsymbol{y}=\mathbb{x}(\tau)$, where $\mathbb{x}:[0,2\pi)\to\Gamma$ denotes the parametrization of the smooth curve $\Gamma$, and let $\phi=\varphi\circ\mathbb{x}$ where $\varphi\in H^{s}(\Gamma)$, $s>5/2.$ Since for $s>5/2$ the density $\phi$ is a twice-continuously differentiable $2\pi$-periodic function on $[0,2\pi)$ [33, Lemma 5.3.3], we have that

[TABLE]

with $\tilde{p}_{0}(\tau|t)=p_{0}(\mathbb{x}(\tau)|\mathbb{x}(t))$ and $\tilde{p}_{1}(\tau|t)=p_{1}(\mathbb{x}(\tau)|\mathbb{x}(t)),$ is also a twice-continuously differentiable 2$\pi$-periodic function. Therefore, expansing $\rho_{D}(\tau|t)$ as a Taylor series around $\tau=t$, we get

[TABLE]

where we have utilized the identities

[TABLE]

that follow from the point conditions in (11) and (12) satisfied by $p$ and $q$, respectively.

Similarly, using the identities

[TABLE]

where $\partial_{n}\tilde{p}_{0}(\tau|t)=\partial_{n}p_{0}(\mathbb{x}(\tau)|\mathbb{x}(t))$ and $\partial_{n}\tilde{p}_{1}(\tau|t)=\partial_{n}p_{1}(\mathbb{x}(\tau)|\mathbb{x}(t)),$ it can be shown that the function

[TABLE]

satisfies

[TABLE]

Therefore, finally, the identities in (17) follow from the fact that $O(|\mathbb{x}(t)-\mathbb{x}(\tau)|^{2})=O(|\tau-t|^{2})$. The proof is now complete. ∎

In order to illustrate the result of Lemma 3.1 we present Figure 1 which displays the functions $\rho_{D}(\tau|t)$ and $\rho_{S}(\tau|t)$ defined in (18) and (20), respectively, obtained by application of $R_{D}$ and $R_{S}$ to a certain smooth density function $\varphi$. As can be observe in this figure, both functions $\rho_{D}(\tau|t)$ and $\rho_{S}(\tau|t)$ vanish quadratically along the line $\tau=t$.

We are now in position to study the singular character of integral operator $K\circ R_{D}-S\circ R_{S}$ in the SD-CFIE (15). As in the proof of Lemma 3.1 we utilize the parametrization $\mathbb{x}$ of the smooth curve $\Gamma$ to define the $2\pi$-periodic function $\phi(t)=\varphi(\mathbb{x}(t))$ as well as the $2\pi$-biperiodic functions $\rho_{D}(\tau|t)$ and $\rho_{S}(\tau|t)$ given in (18) and (20). Using these notations and letting $R=R(t,\tau)=|\mathbb{x}(t)-\mathbb{x}(\tau)|$ and $\mathbb{n}(\tau)=(-x^{\prime}_{2}(\tau),x^{\prime}_{1}(\tau))/|\mathbb{x}^{\prime}(\tau)|$ we have that $v(t)=(K\circ R_{D}-S\circ R_{S})[\varphi](\mathbb{x}(t))$ can be expressed as

[TABLE]

in terms of the weakly-singular kernels

[TABLE]

where the functions

[TABLE]

[TABLE]

can be properly defined at $\tau=t$ so that they are in fact $2\pi$-biperiodic analytic functions [12].

It thus follows from (23), Lemma 3.1 and the fact that $L(t,\tau)=O(|t-\tau|\log|t-\tau|)$, that the integrands in (22) satisfy

[TABLE]

for $\varphi\in H^{s}(\Gamma)$ for $s>5/2$. As it turns out, it can be easily shown that for $s>5/2$ the integrands in (22) are indeed continuously differentiable $2\pi$-biperiodic functions (in $t$ and $\tau$).

To illustrate the smoothness of $L(t,\tau)\rho_{D}(\tau|t)$ and $M(t,\tau)\rho_{S}(\tau|t)$ along the line $\tau=t$ we present Figure 2 which displays both functions for a given smooth density function $\varphi$.

We finish this section with Theorem 3.2 whose proof follows directly from the discussion above:

Theorem 3.2**.**

Let $S$ and $K$ be the single- and double-layer operators defined in (5a) and (5b), respectively, and $R_{D}$ and $R_{S}$ be the smoothing operators defined in (16). Then, the identity

[TABLE]

holds true on $H^{s}(\Gamma)$ for all $s>3/2$. Therefore, in particular, the SD-CFIE (15) is uniquely solvable on $H^{s}(\Gamma)$, $s>3/2$ for all wavenumbers $k\in\mathbb{C}$, $\mathrm{Re}\,k>0$, $\mathrm{Im}\,k\geq 0$, provided $u^{\mathrm{inc}}|_{\Gamma}\in H^{s}(\Gamma)$ and $\eta>0$.

3.2 Neumann problem

We now proceed to derive the smoothed integral equation for exterior Neumann problem (2). Evaluating the normal derivative of the potential (14) on $\Gamma$ and using, once again, the jump conditions (4), we obtain

[TABLE]

where the integral operators $N$ and $K^{\prime}$ are defined in (5c) and (5d), respectively. Therefore, selecting $\boldsymbol{x}_{0}=\boldsymbol{x}$ and using the fact that $\partial_{n}p_{0}(\boldsymbol{x}|\boldsymbol{x})=i\eta$ and $\partial_{n}p_{1}(\boldsymbol{x}|\boldsymbol{x})=0$, we obtain the Smoothed Neumann Combined Field Integral Equation (SN-CFIE)

[TABLE]

for the unknown density function $\varphi$.

Let us now examine the singular character of the integral operator $N\circ R_{D}-K^{\prime}\circ R_{S}$ in the SN-CFIE (25). Using the notations $\rho_{D}(\tau|t)$ and $\rho_{S}(\tau|t)$ introduced in (18) and (20), we have that $v(t)=(N\circ R_{D}-K^{\prime}\circ R_{S})[\varphi](\mathbb{x}(t))$ can be expressed as

[TABLE]

in terms of the integral kernels

[TABLE]

where the functions

[TABLE]

can be properly defined at $\tau=t$ so that they are $2\pi$-biperiodic analytic functions [12].

Therefore, from (27), Lemma 3.1 and the fact that $W(t,\tau)=O(|t-\tau|\log|t-\tau|)$, we obtain that

[TABLE]

as $\tau\to t$ for $\varphi\in H^{s}(\Gamma)$ for $s>5/2$.

As in the case of the SD-CFIE, it can be easily shown that for $s>5/2$ the integrands in (26) are continuously differentiable $2\pi$-biperiodic functions.

Note that, unfortunately, the diagonal values of the integrand in (26), i.e., the limit values of the integrand as $\tau\to t$, depends on the second derivative of $\phi(t)=\varphi(\mathbb{x}(t))$ at $\tau=t$. More precisely, we have

[TABLE]

where, as in the proof of Lemma 3.1, we have used the notations $\tilde{p}_{0}(\tau|t)=p_{0}(\mathbb{x}(\tau)|\mathbb{x}(t))$ and $\tilde{p}_{0}(\tau|t)=p_{1}(\mathbb{x}(\tau)|\mathbb{x}(t))$.

To illustrate the smoothness of the integrand in the expression (26) we present Figures 3 which displays the functions $H(t,\tau)\rho_{D}(\tau|t)$ and $W(t,\tau)\rho_{S}(\tau|t)$.

We end this section with the following assertion that follows directly from the discussion above:

Theorem 3.3**.**

Let $N$ and $K^{\prime}$ be the hypersingular and adjoint double-layer operators defined in (5a) and (5b), respectively, and $R_{D}$ and $R_{S}$ be the smoothing operators defined in (16). Then, the identity

[TABLE]

holds true on $H^{s}(\Gamma)$ for all $s>3/2$. Therefore, in particular, the SN-CFIE (15) admits at most one solution $\varphi\in H^{s}(\Gamma)$, $s>3/2$, for all wavenumbers $k\in\mathbb{C}$, $\mathrm{Re}\,k>0$, $\mathrm{Im}\,k\geq 0$, provided $\partial_{n}u^{\mathrm{inc}}\in H^{s-1}(\Gamma)$ and $\eta>0$.

Remark 3.4**.**

The hypersingular operator can be expressed as

[TABLE]

for all $\varphi\in H^{s}(\Gamma)$, $s>5/2$, and all $\boldsymbol{x}\in\Gamma$, in terms of the smoothing functions $p_{0}$ and $p_{1}$ introduced (11) and (12). It thus follows from the Lemma 3.1 and the smoothness of the kernel $W$ of the adjoint double-layer operator, that the hypersingular operator can be expressed in terms of integrals of continuous functions.

3.3 Evaluation of combined potentials close to the boundary

Once the density function $\varphi$ has been retrieved by solving the integral equation (15) or (25), depending on the boundary condition, the desired solution of the corresponding boundary value problem (1) or (2) is given in terms of the combined potential (8) that can be evaluated everywhere in $\mathbb{R}^{2}\setminus\overline{\Omega}$. From the definition of the single- and double-layer potentials (3) it is clear that for any given target point $\boldsymbol{r}\in\mathbb{R}^{2}\setminus\overline{\Omega}$ the integrands involved in the definition of the combined potential (8) are, in principle, smooth functions of the source/integration point $\boldsymbol{y}\in\Gamma$. In practice, however, numerical issues arise when the observation point $\boldsymbol{r}$ lies “near” the boundary (see discussion on the 5$h$-rule in [3, Remark 6]). In this case the boundary integrands in (8) are still smooth functions of $\boldsymbol{y}\in\Gamma$, but since both the Green function $i/4H_{0}^{(1)}(k|\boldsymbol{r}-\boldsymbol{y}|)$ and its normal derivative blow up as $O(\log|\boldsymbol{r}-\boldsymbol{y}|)$ and $O(|\boldsymbol{r}-\boldsymbol{y}|^{-1})$ respectively, when $\boldsymbol{r}\to\boldsymbol{y}\in\Gamma$, large numbers of quadrature points are needed to properly resolve the nearly singular character of the boundary integrands.

In this section we describe how to utilize the smoothing operators $R_{D}$ and $R_{S}$ introduced in Section 3 above, to substantially mitigate the errors produced by the naive numerical approximation of nearly singular integrals that arise from evaluation of the combined potentials close to the boundary.

Selecting $\boldsymbol{x}_{0}=\bar{\boldsymbol{r}}$ in (14), where $\bar{\boldsymbol{r}}\in\Gamma$ is such that $|\boldsymbol{r}-\bar{\boldsymbol{r}}|=\min_{\boldsymbol{x}\in\Gamma}|\boldsymbol{r}-\boldsymbol{x}|$, we obtain that the combined potential (8) can be expressed as

[TABLE]

It thus follows from the properties of the operators $R_{D}$ and $R_{S}$ established in Lemma 3.1 that for a density function $\varphi\in H^{s}(\Gamma),$ $s>5/2$, the integrands in the expression on the right-hand-side of (30) satisfy

[TABLE]

as $\boldsymbol{y},\boldsymbol{r}\to\bar{\boldsymbol{r}}$, $\boldsymbol{r}\in\mathbb{R}^{2}$, $\bar{\boldsymbol{r}},\boldsymbol{y}\in\Gamma$.

It is demonstrated in Section 5 through numerical examples that use of the smoothed potential (30) instead of the combined potential (3) improves considerably the numerical accuracy of the fields at target points near the boundary.

3.4 Close obstacles

An issue similar to the one described above in Section 3.3 arises in scattering configurations involving two or more obstacles that are close to each other.

Without loss of generality we let $\Omega$ be composed by two disjoint obstacles $\Omega_{1}$ and $\Omega_{2}$ with smooth boundaries $\Gamma_{1}$ and $\Gamma_{2}$ ($\Gamma=\Gamma_{1}\cup\Gamma_{2}$). Clearly, evaluation of any of the integral operators (5) on the curve $\Gamma_{1}$ entails integration on $\Gamma_{2}$ of a certain density function multiplied by the Green function or one of its normal derivatives with respect to the target point $\boldsymbol{x}\in\Gamma_{1}$, the source point $\boldsymbol{y}\in\Gamma_{2}$, or both. If $\Gamma_{1}$ and $\Gamma_{2}$ are close to each other, that is, there are source points $\boldsymbol{y}\in\Gamma_{2}$ that are close to target points $\boldsymbol{x}\in\Gamma_{1}$, the relevant integrands in the (5) become nearly singular as they portrait singularities of the form $O(\log|\boldsymbol{x}-\boldsymbol{y}|)$ in the case of the single-layer, $O(|\boldsymbol{x}-\boldsymbol{y}|^{-1})$ in the case of the double-layer and adjoint double layer operators, and $O(|\boldsymbol{x}-\boldsymbol{y}|^{-2})$ in the case of the hypersingular operator.

Following the ideas presented in Section 3.3 and letting $\varphi|_{\Gamma_{1}}=\varphi_{1}$ and $\varphi|_{\Gamma_{2}}=\varphi_{2}$ denote the restriction of the density function to each one of the curves, we have that the combined field integral operators at a point $\boldsymbol{x}\in\Gamma_{i}$, $i=1,2$, can be expressed as

[TABLE]

in the case of the D-CFIE, and

[TABLE]

in the case of the N-CFIE, where the point $\bar{\boldsymbol{x}}_{j}\in\Gamma_{j}$, $j=1,2$, $j\neq i$, is such that $|\boldsymbol{x}-\bar{\boldsymbol{x}}_{j}|=\min_{\boldsymbol{y}\in\Gamma_{j}}|\boldsymbol{x}-\boldsymbol{y}|$.

Therefore, it follows from Lemma 3.1 that all the integrands on the right hand side of (31) and (32) remain bounded regardless the distance between the curves $\Gamma_{1}$ and $\Gamma_{2}$. In fact, the integrands on $\Gamma_{j}$ in (31) satisfy

[TABLE]

while the integrands on $\Gamma_{j}$ in (32) satisfy

[TABLE]

as $\boldsymbol{y},\boldsymbol{x}\to\bar{\boldsymbol{x}}_{j}$.

4 Smoothing functions

This section presents explicit expressions for the smoothing functions $p_{0}(\boldsymbol{r}|\boldsymbol{x}_{0})$ and $p_{1}(\boldsymbol{r}|\boldsymbol{x}_{0})$ introduced in (11) and (12), respectively, in terms of linear combinations of planes waves (LCPW) of the form $\operatorname{e}^{ik\boldsymbol{d}\cdot(\boldsymbol{r}-\boldsymbol{x}_{0})}$, where $\boldsymbol{d}\in\mathbb{R}^{2}$ is a constant unit vector and where $\boldsymbol{x}_{0}$ is a given point on the curve $\Gamma$. The direction of propagation $\boldsymbol{d}$ of each one of the plane waves will be expressed as a linear combination of the vectors $\boldsymbol{n}_{0}$ and $\boldsymbol{\tau}_{0}$ which denote the unit normal and unit tangent vectors at $\boldsymbol{x}_{0}\in\Gamma$, respectively. More precisely, using the curve parametrization $\mathbb{x}:[0,2\pi)\to\Gamma$ we have $\boldsymbol{x}_{0}=\mathbb{x}(t_{0})$, $\boldsymbol{n}_{0}=(-x_{2}^{\prime}(t_{0}),x_{1}^{\prime}(t_{0}))/|\mathbb{x}^{\prime}(t_{0})|$ and $\boldsymbol{\tau}_{0}=\mathbb{x}^{\prime}(t_{0})/|\mathbb{x}^{\prime}(t_{0})|$ for some $t_{0}\in[0,2\pi)$. Here we note that any LCPW is indeed a smooth homogeneous solution of the Helmholtz equation in all of $\mathbb{R}^{2}$.

We begin then by constructing a LCPW $f_{0}$ that satisfies the point conditions

[TABLE]

From the law of reflections we have that a LCPF of the form $f_{0}(\boldsymbol{r}|\boldsymbol{x}_{0})=c\{\operatorname{e}^{ik\boldsymbol{n}_{0}\cdot(\boldsymbol{r}-\boldsymbol{x}_{0})}+\operatorname{e}^{-ik\boldsymbol{n}_{0}\cdot(\boldsymbol{r}-\boldsymbol{x}_{0})}\}$ satisfies the homogeneous Neumann boundary condition $\partial_{n}f_{0}=0$ on the line tangent to $\Gamma$ at the point $\boldsymbol{x}_{0}$, and, furthermore, it remains constant along the tangent direction $\boldsymbol{\tau}_{0}$. Therefore, enforcing the condition $f_{0}(\boldsymbol{x}_{0}|\boldsymbol{x}_{0})=1$ we obtain that the elementary function $f_{0}(\boldsymbol{r}|\boldsymbol{x}_{0})=\cos(k\boldsymbol{n}_{0}\cdot(\boldsymbol{r}-\boldsymbol{x}_{0}))$ satisfies all the required point conditions (33).

Similarly, considering now a LCPW $g_{0}$ that satisfies the homogeneous Dirichlet boundary condition $g_{0}=0$ on the line tangent to $\Gamma$ at the point $\boldsymbol{x}_{0}$, we find that $g_{0}(\boldsymbol{r}|\boldsymbol{x}_{0})=\sin(k\boldsymbol{n}_{0}\cdot(\boldsymbol{r}-\boldsymbol{x}_{0}))/k$ fulfills all the point conditions

[TABLE]

Combining $f_{0}$ and $g_{0}$ we thus obtain the expression

[TABLE]

for the smoothing function (11).

In order to construct the smoothing function $p_{1}$ that satisfies the conditions in (12), on the other hand, we consider the LCPW that produces the function

[TABLE]

which satisfies

[TABLE]

Therefore, in order to construct $p_{1}$ it suffices to provide a LCPW $f_{1}$ such that

[TABLE]

In order to do so, we first consider the LCPW that produces the function

[TABLE]

which satisfies the point conditions

[TABLE]

Therefore, clearly, a LCPW satisfying all the points conditions listed in (35) is given by

[TABLE]

Finally, combining $f_{1}$ and $g_{1}$ we obtain

[TABLE]

5 Numerical examples

This section presents numerical examples that illustrate the properties of the smoothed CFIEs (15) and (25) and as well as the smoothed potential (30).

5.1 Nyström discretizations

Three different Nyström methods for the numerical solution of the smoothed CFIEs (15) and (25) are briefly reviewed in this section. In order to obtain Nyström discretizations we need to provide quadrature rules for the numerical evaluation of the smoothed integral operators (22) and (26), which are here expressed as

[TABLE]

where the kernels $A$ and $B$ correspond to $A=L$ and $B=M$ in the case of the SD-CFIE (22), and $A=H$ and $B=W$ in the case of the SN-CFIE (26). We recall that the functions $\rho_{D}$ and $\rho_{S}$, introduced in (18) and  (20), respectively, are given in terms of the density function $\phi(t)=\varphi(\mathbb{x}(t))$.

The first and simplest Nyström method considered is based on the direct use of the classical trapezoidal rule (TR), which applied to the integral (37) yields

[TABLE]

where the quadrature weights and quadrature points are given by $w_{j}=h=\pi/n$ and $t_{j}=jh$, $j=0,\ldots,2n-1$, respectively. A linear system of equations for the approximate values of the density function $\phi_{j}\approx\phi(t_{j})=\varphi(\mathbb{x}(t_{j}))$, $j=0,\ldots,2n-1,$ is obtained by equating the right-hand-side of (38) to the corresponding Dirichlet $f_{D}(t)=-u^{\mathrm{inc}}(\mathbb{x}(t))$ or Neumann $f_{N}(t)=-\partial_{n}u^{\mathrm{inc}}(\mathbb{x}(t))$ data at the quadrature points $t=t_{i}$, $i=0,\ldots,2n-1$. Here, the first and second order derivatives of the unknown function $\phi(t)=\varphi(\mathbb{x}(t))$ at $t=t_{i}$—which are needed for evaluation of

[TABLE]

—are approximated by finite differences of $\phi_{j}$, $j=0,\ldots,2n-1$. Note that although the functions $A(t_{i},\tau)\rho_{D}(\tau|t_{i})$ and $B(t_{i},\tau)\rho_{S}(\tau|t_{i})$ in (38) are continuously differentiable, they still have a polylogarithmic singularity of the form $|t_{i}-\tau|^{2}\log|t_{i}-\tau|$ at $\tau=t_{i}$. As expected, this mild singularity limits to $O(h^{3})$ the overall accuracy of the trapezoidal rule and associated Nyström-method solution $\phi_{j}$, $j=0,\ldots,2n-1$, which would otherwise exhibit spectral (superalgebraic) accuracy if the integrands were smooth periodic functions.

To achieve higher order accuracy, of course, Nyström methods that properly handle logarithmic singularities can be utilized instead of the classical trapezoidal rule. Following [12] then, we consider the spectrally accurate Nyström method of Martensen and Kussmaul (MK) which is suited for the discretization of the smoothed integral equations (15) and (25) by virtue of the fact that the integrand in (37) can be expressed as

[TABLE]

where $C_{1}$ and $C_{2}$ are $2\pi$-biperiodic analytic functions. (The explicit form of the function $C_{1}$ and $C_{2}$ can be easily obtained from identities (23) and (27).) It thus follows from [12] that the quadrature rule

[TABLE]

with quadrature weights given by

[TABLE]

provides a spectrally accurate approximation of the smoothed integral operators at the quadrature points $t_{i}=ih$, $i=0,\ldots,2n-1$. Equating the right-hand-side of (41) to the corresponding boundary data $f_{D}$ or $f_{N}$ at the quadrature points $t_{i}$, $i=0,\ldots,2n-1$, we obtain a linear system of equations for the quantities $\phi_{j}$, $j=0,\ldots,2n-1,$ which approximate the desired unknown values $\phi(t_{j}),$ $j=0,\ldots,2n-1$, respectively. The spectral accuracy of the Nyström solution $\phi_{j}$, $j=0,\ldots,2n-1$, is achieved by approximating $\phi^{\prime}(t_{i})$ and $\phi^{\prime\prime}(t_{i})$ in (39c) using FFT-based differentiation of the periodic sequence $\phi_{j}$, $j=0,\ldots,2n-1$.

Alternatively, the spectrally accurate quadrature rule for integral kernels with $O(|t-\tau|^{2}\log|t-\tau|)$ singularities introduced in [6] could be utilized to produce high-order Nyström discretizations of the smoothed CFIEs.

Finally, following [19] we also consider the high-order Nyström method based on Kapur-Rokhlin quadrature rules. Utilizing the same quadrature points $t_{j}=j\pi/n$, $j=0,\ldots,2n-1$, this approach has an advantage over the aforementioned methods in that it does not entail evaluation of the rather involved diagonal term (39c)—as by construction the quadrature weight corresponding to the singular point $\tau=t$ is identically zero (see [19] for details).

5.2 Smooth obstacles

In our first numerical experiment we compare the accuracy of the various Nyström discretizations for the solution of the smoothed CFIEs (15) and (25) for the exterior boundary value problems (1) and (2), respectively, that result from the scattering of an incident plane-wave $u^{\mathrm{inc}}(\boldsymbol{x})=\operatorname{e}^{ik\boldsymbol{d}\cdot\boldsymbol{x}}$, $\boldsymbol{d}=(\cos\pi/8,\sin\pi/8)$, that impinges on a kite-shaped obstacle with boundary $\Gamma=\{(\cos t+0.65(\cos 2t-1),1.5\sin t)\in\mathbb{R}^{2},t\in[0,2\pi)\}.$ The parameter value $\eta=k$ is utilized in all the numerical examples considered in this paper.

Figure 4 displays the maximum relative errors

[TABLE]

in the approximate far-field pattern $\tilde{u}_{\infty}$ obtained from the numerical solution (15) and (25) for various grid sizes $h=\pi/n$. The numerical errors displayed in Figure 4 are measured with respect to a reference highly accurate far-field pattern $u_{\infty}$. Here we recall that the far-field pattern corresponding to an integral equation solution $\varphi$ is given by [12]

[TABLE]

As can be observed in Figure 4, the three Nyström discretizations yield the expected convergence rates which are in turn determined by the order of the associated quadrature rules. Interestingly, it can be noted that for a wide range of grid sizes $h$ the low-order TR method produces more accurate results than the high-order Nyström method based on the Kapur-Rokhlin quadrature rules or orders six (KR-6) and ten (KR-10).

Table 1, on the other hand, displays the number of GMRES iterations required to solve the linear systems resulting from the three Nyström discretizations of the smoothed and non-smoothed integral equations. These results show that the numbers of GMRES iterations required to solve the smoothed integral equations (SD-CFIE (15) and SN-CFIE (25)) do not differ considerably from those required to solve the non-smoothed integral equations (D-CFIE (9) and N-CFIE (10)—these are, of course, expected results in view of the operator identities established in Theorems 3.2 and 3.3. Here we point out that the inordinate large number of GMRES iteration required by the KR-10 method makes it in practice not amenable to iterative linear algebra solvers [19].

Our second numerical experiment illustrates the properties of the smoothed combined potential (30) when evaluated at target points close to boundaries. In this experiment we consider the exterior Dirichlet problem (1) with boundary data $u^{\mathrm{inc}}(\boldsymbol{x})=-H_{0}^{(1)}(k|\boldsymbol{x}|)$ on $\Gamma$ so that the exact solution of (1) is $u^{\mathrm{exact}}(\boldsymbol{x})=H_{0}^{(1)}(k|\boldsymbol{x}|)$ in $\mathbb{R}^{2}\setminus\overline{\Omega}$. Four different approximate solutions $\tilde{u}$ of (1) are obtained by solving the smoothed integral equation (15) by means of the MK and TR methods, and then utilizing the resulting density functions to produce the near fields via the combined and smoothed-combined potentials given in (8) and (30), respectively. Figure 5 displays the logarithm (in base 10) of the absolute errors $|\tilde{u}(\boldsymbol{x})-u^{\mathrm{exact}}(\boldsymbol{x})|$ within a region near the boundary $\Gamma$ of the kite-shaped obstacle. This figure shows that, for the given discretization level, the smoothed potential produces near fields that are at least three digits more accurate than those obtained through the non-smoothed combined potential, in both MK and TR cases.

In the final numerical experiment of this section we illustrate the advantages of the smoothed operators (31) and (32), presented in Section 3.4, for the solution problems involving obstacles that are close to each other. This experiment considers a geometric configuration consisting of two kite-shaped obstacles with boundaries $\Gamma_{1}=\{(\cos t+0.65(\cos 2t-1)-1-d/2,1.5\sin t)\in\mathbb{R}^{2},t\in[0,2\pi)\}$ and $\Gamma_{2}=\{(-\cos t-0.65(\cos 2t-1)+1+d/2,1.5\sin t)\in\mathbb{R}^{2},t\in[0,2\pi)\}$ separated by a distance $d>0$. In order to assess the numerical errors we construct exact solutions to both Dirichlet (1) and Neumann (2) problems by imposing the boundary conditions $u^{\mathrm{inc}}(\boldsymbol{x})=-(H_{0}^{(1)}(k|\boldsymbol{x}-\boldsymbol{x}_{1}|)+H_{0}^{(1)}(k|\boldsymbol{x}-\boldsymbol{x}_{2}|))$ and $\partial_{n}u^{\mathrm{inc}}=-\partial_{n}(H_{0}^{(1)}(k|\boldsymbol{x}-\boldsymbol{x}_{1}|)+H_{0}^{(1)}(k|\boldsymbol{x}-\boldsymbol{x}_{2}|))$ on $\Gamma_{1}\cup\Gamma_{2}$, respectively, where $\boldsymbol{x}_{1}=-(1+d/2,0)\in\Gamma_{1}$ and $\boldsymbol{x}_{2}=(1+d/2,0)\in\Gamma_{2}$ . The exact solution to both problems (1) and (2) is then given by $u_{\mathrm{exact}}(\boldsymbol{x})=H_{0}^{(1)}(k|\boldsymbol{x}-\boldsymbol{x}_{1}|)+H_{0}^{(1)}(k|\boldsymbol{x}-\boldsymbol{x}_{2}|)$.

Figure 6 displays the far-field errors resulting from application of the MK and TR methods to the Dirchlet and Neumann integral equations with smoothed integral operators (31) and (32), respectively, and the far-field errors resulting from application of the MK method to the non-smoothed integral equations D-CFIE (9) and N-CFIE (10). Figure 6a, in particular, displays the errors for various separations $d>0$ and a fixed grid size $h=\pi/32$. This figure reveals the pronounced accuracy deterioration of the MK solutions (blue and red curves) as $d\to 0$ in both Dirichlet (left) and Neumann (right) cases. This deterioration is significantly milder in the case of the smoothed integral equations (red curves). The accuracy of the TR solutions (yellow curves), in turn, does not seem to deteriorate $d\to 0$. Figure 6b, on the other hand, displays the far-field errors for a fixed separation $d=10^{-5}$ and various grid sizes $h>0$. This figure demonstrates that the use the smoothed integral equation (red and yellow curves) improves significantly the accuracy of the far-field in the Neumann case (right). In the Dirichlet case (left), however, the far-fields resulting from smoothed and non-smoothed integral equations exhibit similar convergence rates and accuracies.

Figures 7 and 8, finally, show the logarithm (in base 10) of the absolute errors $|\tilde{u}(\boldsymbol{x})-u^{\mathrm{exact}}(\boldsymbol{x})|$ for the Dirichlet and Neumann problems, respectively, where $\tilde{u}$ denotes any of the approximate solutions obtained by means of the smoothed and non-smoothed integral equations using the smoothed and non-smoothed combined potentials.

5.3 Obstacles with corners

In this section we extend the smoothing procedure to problems of scattering involving obstacles with corners. For presentation simplicity we assume that the boundary $\Gamma$, which is parametrized by $\mathbb{x}:[0,2\pi]\to\Gamma$, has only one corner at the point $\mathbb{x}(0)=\mathbb{x}(2\pi)$.

In order to deal with the corner singularity of the integral equations solutions (cf. [10, 18, 38]) we introduce graded meshes generated by use of the change of variable $t=w(s)$, where letting $p\geq 2$ and

[TABLE]

the function $w:[0,2\pi]\to[0,2\pi]$ is given by

[TABLE]

Note that $w$ is a smooth and monotonically increasing function on the interval $[0,2\pi]$ and its derivatives vanish algebraically at $s=0$ and $s=2\pi$, i.e., $w^{(q)}(0)=w^{(q)}(2\pi)=0$ for $1\leq q\leq p-1$. This change of variable, which was originally introduced by Kress [23], has been extensively utilized to produce high-order Nyström discretizations of boundary integral equations involving domains with corners [1, 15, 37].

Using this transformation then, the integral (37) is expressed as

[TABLE]

Since the function $w^{\prime}(\sigma)$ in (43) vanishes algebraically at $\sigma=0$ and $\sigma=2\pi$, the integrand in (43) can be regarded as a $2\pi$-periodic function which exhibits a certain degree of smoothness around the endpoints $\sigma=0$ and $\sigma=2\pi$ that can be controlled by the parameter $p$. Therefore, for $p$ sufficiently large the smoothed integral equations (15) and (25) resulting from use of this change of variable can be discretized using any of the Nyström methods discussed in the previous section.

Utilizing the quadrature points $s_{j}=\pi/n(j+1/2)$, $j=0,\ldots,2n-1$—which do not include the endpoints 0 and $2\pi$—we obtain a linear system for the approximate the values $\psi_{j}$, $j=0,\ldots,2n-1,$ of the unknown function $\psi(s)=\phi(w(s))$ at the quadrature points $s=s_{j}$, $j=0,\ldots,2n-1$, respectively.

Clearly, the Nyström discretizations of the smoothed integral equations presented in the previous section require approximate expressions for

[TABLE]

and

[TABLE]

in terms of the discrete values $\psi_{j}$, $j=0,\ldots,2n-1$. (Note that $w^{\prime}(s_{i})\neq 0$ for all $i=0,\ldots,2n-1$.) These expressions are here obtained by approximating $\psi^{\prime}(s_{i})$ using finite differences of $\psi_{j}$ in the case of the low-order TR method, and using FFT-based differentiation of $\psi_{j}$ in the case of the higher order MK and KR methods. High accuracy is achieved in the latter case by utilizing the identity

[TABLE]

which allows the desired quantity $\psi^{\prime}(s_{i})$ to be approximated using FFF-based differentiation of the sequences $w^{\prime}(s_{j})\psi_{j}$ and $w^{\prime}(s_{j})$, $j=0,\ldots,2n-1$, which are assumed periodic. We proceed similarly to approximate the second-order derivative $\psi^{\prime\prime}(s_{i})$ that is needed for evaluation of the diagonal term $H(w(s_{i}),w(s_{i}))\rho_{D}(w(s_{i})|w(s_{i}))$.

Figure 9 display the far-field errors in the solution of the problem of scattering of a horizontal plane-wave $u^{\mathrm{inc}}(\boldsymbol{x})=\operatorname{e}^{ik(x\cos\alpha+y\sin\alpha)}$, $\alpha=0$, $k=4$, by a drop-shaped obstacle parametrized by $\mathbb{x}(t)=(2\sin\frac{t}{2},-\sin t)$ [23] that features a convex corner (Figures 9a and 9b), and by a boomerang-shaped obstacle parametrized by $\mathbb{x}(t)=(-\frac{2}{3}\sin\frac{3t}{2},-\sin t)$ [1] that features a concave (reentry) corner (Figures 9c and 9d). The smoothed integral equations (SD-CFIE (15) and SN-CFIE (25)) are discretized using the TR, MK and KR-10 Nyström methods described in the previous section, using the change of variable described above with the parameter value $p=4$. The error curves displayed in these figures demonstrate the accuracy of the three Nyström method considered. The numbers of the GMRES iterations required for the solution of the resulting linear systems are comparable to those reported in Table 1.

Figure 10 (resp. Figure 11), finally, displays the real part of total field (incident field plus scattered field) for the solution of the problem of scattering of a horizontal high-frequency ($k=64$) plane-wave by a drop-shaped (resp. boomerang-shaped) obstacle with Dirichlet and Nuemann boundary conditions. The integral equation solutions were obtained via the MK Nyström method and the near-field were produced through the smoothed potential (30).

Acknowledgement

The author would like to thank Prof. Catalin Turc with the Department of Mathematics of the New Jersey Institute of Technology, for interesting discussion and useful comments on the subject of this manuscript.