High order-accurate solution of scattering integral equations
with unbounded solutions at corners

Using the single-layer representation

$$u(r)=\int_{\Gamma}{G_{k}(r,r^{\prime})\phi d\ell^{\prime}},\quad r\in\Omega^{c},$$

(3)

where

$$G_{k}(r,r^{\prime})=\frac{i}{4}H_{0}^{1}(k|r-r^{\prime}|)$$

(4)

denotes the outgoing Helmholtz free-space Green’s function, the imposition of the Neumann boundary condition (2) results in the Magnetic Field Integral Equation (MFIE) [16]

$$-\frac{1}{2}\phi(r)+\frac{\partial}{\partial n(r)}\int_{\Gamma}G_{k}(r,r^{\prime})\phi(r^{\prime})d\ell(r^{\prime})=-\frac{\partial u^{\mathrm{inc}}}{\partial n(r)},\quad r\in\Gamma.$$

(5)

Once a solution $\phi$ of this equation has been produced and substituted into (3), the solution $u$ of (1)-(2) can be obtained by straightforward numerical quadrature.

The MFIE (5) is not uniquely solvable at a scatterer-dependent countable sequence of frequencies $k=k_{1},k_{2},\dots$, see e.g. [16], around each one of which the associated numerical implementations become ill-conditioned and inaccurate. This problem can be eliminated by utilizing instead a certain combined field integral equation formulation, which is considered in Section 2.2. Unfortunately, however, that formulation incorporates a hypersingular operator whose practical implementation presents certain difficulties, most notably in terms of the associated conditioning characteristics and large iteration numbers it requires if iterative solvers are used for its solution. Following [11], Section 3.1 introduces a second version of the combined equation, where the hypersingular operator is regularized using a smoothing operator. This approach resolves the difficulties associated with the hypersingular nature of the problem, at least for smooth surfaces $\Gamma$. However, new complications emerge when this formulation is applied to domains with corners, since, as is the case for the regular CFIE formulation, its solutions $\phi$ are unbounded at corner points. An additional reformulation of the problem, that introduces the new corner-regularized unknown and associated equation, which is one of the key contributions in this paper, is then presented in Section 3.2. Sections 3.3 and 3.4 detail additional components of the proposed solver designed to address finite precision issues associated with the use of fine singularity-resolving meshes and fine-scale pre-computation problems. The remainder of the present Section 2, in turn, introduces the aforementioned smooth-surface uniquely solvable CFIE equation and associated parametrization and differentiation methods, which are needed for reference in the subsequent Section 3.

Using now the representation

$$u=-i\eta\int_{\Gamma}{G_{k}(r,r^{\prime})\phi d\ell^{\prime}}+\int_{\Gamma}{\frac{\partial G_{k}(r,r^{\prime})}{\partial{n}(r^{\prime})}\phi d\ell^{\prime}},\quad r\in\Omega^{c},$$

(6)

instead of (3) (see [11] and references therein), enforcement of the Neumann boundary conditions (2) yields the Combined Field Integral Equation (CFIE) formulation

$$\frac{i\eta}{2}\phi(r)-i\eta\int_{\Gamma}\frac{\partial G_{k}(r,r^{\prime})}{\partial n(r)}\phi(r^{\prime})d\ell(r^{\prime})+\frac{\partial}{\partial n(r)}\int_{\Gamma}\frac{\partial G_{k}(r,r^{\prime})}{\partial n(r^{\prime})}\phi(r^{\prime})d\ell(r^{\prime})=-\frac{\partial u^{\mathrm{inc}}}{\partial n(r)},$$

(7)

where $\eta\neq 0$ is a coupling parameter. Following  [11], throughout this paper the value $\eta=1$ is used. In equation  (7) and throughout this paper, $\frac{\partial}{\partial n(r)}$ (resp. $\frac{\partial}{\partial n(r^{\prime})}$) denotes differentiation with respect to $r$ in the direction of $n(r)$ (resp. differentiation with respect to $r^{\prime}$ in the direction of $n(r^{\prime})$). The last operator on the left-hand side of (7), namely, the normal derivative of the double-layer operator, is a hypersingular operator—whose computation by direct implementation of normal derivatives can be rather challenging. Fortunately, however, the hypersingular operator may be expressed in the form

$$\frac{\partial}{\partial n(r)}\int_{\Gamma}\frac{\partial G_{k}(r,r^{\prime})}{\partial n(r^{\prime})}\phi(r^{\prime})d\ell(r^{\prime})={k}^{2}\int_{\Gamma}G_{k}(r,r^{\prime})(n(r)\cdot n(r^{\prime}))\phi(r^{\prime})\,d\ell(r^{\prime})+\partial_{\tau(r)}\int_{\Gamma}G_{k}(r,r^{\prime})\partial_{\tau(r^{\prime})}\phi(r^{\prime})\,d\ell(r^{\prime});$$

(8)

(see [11] and references therein), where $\partial_{\tau(r)}$ and $\partial_{\tau(r^{\prime})}$ denote the operators of differentiation with respect to the unit tangent vectors $\tau(r)$ and $\tau(r^{\prime})$ at the respective “target” and “source” points $r$ and $r^{\prime}$ on $\Gamma$; clearly equation (8) does not require evaluation of derivatives in directions other than those tangential to $\Gamma$, and can therefore be computed effectively by means of the discretization methods used in this paper. (Note that the identity (8) holds for closed curves, provided the tangential vector $\tau(r)$ is defined in a consistent manner throughout the scatterer’s boundary, e.g., in counterclockwise fashion.) In sum, defining the operators

$$T_{1}[\phi](r)=\int_{\Gamma}\frac{\partial G_{k}(r,r^{\prime})}{\partial n(r)}\phi(r^{\prime})d\ell(r^{\prime}),$$

(9)

$$T_{2}[\phi](r)=\int_{\Gamma}G_{k}(r,r^{\prime})(n(r)\cdot n(r^{\prime}))\phi(r^{\prime})\,d\ell(r^{\prime})\quad\mbox{and}$$

(10)

$$T_{3}[\phi](r)=\partial_{\tau(r)}\int_{\Gamma}G_{k}(r,r^{\prime})\partial_{\tau(r^{\prime})}\phi(r^{\prime})\,d\ell(r^{\prime}),$$

(11)

the CFIE equation (7) takes the form

$$\frac{i\eta}{2}\phi(r)-i\eta T_{1}[\phi](r)+k^{2}T_{2}[\phi](r)+T_{3}[\phi](r)=-\frac{\partial u^{\mathrm{inc}}}{\partial n(r)}.$$

(12)

Note that the integrals required by the operators $T_{1}$ and $T_{2}$ may be expressed in the forms

$$\mathcal{I}[H,f](r)=\int_{\Gamma}H(r,r^{\prime})f(r^{\prime})d\ell(r^{\prime}),$$

(13)

while $T_{3}$ may be presented in the form

$$\partial_{\tau(r)}\mathcal{I}[H,f](r)=\partial_{\tau(r)}\int_{\Gamma}H(r,r^{\prime})f(r^{\prime})d\ell(r^{\prime}),$$

(14)

for certain kernels $H$ and functions $f$: we have $H(r,r^{\prime})=\frac{\partial G_{k}(r,r^{\prime})}{\partial n(r)}$ and $f(r^{\prime})=\phi(r^{\prime})$ for the operator $T_{1}$; $H(r,r^{\prime})=G_{k}(r,r^{\prime})(n(r)\cdot n(r^{\prime}))$ and $f(r^{\prime})=\phi(r^{\prime})$ for the operator $T_{2}$; and $H(r,r^{\prime})=G_{k}(r,r^{\prime})$ and $f(r^{\prime})=\partial_{\tau(r^{\prime})}\phi(r^{\prime})$ for the opeator $T_{3}$. Clearly, the operator $T_{3}$ requires an additional differentiation of the unknown $\phi$—which, as discussed in. [11], impacts significantly upon the accuracy and conditioning of the resulting numerical implementation.

In preparation for the introduction of the main operator- and corner-regularization methods in Section 3, the following two sections describe, in the context of smooth surfaces, the particular parametrization and discretization methods used, and the associated methods employed for differentiation of surface variables and evaluation of Green function-based integral operators.

Throughout this paper the discretizations of integral operators such as those considered in the previous sections and elsewhere are implemented on the basis of of the Chebyshev-based Rectangular Polar (RP) Nyström methods introduced in [12, 28]. This section briefly reviews the RP method as it applies in the context of smooth surfaces and for operators that do not contain tangential differentiation (e.g., the operators (9) and (10) above); additions presented in Section 2.4 then enable the discretization of the operator (11) with similar efficiency and accuracy.

The RP method relies on use of a partition of curve (surface in 3D) $\Gamma$ into a set a set of $M$ non-overlapping patches $\Gamma^{q}$

$$\Gamma=\bigcup_{q=1}^{M}\Gamma^{q},$$

(15)

as depicted in Figure 1. Without loss of generality it is assumed that each curve $\Gamma^{q}$ is parametrized by a vector function with scalar argument

$$r_{q}:[0,1]\to\Gamma^{q}\in\mathbb{R}^{2}.$$

(16)

(Note that, while previous implementations of the RP method use patch parametrizations defined in the interval $[-1,1]$—given the reliance of the RP method on Chebyshev expansions—, throughout this paper the parameter segment $[0,1]$ is utilized instead, on account of the “zero-centering” technique introduced per Sections 3.2 and 3.3, which requires the use of such a parametrization interval.)

On the basis of the parametrizations $r_{q}$, a function $f(r)$ defined on $\Gamma^{q}$ corresponds to the function $f(r_{q}(s))$ defined for $s\in[0,1]$, which can then be discretized by means of the $Q$-point Chebyshev grid given by [8]

$$s_{j}=\frac{\zeta_{j}+1}{2}\quad j=0,\ldots,Q-1;\quad\omega_{j}=\frac{W_{j}}{2},$$

(17)

where

$$\zeta_{j}=\cos\left(\pi\frac{2j+1}{2Q}\right),\quad j=0,\ldots,Q-1;\quad W_{j}=\frac{2}{Q}\left(1-2\sum_{k=1}^{Q/2}\frac{1}{4k^{2}-1}\cos\left(k\pi\frac{2j+1}{Q}\right)\right).$$

(18)

Clearly, using $Q$ points per patch for a total $M$, the total number of discretization points used is

$$N=MQ.$$

(19)

The implementation of the integral operators at hand requires consideration of the line element

$$d\ell_{q^{\prime}}=L_{q^{\prime}}(s^{\prime})ds,\quad\mbox{where}\quad L_{q^{\prime}}(s^{\prime})\mbox{ = }\sqrt{\left(\frac{dx_{q^{\prime}}}{ds}(s^{\prime})\right)^{2}+\left(\frac{dy_{q^{\prime}}}{ds}(s^{\prime})\right)^{2}},$$

(20)

as well as the unit tangential and normal basis vectors which, at a point $r_{q^{\prime}}(s^{\prime})=\left(x_{q^{\prime}}(s^{\prime}),y_{q^{\prime}}(s^{\prime})\right)\in\Gamma^{q^{\prime}}$, are given by

$${\tau}(r_{q^{\prime}}(s^{\prime}))=\left(\frac{dx_{q^{\prime}}}{ds}(s^{\prime}),\frac{dy_{q^{\prime}}}{ds}(s^{\prime})\right)\bigg{/}L_{q^{\prime}}(s^{\prime}),\quad n(r_{q^{\prime}}(s^{\prime}))=\left(\frac{dx_{q^{\prime}}}{ds}(s^{\prime}),\frac{dy_{q^{\prime}}}{ds}(s^{\prime})\right)\bigg{/}L_{q^{\prime}}(s^{\prime})$$

(21)

In order to achieve high computational efficiency and accuracy the RP quadrature algorithm proceeds by separately treating singular, near-singular, and regular interactions. In view of (13) and (15), the integral $\mathcal{I}[H,f]$ may be expressed in the form

$$\mathcal{I}[H,f](r)=\sum_{{q^{\prime}}=1}^{M}\mathcal{I}^{q^{\prime}}[H,f](r),$$

(22)

where

$$\mathcal{I}^{q^{\prime}}[H,f](r)=\int_{\Gamma^{q^{\prime}}}H(r,r^{\prime})f(r^{\prime})d\ell(r^{\prime})=\int_{0}^{1}H(r,r_{q^{\prime}}(s^{\prime}))f(r_{q^{\prime}}(s^{\prime}))L_{q^{\prime}}(s^{\prime})ds^{\prime}.\\$$

(23)

For target points $r$ away from the source patch $\Gamma^{q^{\prime}}$ (as determined by a certain proximity parameter $\delta^{\mathrm{prox}}>0$) the kernels of the integral operators are smooth over the source patch, and, at least in the smooth-surface case considered in this section, so are the density functions $f$; thus, for such far interactions (on smooth domains) the method directly utilizes Fejér’s first quadrature rule for integration, and the approximation

$$\mathcal{I}^{{q^{\prime}}}[H,f](r)\approx\sum_{j=0}^{N-1}H(r,r_{{q^{\prime}}}(s_{j}))f(r_{{q^{\prime}}}(s_{j}))L_{{q^{\prime}}}(s_{j})\omega_{j},\quad\mbox{for }\quad r\quad\mbox{away from}\quad\Gamma^{{q^{\prime}}}$$

(24)

results, where $s_{j}$ and $\omega_{j}$ denote the integration points and weights introduced in equation (17).

When the target point $r$ is close to or within the source patch, on the other hand, the kernels of the integral operators are near singular or singular, respectively, and accurate evaluation of the corresponding integrals requires use of specialized integration techniques. To produce such integrals, the RP method relies on the Chebyshev approximation

$$f(r_{q^{\prime}}(s^{\prime}))\approx\sum_{m=0}^{N-1}a_{m}^{q^{\prime}}T_{m}(2s^{\prime}-1)$$

(25)

which can be obtained for the necessary functions $f$ (equal either to $\phi$ or its tangential derivative, see equation (14)) from the available values of the density $\phi$ on the Chebyshev mesh. Using this expansion, the operator values in singular and near-singular cases can be efficiently and accurately computed on the basis of the numerically precomputed integration weights

$$\beta_{m}^{q^{\prime}}[H](r)=\int_{0}^{1}H(r,r_{q^{\prime}}(s^{\prime}))\,L_{q^{\prime}}(s^{\prime})\,T_{m}(2s^{\prime}-1)\,ds^{\prime},$$

(26)

These weights can be precomputed by integrating the kernel against each Chebyshev polynomial (e.g., as described in [12, 28]) and stored for repeated use. Using these weights we obtain the accurate approximation

$$\mathcal{I}^{q^{\prime}}[H,f](r)\approx\sum_{m=0}^{N-1}a_{m}^{q^{\prime}}\beta_{m}^{q^{\prime}}[H](r),\quad\mbox{for}\quad r\quad\mbox{on or near}\quad\Gamma^{q^{\prime}}.$$

(27)

This approach applies to each one of the kernels $H$ listed near the end of Section 2.2 and, together with the aforementioned smooth integration methods for $r$ far from $\Gamma^{q^{\prime}}$, it enables the necessary integrals over $\Gamma$ to be computed efficiently via equation (22). Note that, per the discussion above in this section, precomputations are necessary for target discretization points $r=r_{q^{\prime}}(s_{j})$ on the source patch $\Gamma^{q^{\prime}}$, as well as for near singular target discretization points $r=r_{q^{\prime}}(s_{j})$ on the nearby patches $\Gamma^{q^{\prime}}$ for which the distance between the discretization points $r=r_{q^{\prime}}(s_{j})$ and $\Gamma_{q^{\prime}}$ is less than $\delta^{\mathrm{prox}}$.

The evaluation of the operator $T_{3}$ in (11), which contains two tangential derivatives, is discussed in what follows. In view of (15) and writing

$$\left.\partial_{\tau(r)}\int_{\Gamma^{q^{\prime}}}G_{k}(r,r^{\prime})\partial_{\tau(r^{\prime})}\phi(r^{\prime})\,d\ell(r^{\prime})\right|_{r=r_{q}(s)}=\frac{1}{L_{q}(s)}\frac{\partial}{\partial s}\int_{0}^{1}G_{k}(r_{q}(s),r_{q^{\prime}}(s^{\prime}))\frac{\partial\phi(r_{q^{\prime}}(s^{\prime}))}{\partial s^{\prime}}ds^{\prime},$$

(28)

(since $\partial_{\tau(r_{q}(s))}=\frac{1}{L_{q}(s)}\partial/\partial s$ and $\partial_{\tau(r_{q}^{\prime}(s^{\prime}))}=\frac{1}{L_{q^{\prime}}(s^{\prime})}\partial/\partial s^{\prime}$), it suffices to evaluate the right-hand side expression in this equation for all $q$ and $q^{\prime}$ with $1\leq q,q^{\prime}\leq M$ and then sum over $q^{\prime}$ for each $q$. To do this we leverage the Chebyshev structure inherent in the RP algorithm described in Section 2.3 together with the numerically stable Chebyshev differentiation approach [8] (which is briefly reviewed below), all of which results in the following procedure.

1. 1.

   For each $q^{\prime}$ evaluate the derivative $\left.\frac{\partial\phi(r_{q^{\prime}}(s^{\prime}))}{\partial s^{\prime}}\right|_{s^{\prime}=s_{j}}$ for $j=0,\dots,Q-1$ using Chebyshev differentiation.

2. 2.

   For each $q$ and each $q^{\prime}$ evaluate the integral $S^{q^{\prime}}(r_{q}(s))=\int_{0}^{1}G_{k}(r_{q}(s),r_{q^{\prime}}(s^{\prime}))\frac{\partial\phi(r_{q^{\prime}}(s^{\prime}))}{\partial s^{\prime}}ds^{\prime}$ by means of the RP methods described in Section 2.3 with $H(r,r^{\prime})=G_{k}(r,r^{\prime})$ and on the basis of the quantity $f(r_{q^{\prime}}(s_{j}))=\left.\frac{\partial\phi(r_{q^{\prime}}(s^{\prime}))}{\partial s^{\prime}}\right|_{s^{\prime}=s_{j}}$ obtained per point 1, and then obtain the quantity $S(r_{q}(s))=\sum_{q=1}^{M}S^{q^{\prime}}(r_{q}(s))$.

3. 3.

   Calculate $T_{3}[\phi](r_{q}(s))=\frac{1}{L_{q}(s)}\frac{\partial S(r_{q}(s))}{\partial s}$ using once again Chebyshev differentiation.

The aforementioned Chebyshev-based differentiation algorithm [8] numerically obtains the derivative of a given function $g(s)$ defined in the interval $[0,1]$ as follows:

• (a)

  Obtain the Chebyshev expansion $g(s)=\sum_{i=0}^{Q-1}c_{i}T_{i}(2s-1)$ of order $Q-1$ (throughout this paper we have used the value $Q=10$ for differentiation purposes);

• (b)

  Obtain the derivative coefficients $c^{\prime}_{i}$ by means of the recurrence relation

  $$c_{i-1}^{\prime}=c_{i+1}^{\prime}+2ic_{i}\quad(i=Q-1,Q-2,\ldots,1);\quad c_{Q}^{\prime}=c_{Q-1}^{\prime}=0,$$

• (c)

  Evaluate the desired derivative values via the inverse Chebyshev transform $g^{\prime}(s)=\sum_{i=0}^{Q-1}c^{\prime}_{i}T_{i}(2s-1)$.

Computational experiments presented in [11] illustrate the fact that, while the CFIE formulation is uniquely solvable for all wavenumbers, solving it accurately can give rise to certain difficulties. The presence of two tangential derivatives in the hypersingular operator (8) associated with the formulation (7) induces arbitrarily large eigenvalues as the discretizations are refined. As a result, the conditioning of the system suffers and the numbers of iterations for convergence to a required tolerance increases significantly. As shown in [11] these difficulties can be effectively resolved, for smooth surfaces $\Gamma$, by representing the scattered fields in the alternative form

$$u=-i\eta\int_{\Gamma}{G_{k}(r,r^{\prime})\phi(r^{\prime})d\ell^{\prime}}+\int_{\Gamma}{\frac{\partial G_{k}(r,r^{\prime})}{\partial{n}(r^{\prime})}R[\phi(r^{\prime})]d\ell^{\prime}}$$

(29)

where

$$R[\phi](r^{\prime})=\int_{\Gamma}{G_{K}(r^{\prime},r^{\prime\prime})\phi(r^{\prime\prime})d\ell(r^{\prime\prime})}$$

(30)

is a regularizing operator. The wavenumber $K$ can take any real or imaginary values and does not necessarily need to depend on the wavenumeber of the problem $k$, but [11] recommends the use of the value

$$K=ik.$$

Enforcement of the Neumann condition (2) then gives rise to the (Operator) Regularized Combined Field Integral Equation (CFIE-R)

$$\frac{i\eta}{2}\phi(r)-i\eta\frac{\partial}{\partial n(r)}\int_{\Gamma}G_{k}(r,r^{\prime})\phi(r^{\prime})\,d\ell(r^{\prime})+k^{2}\int_{\Gamma}G_{k}(r,r^{\prime})(n(r)\cdot n(r^{\prime}))\left(\int_{\Gamma}{G_{K}(r^{\prime},r^{\prime\prime})\phi(r^{\prime\prime})d\ell(r^{\prime\prime})}\right)\,d\ell(r^{\prime})\\ +\frac{\partial}{\partial\tau(r)}\int_{\Gamma}G_{k}(r,r^{\prime})\frac{\partial}{\partial\tau(r^{\prime})}\left(\int_{\Gamma}{G_{K}(r^{\prime},r^{\prime\prime})\phi(r^{\prime\prime})d\ell(r^{\prime\prime})}\right)\ d\ell(r^{\prime})=-\frac{\partial u^{\mathrm{inc}}(r)}{\partial n(r)}.$$

(31)

As illustrated in [11], at least for smooth surfaces, this formulation leads to well conditioned numerical methods which additionally, if used in conjunction with iterative solvers, require reduced iteration numbers. But the equation suffers from severe accuracy and conditioning problems for surfaces $\Gamma$ containing corners. A “Corner-Regularized” discretization is introduced in Section 3.2; unlike the operator regularization introduced in the present section, the corner regularization is introduced on the basis of the parameterizations of the various smooth curves that make up the piecewise-smooth scatterer $\Gamma$. Since the RP discretization method [12, 28] utilized in this paper otherwise relies on use of certain parametrizations and discretizations of $\Gamma$ in terms of cosine transforms and Chebyshev polynomials, for conciseness Section 3.2 introduces the corner-regularization method on the basis of such parametrizations. We should note, however, that, (i) The corner regularization method presented in Section 3.2 could also be implemented in connection with other Nyström or Boundary-Element discretizations; and, (ii) It is expected that, in conjunction with the electromagnetic operator-regularized equations presented in [10] (3D Maxwell analogous of the ones described in Section 2.2), the methods inherent in the corner-regularized version of of the RP method presented in Section 3.2 should provide an excellent basis for a 3D-Maxwell corner-and-edge regularization methodology.

The proposed corner-regularization approach partitions the domain into a set of $M$ non-overlapping patches, parameterized by the functions $r_{q}=(x_{q}(s),y_{q}(s))$, $q=1,2,\ldots,M$, each one of which contains either one corner or no corners, and it is assumed that, by construction, if the $q$-th patch contains a corner then the corner point equals $r_{q}(0)$—that is, the corner point corresponds to the parameter value $s=0$—a feature that is utilized to eliminate catastrophic cancellation errors in the context of Section 3.3. Note that, in view of this prescription, the tangential vectors obtained by direct differentiation of the parametrization $r_{q}(s)$ on some of the patches $\Gamma^{q}$ need to be flipped (multiplied by $-1$) in order to obtain a consistent prescription (say, counterclockwise) of the tangential vector field around $\Gamma$.

The regularization method used seeks to yield an integral formulation whose unknown density is bounded and sufficiently smooth at corners. To do this, the algorithm relies on the use of a 1-1 change of variables (CoV) $s=s(\theta)$, $s:[0,1]\to[0,1]$ for each patch containing a corner [16, 31, 25]. The corner CoV gives rise to an integration Jacobian that vanishes at the corner together with its derivatives up to a prescribed order, and it maps a uniform mesh in the $\theta$ variable to a mesh in the $s$ variable that is graded toward the corner. The CoV Jacobians, which appear in the integrand multiplying the unknown density, are such that the product of Jacobian and density has a number of bounded derivatives at the corner point—thus resulting in a “regularized” quantity which can be approximated closely by the Chebyshev expansions used in the RP method.

Regularizing CoVs have previously been used for regularization of bounded (Lipschitz continuous) integral densities [16, 31, 25]; the regularization strategy embodied in equations (41) and (44) below, which is applicable to integral equations whose density solutions are unbounded, is one of the main enabling elements introduced in this paper. Other essential algorithmic elements in the proposed approach are introduced in Sections 3.3 and 3.4.

Per equation (35) below the changes of variables utilized in this paper are based on use of the CoV function [16, 31]

$$w(\theta)=2\pi\frac{[v(\theta)]^{p}}{[v(\theta)]^{p}+[v(2\pi-\theta)]^{p}},\quad 0\leq\theta\leq 2\pi,$$

(32)

which is depicted in Figure 2, where

$$v(\theta)=\left(\frac{1}{p}-\frac{1}{2}\right)\left(\frac{\pi-\theta}{\pi}\right)^{3}+\frac{1}{p}\left(\frac{\theta-\pi}{\pi}\right)+\frac{1}{2}.$$

The derivatives of the function $w(\theta)$ of orders less than or equal to $p-1$ vanish at $\theta=0$ and $\theta=2\pi$—so that, in particular, we have

$$w(\theta)\sim\theta^{p}\quad\mbox{as}\quad\theta\to 0.$$

(33)

While simpler changes of variables such as the monomial CoV [25]

$$s(\theta)=\theta^{p},$$

(34)

$p-1$ of whose derivatives vanish at $\theta=0$ can also be effective, the CoV (32), which is used throughout this paper (with exception of the comparison results presented in Figure 12 right), gives rise to significantly better performance than the simpler version (34) on account of the rather evenly split distribution it enjoys of dicretization points near to and away from the endpoints of the interval $0\leq\theta\leq 2\pi$ [16]. Results based on the monomial change of variables (34) presented in Figure 12 right provide an indication of the improvements that result from the use of the CoV (32).

Recalling that, as indicated at the beginning of this section, corners are assumed to occur only at $s=0$ points, on any one of the patches $\Gamma^{q}$ the CoV is induced by the change of variables (32) in accordance with the expression

$$\left\{\begin{aligned} &s_{q}(\theta)=\theta,\quad\mbox{if $\Gamma^{q}$ has no corners},\\ &s_{q}(\theta)=\frac{1}{\pi}w(\pi\theta),\quad\mbox{if $\Gamma^{q}$ has a corner at $\theta=0$.}\end{aligned}\right.$$

(35)

Note that the first $p-1$ derivatives of $s(\theta)$ vanish at $\theta=0$ for patches containing a corner and $\frac{ds}{d\theta}=1$ for patches not containing a corner.

Incorporating the CoV (35) and defining

$$\widetilde{r}_{q}(\theta)=r_{q}(s_{q}(\theta))\quad\mbox{and}\quad\widetilde{L}_{q^{\prime}}(\theta^{\prime})=L_{q^{\prime}}(s(\theta^{\prime}))\frac{ds_{q^{\prime}}(\theta^{\prime})}{d\theta^{\prime}},$$

(36)

the adjoint double-layer integral in (31) (that is, the first integral on the left-hand side in that equation) with $r=\widetilde{r}_{q}(\theta)$ becomes

$$\int_{0}^{1}\frac{\partial G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime})}{\partial n(\widetilde{r}_{q}(\theta))}\phi_{q^{\prime}}(\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\widetilde{L}_{q^{\prime}}(\theta^{\prime})\,d\theta^{\prime}.$$

(37)

(It is important to recall that, in view of the definitions and assumptions introduced above in this section, for patches $\Gamma^{q^{\prime}}$ containing a corner $C$, the corner point is necessarily given by $C=\widetilde{r}_{q}(0)$.) The regularization operator (30), in turn, takes the form

$$R[\phi](\widetilde{r}_{q^{\prime}}(\theta^{\prime}))=\sum_{q^{\prime\prime}=1}^{M}\int_{0}^{1}G_{ik}(\widetilde{r}_{q^{\prime}}(\theta^{\prime})),\widetilde{r}_{q^{\prime\prime}}(\theta^{\prime\prime}))\phi_{q^{\prime\prime}}(\widetilde{r}_{q^{\prime\prime}}(\theta^{\prime\prime}))\widetilde{L}_{q^{\prime\prime}}(\theta^{\prime\prime})\,d\theta^{\prime\prime}.$$

(38)

Similarly, for the last two integrals in (31) we obtain the expressions

$$\int_{0}^{1}G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime})(n(\widetilde{r}_{q}(\theta))\cdot n(\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\widetilde{L}_{q^{\prime}}(\theta^{\prime})R[\phi](\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\,d\theta^{\prime}$$

(39)

and

$$\frac{1}{\widetilde{L}_{q}(\theta)}\frac{\partial}{\partial\theta}\int_{0}^{1}G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\frac{\partial}{\partial\theta^{\prime}}R[\phi](\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\,d\theta^{\prime},$$

(40)

respectively.

In all, incorporating the CoV (35) in (31), using the expressions (37) through (40), and multiplying the result by $\widetilde{L}_{q}(\theta)$, we obtain

$$\frac{i\eta}{2}\phi(\widetilde{r}_{q}(\theta))\widetilde{L}_{q}(\theta)-i\eta\widetilde{L}_{q}(\theta)\sum_{q^{\prime}=1}^{M}\int_{0}^{1}\frac{dG_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime})}{\partial n(\widetilde{r}_{q}(\theta))}\phi(\widetilde{r}_{q^{\prime}}(\theta^{\prime})\widetilde{L}_{q^{\prime}}(\theta^{\prime})\,d\theta^{\prime}\\ +{k}^{2}\widetilde{L}_{q}(\theta)\sum_{q^{\prime}=1}^{M}\int_{0}^{1}G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime}))(n(\widetilde{r}_{q}(\theta))\cdot n(\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\widetilde{L}_{q^{\prime}}(\theta^{\prime})R[\phi](\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\,d\theta^{\prime}\\ +\frac{\partial}{\partial\theta}\sum_{q^{\prime}=1}^{M}\int_{0}^{1}G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\frac{\partial}{\partial\theta^{\prime}}R[\phi](\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\,d\theta^{\prime}=-\widetilde{L}_{q}(\theta)\frac{\partial u^{\text{inc}}(\widetilde{r}_{q}(\theta))}{\partial n(\widetilde{r}_{q}(\theta))}.$$

(41)

Thus, defining the new unknown

$$\psi_{q}(\theta)=\phi(s_{q}(\theta))\widetilde{L}_{q}(\theta),$$

(42)

and expressing the regularization operator (38) in the form

$$\widetilde{R}[\psi_{q}](r)=\sum_{q^{\prime}=1}^{M}\int_{0}^{1}G_{ik}(r,\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\,\psi_{q}(\theta^{\prime})\,d\theta^{\prime},$$

(43)

we finally obtain the change-of-unknown, Corner-Regularized version

$$\frac{i\eta}{2}\psi_{q}(\theta)-i\eta\widetilde{L}_{q}(\theta)\sum_{q^{\prime}=1}^{M}\int_{0}^{1}\frac{\partial G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime})}{\partial n(\widetilde{r}_{q}(\theta))}\psi_{q}(\theta^{\prime})\,d\theta^{\prime}\\ +{k}^{2}\widetilde{L}_{q}(\theta)\sum_{q^{\prime}=1}^{M}\int_{0}^{1}G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime}))(n(\widetilde{r}_{q}(\theta))\cdot n(\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\widetilde{L}_{q^{\prime}}(\theta^{\prime})R[\psi](\widetilde{r}_{q^{\prime}}(\theta))\,d\theta^{\prime}\\ +\frac{\partial}{\partial\theta}\sum_{q^{\prime}=1}^{M}\int_{0}^{1}G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\frac{\partial}{\partial\theta^{\prime}}R[\psi](\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\,d\theta^{\prime}=-\frac{\partial u^{\text{inc}}}{\partial n(\widetilde{r}_{q}(\theta))}\widetilde{L}_{q}(\theta)$$

(44)

of the Operator Regularized R-CFIE equation (31). It is important to emphasize that $\psi(\theta)$ is a much smoother function than $\phi$ (and, thus, significantly easier to approximate by polynomials) on account of the multiplicative line element $\widetilde{L}_{q}(\theta)$ factor contained in $\psi$, whose value together with that of $p-1$ of its derivatives, vanishes at the corner point $\theta=0$ where the field singularity is located. This is another key aspect of the proposed methodology for scattering by domains with corners.

For the purpose of comparison in the numerical results section we additionally consider the following change-of-unknown, corner-regularized version of the MFIE equation (5):

$$-\frac{1}{2}\psi_{q}(\theta)+\widetilde{L}_{q}(\theta)\sum_{q^{\prime}=1}^{M}\int_{0}^{1}\frac{\partial G_{k}(\widetilde{r}_{q}(\theta),\widetilde{r}_{q^{\prime}}(\theta^{\prime})}{\partial n(\widetilde{r}_{q}(\theta))}\psi_{q^{\prime}}(\theta^{\prime})\,d\theta=-\frac{\partial u^{\text{inc}}}{\partial n(\widetilde{r}_{q}(\theta))}\widetilde{L}_{q}(\theta).$$

(45)

Note that in light of the introduction of the new unknown $\psi(\theta)$, which includes a the $\widetilde{L}_{q}(\theta)$ factor, the rectangular polar discretization method needs to be modified as follows:

1. 1.

   Substitute the Fejér evaluation of the far interactions (24) by the corresponding Fejér approximation

   $$\mathcal{I}^{q^{\prime}}[H,\widetilde{f}](r)\approx\sum_{j=0}^{N-1}H(r,\widetilde{r}_{q^{\prime}}(\theta_{j}))\widetilde{f}((\theta_{j}))\omega_{j},\quad\mbox{for }\quad r\quad\mbox{away from}\quad\Gamma^{q^{\prime}}.$$

   (46)

2. 2.

   Substitute (25) by the corresponding expansion

   $$\widetilde{f}_{q^{\prime}}(\theta^{\prime})\approx\sum_{m=0}^{N-1}\widetilde{a}_{m}^{q^{\prime}}T_{m}(2\theta^{\prime}-1)$$

   (47)

   (where, in analogy to Section 2.2, here we set either $\widetilde{f}_{q}(\theta^{\prime})=\psi_{q}(\theta^{\prime})$ or $\widetilde{f}_{q}(\theta^{\prime})=\frac{\partial\psi_{q}}{\partial\theta}(\theta^{\prime})$).

3. 3.

   Replace the precomputation integrals (26) by

   $$\widetilde{\beta}_{m}^{q^{\prime}}[H](r)=\int_{0}^{1}H(r,\widetilde{r}_{q^{\prime}}(\theta^{\prime}))\,\,T_{m}(2\theta^{\prime}-1)\,d\theta^{\prime},$$

   (48)

   for the kernels $H$ listed near the end of Section 2.2, and for target discretization points $r=\widetilde{r}_{q^{\prime}}(\theta_{j})$ on the source patch $\Gamma^{q^{\prime}}$, as well as for near singular target discretization points $r=\widetilde{r}_{q}(\theta_{j})$ on nearby patches $\Gamma^{q}$. (Note that, in particular, the $\widetilde{L}_{q^{\prime}}(\theta^{\prime})$ factor, which is contained in the new unknown $\psi(\theta^{\prime})$, does not appear in the integrand of equation (48).) In view of the singularities introduced by the corner geometries under consideration, certain significant variations of the rectangular polar method [12, 28] are needed, which are described in Sections 3.3 and 3.4.

4. 4.

   Using the precomputed weights (48) evaluate the integral

   $$\mathcal{I}^{q^{\prime}}[H,\widetilde{f}](r)=\sum_{m=0}^{N-1}\widetilde{a}_{m}^{q^{\prime}}\widetilde{\beta}_{m}^{q^{\prime}}[H](r)$$

   (49)

   in terms of the coefficients $\widetilde{a}^{q^{\prime}}_{m}$ in (47), and evaluate the needed overall integrals by a sum analogous to (22). The tangential differentiations needed for the evaluation of the operator $T_{3}$ in (14) are produced by a procedure analogous to the one described in Section 2.4 applied in the $\theta^{\prime}$ variable (instead of the $s^{\prime}$ variable used in that section).