# A Galerkin BEM for high-frequency scattering problems based on frequency-dependent changes of variables

A Galerkin BEM for high-frequency scattering problems based on frequency-dependent changes of... Abstract In this paper we develop a new class of semidiscrete Galerkin boundary element methods for the solution of two-dimensional exterior single-scattering problems. Our approach is based upon construction of Galerkin approximation spaces confined to the asymptotic behavior of the solution through a certain direct sum of appropriate function spaces weighted by the oscillations in the incident field of radiation. Specifically, the function spaces in the illuminated/shadow regions and the shadow boundaries are simply algebraic polynomials, whereas those in the transition regions are generated utilizing novel, yet simple, frequency-dependent changes of variables perfectly matched with the boundary layers of the amplitude in these regions. We rigorously verify for compact, smooth and strictly convex obstacles that, with increasing wavenumber k, these methods require, for any ε > 0, only an $$\mathscr{O}\left ( k^{\varepsilon } \right )$$ increase in the number of degrees of freedom (DoF) to maintain any given accuracy independent on frequency. These theoretical results are confirmed by numerical tests that show that, for sufficiently large DoF, the error tends to decrease with increasing wavenumber k. 1. Introduction High-frequency scattering problems continue to be of immense interest in present day computational science. Indeed, over the course of last two decades, very efficient and effective algorithms have been devised for the numerical solution of scattering problems based on variational (Hesthaven & Warburton, 2004; Davies et al., 2009; Boffi, 2010) and integral equation (Bruno & Kunyansky, 2001; Banjai & Hackbusch, 2008; Tong & Chew, 2010; Bruno et al., 2013) formulations. In exterior scattering simulations methods that rest on variational formulations naturally demand the design and implementation of efficient nonreflecting boundary conditions (Engquist & Majda, 1977; Givoli, 2004; Grote & Sim, 2011) to effectively represent the radiation condition at infinity. On the other hand, solvers based on integral equation formulations (cf. the survey article Chandler-Wilde et al. (2012) on the numerical-asymptotic boundary integral equation methods) readily encode the radiation condition into the equation by choosing an outgoing fundamental solution. Moreover, for surface scattering simulations as considered in this paper, they enable phase extraction, that takes a particularly simple form in single-scattering configurations, and this turns the problem into the estimation of an amplitude whose oscillations are essentially independent of frequency. In this paper we develop a new class of semidiscrete Galerkin boundary element methods (BEMs) for the solution of two-dimensional scattering problems in the exterior of compact, smooth and convex obstacles having everywhere positive curvature. The approach taken here is based on the completely rigorous asymptotic expansions of the aforementioned amplitude developed by Melrose & Taylor (1985) (for their alternative forms, see Lazergui & Boubendir (2016)). These expansions imply the existence of boundary layers around the shadow boundaries (where the rays are tangential to the boundary) and the Galerkin approximation spaces are constructed so as to resolve these layers through frequency-dependent changes of variables. From a theoretical point of view, the resulting methodologies require, for any ε > 0, only an $$\mathscr{O}\left ( k^{\varepsilon } \right )$$ increase in the number of degrees of freedom (DoF) to maintain any given accuracy independent on frequency. These theoretical findings are supported by numerical implementations that show that the method achieves similar accuracy even when frequency-independent numbers of DoF are used. From an application point of view, the tests also show that the algorithms proposed herein are applicable to the class of nonconvex smooth scatterers and incidence directions that do not allow for multiple reflections and more than two shadow boundaries. Indeed, hybrid integral equation methodologies incorporating the asymptotic characteristics of the unknown into the solution strategy have now become the usual practice in the field. The first attempt in this direction is due to Abboud et al. (1994, 1995) where, considering the impedance boundary condition, an h-version BEM was utilized in conjunction with the method of stationary phase for the evaluation of highly oscillatory integrals (see Darrigrand (2002) for a fast multipole variant, and Ganesh & Hawkins (2011) for a fully discrete three-dimensional version). More relevant to our work is the Nyström method proposed by Bruno et al. (2004) for the solution of sound-soft scattering problems in the exterior of smooth convex obstacles. The method therein displays the capability of delivering solutions within any prescribed accuracy in frequency-independent computational times. While, in this approach, the boundary layers of the slowly varying amplitude around the shadow boundaries are resolved through a cubic root change of variables, the associated highly oscillatory integrals are evaluated to high order utilizing novel extensions of the method of stationary phase (for a three-dimensional variant of this approach we refer to Bruno & Geuzaine (2007)). The algorithm in Bruno et al. (2004) has had a great impact in the computational scattering community and, following the basic principles therein, a number of alternative single-scattering solvers have been developed. Giladi (2007) has used a collocation method that integrates Keller’s geometrical theory of diffraction (GTD) to account for creeping rays in the shadow region. Huybrechs and Vandewalle’s collocation method (Huybrechs & Vandewalle, 2007) has utilized the numerical steepest descent method in evaluating highly oscillatory integrals and additional collocation points around shadow boundaries to obtain sparse discretizations. The first rigorous numerical analysis relating to a p-version boundary element implementation of these approaches, due to Domínguez et al. (2007), proved that an increase of $$\mathscr{O}(k^{1/9})$$ in the number of DoF is sufficient to preserve a certain accuracy as $$k \to \infty$$. Parallel with the schemes relating to smooth convex obstacles, Galerkin BEMs based on phase extraction have also been developed for half-planes and convex polygons where the number of DoF is either fixed or must increase in proportion to $$\log k$$ to fix the error with increasing wavenumber k. We refer to the survey article Chandler-Wilde et al. (2012) for an extended review of these procedures; for more recent hp BEMs for convex polygons see Hewett et al. (2013), and for screens and apertures see Hewett et al. (2015). As for multiple scattering problems we refer to Bruno et al. (2005) for an extension of the algorithm in Bruno et al. (2004) to a finite collection of convex obstacles (see also Ecevit & Reitich (2009) and Anand et al. (2010) for a rigorous analysis of this approach in two- and three-dimensional settings, respectively, and Boubendir et al. (2016) for the acceleration of this procedure through use of dynamical Krylov subspaces and Kirchhoff approximations), and to Chandler-Wilde et al. (2015) for a class of nonconvex polygons. Concerning compact, smooth and strictly convex obstacles considered in this paper, let us mention that the methods developed in Bruno et al. (2004), Huybrechs & Vandewalle (2007) and Domínguez et al. (2007) remain asymptotic due to approximation of the solution by zero in the deep shadow region. In our recent frequency-adapted Galerkin BEMs (Ecevit & Özen, 2017), on the other hand, we utilized approximation spaces also in this region for two reasons. First, based on an optimal use of the wavenumber-dependent derivative estimates of the amplitude given by Domínguez et al. (2007, Theorem 5.4), this allowed us to rigorously eliminate the asymptotic terms appearing in the error analysis due to approximation of the solution by zero in the deep shadow region. Secondly, as displayed through numerical tests, use of these additional approximation spaces significantly reduces the condition numbers of the Galerkin linear systems and renders them virtually independent on frequency. This, in return, improves the numerical stability of the method and allows the use of higher DoF within a limited numerical precision which results in improved approximations. Let us mention that our approach in Ecevit & Özen (2017) was based on utilization of appropriate integral equation formulations of the scattering problem and design of Galerkin approximation spaces, as the direct sum of algebraic or trigonometric polynomials weighted by the oscillations in the incident field of radiation. For any given $$m \in \mathbb{N}$$, the number of direct summands was chosen as 4m by assigning one approximation space to each of the illuminated and deep shadow regions, and one to each of the two shadow boundaries, and m − 1 to each one of the four transition regions. As we have shown in Ecevit & Özen (2017), from a theoretical perspective, m had to increase as $$\mathscr{O}(\log k)$$ in order to obtain optimal error bounds, and that these methods can be tuned to demand an increase of $$\mathscr{O}(k^{\varepsilon })$$ (for any ε > 0) in the number of DoF to maintain a prescribed accuracy independent on frequency. For smooth strictly convex scatterers, this is the best theoretical result available in the literature. We note, however, that the numerical results in Ecevit & Özen (2017) were restricted to m = 1 and m = 2, since utilization of larger values of m is simply impractical from an implementation point of view. In the new Galerkin BEMs developed herein we continue to use approximation spaces in the deep shadow region, since this allows the elimination of the asymptotic terms in the error analysis, and their utilization significantly reduces the condition numbers of the Galerkin linear systems. However, we treat the transition regions differently as we utilize novel frequency-dependent changes of variables perfectly matched with the asymptotic behavior of the solution in the transition regions. We thereby eliminate the requirement of increasing the number of direct summands defining the Galerkin approximation spaces with increasing wavenumber k. While this makes the new scheme we propose far easier to implement than that in Ecevit & Özen (2017), from a theoretical point of view, an increase of $$\mathscr{O}(k^{\varepsilon })$$ (for any ε > 0) is still sufficient to fix the approximation error with increasing k. The numerical tests show that, for sufficiently large DoF, the new scheme outperforms the method in Ecevit & Özen (2017) and the relative L2 error tends to decrease with increasing wavenumber k. Let us note that, in fact, accurately approximating the solution in the (deep) shadow region involves significant numerical difficulties since the solution there decays exponentially with increasing wavenumber and geodesic distance to the shadow boundaries. Moreover, its oscillations are further amplified by a phase term that depends on the curvature and geodesic distance to the shadow boundaries (see Remark 3 for a more precise explanation; also cf. Asheim & Huybrechs (2014) and the references therein). For this reason, we test the performance of the new Galerkin BEMs in the shadow region in terms of pointwise relative errors and, in this connection, discuss the decaying behavior of the solution in the shadow region depending on the wavenumber and the radius of curvature. Further, we compare the numerical results with those of Asheim & Huybrechs (2014), where they have developed adequate phase extraction techniques (based on suitable combinations of GTD and Melrose--Taylor asymptotics) more appropriate for approximating the solution in the shadow region. Incidentally, it is important to note that, just as the methods in Domínguez et al. (2007) and Ecevit & Özen (2017), the methods developed here are semidiscrete as they transform the solution of the scattering problem to the evaluation of highly-oscillatory double integrals. Efficient evaluation of integrals arising in connection with the method in Domínguez et al. (2007) is discussed in both theoretical and practical terms by Kim (2012). With some additional effort, the approach therein can be extended to cover the integrals associated with the method in Ecevit & Özen (2017). For the new methods we propose here, the structure of these highly-oscillatory double integrals are slightly different, but they can be treated similarly. The actual theoretical and implementation details are left for future work. The paper is organized as follows. In §2, we describe the exterior sound-soft scattering problem along with the relevant integral equations and associated Galerkin formulations. In §3, we introduce the new Galerkin schemes for high-frequency single-scattering problems, and state the associated convergence theorem for compact, smooth, and strictly convex obstacles which constitutes the main result of the paper. To allow a direct comparison, in the same section, we also present a more general version of our algorithm in Ecevit & Özen (2017) along with the corresponding approximation properties. The proof of the main result is given in §4. In §5.1 we present the highly-oscillatory double integrals arising in connection with the methods developed herein and discuss implementation details. §5.2 is reserved for numerical tests that provide a comparison of our methods developed herein and Ecevit & Özen (2017). Finally, in §5.3, we examine approximations in the shadow region particularly in contrast with those in Asheim & Huybrechs (2014). 2. The scattering problem and Galerkin formulation The two-dimensional scattering problem we consider in this manuscript is related with the determination of the scattered field u that satisfies the Helmholtz equation \begin{align*} \left( \varDelta + k^{2} \right) u = 0 \end{align*} in the exterior of a compact smooth strictly convex obstacle K, the Sommerfeld radiation condition at infinity that amounts to requiring \begin{align*} \lim_{\left| x \right| \to \infty} \left| x \right|{}^{1/2} \left[ \partial_{\left| x \right|} - i k \right] u = 0 \end{align*} uniformly for all directions $$x/\left | x \right |$$, and the sound-soft boundary condition \begin{align*} u = - u^{\textrm{inc}} \end{align*} for a plane-wave incidence $$u^{\textrm{inc}} \left ( x \right ) = e^{ik \alpha \cdot x}$$ with direction α ($$\left | \alpha \right | =1$$) impinging on K. As is well known, the scattered field u can be reconstructed by means of either the direct or the indirect approach (Colton & Kress, 1992); for discussions on the direct and indirect approaches we refer to Bruno et al. (2004, pp. 631–632) and Chandler-Wilde et al. (2012, p. 132). As in the previous attempts aimed at frequency-independent simulations (Bruno et al., 2004; Domínguez et al., 2007; Giladi, 2007; Huybrechs & Vandewalle, 2007; Ecevit & Özen, 2017;), here we favor the former wherein the associated (unknown) surface density is the normal derivative of the total field (known as the surface current in electromagnetism) $$\eta = \partial _{\nu } \left ( u+ u^{\textrm{inc}} \right )$$ on ∂K. Once η is available, the scattered field can be recovered through the single-layer potential \begin{align*} u(x) = - \int_{\partial K} \varPhi(x,y) \, \eta(y) \, \textrm{d}s(y), \end{align*} where \begin{align*} \varPhi(x,y) = \dfrac{i}{4} \ H_{0}^{(1)}(k|x-y|) \end{align*} is the fundamental solution of the Helmholtz equation, and $$H_{0}^{(1)}$$ is the Hankel function of the first kind and order zero. The new unknown, namely η, can be expressed as the unique solution of a variety of integral equations (Colton & Kress, 1992; Spence et al., 2011) taking the form of an operator equation \begin{align} \mathscr{R}_{k} \, \eta = f_{k} \end{align} (2.1) in $$L^{2} \left ( \partial K \right )$$. The solution of (2.1) corresponds exactly to that of \begin{align} B_{k}(\mu,\eta) = F_{k}(\mu), \quad \textrm{for all } \mu \in L^{2}(\partial K), \end{align} (2.2) where the sesquilinear form Bk and bounded linear functional Fk are defined by \begin{align} B_{k}(\mu,\eta) = \langle \mu, \mathscr{R}_{k} \, \eta \rangle_{L^{2}(\partial K)} \quad \textrm{and} \quad F_{k}(\mu) = \langle \mu, f_{k} \rangle_{L^{2}(\partial K)}. \end{align} (2.3) Equation (2.2), in turn, is amenable to a treatment by the Galerkin method wherein one determines the Galerkin solution$$\hat{\eta }$$ approximating the exact solution η in a given finite dimensional Galerkin subspace$$\hat{X}_{k}$$ requiring that \begin{align} B_{k}(\hat{\mu},\hat{\eta}) = F_{k}(\hat{\mu}), \quad \textrm{for all } \hat{\mu} \in \hat{X}_{k}. \end{align} (2.4) Further, provided the sesquilinear form Bk is continuous with a continuity constant Ck and strictly coercive with a coercivity constant ck so that \begin{align*} \left| B_{k}(\mu,\eta) \right| \le C_{k} \, \Vert \mu \Vert \Vert \eta \Vert \qquad \textrm{and} \qquad \mathscr{R}e \ B_{k}(\mu,\mu) \ge c_{k} \, \Vert \mu \Vert^{2} \end{align*} for all μ, η ∈ L2(∂K), equation (2.4) is uniquely solvable and C$${\acute{\textrm{e}}}$$a's lemma implies \begin{align} \Vert \eta -\hat{\eta} \Vert \leq \dfrac{C_{k}}{c_{k}} \inf_{\hat{\mu} \in \hat{X}_{k}} \Vert \eta- \hat{\mu} \Vert. \end{align} (2.5) Remark 2.1 Among the aforementioned integral equations, in this connection, combined field integral equation (CFIE) and star-combined integral equation (SCIE) step forward as the continuity and coercivity properties of the associated sesquilinear forms are well-understood. More precisely, the sesquilinear form associated with CFIE is continuous for Lipschitz domains for all k > 0 with $$C_{k} = \mathscr{O} \left ( k^{1/2} \right )$$ as $$k \to \infty$$ (Chandler-Wilde et al., 2009, Theorem 3.6), and coercive (for k ≫ 1) for strictly convex domains having piecewise analytic C3 boundaries with everywhere positive curvature such that, for any 0 < δ < 1/2, ck ≥ (1/2) − δ for all sufficiently large k (Spence et al., 2015, Theorem 1.2). On the other hand, the sesquilinear form corresponding to SCIE is both continuous and coercive (for k > 0) for star-shaped Lipschitz domains with $$C_{k} = \mathscr{O} \left ( k^{1/2} \right )$$ as $$k \to \infty$$ and ck independent of k (Spence et al., 2011). Moreover, recently obtained bounds, based on the results in Han & Tacy (2015), on the single- and double-layer potentials (see Galkowski (2014, Theorems 4.29 and 4.32) and Galkowski et al. (2016, Theorem 1.4) for CFIE, and Galkowski (2014, Theorems 4.29 and 4.32) and Galkowski et al. (2016, Theorem 1.9, Remark 4.2) for SCIE) show that, for strictly convex domains having everywhere positive curvature, $$C_{k} = \mathscr{O} \left ( k^{1/3} \right )$$ as $$k \to \infty$$. 3. Galerkin approximation spaces based on frequency-dependent changes of variables The developments in this section are independent of the integral equation used as they relate, specifically, to the construction of Galerkin approximation spaces $$\hat{X}_{k}$$ whose dimension should increase only as $$\log k$$ with increasing wavenumber k to ensure that the relative error \begin{align*} \inf_{\hat{\mu} \in \hat{X}_{k}} \dfrac{\Vert \eta- \hat{\mu} \Vert}{\Vert \eta \Vert} \end{align*} in connection with the infimum on the right-hand side of (2.5) is independent of k. Considering a smooth convex obstacle K illuminated by a plane-wave uinc(x) = eikα⋅x, our approach is based on phase extraction \begin{align*} \eta \left( x \right) = e^{ik \alpha \cdot x} \, \eta^{\textrm{slow}} \left( x \right), \qquad x \in \partial K, \end{align*} and design of approximation spaces adopted to the asymptotic behavior of the amplitude ηslow that was initially characterized by Melrose & Taylor (1985) around the shadow boundaries which we have later generalized to the entire boundary (Ecevit & Reitich, 2009). Remark 3.1 The next theorem is included to highlight the difficulties associated with the construction of Galerkin approximation spaces, and it is not used in any of the proofs. In particular, the definition of Hörmander classes and a detailed understanding of asymptotic expansions are not necessary to follow the developments in the rest of the paper (the interested reader may consult to Section 2.2 in Ecevit & Reitich (2009) for a short exposition on Hörmander classes and asymptotic expansions). What is needed for the analysis in this paper is the wavenumber explicit derivative estimates of ηslow as given in Theorem 4.2 below which is based on Domínguez et al. (2007, Theorem 5.4). Theorem 3.1 (Ecevit & Reitich, 2009, Corollary 2.1) Let $$K \subset \mathbb{R}^{2}$$ be a compact, strictly convex set with smooth boundary ∂K. Then ηslow = ηslow(x, k) belongs to the Hörmander class $$S^{1}_{2/3,1/3} \left ( \partial K \times \left ( 0, \infty \right ) \right )$$ and admits an asymptotic expansion \begin{align} \eta^{\textrm{slow}}(x,k) \sim \sum_{p,q \ge 0} a_{p,q}\left( x,k \right) \end{align} (3.1) with \begin{align*} a_{p,q}(x,k) = k^{2/3-2p/3-q} \, b_{p,q}(x) \, \varPsi^{(p)}(k^{1/3}Z(x)), \end{align*} where bp, q and Ψ are complex-valued $$C^{\infty }$$ functions and Z is a real-valued $$C^{\infty }$$ function that is positive on the illuminated region ∂KIL = {x ∈ ∂K : α ⋅ ν(x) < 0}, negative on the shadow region ∂KSR = {x ∈ ∂K : α ⋅ ν(x) > 0} and vanishes precisely to the first order on the shadow boundaries ∂KSB = {x ∈ ∂K : α ⋅ ν(x) = 0}. Moreover, the function Ψ admits the asymptotic expansion \begin{align*} \varPsi(\tau) \sim \sum_{j=0}^{\infty} c_{j} \tau^{1-3j} \qquad \textrm{as } \tau \to \infty, \end{align*} and Ψ is rapidly decreasing in the sense of Schwartz as $$\tau \to -\infty$$. Theorem 3.1 clearly displays the challenges associated with the construction of approximation spaces adapted to the asymptotic behavior of ηslow as it shows that while ηslow admits a classical asymptotic expansion in the illuminated region and rapidly decays in the shadow region, it possesses boundary layers in the $$\mathscr{O}(k^{-1/3})$$ neighborhoods of the shadow boundaries (see Fig. 1). We overcome this difficulty by constructing approximation spaces improving upon our approach in Ecevit & Özen (2017), and better adapted to the changes in frequency through use of wavenumber dependent changes of variables that resolve the aforementioned boundary layers, and that provides a smooth transition from the shadow boundaries into the illuminated and shadow regions. The resulting schemes, when compared with our algorithms in Ecevit & Özen (2017), are easier to implement since the Galerkin spaces are represented as the direct sum of a fixed number of approximation spaces rather than a number that should increase in proportion to $$\log k$$ as in Ecevit & Özen (2017). Moreover, as we explain, they display better approximation properties from a theoretical perspective since they provide an improvement of order $$\sqrt{\log k}$$ in the error for the same order of DoF as $$k\to \infty$$. Fig. 1. View largeDownload slide The unit circle illuminated from the left and the associated boundary layers of the solutions. Fig. 1. View largeDownload slide The unit circle illuminated from the left and the associated boundary layers of the solutions. Remark 3.2 ηslow is in fact both exponentially small and highly oscillatory in the deep shadow region which is not directly apparent from the rapid decay of Ψ(τ) as $$\tau \to -\infty$$ mentioned in Theorem 3.1. Indeed, there exist constants c0≠0 and β > 0 such that \begin{align*} \varPsi(\tau) = c_{0} \, e^{-i\tau^{3}/3-i\tau\alpha_{1}} (1+ \mathscr{O}(e^{\beta\tau})), \qquad \textrm{as } \tau \to -\infty, \end{align*} where $$\alpha _{1} = e^{-2i \pi \nu _{1}/3}$$ and ν1 < 0 is the right-most root of the airy function Ai (Domínguez et al., 2007, Theorem 5.1). Therefore, in the shadow region, the leading term in (3.1) is asymptotically given by \begin{align} c_{0} \, k^{2/3} \, b_{0,0}(x) \, e^{ik(Z(x))^{3}/3-ik^{1/3}Z(x)\alpha_{1}} \end{align} (3.2) which justifies the exponential decay and highly oscillatory behavior of ηslow in the deep shadow region. On the other hand, to a leading order approximation, Z(x) can be expressed as a constant multiple of the integral of κ2/3 (κ being the curvature) between x and the nearest shadow boundary point (Asheim & Huybrechs, 2014). Accordingly, equation (3.2) also justifies the statement in the introduction, namely that in the shadow region the solution decays exponentially with increasing wavenumber and geodesic distance to the shadow boundaries, and its oscillations are further amplified by a phase term that depends on the curvature and geodesic distance to the shadow boundaries. These observations motivate one to either approximate η in the deep shadow region simply by zero (which would clearly add asymptotic terms in the error analysis as in Domínguez et al. (2007)), or construct the Galerkin approximation spaces so as to reflect the exponentially decaying and highly oscillatory behavior of ηslow in the deep shadow region. While the Galerkin approximation spaces we propose below do not approximate η by zero in the deep shadow region, they are not designed to capture the oscillations in (3.2) at large k either. Indeed they do not need to since, as we shall see, the exponential damping implies that the contribution from the deep shadow region is provably negligible. To describe our approximation spaces, we choose γ to be the L−periodic arc length parameterization of ∂K in the counterclockwise orientation with $$\alpha \cdot \nu \left ( \gamma (0) \right ) = 1$$. This ensures that if 0 < t1 < t2 < L are the points corresponding to the shadow boundaries \begin{align*} \gamma \left( \left\{ t_{1},t_{2} \right\} \right) = \partial K^{SB}, \end{align*} then the illuminated and shadow regions are given by (see Fig. 1(a)) \begin{align*} \gamma \left( \left( t_{1},t_{2} \right) \right) = \partial K^{IL} \qquad \textrm{and} \qquad \gamma \left( \left( t_{2},t_{1}+L \right) \right) = \partial K^{SR}. \end{align*} For k > 1, we introduce two types of approximation spaces confined to the regions depicted in Fig. 2(a)–(b). In both cases, we define the illuminated transition and shadow transition intervals as \begin{align*} I_{IT_{1}} & = [t_{1} + \xi_{1} k^{-1/3}, t_{1} + \xi_{1}^{\prime}] = [a_{1},b_{1}], \\ I_{IT_{2}} & = [t_{2} - \xi_{2}^{\prime}, t_{2} - \xi_{2} k^{-1/3} ] = [a_{2},b_{2}], \\ I_{ST_{1}} & = [t_{1} - \zeta_{1}^{\prime}, t_{1} - \zeta_{1} k^{-1/3} ] = [a_{3},b_{3}], \\ I_{ST_{2}} & = [t_{2} + \zeta_{2} k^{-1/3}, t_{2} + \zeta_{2}^{\prime}] = [a_{4},b_{4}], \end{align*} and the shadow boundary intervals as \begin{align*} I_{sb:{1}} & = [t_{1} - \zeta_{1} k^{-1/3}, t_{1} + \xi_{1} k^{-1/3}] = [a_{5},b_{5}], \\ I_{sb:{2}} & = [t_{2} - \xi_{2} k^{-1/3}, t_{2} + \zeta_{2} k^{-1/3} ] = [a_{6},b_{6}], \end{align*} where the parameters $$\xi _{j}, \xi _{j}^{\prime }, \zeta _{j},\zeta _{j}^{\prime }>0$$ (j = 1, 2) are chosen so that \begin{align*} t_{1} + \xi_{1} \le t_{1} + \xi_{1}^{\prime} \overset{(A)}{\le} t_{2} - \xi_{2}^{\prime} \le t_{2} - \xi_{2}, \end{align*} and \begin{align*} t_{2} + \zeta_{2} \le t_{2} + \zeta_{2}^{\prime} \overset{(B)}{\le} L + t_{1} - \zeta_{1}^{\prime} \le L + t_{1} - \zeta_{1}. \end{align*} Fig. 2. View largeDownload slide Regions on the boundary of the unit circle. Fig. 2. View largeDownload slide Regions on the boundary of the unit circle. Note specifically that the regions in Fig. 2(a) correspond to equalities in (A) and (B), and those in Fig. 2(b) to strict inequalities. In the latter case, we define the illuminated and deep shadow intervals as \begin{align*} I_{IL} & = [t_{1} + \xi_{1}^{\prime}, t_{2} - \xi_{2}^{\prime}] = [a_{7},b_{7}], \\ I_{DS} & = [ t_{2} + \zeta_{2}^{\prime}, L + t_{1} - \zeta_{1}^{\prime}] = [a_{8},b_{8}]. \end{align*} With these choices, given $$\mathbf{d} = \left ( d_{1}, \ldots , d_{J} \right ) \in \mathbb{Z}_{+}^{J}$$ where $$\mathbb{Z}_{+} = \mathbb{N} \cup \{ 0 \}$$ (with J = 6 and J = 8 for Fig. 2(a) and (b), respectively), we define our (|d| + J)−dimensional Galerkin approximation spaces based on algebraic polynomials and frequency dependent changes of variables as \begin{align} \ \end{align} (3.3) Here $$_{\mathscr{I}_{j}}$$ is the characteristic function of $$\mathscr{I}_{j} = [a_{j},b_{j}]$$, and \begin{align*} \mathscr{A}_{d_{j}}^{\mathscr{C}} = \begin{cases}{cl} \mathbb{P}_{d_{j}} \circ \phi^{-1}, & \mbox{if } \mathscr{I}_{j} \mbox{ is a transition region}, \\[0.5 em] \mathbb{P}_{d_{j}}, & \textrm{otherwise}, \end{cases} \end{align*} where $$\mathbb{P}_{d_{j}}$$ is the vector space of algebraic polynomials of degree at most dj, and ϕ is the change of variables on the transition intervals given explicitly as \begin{align*} \phi \left( s \right) = \begin{cases} t_{1} + \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{IT_{1}},\\[0.3 em] t_{2} - \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{IT_{2}},\\[0.3 em] t_{1} - \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{ST_{1}},\\[0.3 em] t_{2} + \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{ST_{2}}, \end{cases} \end{align*} wherein φ is the affine map \begin{align*} \varphi (s) = \begin{cases} \xi_{1} + \left( \xi_{1}^{\prime}-\xi_{1} \right) \dfrac{s-a_{1}}{b_{1}-a_{1}}, & s \in I_{IT_{1}},\\[1 em] \xi_{2}^{\prime} + \left( \xi_{2}-\xi_{2}^{\prime} \right) \dfrac{s-a_{2}}{b_{2}-a_{2}}, & s \in I_{IT_{2}},\\[1 em] \zeta_{1}^{\prime} + \left( \zeta_{1}-\zeta_{1}^{\prime} \right) \dfrac{s-a_{3}}{b_{3}-a_{3}}, & s \in I_{ST_{1}},\\[1 em] \zeta_{2} + \left( \zeta_{2}^{\prime}-\zeta_{2} \right) \dfrac{s-a_{4}}{b_{4}-a_{4}}, & s \in I_{ST_{2}}, \end{cases} \end{align*} and ψ is the linear map \begin{align*} \psi (s) = -\dfrac{1}{3} \begin{cases} \dfrac{b_{1}-s}{b_{1}-a_{1}}, & s \in I_{IT_{1}},\\[1 em] \dfrac{s-a_{2}}{b_{2}-a_{2}}, & s \in I_{IT_{2}},\\[1 em] \dfrac{s-a_{3}}{b_{3}-a_{3}}, & s \in I_{ST_{1}},\\[1 em] \dfrac{b_{4}-s}{b_{4}-a_{4}}, & s \in I_{ST_{2}}. \end{cases} \end{align*} The change of variables ϕ is constructed so that, for j = 1, 2, 3, 4, the map ϕ : [aj, bj] → [aj, bj] is an orientation preserving diffeomorphism and the exponent ψ of k increases linearly from − 1/3 to 0 as we move away from the shadow boundaries into the illuminated or shadow regions. While on the one hand this choice guarantees that the DoF assigned to the $$\mathscr{O}(k^{-1/3})$$ neighborhoods of the shadow boundaries remains fixed with increasing wavenumber k, and on the other hand it also ensures that the approximation spaces are perfectly adapted to the asymptotic behavior of the solution. Remark 3.3 Note that, by construction, $$\gamma ( \cup _{j=1}^{J} \mathscr{I}_{j} ) = \partial K$$ and the intervals $$\mathscr{I}_{j}$$ intersect either trivially or only at an end point. Therefore, we can clearly identify $$L^{2}\left ( \partial K \right )$$ and $$L^{2} ( \cup _{j=1}^{J} \mathscr{I}_{j} )$$ through the (L−periodic arc length) parametrization γ. This will be the convention we shall follow without any further reference in the rest of the paper. In particular, we shall write η(s, k), ηslow(s, k), ν(s), etc. rather that η(γ(s), k), ηslow(γ(s), k), ν(γ(s)), etc. With the above definitions, the Galerkin formulation of problem (2.1) is equivalent to finding the unique $$\hat{\eta } \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ such that \begin{align} B_{k}(\hat{\mu},\hat{\eta}) = F_{k}(\hat{\mu}), \quad \textrm{for all } \hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}. \end{align} (3.4) While, in this connection, the following theorem constitutes the main result of the paper and provides the approximation properties of the solution of equation (3.4), its proof is deferred to the next section. Remark 3.4 In all the estimates that follow, the notation $$A \lesssim _{n}B$$ (respectively $$A \lesssim _{n,k_{0}} B$$) will mean that 0 ≤ A ≤ cB for all sufficiently large values of k (respectively for all k ≥ k0), where c is a positive constant independent of k and, where relevant, of s, but its value may be different for different values of n. Note that c depends on the geometry of the scatterer, the direction of incidence α and, where relevant, the parameters ξj, ξj′ etc. used to construct the Galerkin approximation spaces. Theorem 3.2 Given k0 > 1 and k ≥ k0, suppose that the sesquilinear form Bk in (2.3) associated with the integral operator $$\mathscr{R}_{k}$$ in (2.1) is continuous with a continuity constant Ck and coercive with a coercivity constant ck. Then, for all nj ∈{0, …, dj + 1} (j = 1, …, J), we have \begin{align} \Vert \eta - \hat{\eta} \Vert_{L^{2}(\partial K)} \lesssim_{n_{1},\ldots,n_{J},k_{0}} \dfrac{C_{k}}{c_{k}} \, k \sum_{j=1}^{J} \dfrac{E(k,j)}{\left( d_{j} \right)^{n_{j}}} \end{align} (3.5) for the Galerkin solution $$\hat{\eta }$$ to (3.4), where \begin{align*} E(k,j) = \begin{cases}{ll} \left( \log k \right)^{n_{j}+1/2}, & j =1,2,3,4 \ (\textrm{transition regions}), \\[0.5em] k^{-1/6},& j =5,6 \, \qquad (\textrm{shadow boundaries}); \end{cases} \end{align*} if J = 8, then \begin{align*} E(k,j) = 1, \qquad j = 7,8 \quad (\textrm{illuminated and shadow regions}). \end{align*} Since the order of magnitude of $$\Vert \eta \Vert _{L^{2}(\partial K)}$$ is k as $$k \to \infty$$ (see e.g. equation (1.15) in Melrose & Taylor (1985) or equation (3.43) in Ecevit & Reitich (2009)), if we assign the same polynomial degree to each interval $$\mathscr{I}_{j}$$, we obtain the following estimate for the relative error. Corollary 3.1 Under the assumptions of Theorem 3.2, if the same polynomial degree d = d1 = … = dJ is used on each interval, then for all n ∈ {0, …, d + 1}, there holds \begin{align} \dfrac{\Vert \eta - \hat{\eta} \Vert_{_{L^{2}(\partial K)}}}{\Vert \eta \Vert_{_{L^{2}(\partial K)}}} \lesssim_{n,k_{0}} \dfrac{C_{k}}{c_{k}} \dfrac{\left( \log k \right)^{n+1/2}}{d^{n}} \end{align} (3.6) for the Galerkin solution $$\hat{\eta }$$ to (3.4). Proof. As was shown in Ecevit & Reitich (2009, Theorem 3.3; see equation (3.43) therein), given a compact subset S of the illuminated region ∂KIL and a wavenumber k0 > 0, there exists a constant C1 depending on k0, α, S, and ∂K such that \begin{align*} \left|\eta^{\textrm{slow}}(s,k) - 2ik \, \alpha \cdot \nu(s)\right| \le C_{1} \end{align*} for all k ≥ k0 and all x ∈ S. Thus, for some constant C2 depending on k0, α and ∂K, we have $$0< C_{2} k \le \Vert \eta \Vert _{L^{2}(\partial K)}$$ for all sufficiently large k. Since $$\Vert \eta \Vert _{L^{2}(\partial K)}$$ depends continuously on k and is never zero for k > 0, by modifying the constant C2, if necessary, we obtain $$0< C_{2} k \le \Vert \eta \Vert _{L^{2}(\partial K)}$$ for all k ≥ k0. Therefore, we can replace k appearing on the right-hand side of (3.5) by $$\Vert \eta \Vert _{L^{2}(\partial K)}$$. Finally, taking d = d1 = … = dJ and n = n1 = … = nJ prints (3.6). Theorem 3.2 and Corollary 3.1 display the convergence characteristics of the Galerkin approximation spaces based on frequency-dependent changes of variables. In order to allow for a direct comparison, here we present a more flexible version of our algorithm in Ecevit & Özen (2017) together with the associated convergence results. To this end, given $$m \in \mathbb{N}$$, 0 ≤ εm < εm−1 < ⋯ < ε1 < 1/3, and constants ξ1, ξ2, ζ1, ζ2 > 0 with t1 − ξ1 < t2 − ξ2 and t2 + ζ2 < L + t1 − ζ1, we define the associated illuminated transition and shadow transition intervals as \begin{align*} I_{IT_{1}} & = [t_{1} + \xi_{1} k^{-1/3+\varepsilon_{m}}, t_{1} + \xi_{1} k^{-1/3+\varepsilon_{1}}], \\ I_{IT_{2}} & = [t_{2} - \xi_{2} k^{-1/3+\varepsilon_{1}}, t_{2} - \xi_{2} k^{-1/3+\varepsilon_{m}}], \\ I_{ST_{1}} & = [t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{1}}, t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{m}}], \\ I_{ST_{2}} & = [t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{m}}, t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{1}}], \end{align*} and, rather than utilizing a change of variables in the transition intervals, we divide each one into m − 1 subintervals by setting \begin{align*} I_{IT_{1}}^{j} & = [t_{1} + \xi_{1} k^{-1/3+\varepsilon_{j+1}}, t_{1} + \xi_{1} k^{-1/3+\varepsilon_{j}}], \\ I_{IT_{2}}^{j} & = [t_{2} - \xi_{2} k^{-1/3+\varepsilon_{j}}, t_{2} - \xi_{2} k^{-1/3+\varepsilon_{j+1}}], \\ I_{ST_{1}}^{j} & = [t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{j}}, t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{j+1}}], \\ I_{ST_{2}}^{j} & = [t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{j+1}}, t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{j}}], \end{align*} for j = 1, …, m − 1. In addition to these 4m − 4 transition intervals, we further define the shadow boundary intervals as \begin{align*} I_{sb:{1}} & = [t_{1} - \zeta_{1} k^{-1/3 +\varepsilon_{m}}, t_{1} + \xi_{1} k^{-1/3+\varepsilon_{m}}], \\ I_{sb:{2}} & = [t_{2} - \xi_{1} k^{-1/3 +\varepsilon_{m}}, t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{m}}], \end{align*} and the illuminated region and shadow region intervals as \begin{align*} I_{IL} & = [t_{1} + \xi_{1}k^{-1/3+\varepsilon_{1}}, t_{2} - \xi_{2}k^{-1/3+\varepsilon_{1}}], \\ I_{DS} & = [t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{1}}, L+t_{1}-\zeta_{1} k^{-1/3+\varepsilon_{1}}]. \end{align*} These give rise to a total of 4m intervals which we shall denote as $$\mathscr{I}_{j}$$ (j = 1, …, 4m). Reasoning as before, we identify $$L^{2} \left ( \partial K \right )$$ with $$L^{2} ( \cup _{j=1}^{4m} \mathscr{I}_{j} )$$. Given $$\mathbf{d} = \left (d_{1},\ldots ,d_{4m} \right ) \in \mathbb{Z}_{+}^{4m}$$, we now define the $$( \left | \mathbf{d} \right | + 4m)-$$dimensional Galerkin approximation space based on algebraic polynomials as \begin{align} \ \end{align} (3.7) The associated Galerkin formulation of (2.1) is to find the unique $$\hat{\eta } \in \mathscr{A}_{\mathbf{d}}$$ such that \begin{align} B_{k}(\hat{\mu},\hat{\eta}) = F_{k}(\hat{\mu}), \quad \textrm{for all } \hat{\mu} \in \mathscr{A}_{\mathbf{d}}. \end{align} (3.8) Incidentally, the Galerkin approximation spaces defined by equation (3.7) provide a more flexible version of those in Ecevit & Özen (2017), since here we allow for different ξ and ζ values rather than the values ξ1 = ξ2 and ζ1 = ζ2 we used in Ecevit & Özen (2017). This clearly renders the new approximation spaces better adapted to different geometries. Theorem 3.3 Suppose that k is sufficiently large and the sesquilinear form Bk in (2.3) associated with the integral operator $$\mathscr{R}_{k}$$ in (2.1) is continuous with a continuity constant Ck and coercive with a coercivity constant ck. Then, for all nj ∈{0, …, dj + 1} (j = 1, …, 4m), we have \begin{align*} \Vert \eta - \hat{\eta} \Vert_{L^{2}(\partial K)} \lesssim_{n_{1},\ldots,n_{4m}} \dfrac{C_{k}}{c_{k}} \, k \, \sum_{j = 1}^{4m} \dfrac{1+E(k,j)}{\left( d_{j} \right)^{n_{j}}} \end{align*} for the Galerkin solution $$\hat{\eta }$$ to (3.8), where \begin{align*} E(k,j) = \begin{cases} k^{-(1+3\varepsilon_{r+1})/2} \left( k^{(\epsilon_{r}-\varepsilon_{r+1})/2} \right)^{n_{j}}\!\!\!\!\!, & j=1,\ldots,4m-4 \, \quad (\textrm{transition regions}), \\[0.5 em] k^{-1/2} \left( k^{\varepsilon_{m}} \right)^{n_{j}}, & j=4m-3,4m-2 \quad (\textrm{shadow boundaries}), \\[0.5 em] k^{-(1+3\varepsilon_{1})/2} \left( k^{(1/3-\varepsilon_{1})/2} \right)^{n_{j}}, & j=4m-1,4m \qquad (\text{illuminated \& shadow reg.}). \end{cases} \end{align*} The proof of Theorem 3.3 is similar to that of Theorem 1 in Ecevit & Özen (2017) and is therefore skipped. On the other hand, exactly as Theorem 1 therein implies Corollary 1 in Ecevit & Özen (2017), Theorem 3.3 above yields the following result. Corollary 3.2 Under the assumptions of Theorem 3.3, if εj are chosen as \begin{align*} \varepsilon_{j} = \dfrac{1}{3} \, \dfrac{2m-2j+1}{2m+1}, \qquad j = 1,\ldots,m, \end{align*} and the same polynomial degree d = d1 = … = d4m is used on each interval, then for all n ∈{0, …, d + 1}, there holds \begin{align} \dfrac{\Vert \eta - \hat{\eta} \Vert_{_{L^{2}(\partial K)}}}{\Vert \eta \Vert_{_{L^{2}(\partial K)}}} \lesssim_{n} \dfrac{C_{k}}{c_{k}} \, m \, \dfrac{1+ k^{-\frac{1}{2}} \left( k^{\frac{1}{6m+3}} \right)^{n}}{d^{n}} \end{align} (3.9) for the Galerkin solution $$\hat{\eta }$$ to (3.8). Remark 3.5 Let us say that ‘f increases in proportion to g’ to mean that f and g are positive functions of k and f/g is bounded from above and below by positive constants independent on k. As mentioned in Remark 2.1, $$C_{k}/c_{k} = \mathscr{O}(k^{\delta })$$ as $$k \to \infty$$ with δ = 1/3 for CFIE and SCIE. In this case, for the estimate in Corollary 3.2, if m increases in proportion to $$\log k$$—which implies that $$k^{1/(6m+3)} = \exp (\log k/(6m+3))$$ is bounded independently of k—and the local polynomial degree d increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{\varepsilon _{2}}$$ with ε1, ε2 > 0 so that the total number of DoF, namely 4m(d + 1), increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{1+\varepsilon _{2}}$$, then the right-hand side of (3.9) is $$\mathscr{O}(k^{\delta -n\varepsilon _{1}} \, (\log k)^{1-n\varepsilon _{2}}$$) as $$k \to \infty$$. For the estimate (3.6) in Corollary 3.1, if the local polynomial degree d increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{1+\varepsilon _{2}}$$ which implies that the total number of DoF, namely J(d + 1) with J = 6 or J = 8, increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{1+\varepsilon _{2}}$$ as above, then the right-hand side of (3.6) is $$\mathscr{O}(k^{\delta - n\varepsilon _{1}} \, (\log k)^{1/2-n\varepsilon _{2}})$$ as $$k \to \infty$$. In this sense, when compared with the frequency-adapted Galerkin BEM, the method based on frequency-dependent changes of variables provides an improvement of $$\sqrt{\log k}$$ in numerical accuracy as $$k \to \infty$$. Since n here can be arbitrarily large and $$k^{\delta _{1} - n\varepsilon _{1}} \, (\log k)^{\delta _{2}-n\varepsilon _{2}} = \mathscr{O}(1)$$ as $$k \to \infty$$ provided that either ε1 ≥ δ1/n and ε2 ≥ δ2/n or ε1 > δ1/n, it follows that, for any ε > 0, increasing the total number of DoF associated with either of the Galerkin schemes as $$\mathscr{O}(k^{\varepsilon })$$ is sufficient to obtain any prescribed accuracy independent on frequency. 4. Error analysis In this section we present the proof of Theorem 3.2. In light of inequality (2.5), it is sufficient to estimate \begin{align*} \inf_{\hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}} \Vert \eta- \hat{\mu} \Vert_{L^{2}(\partial K)}. \end{align*} To this end, we make use of the following classical result from approximation theory. Theorem 4.1 (Best approximation by algebraic polynomials (Schwab, 1998)) Given an interval $$\mathscr{I} = (a,b)$$ and $$n \in \mathbb{Z}_{+}$$, introduce the seminorms (for suitable f) by \begin{align} \left|\, f \right|{}_{n,\mathscr{I}} = \left[{\int_{a}^{b}} \left| D^{n} f(s) \right|{}^{2} \left( s-a \right)^{n} \left( b-s \right)^{n} \textrm{d}s \right]^{1/2} . \end{align} (4.1) Then, for all n ∈{0, …, d + 1}, there holds \begin{align*} \inf_{p \in \mathbb{P}_{d}} \Vert\, f-p \Vert \lesssim_{n} \left|\, f \right|{}_{n,\mathscr{I}} \, d^{-n}. \end{align*} Use of Theorem 4.1, in turn, requires the knowledge of derivative estimates of ηslow which we present next. In the rest of the paper we use the convention that an empty sum is zero. Theorem 4.2 Given k0 > 0, there holds \begin{align} \left| {D_{s}^{n}} \eta^{\textrm{slow}}(s,k) \right| \lesssim_{n,k_{0}} k + \sum\limits_{m=4}^{n+2} \left( k^{-1/3} + \left| w(s) \right| \right)^{-m} \end{align} (4.2) for all $$n \in \mathbb{Z}_{+}$$ and all k ≥ k0. Here w(s) = (s − t1)(t2 − s). Proof. The same estimate is shown to hold for all sufficiently large k in Domínguez et al. (2007, Theorem 5.4). Since $${D_{s}^{n}} \eta ^{\textrm{slow}}(s,k)$$ depends continuously on s and k, the result follows. Remark 4.1 Let us note that estimate (4.2) is valid over the entire boundary, but it is not sharp in the shadow region well away from the shadow boundary points, where the solution decays exponentially with increasing wavenumber k and the geodesic distance to the shadow boundary points (see Domínguez et al. (2007), Asheim & Huybrechs (2014), and the references therein). Nevertheless, as in Ecevit & Özen (2017), it allows us to eliminate, through utilization of additional approximation spaces in the shadow region, the asymptotic terms appearing in the error analysis due to approximation by zero in the deep shadow region (compare the estimates in Theorems 3.2 and 3.3 with Domínguez et al. (2007, Theorem 6.7)). Further, as we discussed in detail in Ecevit & Özen (2017), the benefit of this utilization is directly apparent from the significantly reduced condition numbers of the Galerkin linear systems. We continue with the derivation of estimates on the derivatives of the change of variables ϕ on the transition interval $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4). Proposition 4.3 Given k0 > 1, there holds \begin{align*} \left| {D^{n}_{s}} \phi \right| \lesssim_{n,k_{0}} \left( \log k \right)^{n} k^{\psi} \quad \textrm{on } \mathscr{I}_{j} \quad (\ j=1,2,3,4), \end{align*} for all $$n \in \mathbb{N}$$ and all k ≥ k0. Proof. Since the proof is similar for j = 1, 2, 3, 4, we concentrate on the case j = 1. Now since φ′′ = ψ′′ = 0, direct computations entail \begin{align} {D^{n}_{s}} \phi = \sum_{j=0}^{n} \binom{n}{j} \, D^{n-j}_{s} \varphi \ {D^{j}_{s}} k^{\psi} = \varphi \ {D^{n}_{s}} k^{\psi} + n \ {D^{1}_{s}} \varphi \ D^{n-1}_{s} k^{\psi}, \qquad n \ge 1, \end{align} (4.3) and \begin{align} {D^{n}_{s}} k^{\psi} = \left( {D^{1}_{s}} \psi \right)^{n} \left( \log k \right)^{n} k^{\psi}, \qquad n \ge 0. \end{align} (4.4) Using (4.4) in (4.3), we obtain \begin{align} {D^{n}_{s}} \phi = \left( \varphi \ {D^{1}_{s}} \psi \ \log k + n \ {D^{1}_{s}} \varphi \right) \left( {D^{1}_{s}} \psi \right)^{n-1} \left( \log k \right)^{n-1} k^{\psi}, \qquad n \ge 1. \end{align} (4.5) Since, for k ≥ k0 > 1, \begin{align*} \dfrac{1}{\left| b_{1}-a_{1} \right|} = \dfrac{1}{\xi_{1}^{\prime} - \xi_{1}k^{-1/3}} \le \dfrac{1}{\xi_{1}^{\prime} - \xi_{1}k_{0}^{-1/3}} \lesssim_{k_{0}} 1 \end{align*} it follows for $$s \in \mathscr{I}_{1} = I_{IT_{1}}$$ that \begin{align*} \left| {D^{1}_{s}} \psi \right| = \dfrac{1}{3} \, \dfrac{1}{b_{1}-a_{1}} \lesssim_{k_{0}} 1 \qquad \textrm{and} \qquad \left| {D^{1}_{s}} \varphi \right| = \dfrac{\xi_{1}^{\prime} - \xi_{1}}{b_{1}-a_{1}} \lesssim_{k_{0}} 1 \end{align*} and, clearly, $$\xi _{1} \le \varphi \le \xi _{1}^{\prime }$$. Use of these inequalities in (4.5) yields the desired result. Next we combine Theorem 4.2 and Proposition 4.3 to derive estimates on the derivatives of the composition ηslow ∘ ϕ. Proposition 4.4 Given k0 > 0, there holds \begin{align*} \left| {D_{s}^{n}} (\eta^{slow} \circ \phi) \right| \lesssim_{n,k_{0}} k \left( \log k \right)^{n} \quad \textrm{on } \mathscr{I}_{j} \quad \left( j=1,2,3,4 \right)\!, \end{align*} for all $$n \in \mathbb{N}$$ and all k ≥ k0. Proof. We fix k ≥ k0 > 1 and the interval $$\mathscr{I}_{j}$$ (j = 1, …, 4). Faá Di Bruno’s formula for the derivatives of a composition states \begin{align*} D^{n} \left(\, f \circ g \right) \left( t \right) = \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} (D^{m} f) (g(t)) \prod_{\ell=1}^{n} \dfrac{\ell}{m_{\ell}!} \left( \dfrac{D^{\ell}g(t)}{\ell!} \right)^{m_{\ell}}, \end{align*} where $$\mathscr{F}_{n} = \{ (m_{1},\ldots ,m_{n}) \in \mathbb{Z}^{n}_{+} : n = \sum _{\ell =1}^{n} \ell m_{\ell } \}$$; here $$m = \sum _{\ell =1}^{n} m_{\ell }$$. This yields \begin{align*} \left| {D^{n}_{s}} \left( \eta^{\textrm{slow}} \circ \phi \right) \right| \lesssim_{n} \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| \prod_{\ell=1}^{n} \left| D^{\ell}_{s} \phi \right|{}^{m_{\ell}} \end{align*} so that an appeal to Proposition 4.3 entails \begin{align*} \left| {D^{n}_{s}} \left( \eta^{\textrm{slow}} \circ \phi \right) \right| \lesssim_{n,k_{0}} \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| \prod_{\ell=1}^{n} \left( \left( \log k \right)^{\ell} k^{\psi} \right)^{m_{\ell}}. \end{align*} Since $$n = \sum _{\ell =1}^{n} \ell m_{\ell }$$ and $$m = \sum _{\ell =1}^{n} m_{\ell }$$, we therefore obtain \begin{align*} \left| {D^{n}_{s}} \left( \eta^{\textrm{slow}} \circ \phi \right) \right| & \lesssim_{n,k_{0}} \left( \log k \right)^{n} \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| k^{m\psi} \\ & \lesssim_{n,k_{0}} \left( \log k \right)^{n} \sum_{m=0}^{n} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| k^{m\psi}. \end{align*} It is hence sufficient to show, for $$m \in \mathbb{Z}_{+}$$, that \begin{align*} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| k^{m\psi} \lesssim_{m,k_{0}} k. \end{align*} To this end, we note that if 0 ≤ ℓ ≤ m, then \begin{align*} \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-\ell} & = \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-m} \left( k^{-1/3} + \left| w (\phi) \right| \right)^{ m-\ell} \\ & \leq \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-m} \left( k_{0}^{-1/3} + L^{2} \right)^{m-\ell} \\ &\lesssim_{m,k_{0}} \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-m} \end{align*} so that an appeal to Theorem 4.2 yields \begin{align*} \left| {D_{s}^{m}} \eta^{\textrm{slow}}(\phi) \right| k^{m\psi} & \lesssim_{m,k_{0}} \left[ k + \sum_{\ell=4}^{m+2} \left( k^{-1/3} + \left| w(\phi) \right| \right)^{-\ell} \right] k^{m\psi} \\ & \lesssim_{m,k_{0}} \left[ k + \left( k^{-1/3} + \left| w(\phi) \right| \right)^{-\left( m+2 \right)} \right] k^{m\psi} \\ & \lesssim_{m,k_{0}} k + \left( \dfrac{k^{\psi}}{k^{-1/3} + \left| w(\phi) \right|} \right)^{m} \left( k^{-1/3} + \left| w(\phi) \right| \right)^{-2} \\ & \lesssim_{m,k_{0}} k + \left( \dfrac{k^{\psi}}{\left| w(\phi) \right|} \right)^{m} k^{2/3}, \end{align*} where, in the third inequality, we used that ψ ≤ 0. Thus, it is now enough to show that the quotient $$k^{\psi }/\left | w(\phi ) \right |$$ is bounded by a constant independent of k. This estimation is similar on each of the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4) and we focus on $$\mathscr{I}_{1}$$. Indeed, on $$\mathscr{I}_{1}=I_{IT_{1}}$$, we have ϕ − t1 = φkψ, φ ≥ ξ1 > 0 and $$t_{2} - \phi \ge t_{2} - \left ( t_{1} + \xi _{1}^{\prime } \right ) \ge \xi _{2}^{\prime }> 0$$  so that \begin{align*} \dfrac{k^{\psi}}{\left| w(\phi) \right|} = \dfrac{k^{\psi}}{\left( \phi -t_{1} \right) \left( t_{2} - \phi \right)} = \dfrac{k^{\psi}}{\varphi k^{\psi} \left( t_{2} - \phi \right)} \le \dfrac{1}{\xi_{1} \, \xi_{2}^{\prime}}. \end{align*} This finishes the proof. Next, we estimate the seminorms (4.1) for the composition ηslow ∘ ϕ on the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4). Corollary 4.1 On the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4), given k0 > 1, there holds \begin{align*} \left| \eta^{\textrm{slow}} \circ \phi \right|{}_{n,\mathscr{I}_{j}} \lesssim_{n,k_{0}} k \left( \log k \right)^{n} \end{align*} for all $$n \in \mathbb{N}$$ and all k ≥ k0. Proof. On account of Proposition 4.4, we estimate for j = 1, 2, 3, 4 \begin{align*} \left| \eta^{\textrm{slow}} \circ \phi \right|{}^{2}_{n,\mathscr{I}_{j}} & = \int_{a_{j}}^{b_{j}} \left| {D^{n}_{s}}\left( \eta^{\textrm{slow}} \circ \phi \right) (s) \right|{}^{2} \left( s - a_{j} \right)^{n} \left( b_{j} -s \right)^{n} \textrm{d}s \\ & \lesssim_{n,k_{0}} k^{2} \left( \log k \right)^{2n} \int_{a_{j}}^{b_{j}} \left( s - a_{j} \right)^{n} \left( b_{j} -s \right)^{n} \textrm{d}s \\ & \lesssim_{n,k_{0}} k^{2} \left( \log k \right)^{2n}, \end{align*} where we used that 0 < bj − aj < L. Thus, the result. We are now ready to prove Theorem 3.2. Proof. (of Theorem 3.2): Céa’s lemma (cf. inequality (2.5)) implies \begin{align} \Vert \eta -\hat{\eta} \Vert_{L^{2} \left( \partial K \right)} \leq \dfrac{C_{k}}{c_{k}} \inf_{\hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}} \Vert \eta- \hat{\mu} \Vert_{L^{2} \left( \partial K \right)} \end{align} (4.6) for the unique solution $$\hat{\eta }$$ of the Galerkin formulation (3.4). As we identify $$L^{2} \left ( \partial K \right )$$ with $$L^{2} ( \cup _{j=1}^{J} \mathscr{I}_{j} )$$ through the L−periodic arc length parameterization γ, there holds \begin{align*} \Vert \eta- \hat{\mu} \Vert_{L^{2}(\partial K)} = \Vert \eta- \hat{\mu} \Vert_{L^{2} \left( \cup_{j=1}^{J} \mathscr{I}_{j} \right)} \le \sum_{j=1}^{J} \Vert \eta- \hat{\mu} \Vert_{L^{2} \left( \mathscr{I}_{j} \right)} \end{align*} for any $$\hat{\mu } \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$. Accordingly, the very definition of Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ entails \begin{align} \inf_{\hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}} \Vert \eta- \hat{\mu} \Vert_{L^{2}(\partial K)} \le \sum_{j=1}^{4} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} - p \circ \phi^{-1} \Vert_{L^{2}(\mathscr{I}_{j})} + \sum_{j=5}^{J} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} - p \Vert_{L^{2}(\mathscr{I}_{j})}. \end{align} (4.7) On the other hand, utilizing the change of variables ϕ on the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4), for any $$p \in \mathbb{P}_{d_{j}}$$, we have \begin{align*} \Vert \eta^{\textrm{slow}} - p \circ \phi^{-1} \Vert^{2}_{L^{2}(\mathscr{I}_{j})} & = \int_{a_{j}}^{b_{j}} \left| \left( \eta^{\textrm{slow}} - p \circ \phi^{-1}\right)(s) \right|{}^{2} \textrm{d}s \\ & = \int_{a_{j}}^{b_{j}} \left| \left( \eta^{\textrm{slow}} \circ \phi - p \right)(s) \right|{}^{2} {D^{1}_{s}}\phi(s) \, \textrm{d}s \\ & \lesssim_{k_{0}} \log k \ \Vert \eta^{\textrm{slow}} \circ \phi - p \Vert^{2}_{L^{2}(\mathscr{I}_{j})}, \end{align*} where we used Proposition 4.3 in conjunction with the fact that kψ < 1. Combining this last estimate with (4.6) and (4.7), we deduce \begin{align*} \Vert \eta -\hat{\eta} \Vert_{L^{2} \left( \partial K \right)} \lesssim_{k_{0}} \dfrac{C_{k}}{c_{k}} \left\{ \sum_{j=1}^{4} \left( \log k \right)^{1/2} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} \circ \phi - p \Vert_{L^{2}(\mathscr{I}_{j})} + \sum_{j=5}^{J} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} - p \Vert_{L^{2}(\mathscr{I}_{j})} \right\} \end{align*} and this, on account of Theorem 4.1, implies \begin{align*} \Vert \eta -\hat{\eta} \Vert_{L^{2} \left( \partial K \right)} \lesssim_{n_{1},\ldots,n_{J},k_{0}} \dfrac{C_{k}}{c_{k}} \left\{ \sum_{j=1}^{4} \left( \log k \right)^{1/2} \left| \eta^{\textrm{slow}} \circ \phi \right|{}_{n_{j},\mathscr{I}_{j}} d_{j}^{-n_{j}} + \sum_{j=5}^{J} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}} d_{j}^{-n_{j}} \right\}\!\!. \end{align*} Therefore, to complete the proof, it suffices to show that \begin{align} \left| \eta^{\textrm{slow}} \circ \phi \right|{}_{n_{j},\mathscr{I}_{j}} \lesssim_{n_{j},k_{0}} k \left( \log k \right)^{n_{j}}, \qquad j = 1,2,3,4, \end{align} (4.8) and \begin{align} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}} \lesssim_{n_{j},k_{0}} k \, k^{-1/6}, \qquad j = 5,6, \end{align} (4.9) and (if J = 8) \begin{align} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}} \lesssim_{n_{j},k_{0}} k, \qquad j = 7,8. \end{align} (4.10) While the estimates in (4.8) are given by Corollary 4.1, for the shadow boundary intervals $$\mathscr{I}_{j}$$ (j = 5, 6) we use Theorem 4.2 to deduce \begin{align*} \left| D_{s}^{n_{j}} \eta^{\textrm{slow}}(s,k) \right| \lesssim_{n_{j},k_{0}} k + \sum_{m=4}^{n_{j}+2} \left( k^{-1/3} + \left| w(s) \right| \right)^{-m} \lesssim_{n_{j},k_{0}} k + k^{(n_{j}+2)/3}; \end{align*} this implies \begin{align*} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}}^{2} & = \int_{a_{j}}^{b_{j}} \left| D^{n_{j}}_{s} \eta^{\textrm{slow}}(s) \right|{}^{2} \left( s -a_{j} \right)^{n_{j}} \left( b_{j}-s \right)^{n_{j}} \textrm{d}s \\ & \lesssim_{n_{j},k_{0}} \left( k + k^{(n_{j}+2)/3} \right)^{2} \left( b_{j} -a_{j} \right)^{2n_{j}+1} \\ & \lesssim_{n_{j},k_{0}} \left( k + k^{(n_{j}+2)/3} \right)^{2} \left( k^{-1/3} \right)^{2n_{j}+1} \\ & \lesssim_{n_{j},k_{0}} \left( k \, k^{-1/6} \right)^{2}, \end{align*} which justifies the estimates in (4.9). This completes the proof when J = 6. When J = 8, for the illuminated and shadow region intervals $$\mathscr{I}_{j}$$ (j = 7, 8), we use Theorem 4.2 to estimate \begin{align*} \left| D_{s}^{n_{j}} \eta^{\textrm{slow}}(s,k) \right| \lesssim_{n_{j},k_{0}} k + \sum_{m=4}^{n_{j}+2} \left( k^{-1/3} + \left| w(s) \right| \right)^{-m} \lesssim_{n_{j},k_{0}} k + \sum_{m=4}^{n_{j}+2} \left| w(s) \right|{}^{-m} \lesssim_{n_{j},k_{0}} k \end{align*} so that \begin{align*} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}}^{2} & = \int_{a_{j}}^{b_{j}} \left| D^{n_{j}}_{s} \eta^{\textrm{slow}}(s) \right|{}^{2} \left( s -a_{j} \right)^{n_{j}} \left( b_{j}-s \right)^{n_{j}} \textrm{d}s \\ & \lesssim_{n_{j},k_{0}} k^{2} \left( b_{j} -a_{j} \right)^{2n_{j}+1} \\ & \lesssim_{n_{j},k_{0}} k^{2}, \end{align*} which verifies (4.10). This finishes the proof. 5. Numerical tests In this section, we first discuss the implementation details related with the Galerkin BEMs based on frequency-dependent changes of variables. Next we present numerical tests exhibiting their performance in comparison with the frequency-adapted Galerkin BEMs. Finally, we discuss the accuracy of numerical solutions in the shadow region and compare the results with those in Asheim & Huybrechs (2014). 5.1. Implementation details As we mentioned, the developments central to this paper are independent of the integral equation used. In order to allow a simple performance comparison with the aforementioned algorithms, we base our numerical implementations on the CFIE wherein the integral operator and the right-hand side are given by \begin{align*} \mathscr{R}_{k} = \dfrac{1}{2} \, I + \mathscr{D} - ik \mathscr{S} \quad \textrm{and} \quad f_{k} = \dfrac{\partial u^{\textrm{inc}}}{\partial \nu} - ik u^{\textrm{inc}}. \end{align*} Here, $$\mathscr{S}$$ is the acoustic single-layer integral operator and $$\mathscr{D}$$ is its normal derivative, and they are defined by \begin{align*} \mathscr{S} \eta(x) =&\, \int_{\partial K} \varPhi(x,y) \, \eta(y) \, \textrm{d}s(y), \qquad x \in \partial K,\\ \mathscr{D} \eta(x) =&\, \int_{\partial K} \dfrac{ \partial \varPhi(x,y)}{\partial \nu(x)} \, \eta(y) \, \textrm{d}s(y), \qquad x \in \partial K, \end{align*} where ν(x) is the outward unit normal to ∂K. In this case, the conjugates of the entries of the Galerkin matrices arising in connection with the approximation spaces $$\mathscr{A}_{d}$$ or $$\mathscr{A}_{d}^{\mathscr{C}}$$ are of the form \begin{align} \overline{B_{k}(\hat{\mu}_{j},\hat{\mu}_{j^{\prime}})} = \dfrac{1}{2} \int_{\mathscr{I}_{j} \cap \mathscr{I}_{j^{\prime}}} \overline{\hat{\mu}_{j}(s)} \, \hat{\mu}_{j^{\prime}}(s) \, \textrm{d}s + \iint_{\mathscr{I}_{j} \times \mathscr{I}_{j^{\prime}}} F(s,t,k) \, \overline{\hat{\mu}_{j}(s)} \, \hat{\mu}_{j^{\prime}}(t) \, \textrm{d}t \, \textrm{d}s, \end{align} (5.1) where $$\hat{\mu }_{j}$$ and $$\hat{\mu }_{j^{\prime }}$$ are elements of $$\mathscr{A}_{d}$$ or $$\mathscr{A}_{d}^{\mathscr{C}}$$ supported on $$\mathscr{I}_{j}$$ and $$\mathscr{I}_{j^{\prime }}$$, and F(s, t, k) = ∂ν(s)Φ(γ(s), γ(t)) − ikΦ(γ(s), γ(t)). Note that $$d H^{(1)}_{0}(z)/dz = - H^{(1)}_{1}(z)$$ and \begin{align*} H_{r}^{(1)}(z) = e^{iz} \left\{ \sum_{\ell=0}^{N} \dfrac{c_{r,\ell}}{z^{\ell+1/2}} + \theta_{r,N}(z) \, \dfrac{c_{r,N+1}}{z^{N+3/2}}\right\}\!\!,\qquad r,N \ge 0 \end{align*} with $$\left | \theta _{r,N}(z) \right | < 1$$ for all z ≫ 1 for some constants cr, ℓ (Gradshteyn & Ryzhik, 2000). Therefore, setting \begin{align*} G(s,t,k) = e^{-ik|\gamma(s)-\gamma(t)|} F(s,t,k) \qquad \textrm{and} \qquad \varPsi(s,t) = \alpha \cdot \left(\gamma(t)-\gamma(s)\right)+\left|\gamma(s)-\gamma(t)\right| \end{align*} and using p and q to denote polynomials, the right-hand side of (5.1) takes on the generic forms \begin{align} \dfrac{1}{2} \int_{\mathscr{I}_{j} \cap \mathscr{I}_{j^{\prime}}} p(s) \, q(s) \, \textrm{d}s + \iint_{\mathscr{I}_{j} \times \mathscr{I}_{j^{\prime}}} e^{ik \, \varPsi(s,t)} \, G(s,t,k) \, p(s) \, q(t) \, \textrm{d}t \, \textrm{d}s \end{align} (5.2) for $$\mathscr{A}_{d}$$, and \begin{align} \dfrac{1}{2} \int_{\mathscr{I}_{j} \cap \mathscr{I}_{j^{\prime}}} p\left(\phi_{j}^{-1}(s)\right) \, q\left(\phi_{j^{\prime}}^{-1}(s)\right) \, \textrm{d}s + \iint_{\mathscr{I}_{j} \times \mathscr{I}_{j^{\prime}}} e^{ik \, \varPsi(s,t)} \, G(s,t,k) \, p\left(\phi_{j}^{-1}(s)\right) \, q\left(\phi_{j^{\prime}}^{-1}(t)\right) \, \textrm{d}t \, \textrm{d}s \end{align} (5.3) for $$\mathscr{A}^{\mathscr{C}}_{d}$$. Here ϕj = ϕ if $$\mathscr{I}_{j}$$ is an interval of change of variables, and it is the identity map otherwise (the same comment applies to $$\phi _{j^{\prime }}$$). While the one-dimensional integrals in (5.2) and (5.3) can be evaluated through any suitable quadrature rule, the highly-oscillatory double integrals therein require special treatment. Indeed, Kim (2012) provides a detailed theoretical discussion on the efficient evaluation of the double integrals in (5.2) by analyzing the regularity of the function G and the stationary points of the phase function Ψ for compact, smooth, and strictly convex obstacles. The discussions there are confined to the method in Domínguez et al. (2007) and thus limited to the use of three regions only, but provide an excellent source for the extension and numerical analysis of the methodologies discussed there to the more general case here involving 4m regions. Concerning the double integrals in (5.3) one would be inclined to substitute $$u = \phi _{j}^{-1}(s)$$ and $$v = \phi _{j^{\prime }}^{-1}(t)$$, but this transforms the phase functions into $$\varPsi (\phi _{j}(u),\phi _{j^{\prime \prime }}(v))$$. In this case the associated double integrals cannot be evaluated efficiently through standard quadrature rules for highly oscillatory integrals since they are based on the assumption that the phase functions do not depend on k (here ϕj and $$\phi _{j^{\prime }}$$ depend on k if they correspond to a transition interval). It is therefore reasonable to invert these functions so as to preserve the form of the phase functions. In this connection, observe that if $$\mathscr{I}_{j} = I_{\textrm{IT}_{1}}$$, then ϕj can be expressed as \begin{align} \phi_{j}(s) = t_{1} + \left(\xi_{1} + \left(\xi^{\prime}_{1}-\xi_{1}\right) s^{\prime}\right) \, k^{(s^{\prime}-1)/3}, \qquad s \in [a_{1},b_{1}] = \left[t_{1} + \xi_{1} k^{-1/3}, t_{1} + \xi_{1}^{\prime}\right] \end{align} (5.4) with $$0 \le s^{\prime } = \frac{s-a_{1}}{b_{1}-a_{1}} \le 1$$, and ϕj admits similar representations in the other transition intervals. It follows from (5.4) that ϕj does not oscillate with increasing k and can be inverted numerically. Here we suggest a more direct approach observing first that when $$\xi _{1} = \xi ^{\prime }_{1}$$ \begin{align*} \phi_{j}^{-1}(s) = a_{1} + (b_{1}-a_{1}) \left( 3\dfrac{\log(s-t_{1})-\log \xi_{1}}{\log k} + 1 \right)\!. \end{align*} Otherwise, $$\xi _{1} < \xi ^{\prime }_{1}$$ and equation (5.4) can be rewritten as \begin{align} \dfrac{1}{3} \dfrac{\phi_{j}(s)-t_{1}}{\xi^{\prime}_{1}-\xi_{1}} = \left( \dfrac{1}{3} \dfrac{\xi^{\prime}_{1}}{\xi^{\prime}_{1}-\xi_{1}} + u \right) k^{u}, \qquad s \in [a_{1},b_{1}] \end{align} (5.5) with $$-\frac{1}{3}\le u = \frac{s^{\prime }-1}{3} \le 0$$. In this case, since equations (5.4) and (5.5) are related by affine transformations, the problem of numerically inverting ϕj is equivalent to that of f(u) = (c + u)ku on the interval [−1/3, 0] where c > 1/3 is a constant. Using the affine transformations v = c + u and $$w = v \log k$$, this reduces to inverting g(w) = wew on $$[(c-1/3) \log k, c \log k]$$. The inverse of g is known as the Lambert or omega function and can be evaluated efficiently (Corless et al., 1996). With this information, the highly oscillatory double integrals in (5.3) can be efficiently evaluated based on extensions of the methodologies in Kim (2012) to J = 6 and J = 8 regions. Further, the explicit forms of the functions ϕj (as in (5.4)) also allow for the extension of the regularity analysis and numerical quadrature discussed therein to cover the highly oscillatory double integrals in (5.3). Leaving these details for future work, here we evaluate the double integrals appearing as the entries of Galerkin matrices using a combination of the Nyström and trapezoidal rule discretizations (Colton & Kress, 1992) for the inner integrals and the trapezoidal rule for the outer integrals. This latter rule is also used in the evaluation of one-dimensional integrals. In order to preserve the high-order approximation properties of these numerical integration rules for smooth and periodic integrands, as in Ecevit & Özen (2017), we additionally utilize a smooth partition of unity confined to the regions on the boundary of the scatterer described in Section 3. In each case, based on our experience in Ecevit & Özen (2017), the number of discretization points is chosen approximately as 10 to 12 points per wave length. An important component of our algorithms relates to the choice of bases for the spaces of algebraic polynomials and, in order to minimize the numerical instabilities arising from the use of high-degree polynomials (Ecevit & Özen, 2017), on any generic interval $$\mathscr{I} = [a,b]$$, we use the bases {ρr : r = 0, …, d} and {(ρ∘ϕ−1)r : r = 0, …, d} for $$\mathbb{P}_{d}$$ and $$\mathbb{P}_{d} \circ \phi ^{-1}$$, respectively, where ρ is the affine function that maps the interval $$\mathscr{I}$$ onto [−1, 1]. In the same vein, the choice of the parameters ξ, ξ′, ζ, ζ′ appearing in the definitions of the intervals $$I_{IT_{j}},I_{ST_{j}},I_{sb:{j}}$$ (j = 1, 2) and IIT, IDS related with the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces is of great importance, since a random choice may result in a loss of accuracy due to poor resolution of the boundary layers in the solution. We therefore optimize these parameters for a small wavenumber (k = 10 for all the implementations here) through a simple iterative procedure, and use these values for all larger wavenumbers. For instance, when J = 6, we take $$\xi _{1}^{\prime } = \xi _{2}^{\prime }$$, $$\zeta _{1}^{\prime }=\zeta _{2}^{\prime }$$ and initially require that $$\alpha \cdot \gamma (\xi _{1}^{\prime }) = -1$$ and $$\alpha \cdot \gamma (\zeta _{1}^{\prime }) = 1$$. We then take ξ1 to be the mid-point between t1 and $$\xi _{1}^{\prime }$$, and change it in small increments until the local error in IT1 ∪ SB1 is minimized. We treat the triplet $$(t_{2},\xi _{2},\xi _{2}^{\prime })$$ similarly. Finally, we fix ξ1 and ξ2, and change $$\xi _{1}^{\prime }=\xi _{2}^{\prime }$$ in small increments until the error in IL ∪ IT1 ∪ IT2 ∪ SB1 ∪ SB2 is minimized. The optimization of ζ parameters is realized similarly. The computed values are taken as the initial guess for the next iterate, and the procedure is repeated with smaller increments, with step size half that of the previous one, until the variations in the global error stabilizes. Following the prescriptions in Ecevit & Özen (2017, §4.2), then we introduce a smooth partition of unity confined to the regions$$I_{IT_{j}},I_{ST_{j}}, I_{sb:{j}}$$(j = 1, 2) and IIT, IDS, and optimize the shapes of hat functions therein using a similar iterative procedure. We follow a similar procedure for the optimizaton of the parameters appearing in connection with the $$\mathscr{A}_{\mathbf{d}}$$ spaces. For further details, we refer to Eruslu (2015). Fig. 3. View largeDownload slide Single-scattering geometries and associated incidence directions α. Fig. 3. View largeDownload slide Single-scattering geometries and associated incidence directions α. 5.2. Performance comparison of $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces Here we present numerical tests comparing the performances of frequency-adapted Galerkin BEMs and Galerkin BEMs based on frequency-dependent changes of variables. For the former, the Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}$$ are generated by taking m = 2 and assigning the same polynomial degree d = 3, 6, 9, 12, 15 to each one of the eight direct summands which give rise to the number of DoF 8(d + 1) = 32, 56, 80, 104, 128. Let us mention in this connection that the numerical tests we presented in Ecevit & Özen (2017) are restricted to the choices m = 1 and m = 2 (with some variations) for the very reason that utilization of larger values of m is simply impractical from an implementation point of view. On the other hand, these tests have displayed that the condition numbers of the Galerkin linear systems corresponding to m = 2 are significantly smaller when compared to m = 1, and this enhances numerical stability with increasing polynomial degrees and thereby allows for more accurate approximations. For the latter, the Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ are generated by taking J = 6 and assigning the same polynomial degree d = 4, 8, 12, 16, 20 to each one of the six direct summands resulting in the number of DoF 6(d + 1) = 30, 54, 78, 102, 126. We consider three different single-scattering geometries (see Fig. 3) consisting of the unit circle, the ellipse with major/minor axes (aligned with the x/y-axes) of 2/1, and the kite shaped obstacle given parametrically as $$\{(\cos t + 0.65 \cos 2t - 0.65, 1.5 \sin t) \, : \, t \in [0,2\pi ]\}$$. The unit circle is the standard test example for single-scattering solvers since circles are the only two-dimensional obstacles for which explicit solutions are available through a straightforward Fourier analysis (see equation (5.7) below). As for the ellipse and kite shaped obstacles, we compare the outcome of our numerical implementations with highly accurate reference solutions obtained by a combination of the Nyström and trapezoidal rule discretizations applied to the CFIE (Colton & Kress, 1992). Let us note that the kite shaped obstacle and the associated direction of incidence is chosen to display that the methodologies discussed in this paper can be applied to more general single-scattering problems involving nonconvex smooth scatterers, provided that there are exactly two shadow boundaries and no multiple reflections. The numerical results in Table 1 display the DoF vs. the relative L2 error \begin{align} \dfrac{\Vert \eta - \hat{\eta} \Vert_{L^{2}(\partial K)}}{\Vert \eta \Vert_{L^{2}(\partial K)}} \end{align} (5.6) related with the approximation spaces $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ for the wavenumbers k = 50, 100, 200, 400, 800. For any fixed value of k, the table depicts for both $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces that the error decreases with increasing DoF. The data corresponding to the unit circle shows that, for fixed DoF, the approximation errors associated with the $$\mathscr{A}_{\mathbf{d}}$$ spaces vary depending on the value of k, and those related with the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces tend to decrease with increasing k. Further, $$\mathscr{A}_{\mathbf{d}}$$ spaces provide better approximations for smaller values of k, and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces perform better for larger values of k (when DoF ≥ 54). In the case of the ellipse, for fixed DoF, the approximation errors corresponding to the $$\mathscr{A}_{\mathbf{d}}$$ spaces increase with increasing k, and those related with the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces decrease with increasing k; and when DOF ≥ 54, the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces perform better for all values of k. In light of these observations, we can therefore infer that, for sufficiently large DoF, the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces are better adopted to the high-frequency structure of the scattering problem for convex obstacles. Finally, in the case of the nonconvex kite shaped scatterer, no reasonable performance comparison can be made except that for large DoF the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces seem to provide better approximations. Table 1 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces for the single-scattering configurations in Fig. 3 and for various different wavenumbers k Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Table 1 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces for the single-scattering configurations in Fig. 3 and for various different wavenumbers k Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Let us finally note that the error term on the right-hand side of equation (3.6) in Corollary 3.1 suggests (excluding the fact that $$C_{k}/c_{k} = \mathscr{O}(k^{1/3})$$ as $$k \to \infty$$ for CFIE) that DoF should increase in proportion to $$\log k$$ with increasing k in order to maintain an error independent on k. Although the data in Table 1 essentially show for the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces that the error decreases with increasing k for sufficiently large DoF, in Table 2 we display DoF vs. the relative L2 error for the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces when $$\textrm{DoF} \approx 80 \times \log _{10}(k/20)$$. Table 2 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces when DoF increases logarithmically with increasing k: DoF $$\approx 80 \times \log _{10}(k/20)$$ Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 View Large Table 2 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces when DoF increases logarithmically with increasing k: DoF $$\approx 80 \times \log _{10}(k/20)$$ Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 View Large 5.3. Approximations in the shadow region Here we focus on approximations in the shadow region, and we present numerical tests exhibiting the performance of Galerkin BEMs based on frequency dependent changes of variables. To this end, in Fig. 4, we consider the unit circle parametrized as $$(\cos \theta ,\sin \theta )$$, 0 ≤ θ ≤ 2π, and illuminated by a plane wave incidence with direction α = (1, 0), and we display the pointwise relative error as a function of the radial angle θ for k = 10, 20, 40, 80, 120. In these experiments, we have constructed the Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ defined in equation (3.3) by using J = 6 regions. We have assigned the polynomial degrees d = 3, 7, 11 to each of the six regions and this resulted in the total numbers of DoF = 6(d + 1) = 24, 48, 72. Fig. 4. View largeDownload slide Pointwise relative error on the boundary of the unit circle for increasing DoF and increasing wavenumbers k. Fig. 4. View largeDownload slide Pointwise relative error on the boundary of the unit circle for increasing DoF and increasing wavenumbers k. The pointwise relative errors displayed in Fig. 4 show that the quality of approximations in the illuminated region, namely the interval (π/2, 3π/2), improves with increasing DoF and it does not essentially depend on the wavenumber k. The approximations in the shadow region also improve with increasing DoF while degrading for increasing values of the wavenumber k. Except for k = 40, there is essentially no difference between the approximations in the intervals [0, θ0] and [θ1, 2π] for DoF = 48 and 72 (here θ0 = 0.86 and θ1 = 2π − θ0 and they are depicted by the vertical dotted lines in these plots). As we discuss below, this behavior is related with the asymptotic behavior of the amplitude ηslow in the shadow region. On the other hand, let us emphasize that the change of variables used to construct the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ Galerkin approximation spaces are designed to resolve the boundary layers around the shadow boundaries (i.e. the regions where the derivatives of ηslow actually blow up with increasing k—as implied by the estimates in Theorem 4.2) by assigning more DoF in these regions when compared to the deep shadow regions with increasing k. Moreover, from a numerical perspective, it is also important to note that the condition numbers of the Galerkin linear systems increase with increasing DoF, e.g. they are about 1e + 3, 3.2e + 5 and 3.2e + 8 (almost independently of the value of k) for DoF = 24, 48 and 72, respectively. Pointwise relative error on the boundary of the unit circle for increasing DoF and increasing wavenumbers k. As for the asymptotic behavior of the amplitude ηslow in the shadow region, its phase in this region is modulated depending on the wavenumber, the curvature, and the geodesic distance to the shadow boundaries. Moreover, while its size increases linearly in the illuminated region with increasing wavenumber k, it decays exponentially in the shadow region. Furthermore, this decay depends additionally on the curvature and the geodesic distance to the shadow boundary points. A detailed discussion on these issues can be found in Asheim & Huybrechs (2014) where, based on suitable combinations of the GTD and Melrose–Taylor asymptotics, the authors design numerical schemes more appropriate for approximations in the shadow region. To illustrate this decaying behavior, let us consider a circular scatterer of radius r parametrized as $$(r\cos \theta ,r\sin \theta )$$, 0 ≤ θ ≤ 2π, and illuminated by a plane wave incidence with direction α = (1, 0). In this case, simple calculations based on Fourier analysis (Bowman et al., 1970) shows that η can be expressed as a function of the radial angle θ in the form \begin{align} \eta(\theta) = - \dfrac{2i}{\pi r} \sum_{n = -\infty}^{\infty} \dfrac{e^{in(\frac{\pi}{2}+\theta)}}{H^{(1)}_{n}(kr)}, \qquad 0 \le \theta \le 2\pi, \end{align} (5.7) where $$H^{(1)}_{n}$$ is the Hankel function of the first kind and order n. In Fig. 5, we display the absolute value of η as a function of the radial angle θ for increasing values of k for circles of radii r = 1 and r = 1/2. As this figure depicts, the region where |η(θ)|≥ 1 shrinks towards the illuminated region (i.e. the interval (π/2, 3π/2)) with increasing k, and its size is smaller for r = 1. Moreover, $$\max |\eta |$$ is approximately 2k for both r = 1 and r = 1/2, but as k increases from 10 to 1280 the values of $$\min |\eta |$$ decreases from 3.2e − 2 to 1.1e − 10 for r = 1 and from 1.7e − 1 to 6.3e − 8 for r = 1/2. Therefore, it can be inferred that accurately approximating the solution in the (deep) shadow region requires the use of significantly higher numerical precision, not only with increasing wavenumber k, but also increasing radius of curvature. Fig. 5. View largeDownload slide |η(θ)| plotted against the radial angle θ for circles of radii r = 1 (left) and r = 1/2 (right) and increasing values of k. Fig. 5. View largeDownload slide |η(θ)| plotted against the radial angle θ for circles of radii r = 1 (left) and r = 1/2 (right) and increasing values of k. With these observations, let us finally compare the numerical results in Fig. 4 above with those presented by Asheim & Huybrechs (2014). Just as in this figure, their implementations are confined to k = 10, 20, 40, 80, 120. Figure 11 therein concerns scattering by a circle of radius r = 1/2 and displays the pointwise relative errors for DoF = 20 (note that, since the number of wave lengths λ = 2π/k on the surface of the unit circle is twice that of the circle of radius r = 1/2, this corresponds to DoF = 40 for the unit circle). There they approximate the solution in $$[0,\theta ^{\prime }_{0}] \cup [\theta ^{\prime }_{1},2\pi ]$$ using the values of the approximate solution at the radial angles $$\theta ^{\prime }_{0} = \pi /5$$ and $$\theta ^{\prime }_{1} = 2\pi -\theta ^{\prime }_{0}$$ (these angles are exactly the vertical dotted lines on the right pane in Fig. 5 above). Evaluating η(θ) through equation (5.7) to machine accuracy for r = 1/2 reveals that $$0.14 \le |\eta (\theta ^{\prime }_{j})| \le 0.74$$ (j = 0, 1) for k = 10, 20, 40, 80, 120. For r = 1, the vertical dotted lines in the left pane in Fig. 5 as well as those in Fig. 4 correspond exactly to θ0 = 0.86 and θ1 = 2π − θ0, and evaluating η(θ) through formula (5.7) we see that, in this case, 0.14 ≤|η(θj)|≤ 0.67 (j = 0, 1) for k = 10, 20, 40, 80, 120. Table 3 $$\min |\eta (\theta )|$$ for circles of radii r = 1 and r = 1/2 and the ellipse for the wavenumbers k = 10, 20, 40, 80, 120 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 Table 3 $$\min |\eta (\theta )|$$ for circles of radii r = 1 and r = 1/2 and the ellipse for the wavenumbers k = 10, 20, 40, 80, 120 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 In light of the discussions in Asheim & Huybrechs (2014), this explains the lack of accuracy in the pointwise relative errors on [0, θ0] ∪ [θ1, 2π] in Fig. 4 above. Let us mention in this connection that the quality of numerical approximations in Fig. 11 therein are virtually independent on frequency. This is, however, not the case in Fig. 14 therein which concerns scattering by an ellipse parametrized as $$(\frac{1}{2}\cos \theta , \sin \theta )$$ that is also rotated by π/6 in the counterclockwise direction and illuminated by a plane wave incidence with direction α = (1, 0). The circumference of this scatterer is approximately 4.97 and they use DoF = 30 (note that, with the same reasoning as above, DoF = 20 for the circle of radius r = 1/2 would correspond to about DoF = 32 for this ellipse). In the illuminated region the error decreases as k increases from 10 to 80, but there are virtually no differences between the errors for k = 80 and 120. In the shadow region, the error degrades and blows up in the deep shadow region just as in Fig. 4 above. Let us finally note that the differences between the sizes of blow up in the pointwise relative errors in the deep shadow regions are related with the minimum values of |η(θ)| depicted in Table 3. 6. Conclusions In this paper we proposed a new class of Galerkin BEMs based on novel changes of variables for the solution of two-dimensional single-scattering problems. The Galerkin approximation spaces, generated in the form of a direct sum of algebraic polynomial spaces weighted by the oscillations in the incident field of radiation, are additionally coupled with novel frequency dependent changes of variables in the transition regions. As we have shown, this construction ensures that the global approximation spaces are perfectly matched with the changes in the boundary layers of the solutions with increasing wavenumber k. Furthermore, for large frequencies, they provide significant savings over their counterparts in Ecevit & Özen (2017) in regards to the total number of DoF necessary to obtain a prescribed accuracy. Acknowledgments We thank the anonymous referees for their constructive comments which significantly improved the presentation in the paper. References Abboud , T. , Nédélec , J.-C. & Zhou , B. ( 1994 ) Méthode des équations intégrales pour les hautes fréquences . C. R. Acad. Sci. Paris Sér. I Math. , 318 , 165 -- 170 . Abboud , T. , Nédélec , J.-C. & Zhou , B. ( 1995 ) Improvements of the integral equation method for high frequency problems . In Proc. of 3rd Int. Conf. on Mathematical Aspects of Wave Propagation Phenomena, SIAM, Philadelphia, PA, USA , 178 -- 187 . Anand , A. , Boubendir , Y. , Ecevit , F. & Reitich , F. ( 2010 ) Analysis of multiple scattering iterations for high-frequency scattering problems. II. The three-dimensional scalar case . Numer. Math. , 114 , 373 -- 427 . Google Scholar CrossRef Search ADS Asheim , A. & Huybrechs , D. ( 2014 ) Extraction of uniformly accurate phase functions across smooth shadow boundaries in high frequency scattering problems . SIAM J. Appl. Math. , 74 , 454 -- 476 . Google Scholar CrossRef Search ADS Banjai , L. & Hackbusch , W. ( 2008 ) Hierarchical matrix techniques for low- and high-frequency Helmholtz problems . IMA J. Numer. Anal. , 28 , 46 -- 79 . Google Scholar CrossRef Search ADS Boffi , D . ( 2010 ) Finite element approximation of eigenvalue problems . Acta Numer ., 19 , 1 -- 120 . Google Scholar CrossRef Search ADS Boubendir , Y. , Ecevit , F. & Reitich , F. ( 2016 ) Acceleration of an iterative method for the evaluation of high-frequency multiple scattering effects. ArXiv e-prints . Bowman , J. J. , Senior , T. B. A. & Uslenghi , P. L. E. ( 1987 ) Electromagnetic and Acoustic Scattering by Simple Shapes . New York : Hemisphere Publishing Corp ., xix + 747 . Bruno , O. P. & Geuzaine , C. A. ( 2007 ) An O(1) integration scheme for three-dimensional surface scattering problems . J. Comput. Appl. Math. , 204 , 463 -- 476 . Google Scholar CrossRef Search ADS Bruno , O. P. & Kunyansky , L. A. ( 2001 ) A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications . J. Comput. Phys. , 169 , 80 -- 110 . Google Scholar CrossRef Search ADS Bruno , O. , Geuzaine , C. & Reitich , F. ( 2005 ) On the O(1) solution of multiple-scattering problems . IEEE Trans. Magn. , 41 , 1488 -- 1491 . Google Scholar CrossRef Search ADS Bruno , O. P. , Geuzaine , C. A. , Monro , J. A. Jr. & Reitich , F. ( 2004 ) Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case . Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. , 362 , 629 -- 645 . Google Scholar CrossRef Search ADS Bruno , O. P. , Domínguez , V. & Sayas , F.-J. ( 2013 ) Convergence analysis of a high-order Nyström integral-equation method for surface scattering problems . Numer. Math. , 124 , 603 -- 645 . Google Scholar CrossRef Search ADS Chandler-Wilde , S. N. , Graham , I. G. , Langdon , S. & Lindner , M. ( 2009 ) Condition number estimates for combined potential boundary integral operators in acoustic scattering . J. Integral Equations Appl. , 21 , 229 -- 279 . Google Scholar CrossRef Search ADS Chandler-Wilde , S. N. , Graham , I. G. , Langdon , S. & Spence , E. A. ( 2012 ) Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering . Acta Numer ., 21 , 89 -- 305 . Google Scholar CrossRef Search ADS Chandler-Wilde , S. N. , Hewett , D. P. , Langdon , S. & Twigger , A. ( 2015 ) A high frequency boundary element method for scattering by a class of nonconvex obstacles . Numer. Math. , 129 , 647 -- 689 . Google Scholar CrossRef Search ADS Colton , D. & Kress , R. ( 1992 ) Inverse Acoustic and Electromagnetic Scattering Theory . Applied Mathematical Sciences, vol. 93. Berlin: Springer, pp. x+305 . Corless , R. M. , Gonnet , G. H. , Hare , D. E. G. , Jeffrey , D. J. & Knuth , D. E. ( 1996 ) On the LambertW function . Adv. Comput. Math. , 5 , 329 -- 359 . Google Scholar CrossRef Search ADS Domínguez , V. , Graham , I. G. & Smyshlyaev , V. P. ( 2007 ) A hybrid numerical-asymptotic boundary integral method for high-frequency acoustic scattering . Numer. Math. , 106 , 471 -- 510 . Google Scholar CrossRef Search ADS Darrigrand , E . ( 2002 ) Coupling of fast multipole method and microlocal discretization for the 3-D Helmholtz equation . J. Comput. Phys. , 181 , 126 -- 154 . Google Scholar CrossRef Search ADS Davies , R. W. , Morgan , K. & Hassan , O. ( 2009 ) A high order hybrid finite element method applied to the solution of electromagnetic wave scattering problems in the time domain . Comput. Mech. , 44 , 321 -- 331 . Google Scholar CrossRef Search ADS Ecevit , F. & Özen , H. Ç. ( 2017 ) Frequency-adapted galerkin boundary element methods for convex scattering problems . Numer. Math. , 135 , 27 -- 71 . Google Scholar CrossRef Search ADS Ecevit , F. & Reitich , F. ( 2009 ) Analysis of multiple scattering iterations for high-frequency scattering problems. I. The two-dimensional case . Numer. Math. , 114 , 271 -- 354 . Google Scholar CrossRef Search ADS Engquist , B. & Majda , A. ( 1977 ) Absorbing boundary conditions for the numerical simulation of waves . Math. Comp. , 31 , 629 -- 651 . Google Scholar CrossRef Search ADS Eruslu , H. H. ( 2015 ) An optimal change of variables scheme for single scattering problems . Master’s Thesis . Istanbul : Boğaziçi University . Galkowski , J. ( 2014 ) Distribution of resonances in scattering by thin barriers . ArXiv e-prints . Galkowski , J. , Müller , E. H. & Spence , E. A. ( 2016 ) Wavenumber-explicit analysis for the Helmholtz h-BEM: error estimates and iteration counts for the Dirichlet problem . ArXiv e-prints . Ganesh , M. & Hawkins , S. C. ( 2011 ) A fully discrete Galerkin method for high frequency exterior acoustic scattering in three dimensions . J. Comput. Phys. , 230 , 104 -- 125 . Google Scholar CrossRef Search ADS Giladi , E . ( 2007 ) Asymptotically derived boundary elements for the Helmholtz equation in high frequencies . J. Comput. Appl. Math. , 198 , 52 -- 74 . Google Scholar CrossRef Search ADS Givoli , D . ( 2004 ) High-order local non-reflecting boundary conditions: a review . Wave Motion , 39 , 319 -- 326 . Google Scholar CrossRef Search ADS Gradshteyn , I. S. & Ryzhik , I. M. ( 2000 ) Table of Integrals, Series, and Products, 6th edn. San Diego, CA: Academic Press Inc., pp. xlvii+1163 . Grote , M. J. & Sim , I. ( 2011 ) Local nonreflecting boundary condition for time-dependent multiple scattering . J. Comput. Phys. , 230 , 3135 -- 3154 . Google Scholar CrossRef Search ADS Han , X. & Tacy , M. ( 2015 ) Sharp norm estimates of layer potentials and operators at high frequency . J. Funct. Anal. , 269 , 2890 -- 2926 . Google Scholar CrossRef Search ADS Hesthaven , J. & Warburton , T. ( 2004 ) High-order accurate methods for time-domain electromagnetics. CMES Comput. Model. Eng. Sci. , 5 , 395 -- 407 . Hewett , D. P. , Langdon , S. & Melenk , J. M. ( 2013 ) A high frequency $hp$ boundary element method for scattering by convex polygons . SIAM J. Numer. Anal. , 51 , 629 -- 653 . Google Scholar CrossRef Search ADS Hewett , D. P. , Langdon , S. & Chandler-Wilde , S. N. ( 2015 ) A frequency-independent boundary element method for scattering by two-dimensional screens and apertures . IMA J. Numer. Anal. , 35 , 1698 -- 1728 . Google Scholar CrossRef Search ADS Huybrechs , D. & Vandewalle , S. ( 2007 ) A sparse discretization for integral equation formulations of high frequency scattering problems . SIAM J. Sci. Comput. , 29 , 2305 -- 2328 . Google Scholar CrossRef Search ADS Kim , T. ( 2012 ) Asymptotic and numerical methods for high-frequency scattering problems. Ph.D. Thesis , University of Bath, UK . Lazergui , S. & Boubendir , Y. ( 2016 ) Asymptotic expansions of the Helmholtz equation solutions using approximations of the Dirichlet to Neumann operator . ArXiv e-prints . Melrose , R. B. & Taylor , M. E. ( 1985 ) Near peak scattering and the corrected Kirchhoff approximation for a convex obstacle . Adv. Math. , 55 , 242 -- 315 . Google Scholar CrossRef Search ADS Schwab , C. ( 1998 ) p- and hp-finite Element Methods: Theory and Applications in Solid and Fluid Mechanics. Numerical Mathematics and Scientific Computation. The Clarendon Press, New York: Oxford University Press, pp. xii+374 . Spence , E. A. , Chandler-Wilde , S. N. , Graham , I. G. & Smyshlyaev , V. P. ( 2011 ) A new frequency-uniform coercive boundary integral equation for acoustic scattering . Comm. Pure Appl. Math. , 64 , 1384 -- 1415 . Google Scholar CrossRef Search ADS Spence , E. A. , Kamotski , I. V. & Smyshlyaev , V. P. ( 2015 ) Coercivity of combined boundary integral equations in high-frequency scattering . Comm. Pure Appl. Math. , 68 , 1587 -- 1639 . Google Scholar CrossRef Search ADS Tong , M. S. & Chew , W. C. ( 2010 ) Multilevel fast multipole acceleration in the Nyström discretization of surface electromagnetic integral equations for composite objects . IEEE Trans. Antennas Propagation , 58 , 3411 -- 3416 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# A Galerkin BEM for high-frequency scattering problems based on frequency-dependent changes of variables

, Volume Advance Article – Feb 5, 2018
31 pages

/lp/ou_press/a-galerkin-bem-for-high-frequency-scattering-problems-based-on-qx1qob3AUQ
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx079
Publisher site
See Article on Publisher Site

### Abstract

Abstract In this paper we develop a new class of semidiscrete Galerkin boundary element methods for the solution of two-dimensional exterior single-scattering problems. Our approach is based upon construction of Galerkin approximation spaces confined to the asymptotic behavior of the solution through a certain direct sum of appropriate function spaces weighted by the oscillations in the incident field of radiation. Specifically, the function spaces in the illuminated/shadow regions and the shadow boundaries are simply algebraic polynomials, whereas those in the transition regions are generated utilizing novel, yet simple, frequency-dependent changes of variables perfectly matched with the boundary layers of the amplitude in these regions. We rigorously verify for compact, smooth and strictly convex obstacles that, with increasing wavenumber k, these methods require, for any ε > 0, only an $$\mathscr{O}\left ( k^{\varepsilon } \right )$$ increase in the number of degrees of freedom (DoF) to maintain any given accuracy independent on frequency. These theoretical results are confirmed by numerical tests that show that, for sufficiently large DoF, the error tends to decrease with increasing wavenumber k. 1. Introduction High-frequency scattering problems continue to be of immense interest in present day computational science. Indeed, over the course of last two decades, very efficient and effective algorithms have been devised for the numerical solution of scattering problems based on variational (Hesthaven & Warburton, 2004; Davies et al., 2009; Boffi, 2010) and integral equation (Bruno & Kunyansky, 2001; Banjai & Hackbusch, 2008; Tong & Chew, 2010; Bruno et al., 2013) formulations. In exterior scattering simulations methods that rest on variational formulations naturally demand the design and implementation of efficient nonreflecting boundary conditions (Engquist & Majda, 1977; Givoli, 2004; Grote & Sim, 2011) to effectively represent the radiation condition at infinity. On the other hand, solvers based on integral equation formulations (cf. the survey article Chandler-Wilde et al. (2012) on the numerical-asymptotic boundary integral equation methods) readily encode the radiation condition into the equation by choosing an outgoing fundamental solution. Moreover, for surface scattering simulations as considered in this paper, they enable phase extraction, that takes a particularly simple form in single-scattering configurations, and this turns the problem into the estimation of an amplitude whose oscillations are essentially independent of frequency. In this paper we develop a new class of semidiscrete Galerkin boundary element methods (BEMs) for the solution of two-dimensional scattering problems in the exterior of compact, smooth and convex obstacles having everywhere positive curvature. The approach taken here is based on the completely rigorous asymptotic expansions of the aforementioned amplitude developed by Melrose & Taylor (1985) (for their alternative forms, see Lazergui & Boubendir (2016)). These expansions imply the existence of boundary layers around the shadow boundaries (where the rays are tangential to the boundary) and the Galerkin approximation spaces are constructed so as to resolve these layers through frequency-dependent changes of variables. From a theoretical point of view, the resulting methodologies require, for any ε > 0, only an $$\mathscr{O}\left ( k^{\varepsilon } \right )$$ increase in the number of degrees of freedom (DoF) to maintain any given accuracy independent on frequency. These theoretical findings are supported by numerical implementations that show that the method achieves similar accuracy even when frequency-independent numbers of DoF are used. From an application point of view, the tests also show that the algorithms proposed herein are applicable to the class of nonconvex smooth scatterers and incidence directions that do not allow for multiple reflections and more than two shadow boundaries. Indeed, hybrid integral equation methodologies incorporating the asymptotic characteristics of the unknown into the solution strategy have now become the usual practice in the field. The first attempt in this direction is due to Abboud et al. (1994, 1995) where, considering the impedance boundary condition, an h-version BEM was utilized in conjunction with the method of stationary phase for the evaluation of highly oscillatory integrals (see Darrigrand (2002) for a fast multipole variant, and Ganesh & Hawkins (2011) for a fully discrete three-dimensional version). More relevant to our work is the Nyström method proposed by Bruno et al. (2004) for the solution of sound-soft scattering problems in the exterior of smooth convex obstacles. The method therein displays the capability of delivering solutions within any prescribed accuracy in frequency-independent computational times. While, in this approach, the boundary layers of the slowly varying amplitude around the shadow boundaries are resolved through a cubic root change of variables, the associated highly oscillatory integrals are evaluated to high order utilizing novel extensions of the method of stationary phase (for a three-dimensional variant of this approach we refer to Bruno & Geuzaine (2007)). The algorithm in Bruno et al. (2004) has had a great impact in the computational scattering community and, following the basic principles therein, a number of alternative single-scattering solvers have been developed. Giladi (2007) has used a collocation method that integrates Keller’s geometrical theory of diffraction (GTD) to account for creeping rays in the shadow region. Huybrechs and Vandewalle’s collocation method (Huybrechs & Vandewalle, 2007) has utilized the numerical steepest descent method in evaluating highly oscillatory integrals and additional collocation points around shadow boundaries to obtain sparse discretizations. The first rigorous numerical analysis relating to a p-version boundary element implementation of these approaches, due to Domínguez et al. (2007), proved that an increase of $$\mathscr{O}(k^{1/9})$$ in the number of DoF is sufficient to preserve a certain accuracy as $$k \to \infty$$. Parallel with the schemes relating to smooth convex obstacles, Galerkin BEMs based on phase extraction have also been developed for half-planes and convex polygons where the number of DoF is either fixed or must increase in proportion to $$\log k$$ to fix the error with increasing wavenumber k. We refer to the survey article Chandler-Wilde et al. (2012) for an extended review of these procedures; for more recent hp BEMs for convex polygons see Hewett et al. (2013), and for screens and apertures see Hewett et al. (2015). As for multiple scattering problems we refer to Bruno et al. (2005) for an extension of the algorithm in Bruno et al. (2004) to a finite collection of convex obstacles (see also Ecevit & Reitich (2009) and Anand et al. (2010) for a rigorous analysis of this approach in two- and three-dimensional settings, respectively, and Boubendir et al. (2016) for the acceleration of this procedure through use of dynamical Krylov subspaces and Kirchhoff approximations), and to Chandler-Wilde et al. (2015) for a class of nonconvex polygons. Concerning compact, smooth and strictly convex obstacles considered in this paper, let us mention that the methods developed in Bruno et al. (2004), Huybrechs & Vandewalle (2007) and Domínguez et al. (2007) remain asymptotic due to approximation of the solution by zero in the deep shadow region. In our recent frequency-adapted Galerkin BEMs (Ecevit & Özen, 2017), on the other hand, we utilized approximation spaces also in this region for two reasons. First, based on an optimal use of the wavenumber-dependent derivative estimates of the amplitude given by Domínguez et al. (2007, Theorem 5.4), this allowed us to rigorously eliminate the asymptotic terms appearing in the error analysis due to approximation of the solution by zero in the deep shadow region. Secondly, as displayed through numerical tests, use of these additional approximation spaces significantly reduces the condition numbers of the Galerkin linear systems and renders them virtually independent on frequency. This, in return, improves the numerical stability of the method and allows the use of higher DoF within a limited numerical precision which results in improved approximations. Let us mention that our approach in Ecevit & Özen (2017) was based on utilization of appropriate integral equation formulations of the scattering problem and design of Galerkin approximation spaces, as the direct sum of algebraic or trigonometric polynomials weighted by the oscillations in the incident field of radiation. For any given $$m \in \mathbb{N}$$, the number of direct summands was chosen as 4m by assigning one approximation space to each of the illuminated and deep shadow regions, and one to each of the two shadow boundaries, and m − 1 to each one of the four transition regions. As we have shown in Ecevit & Özen (2017), from a theoretical perspective, m had to increase as $$\mathscr{O}(\log k)$$ in order to obtain optimal error bounds, and that these methods can be tuned to demand an increase of $$\mathscr{O}(k^{\varepsilon })$$ (for any ε > 0) in the number of DoF to maintain a prescribed accuracy independent on frequency. For smooth strictly convex scatterers, this is the best theoretical result available in the literature. We note, however, that the numerical results in Ecevit & Özen (2017) were restricted to m = 1 and m = 2, since utilization of larger values of m is simply impractical from an implementation point of view. In the new Galerkin BEMs developed herein we continue to use approximation spaces in the deep shadow region, since this allows the elimination of the asymptotic terms in the error analysis, and their utilization significantly reduces the condition numbers of the Galerkin linear systems. However, we treat the transition regions differently as we utilize novel frequency-dependent changes of variables perfectly matched with the asymptotic behavior of the solution in the transition regions. We thereby eliminate the requirement of increasing the number of direct summands defining the Galerkin approximation spaces with increasing wavenumber k. While this makes the new scheme we propose far easier to implement than that in Ecevit & Özen (2017), from a theoretical point of view, an increase of $$\mathscr{O}(k^{\varepsilon })$$ (for any ε > 0) is still sufficient to fix the approximation error with increasing k. The numerical tests show that, for sufficiently large DoF, the new scheme outperforms the method in Ecevit & Özen (2017) and the relative L2 error tends to decrease with increasing wavenumber k. Let us note that, in fact, accurately approximating the solution in the (deep) shadow region involves significant numerical difficulties since the solution there decays exponentially with increasing wavenumber and geodesic distance to the shadow boundaries. Moreover, its oscillations are further amplified by a phase term that depends on the curvature and geodesic distance to the shadow boundaries (see Remark 3 for a more precise explanation; also cf. Asheim & Huybrechs (2014) and the references therein). For this reason, we test the performance of the new Galerkin BEMs in the shadow region in terms of pointwise relative errors and, in this connection, discuss the decaying behavior of the solution in the shadow region depending on the wavenumber and the radius of curvature. Further, we compare the numerical results with those of Asheim & Huybrechs (2014), where they have developed adequate phase extraction techniques (based on suitable combinations of GTD and Melrose--Taylor asymptotics) more appropriate for approximating the solution in the shadow region. Incidentally, it is important to note that, just as the methods in Domínguez et al. (2007) and Ecevit & Özen (2017), the methods developed here are semidiscrete as they transform the solution of the scattering problem to the evaluation of highly-oscillatory double integrals. Efficient evaluation of integrals arising in connection with the method in Domínguez et al. (2007) is discussed in both theoretical and practical terms by Kim (2012). With some additional effort, the approach therein can be extended to cover the integrals associated with the method in Ecevit & Özen (2017). For the new methods we propose here, the structure of these highly-oscillatory double integrals are slightly different, but they can be treated similarly. The actual theoretical and implementation details are left for future work. The paper is organized as follows. In §2, we describe the exterior sound-soft scattering problem along with the relevant integral equations and associated Galerkin formulations. In §3, we introduce the new Galerkin schemes for high-frequency single-scattering problems, and state the associated convergence theorem for compact, smooth, and strictly convex obstacles which constitutes the main result of the paper. To allow a direct comparison, in the same section, we also present a more general version of our algorithm in Ecevit & Özen (2017) along with the corresponding approximation properties. The proof of the main result is given in §4. In §5.1 we present the highly-oscillatory double integrals arising in connection with the methods developed herein and discuss implementation details. §5.2 is reserved for numerical tests that provide a comparison of our methods developed herein and Ecevit & Özen (2017). Finally, in §5.3, we examine approximations in the shadow region particularly in contrast with those in Asheim & Huybrechs (2014). 2. The scattering problem and Galerkin formulation The two-dimensional scattering problem we consider in this manuscript is related with the determination of the scattered field u that satisfies the Helmholtz equation \begin{align*} \left( \varDelta + k^{2} \right) u = 0 \end{align*} in the exterior of a compact smooth strictly convex obstacle K, the Sommerfeld radiation condition at infinity that amounts to requiring \begin{align*} \lim_{\left| x \right| \to \infty} \left| x \right|{}^{1/2} \left[ \partial_{\left| x \right|} - i k \right] u = 0 \end{align*} uniformly for all directions $$x/\left | x \right |$$, and the sound-soft boundary condition \begin{align*} u = - u^{\textrm{inc}} \end{align*} for a plane-wave incidence $$u^{\textrm{inc}} \left ( x \right ) = e^{ik \alpha \cdot x}$$ with direction α ($$\left | \alpha \right | =1$$) impinging on K. As is well known, the scattered field u can be reconstructed by means of either the direct or the indirect approach (Colton & Kress, 1992); for discussions on the direct and indirect approaches we refer to Bruno et al. (2004, pp. 631–632) and Chandler-Wilde et al. (2012, p. 132). As in the previous attempts aimed at frequency-independent simulations (Bruno et al., 2004; Domínguez et al., 2007; Giladi, 2007; Huybrechs & Vandewalle, 2007; Ecevit & Özen, 2017;), here we favor the former wherein the associated (unknown) surface density is the normal derivative of the total field (known as the surface current in electromagnetism) $$\eta = \partial _{\nu } \left ( u+ u^{\textrm{inc}} \right )$$ on ∂K. Once η is available, the scattered field can be recovered through the single-layer potential \begin{align*} u(x) = - \int_{\partial K} \varPhi(x,y) \, \eta(y) \, \textrm{d}s(y), \end{align*} where \begin{align*} \varPhi(x,y) = \dfrac{i}{4} \ H_{0}^{(1)}(k|x-y|) \end{align*} is the fundamental solution of the Helmholtz equation, and $$H_{0}^{(1)}$$ is the Hankel function of the first kind and order zero. The new unknown, namely η, can be expressed as the unique solution of a variety of integral equations (Colton & Kress, 1992; Spence et al., 2011) taking the form of an operator equation \begin{align} \mathscr{R}_{k} \, \eta = f_{k} \end{align} (2.1) in $$L^{2} \left ( \partial K \right )$$. The solution of (2.1) corresponds exactly to that of \begin{align} B_{k}(\mu,\eta) = F_{k}(\mu), \quad \textrm{for all } \mu \in L^{2}(\partial K), \end{align} (2.2) where the sesquilinear form Bk and bounded linear functional Fk are defined by \begin{align} B_{k}(\mu,\eta) = \langle \mu, \mathscr{R}_{k} \, \eta \rangle_{L^{2}(\partial K)} \quad \textrm{and} \quad F_{k}(\mu) = \langle \mu, f_{k} \rangle_{L^{2}(\partial K)}. \end{align} (2.3) Equation (2.2), in turn, is amenable to a treatment by the Galerkin method wherein one determines the Galerkin solution$$\hat{\eta }$$ approximating the exact solution η in a given finite dimensional Galerkin subspace$$\hat{X}_{k}$$ requiring that \begin{align} B_{k}(\hat{\mu},\hat{\eta}) = F_{k}(\hat{\mu}), \quad \textrm{for all } \hat{\mu} \in \hat{X}_{k}. \end{align} (2.4) Further, provided the sesquilinear form Bk is continuous with a continuity constant Ck and strictly coercive with a coercivity constant ck so that \begin{align*} \left| B_{k}(\mu,\eta) \right| \le C_{k} \, \Vert \mu \Vert \Vert \eta \Vert \qquad \textrm{and} \qquad \mathscr{R}e \ B_{k}(\mu,\mu) \ge c_{k} \, \Vert \mu \Vert^{2} \end{align*} for all μ, η ∈ L2(∂K), equation (2.4) is uniquely solvable and C$${\acute{\textrm{e}}}$$a's lemma implies \begin{align} \Vert \eta -\hat{\eta} \Vert \leq \dfrac{C_{k}}{c_{k}} \inf_{\hat{\mu} \in \hat{X}_{k}} \Vert \eta- \hat{\mu} \Vert. \end{align} (2.5) Remark 2.1 Among the aforementioned integral equations, in this connection, combined field integral equation (CFIE) and star-combined integral equation (SCIE) step forward as the continuity and coercivity properties of the associated sesquilinear forms are well-understood. More precisely, the sesquilinear form associated with CFIE is continuous for Lipschitz domains for all k > 0 with $$C_{k} = \mathscr{O} \left ( k^{1/2} \right )$$ as $$k \to \infty$$ (Chandler-Wilde et al., 2009, Theorem 3.6), and coercive (for k ≫ 1) for strictly convex domains having piecewise analytic C3 boundaries with everywhere positive curvature such that, for any 0 < δ < 1/2, ck ≥ (1/2) − δ for all sufficiently large k (Spence et al., 2015, Theorem 1.2). On the other hand, the sesquilinear form corresponding to SCIE is both continuous and coercive (for k > 0) for star-shaped Lipschitz domains with $$C_{k} = \mathscr{O} \left ( k^{1/2} \right )$$ as $$k \to \infty$$ and ck independent of k (Spence et al., 2011). Moreover, recently obtained bounds, based on the results in Han & Tacy (2015), on the single- and double-layer potentials (see Galkowski (2014, Theorems 4.29 and 4.32) and Galkowski et al. (2016, Theorem 1.4) for CFIE, and Galkowski (2014, Theorems 4.29 and 4.32) and Galkowski et al. (2016, Theorem 1.9, Remark 4.2) for SCIE) show that, for strictly convex domains having everywhere positive curvature, $$C_{k} = \mathscr{O} \left ( k^{1/3} \right )$$ as $$k \to \infty$$. 3. Galerkin approximation spaces based on frequency-dependent changes of variables The developments in this section are independent of the integral equation used as they relate, specifically, to the construction of Galerkin approximation spaces $$\hat{X}_{k}$$ whose dimension should increase only as $$\log k$$ with increasing wavenumber k to ensure that the relative error \begin{align*} \inf_{\hat{\mu} \in \hat{X}_{k}} \dfrac{\Vert \eta- \hat{\mu} \Vert}{\Vert \eta \Vert} \end{align*} in connection with the infimum on the right-hand side of (2.5) is independent of k. Considering a smooth convex obstacle K illuminated by a plane-wave uinc(x) = eikα⋅x, our approach is based on phase extraction \begin{align*} \eta \left( x \right) = e^{ik \alpha \cdot x} \, \eta^{\textrm{slow}} \left( x \right), \qquad x \in \partial K, \end{align*} and design of approximation spaces adopted to the asymptotic behavior of the amplitude ηslow that was initially characterized by Melrose & Taylor (1985) around the shadow boundaries which we have later generalized to the entire boundary (Ecevit & Reitich, 2009). Remark 3.1 The next theorem is included to highlight the difficulties associated with the construction of Galerkin approximation spaces, and it is not used in any of the proofs. In particular, the definition of Hörmander classes and a detailed understanding of asymptotic expansions are not necessary to follow the developments in the rest of the paper (the interested reader may consult to Section 2.2 in Ecevit & Reitich (2009) for a short exposition on Hörmander classes and asymptotic expansions). What is needed for the analysis in this paper is the wavenumber explicit derivative estimates of ηslow as given in Theorem 4.2 below which is based on Domínguez et al. (2007, Theorem 5.4). Theorem 3.1 (Ecevit & Reitich, 2009, Corollary 2.1) Let $$K \subset \mathbb{R}^{2}$$ be a compact, strictly convex set with smooth boundary ∂K. Then ηslow = ηslow(x, k) belongs to the Hörmander class $$S^{1}_{2/3,1/3} \left ( \partial K \times \left ( 0, \infty \right ) \right )$$ and admits an asymptotic expansion \begin{align} \eta^{\textrm{slow}}(x,k) \sim \sum_{p,q \ge 0} a_{p,q}\left( x,k \right) \end{align} (3.1) with \begin{align*} a_{p,q}(x,k) = k^{2/3-2p/3-q} \, b_{p,q}(x) \, \varPsi^{(p)}(k^{1/3}Z(x)), \end{align*} where bp, q and Ψ are complex-valued $$C^{\infty }$$ functions and Z is a real-valued $$C^{\infty }$$ function that is positive on the illuminated region ∂KIL = {x ∈ ∂K : α ⋅ ν(x) < 0}, negative on the shadow region ∂KSR = {x ∈ ∂K : α ⋅ ν(x) > 0} and vanishes precisely to the first order on the shadow boundaries ∂KSB = {x ∈ ∂K : α ⋅ ν(x) = 0}. Moreover, the function Ψ admits the asymptotic expansion \begin{align*} \varPsi(\tau) \sim \sum_{j=0}^{\infty} c_{j} \tau^{1-3j} \qquad \textrm{as } \tau \to \infty, \end{align*} and Ψ is rapidly decreasing in the sense of Schwartz as $$\tau \to -\infty$$. Theorem 3.1 clearly displays the challenges associated with the construction of approximation spaces adapted to the asymptotic behavior of ηslow as it shows that while ηslow admits a classical asymptotic expansion in the illuminated region and rapidly decays in the shadow region, it possesses boundary layers in the $$\mathscr{O}(k^{-1/3})$$ neighborhoods of the shadow boundaries (see Fig. 1). We overcome this difficulty by constructing approximation spaces improving upon our approach in Ecevit & Özen (2017), and better adapted to the changes in frequency through use of wavenumber dependent changes of variables that resolve the aforementioned boundary layers, and that provides a smooth transition from the shadow boundaries into the illuminated and shadow regions. The resulting schemes, when compared with our algorithms in Ecevit & Özen (2017), are easier to implement since the Galerkin spaces are represented as the direct sum of a fixed number of approximation spaces rather than a number that should increase in proportion to $$\log k$$ as in Ecevit & Özen (2017). Moreover, as we explain, they display better approximation properties from a theoretical perspective since they provide an improvement of order $$\sqrt{\log k}$$ in the error for the same order of DoF as $$k\to \infty$$. Fig. 1. View largeDownload slide The unit circle illuminated from the left and the associated boundary layers of the solutions. Fig. 1. View largeDownload slide The unit circle illuminated from the left and the associated boundary layers of the solutions. Remark 3.2 ηslow is in fact both exponentially small and highly oscillatory in the deep shadow region which is not directly apparent from the rapid decay of Ψ(τ) as $$\tau \to -\infty$$ mentioned in Theorem 3.1. Indeed, there exist constants c0≠0 and β > 0 such that \begin{align*} \varPsi(\tau) = c_{0} \, e^{-i\tau^{3}/3-i\tau\alpha_{1}} (1+ \mathscr{O}(e^{\beta\tau})), \qquad \textrm{as } \tau \to -\infty, \end{align*} where $$\alpha _{1} = e^{-2i \pi \nu _{1}/3}$$ and ν1 < 0 is the right-most root of the airy function Ai (Domínguez et al., 2007, Theorem 5.1). Therefore, in the shadow region, the leading term in (3.1) is asymptotically given by \begin{align} c_{0} \, k^{2/3} \, b_{0,0}(x) \, e^{ik(Z(x))^{3}/3-ik^{1/3}Z(x)\alpha_{1}} \end{align} (3.2) which justifies the exponential decay and highly oscillatory behavior of ηslow in the deep shadow region. On the other hand, to a leading order approximation, Z(x) can be expressed as a constant multiple of the integral of κ2/3 (κ being the curvature) between x and the nearest shadow boundary point (Asheim & Huybrechs, 2014). Accordingly, equation (3.2) also justifies the statement in the introduction, namely that in the shadow region the solution decays exponentially with increasing wavenumber and geodesic distance to the shadow boundaries, and its oscillations are further amplified by a phase term that depends on the curvature and geodesic distance to the shadow boundaries. These observations motivate one to either approximate η in the deep shadow region simply by zero (which would clearly add asymptotic terms in the error analysis as in Domínguez et al. (2007)), or construct the Galerkin approximation spaces so as to reflect the exponentially decaying and highly oscillatory behavior of ηslow in the deep shadow region. While the Galerkin approximation spaces we propose below do not approximate η by zero in the deep shadow region, they are not designed to capture the oscillations in (3.2) at large k either. Indeed they do not need to since, as we shall see, the exponential damping implies that the contribution from the deep shadow region is provably negligible. To describe our approximation spaces, we choose γ to be the L−periodic arc length parameterization of ∂K in the counterclockwise orientation with $$\alpha \cdot \nu \left ( \gamma (0) \right ) = 1$$. This ensures that if 0 < t1 < t2 < L are the points corresponding to the shadow boundaries \begin{align*} \gamma \left( \left\{ t_{1},t_{2} \right\} \right) = \partial K^{SB}, \end{align*} then the illuminated and shadow regions are given by (see Fig. 1(a)) \begin{align*} \gamma \left( \left( t_{1},t_{2} \right) \right) = \partial K^{IL} \qquad \textrm{and} \qquad \gamma \left( \left( t_{2},t_{1}+L \right) \right) = \partial K^{SR}. \end{align*} For k > 1, we introduce two types of approximation spaces confined to the regions depicted in Fig. 2(a)–(b). In both cases, we define the illuminated transition and shadow transition intervals as \begin{align*} I_{IT_{1}} & = [t_{1} + \xi_{1} k^{-1/3}, t_{1} + \xi_{1}^{\prime}] = [a_{1},b_{1}], \\ I_{IT_{2}} & = [t_{2} - \xi_{2}^{\prime}, t_{2} - \xi_{2} k^{-1/3} ] = [a_{2},b_{2}], \\ I_{ST_{1}} & = [t_{1} - \zeta_{1}^{\prime}, t_{1} - \zeta_{1} k^{-1/3} ] = [a_{3},b_{3}], \\ I_{ST_{2}} & = [t_{2} + \zeta_{2} k^{-1/3}, t_{2} + \zeta_{2}^{\prime}] = [a_{4},b_{4}], \end{align*} and the shadow boundary intervals as \begin{align*} I_{sb:{1}} & = [t_{1} - \zeta_{1} k^{-1/3}, t_{1} + \xi_{1} k^{-1/3}] = [a_{5},b_{5}], \\ I_{sb:{2}} & = [t_{2} - \xi_{2} k^{-1/3}, t_{2} + \zeta_{2} k^{-1/3} ] = [a_{6},b_{6}], \end{align*} where the parameters $$\xi _{j}, \xi _{j}^{\prime }, \zeta _{j},\zeta _{j}^{\prime }>0$$ (j = 1, 2) are chosen so that \begin{align*} t_{1} + \xi_{1} \le t_{1} + \xi_{1}^{\prime} \overset{(A)}{\le} t_{2} - \xi_{2}^{\prime} \le t_{2} - \xi_{2}, \end{align*} and \begin{align*} t_{2} + \zeta_{2} \le t_{2} + \zeta_{2}^{\prime} \overset{(B)}{\le} L + t_{1} - \zeta_{1}^{\prime} \le L + t_{1} - \zeta_{1}. \end{align*} Fig. 2. View largeDownload slide Regions on the boundary of the unit circle. Fig. 2. View largeDownload slide Regions on the boundary of the unit circle. Note specifically that the regions in Fig. 2(a) correspond to equalities in (A) and (B), and those in Fig. 2(b) to strict inequalities. In the latter case, we define the illuminated and deep shadow intervals as \begin{align*} I_{IL} & = [t_{1} + \xi_{1}^{\prime}, t_{2} - \xi_{2}^{\prime}] = [a_{7},b_{7}], \\ I_{DS} & = [ t_{2} + \zeta_{2}^{\prime}, L + t_{1} - \zeta_{1}^{\prime}] = [a_{8},b_{8}]. \end{align*} With these choices, given $$\mathbf{d} = \left ( d_{1}, \ldots , d_{J} \right ) \in \mathbb{Z}_{+}^{J}$$ where $$\mathbb{Z}_{+} = \mathbb{N} \cup \{ 0 \}$$ (with J = 6 and J = 8 for Fig. 2(a) and (b), respectively), we define our (|d| + J)−dimensional Galerkin approximation spaces based on algebraic polynomials and frequency dependent changes of variables as \begin{align} \ \end{align} (3.3) Here $$_{\mathscr{I}_{j}}$$ is the characteristic function of $$\mathscr{I}_{j} = [a_{j},b_{j}]$$, and \begin{align*} \mathscr{A}_{d_{j}}^{\mathscr{C}} = \begin{cases}{cl} \mathbb{P}_{d_{j}} \circ \phi^{-1}, & \mbox{if } \mathscr{I}_{j} \mbox{ is a transition region}, \\[0.5 em] \mathbb{P}_{d_{j}}, & \textrm{otherwise}, \end{cases} \end{align*} where $$\mathbb{P}_{d_{j}}$$ is the vector space of algebraic polynomials of degree at most dj, and ϕ is the change of variables on the transition intervals given explicitly as \begin{align*} \phi \left( s \right) = \begin{cases} t_{1} + \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{IT_{1}},\\[0.3 em] t_{2} - \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{IT_{2}},\\[0.3 em] t_{1} - \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{ST_{1}},\\[0.3 em] t_{2} + \varphi \left( s \right) k^{\psi \left( s \right)}, & s \in I_{ST_{2}}, \end{cases} \end{align*} wherein φ is the affine map \begin{align*} \varphi (s) = \begin{cases} \xi_{1} + \left( \xi_{1}^{\prime}-\xi_{1} \right) \dfrac{s-a_{1}}{b_{1}-a_{1}}, & s \in I_{IT_{1}},\\[1 em] \xi_{2}^{\prime} + \left( \xi_{2}-\xi_{2}^{\prime} \right) \dfrac{s-a_{2}}{b_{2}-a_{2}}, & s \in I_{IT_{2}},\\[1 em] \zeta_{1}^{\prime} + \left( \zeta_{1}-\zeta_{1}^{\prime} \right) \dfrac{s-a_{3}}{b_{3}-a_{3}}, & s \in I_{ST_{1}},\\[1 em] \zeta_{2} + \left( \zeta_{2}^{\prime}-\zeta_{2} \right) \dfrac{s-a_{4}}{b_{4}-a_{4}}, & s \in I_{ST_{2}}, \end{cases} \end{align*} and ψ is the linear map \begin{align*} \psi (s) = -\dfrac{1}{3} \begin{cases} \dfrac{b_{1}-s}{b_{1}-a_{1}}, & s \in I_{IT_{1}},\\[1 em] \dfrac{s-a_{2}}{b_{2}-a_{2}}, & s \in I_{IT_{2}},\\[1 em] \dfrac{s-a_{3}}{b_{3}-a_{3}}, & s \in I_{ST_{1}},\\[1 em] \dfrac{b_{4}-s}{b_{4}-a_{4}}, & s \in I_{ST_{2}}. \end{cases} \end{align*} The change of variables ϕ is constructed so that, for j = 1, 2, 3, 4, the map ϕ : [aj, bj] → [aj, bj] is an orientation preserving diffeomorphism and the exponent ψ of k increases linearly from − 1/3 to 0 as we move away from the shadow boundaries into the illuminated or shadow regions. While on the one hand this choice guarantees that the DoF assigned to the $$\mathscr{O}(k^{-1/3})$$ neighborhoods of the shadow boundaries remains fixed with increasing wavenumber k, and on the other hand it also ensures that the approximation spaces are perfectly adapted to the asymptotic behavior of the solution. Remark 3.3 Note that, by construction, $$\gamma ( \cup _{j=1}^{J} \mathscr{I}_{j} ) = \partial K$$ and the intervals $$\mathscr{I}_{j}$$ intersect either trivially or only at an end point. Therefore, we can clearly identify $$L^{2}\left ( \partial K \right )$$ and $$L^{2} ( \cup _{j=1}^{J} \mathscr{I}_{j} )$$ through the (L−periodic arc length) parametrization γ. This will be the convention we shall follow without any further reference in the rest of the paper. In particular, we shall write η(s, k), ηslow(s, k), ν(s), etc. rather that η(γ(s), k), ηslow(γ(s), k), ν(γ(s)), etc. With the above definitions, the Galerkin formulation of problem (2.1) is equivalent to finding the unique $$\hat{\eta } \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ such that \begin{align} B_{k}(\hat{\mu},\hat{\eta}) = F_{k}(\hat{\mu}), \quad \textrm{for all } \hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}. \end{align} (3.4) While, in this connection, the following theorem constitutes the main result of the paper and provides the approximation properties of the solution of equation (3.4), its proof is deferred to the next section. Remark 3.4 In all the estimates that follow, the notation $$A \lesssim _{n}B$$ (respectively $$A \lesssim _{n,k_{0}} B$$) will mean that 0 ≤ A ≤ cB for all sufficiently large values of k (respectively for all k ≥ k0), where c is a positive constant independent of k and, where relevant, of s, but its value may be different for different values of n. Note that c depends on the geometry of the scatterer, the direction of incidence α and, where relevant, the parameters ξj, ξj′ etc. used to construct the Galerkin approximation spaces. Theorem 3.2 Given k0 > 1 and k ≥ k0, suppose that the sesquilinear form Bk in (2.3) associated with the integral operator $$\mathscr{R}_{k}$$ in (2.1) is continuous with a continuity constant Ck and coercive with a coercivity constant ck. Then, for all nj ∈{0, …, dj + 1} (j = 1, …, J), we have \begin{align} \Vert \eta - \hat{\eta} \Vert_{L^{2}(\partial K)} \lesssim_{n_{1},\ldots,n_{J},k_{0}} \dfrac{C_{k}}{c_{k}} \, k \sum_{j=1}^{J} \dfrac{E(k,j)}{\left( d_{j} \right)^{n_{j}}} \end{align} (3.5) for the Galerkin solution $$\hat{\eta }$$ to (3.4), where \begin{align*} E(k,j) = \begin{cases}{ll} \left( \log k \right)^{n_{j}+1/2}, & j =1,2,3,4 \ (\textrm{transition regions}), \\[0.5em] k^{-1/6},& j =5,6 \, \qquad (\textrm{shadow boundaries}); \end{cases} \end{align*} if J = 8, then \begin{align*} E(k,j) = 1, \qquad j = 7,8 \quad (\textrm{illuminated and shadow regions}). \end{align*} Since the order of magnitude of $$\Vert \eta \Vert _{L^{2}(\partial K)}$$ is k as $$k \to \infty$$ (see e.g. equation (1.15) in Melrose & Taylor (1985) or equation (3.43) in Ecevit & Reitich (2009)), if we assign the same polynomial degree to each interval $$\mathscr{I}_{j}$$, we obtain the following estimate for the relative error. Corollary 3.1 Under the assumptions of Theorem 3.2, if the same polynomial degree d = d1 = … = dJ is used on each interval, then for all n ∈ {0, …, d + 1}, there holds \begin{align} \dfrac{\Vert \eta - \hat{\eta} \Vert_{_{L^{2}(\partial K)}}}{\Vert \eta \Vert_{_{L^{2}(\partial K)}}} \lesssim_{n,k_{0}} \dfrac{C_{k}}{c_{k}} \dfrac{\left( \log k \right)^{n+1/2}}{d^{n}} \end{align} (3.6) for the Galerkin solution $$\hat{\eta }$$ to (3.4). Proof. As was shown in Ecevit & Reitich (2009, Theorem 3.3; see equation (3.43) therein), given a compact subset S of the illuminated region ∂KIL and a wavenumber k0 > 0, there exists a constant C1 depending on k0, α, S, and ∂K such that \begin{align*} \left|\eta^{\textrm{slow}}(s,k) - 2ik \, \alpha \cdot \nu(s)\right| \le C_{1} \end{align*} for all k ≥ k0 and all x ∈ S. Thus, for some constant C2 depending on k0, α and ∂K, we have $$0< C_{2} k \le \Vert \eta \Vert _{L^{2}(\partial K)}$$ for all sufficiently large k. Since $$\Vert \eta \Vert _{L^{2}(\partial K)}$$ depends continuously on k and is never zero for k > 0, by modifying the constant C2, if necessary, we obtain $$0< C_{2} k \le \Vert \eta \Vert _{L^{2}(\partial K)}$$ for all k ≥ k0. Therefore, we can replace k appearing on the right-hand side of (3.5) by $$\Vert \eta \Vert _{L^{2}(\partial K)}$$. Finally, taking d = d1 = … = dJ and n = n1 = … = nJ prints (3.6). Theorem 3.2 and Corollary 3.1 display the convergence characteristics of the Galerkin approximation spaces based on frequency-dependent changes of variables. In order to allow for a direct comparison, here we present a more flexible version of our algorithm in Ecevit & Özen (2017) together with the associated convergence results. To this end, given $$m \in \mathbb{N}$$, 0 ≤ εm < εm−1 < ⋯ < ε1 < 1/3, and constants ξ1, ξ2, ζ1, ζ2 > 0 with t1 − ξ1 < t2 − ξ2 and t2 + ζ2 < L + t1 − ζ1, we define the associated illuminated transition and shadow transition intervals as \begin{align*} I_{IT_{1}} & = [t_{1} + \xi_{1} k^{-1/3+\varepsilon_{m}}, t_{1} + \xi_{1} k^{-1/3+\varepsilon_{1}}], \\ I_{IT_{2}} & = [t_{2} - \xi_{2} k^{-1/3+\varepsilon_{1}}, t_{2} - \xi_{2} k^{-1/3+\varepsilon_{m}}], \\ I_{ST_{1}} & = [t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{1}}, t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{m}}], \\ I_{ST_{2}} & = [t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{m}}, t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{1}}], \end{align*} and, rather than utilizing a change of variables in the transition intervals, we divide each one into m − 1 subintervals by setting \begin{align*} I_{IT_{1}}^{j} & = [t_{1} + \xi_{1} k^{-1/3+\varepsilon_{j+1}}, t_{1} + \xi_{1} k^{-1/3+\varepsilon_{j}}], \\ I_{IT_{2}}^{j} & = [t_{2} - \xi_{2} k^{-1/3+\varepsilon_{j}}, t_{2} - \xi_{2} k^{-1/3+\varepsilon_{j+1}}], \\ I_{ST_{1}}^{j} & = [t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{j}}, t_{1} - \zeta_{1} k^{-1/3+\varepsilon_{j+1}}], \\ I_{ST_{2}}^{j} & = [t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{j+1}}, t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{j}}], \end{align*} for j = 1, …, m − 1. In addition to these 4m − 4 transition intervals, we further define the shadow boundary intervals as \begin{align*} I_{sb:{1}} & = [t_{1} - \zeta_{1} k^{-1/3 +\varepsilon_{m}}, t_{1} + \xi_{1} k^{-1/3+\varepsilon_{m}}], \\ I_{sb:{2}} & = [t_{2} - \xi_{1} k^{-1/3 +\varepsilon_{m}}, t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{m}}], \end{align*} and the illuminated region and shadow region intervals as \begin{align*} I_{IL} & = [t_{1} + \xi_{1}k^{-1/3+\varepsilon_{1}}, t_{2} - \xi_{2}k^{-1/3+\varepsilon_{1}}], \\ I_{DS} & = [t_{2} + \zeta_{2} k^{-1/3+\varepsilon_{1}}, L+t_{1}-\zeta_{1} k^{-1/3+\varepsilon_{1}}]. \end{align*} These give rise to a total of 4m intervals which we shall denote as $$\mathscr{I}_{j}$$ (j = 1, …, 4m). Reasoning as before, we identify $$L^{2} \left ( \partial K \right )$$ with $$L^{2} ( \cup _{j=1}^{4m} \mathscr{I}_{j} )$$. Given $$\mathbf{d} = \left (d_{1},\ldots ,d_{4m} \right ) \in \mathbb{Z}_{+}^{4m}$$, we now define the $$( \left | \mathbf{d} \right | + 4m)-$$dimensional Galerkin approximation space based on algebraic polynomials as \begin{align} \ \end{align} (3.7) The associated Galerkin formulation of (2.1) is to find the unique $$\hat{\eta } \in \mathscr{A}_{\mathbf{d}}$$ such that \begin{align} B_{k}(\hat{\mu},\hat{\eta}) = F_{k}(\hat{\mu}), \quad \textrm{for all } \hat{\mu} \in \mathscr{A}_{\mathbf{d}}. \end{align} (3.8) Incidentally, the Galerkin approximation spaces defined by equation (3.7) provide a more flexible version of those in Ecevit & Özen (2017), since here we allow for different ξ and ζ values rather than the values ξ1 = ξ2 and ζ1 = ζ2 we used in Ecevit & Özen (2017). This clearly renders the new approximation spaces better adapted to different geometries. Theorem 3.3 Suppose that k is sufficiently large and the sesquilinear form Bk in (2.3) associated with the integral operator $$\mathscr{R}_{k}$$ in (2.1) is continuous with a continuity constant Ck and coercive with a coercivity constant ck. Then, for all nj ∈{0, …, dj + 1} (j = 1, …, 4m), we have \begin{align*} \Vert \eta - \hat{\eta} \Vert_{L^{2}(\partial K)} \lesssim_{n_{1},\ldots,n_{4m}} \dfrac{C_{k}}{c_{k}} \, k \, \sum_{j = 1}^{4m} \dfrac{1+E(k,j)}{\left( d_{j} \right)^{n_{j}}} \end{align*} for the Galerkin solution $$\hat{\eta }$$ to (3.8), where \begin{align*} E(k,j) = \begin{cases} k^{-(1+3\varepsilon_{r+1})/2} \left( k^{(\epsilon_{r}-\varepsilon_{r+1})/2} \right)^{n_{j}}\!\!\!\!\!, & j=1,\ldots,4m-4 \, \quad (\textrm{transition regions}), \\[0.5 em] k^{-1/2} \left( k^{\varepsilon_{m}} \right)^{n_{j}}, & j=4m-3,4m-2 \quad (\textrm{shadow boundaries}), \\[0.5 em] k^{-(1+3\varepsilon_{1})/2} \left( k^{(1/3-\varepsilon_{1})/2} \right)^{n_{j}}, & j=4m-1,4m \qquad (\text{illuminated \& shadow reg.}). \end{cases} \end{align*} The proof of Theorem 3.3 is similar to that of Theorem 1 in Ecevit & Özen (2017) and is therefore skipped. On the other hand, exactly as Theorem 1 therein implies Corollary 1 in Ecevit & Özen (2017), Theorem 3.3 above yields the following result. Corollary 3.2 Under the assumptions of Theorem 3.3, if εj are chosen as \begin{align*} \varepsilon_{j} = \dfrac{1}{3} \, \dfrac{2m-2j+1}{2m+1}, \qquad j = 1,\ldots,m, \end{align*} and the same polynomial degree d = d1 = … = d4m is used on each interval, then for all n ∈{0, …, d + 1}, there holds \begin{align} \dfrac{\Vert \eta - \hat{\eta} \Vert_{_{L^{2}(\partial K)}}}{\Vert \eta \Vert_{_{L^{2}(\partial K)}}} \lesssim_{n} \dfrac{C_{k}}{c_{k}} \, m \, \dfrac{1+ k^{-\frac{1}{2}} \left( k^{\frac{1}{6m+3}} \right)^{n}}{d^{n}} \end{align} (3.9) for the Galerkin solution $$\hat{\eta }$$ to (3.8). Remark 3.5 Let us say that ‘f increases in proportion to g’ to mean that f and g are positive functions of k and f/g is bounded from above and below by positive constants independent on k. As mentioned in Remark 2.1, $$C_{k}/c_{k} = \mathscr{O}(k^{\delta })$$ as $$k \to \infty$$ with δ = 1/3 for CFIE and SCIE. In this case, for the estimate in Corollary 3.2, if m increases in proportion to $$\log k$$—which implies that $$k^{1/(6m+3)} = \exp (\log k/(6m+3))$$ is bounded independently of k—and the local polynomial degree d increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{\varepsilon _{2}}$$ with ε1, ε2 > 0 so that the total number of DoF, namely 4m(d + 1), increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{1+\varepsilon _{2}}$$, then the right-hand side of (3.9) is $$\mathscr{O}(k^{\delta -n\varepsilon _{1}} \, (\log k)^{1-n\varepsilon _{2}}$$) as $$k \to \infty$$. For the estimate (3.6) in Corollary 3.1, if the local polynomial degree d increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{1+\varepsilon _{2}}$$ which implies that the total number of DoF, namely J(d + 1) with J = 6 or J = 8, increases in proportion to $$k^{\varepsilon _{1}} (\log k)^{1+\varepsilon _{2}}$$ as above, then the right-hand side of (3.6) is $$\mathscr{O}(k^{\delta - n\varepsilon _{1}} \, (\log k)^{1/2-n\varepsilon _{2}})$$ as $$k \to \infty$$. In this sense, when compared with the frequency-adapted Galerkin BEM, the method based on frequency-dependent changes of variables provides an improvement of $$\sqrt{\log k}$$ in numerical accuracy as $$k \to \infty$$. Since n here can be arbitrarily large and $$k^{\delta _{1} - n\varepsilon _{1}} \, (\log k)^{\delta _{2}-n\varepsilon _{2}} = \mathscr{O}(1)$$ as $$k \to \infty$$ provided that either ε1 ≥ δ1/n and ε2 ≥ δ2/n or ε1 > δ1/n, it follows that, for any ε > 0, increasing the total number of DoF associated with either of the Galerkin schemes as $$\mathscr{O}(k^{\varepsilon })$$ is sufficient to obtain any prescribed accuracy independent on frequency. 4. Error analysis In this section we present the proof of Theorem 3.2. In light of inequality (2.5), it is sufficient to estimate \begin{align*} \inf_{\hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}} \Vert \eta- \hat{\mu} \Vert_{L^{2}(\partial K)}. \end{align*} To this end, we make use of the following classical result from approximation theory. Theorem 4.1 (Best approximation by algebraic polynomials (Schwab, 1998)) Given an interval $$\mathscr{I} = (a,b)$$ and $$n \in \mathbb{Z}_{+}$$, introduce the seminorms (for suitable f) by \begin{align} \left|\, f \right|{}_{n,\mathscr{I}} = \left[{\int_{a}^{b}} \left| D^{n} f(s) \right|{}^{2} \left( s-a \right)^{n} \left( b-s \right)^{n} \textrm{d}s \right]^{1/2} . \end{align} (4.1) Then, for all n ∈{0, …, d + 1}, there holds \begin{align*} \inf_{p \in \mathbb{P}_{d}} \Vert\, f-p \Vert \lesssim_{n} \left|\, f \right|{}_{n,\mathscr{I}} \, d^{-n}. \end{align*} Use of Theorem 4.1, in turn, requires the knowledge of derivative estimates of ηslow which we present next. In the rest of the paper we use the convention that an empty sum is zero. Theorem 4.2 Given k0 > 0, there holds \begin{align} \left| {D_{s}^{n}} \eta^{\textrm{slow}}(s,k) \right| \lesssim_{n,k_{0}} k + \sum\limits_{m=4}^{n+2} \left( k^{-1/3} + \left| w(s) \right| \right)^{-m} \end{align} (4.2) for all $$n \in \mathbb{Z}_{+}$$ and all k ≥ k0. Here w(s) = (s − t1)(t2 − s). Proof. The same estimate is shown to hold for all sufficiently large k in Domínguez et al. (2007, Theorem 5.4). Since $${D_{s}^{n}} \eta ^{\textrm{slow}}(s,k)$$ depends continuously on s and k, the result follows. Remark 4.1 Let us note that estimate (4.2) is valid over the entire boundary, but it is not sharp in the shadow region well away from the shadow boundary points, where the solution decays exponentially with increasing wavenumber k and the geodesic distance to the shadow boundary points (see Domínguez et al. (2007), Asheim & Huybrechs (2014), and the references therein). Nevertheless, as in Ecevit & Özen (2017), it allows us to eliminate, through utilization of additional approximation spaces in the shadow region, the asymptotic terms appearing in the error analysis due to approximation by zero in the deep shadow region (compare the estimates in Theorems 3.2 and 3.3 with Domínguez et al. (2007, Theorem 6.7)). Further, as we discussed in detail in Ecevit & Özen (2017), the benefit of this utilization is directly apparent from the significantly reduced condition numbers of the Galerkin linear systems. We continue with the derivation of estimates on the derivatives of the change of variables ϕ on the transition interval $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4). Proposition 4.3 Given k0 > 1, there holds \begin{align*} \left| {D^{n}_{s}} \phi \right| \lesssim_{n,k_{0}} \left( \log k \right)^{n} k^{\psi} \quad \textrm{on } \mathscr{I}_{j} \quad (\ j=1,2,3,4), \end{align*} for all $$n \in \mathbb{N}$$ and all k ≥ k0. Proof. Since the proof is similar for j = 1, 2, 3, 4, we concentrate on the case j = 1. Now since φ′′ = ψ′′ = 0, direct computations entail \begin{align} {D^{n}_{s}} \phi = \sum_{j=0}^{n} \binom{n}{j} \, D^{n-j}_{s} \varphi \ {D^{j}_{s}} k^{\psi} = \varphi \ {D^{n}_{s}} k^{\psi} + n \ {D^{1}_{s}} \varphi \ D^{n-1}_{s} k^{\psi}, \qquad n \ge 1, \end{align} (4.3) and \begin{align} {D^{n}_{s}} k^{\psi} = \left( {D^{1}_{s}} \psi \right)^{n} \left( \log k \right)^{n} k^{\psi}, \qquad n \ge 0. \end{align} (4.4) Using (4.4) in (4.3), we obtain \begin{align} {D^{n}_{s}} \phi = \left( \varphi \ {D^{1}_{s}} \psi \ \log k + n \ {D^{1}_{s}} \varphi \right) \left( {D^{1}_{s}} \psi \right)^{n-1} \left( \log k \right)^{n-1} k^{\psi}, \qquad n \ge 1. \end{align} (4.5) Since, for k ≥ k0 > 1, \begin{align*} \dfrac{1}{\left| b_{1}-a_{1} \right|} = \dfrac{1}{\xi_{1}^{\prime} - \xi_{1}k^{-1/3}} \le \dfrac{1}{\xi_{1}^{\prime} - \xi_{1}k_{0}^{-1/3}} \lesssim_{k_{0}} 1 \end{align*} it follows for $$s \in \mathscr{I}_{1} = I_{IT_{1}}$$ that \begin{align*} \left| {D^{1}_{s}} \psi \right| = \dfrac{1}{3} \, \dfrac{1}{b_{1}-a_{1}} \lesssim_{k_{0}} 1 \qquad \textrm{and} \qquad \left| {D^{1}_{s}} \varphi \right| = \dfrac{\xi_{1}^{\prime} - \xi_{1}}{b_{1}-a_{1}} \lesssim_{k_{0}} 1 \end{align*} and, clearly, $$\xi _{1} \le \varphi \le \xi _{1}^{\prime }$$. Use of these inequalities in (4.5) yields the desired result. Next we combine Theorem 4.2 and Proposition 4.3 to derive estimates on the derivatives of the composition ηslow ∘ ϕ. Proposition 4.4 Given k0 > 0, there holds \begin{align*} \left| {D_{s}^{n}} (\eta^{slow} \circ \phi) \right| \lesssim_{n,k_{0}} k \left( \log k \right)^{n} \quad \textrm{on } \mathscr{I}_{j} \quad \left( j=1,2,3,4 \right)\!, \end{align*} for all $$n \in \mathbb{N}$$ and all k ≥ k0. Proof. We fix k ≥ k0 > 1 and the interval $$\mathscr{I}_{j}$$ (j = 1, …, 4). Faá Di Bruno’s formula for the derivatives of a composition states \begin{align*} D^{n} \left(\, f \circ g \right) \left( t \right) = \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} (D^{m} f) (g(t)) \prod_{\ell=1}^{n} \dfrac{\ell}{m_{\ell}!} \left( \dfrac{D^{\ell}g(t)}{\ell!} \right)^{m_{\ell}}, \end{align*} where $$\mathscr{F}_{n} = \{ (m_{1},\ldots ,m_{n}) \in \mathbb{Z}^{n}_{+} : n = \sum _{\ell =1}^{n} \ell m_{\ell } \}$$; here $$m = \sum _{\ell =1}^{n} m_{\ell }$$. This yields \begin{align*} \left| {D^{n}_{s}} \left( \eta^{\textrm{slow}} \circ \phi \right) \right| \lesssim_{n} \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| \prod_{\ell=1}^{n} \left| D^{\ell}_{s} \phi \right|{}^{m_{\ell}} \end{align*} so that an appeal to Proposition 4.3 entails \begin{align*} \left| {D^{n}_{s}} \left( \eta^{\textrm{slow}} \circ \phi \right) \right| \lesssim_{n,k_{0}} \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| \prod_{\ell=1}^{n} \left( \left( \log k \right)^{\ell} k^{\psi} \right)^{m_{\ell}}. \end{align*} Since $$n = \sum _{\ell =1}^{n} \ell m_{\ell }$$ and $$m = \sum _{\ell =1}^{n} m_{\ell }$$, we therefore obtain \begin{align*} \left| {D^{n}_{s}} \left( \eta^{\textrm{slow}} \circ \phi \right) \right| & \lesssim_{n,k_{0}} \left( \log k \right)^{n} \sum_{(m_{1},\ldots,m_{n}) \in \mathscr{F}_{n}} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| k^{m\psi} \\ & \lesssim_{n,k_{0}} \left( \log k \right)^{n} \sum_{m=0}^{n} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| k^{m\psi}. \end{align*} It is hence sufficient to show, for $$m \in \mathbb{Z}_{+}$$, that \begin{align*} \left| \left({D^{m}_{s}} \eta^{\textrm{slow}}\right) (\phi) \right| k^{m\psi} \lesssim_{m,k_{0}} k. \end{align*} To this end, we note that if 0 ≤ ℓ ≤ m, then \begin{align*} \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-\ell} & = \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-m} \left( k^{-1/3} + \left| w (\phi) \right| \right)^{ m-\ell} \\ & \leq \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-m} \left( k_{0}^{-1/3} + L^{2} \right)^{m-\ell} \\ &\lesssim_{m,k_{0}} \left( k^{-1/3} + \left| w (\phi) \right| \right)^{-m} \end{align*} so that an appeal to Theorem 4.2 yields \begin{align*} \left| {D_{s}^{m}} \eta^{\textrm{slow}}(\phi) \right| k^{m\psi} & \lesssim_{m,k_{0}} \left[ k + \sum_{\ell=4}^{m+2} \left( k^{-1/3} + \left| w(\phi) \right| \right)^{-\ell} \right] k^{m\psi} \\ & \lesssim_{m,k_{0}} \left[ k + \left( k^{-1/3} + \left| w(\phi) \right| \right)^{-\left( m+2 \right)} \right] k^{m\psi} \\ & \lesssim_{m,k_{0}} k + \left( \dfrac{k^{\psi}}{k^{-1/3} + \left| w(\phi) \right|} \right)^{m} \left( k^{-1/3} + \left| w(\phi) \right| \right)^{-2} \\ & \lesssim_{m,k_{0}} k + \left( \dfrac{k^{\psi}}{\left| w(\phi) \right|} \right)^{m} k^{2/3}, \end{align*} where, in the third inequality, we used that ψ ≤ 0. Thus, it is now enough to show that the quotient $$k^{\psi }/\left | w(\phi ) \right |$$ is bounded by a constant independent of k. This estimation is similar on each of the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4) and we focus on $$\mathscr{I}_{1}$$. Indeed, on $$\mathscr{I}_{1}=I_{IT_{1}}$$, we have ϕ − t1 = φkψ, φ ≥ ξ1 > 0 and $$t_{2} - \phi \ge t_{2} - \left ( t_{1} + \xi _{1}^{\prime } \right ) \ge \xi _{2}^{\prime }> 0$$  so that \begin{align*} \dfrac{k^{\psi}}{\left| w(\phi) \right|} = \dfrac{k^{\psi}}{\left( \phi -t_{1} \right) \left( t_{2} - \phi \right)} = \dfrac{k^{\psi}}{\varphi k^{\psi} \left( t_{2} - \phi \right)} \le \dfrac{1}{\xi_{1} \, \xi_{2}^{\prime}}. \end{align*} This finishes the proof. Next, we estimate the seminorms (4.1) for the composition ηslow ∘ ϕ on the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4). Corollary 4.1 On the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4), given k0 > 1, there holds \begin{align*} \left| \eta^{\textrm{slow}} \circ \phi \right|{}_{n,\mathscr{I}_{j}} \lesssim_{n,k_{0}} k \left( \log k \right)^{n} \end{align*} for all $$n \in \mathbb{N}$$ and all k ≥ k0. Proof. On account of Proposition 4.4, we estimate for j = 1, 2, 3, 4 \begin{align*} \left| \eta^{\textrm{slow}} \circ \phi \right|{}^{2}_{n,\mathscr{I}_{j}} & = \int_{a_{j}}^{b_{j}} \left| {D^{n}_{s}}\left( \eta^{\textrm{slow}} \circ \phi \right) (s) \right|{}^{2} \left( s - a_{j} \right)^{n} \left( b_{j} -s \right)^{n} \textrm{d}s \\ & \lesssim_{n,k_{0}} k^{2} \left( \log k \right)^{2n} \int_{a_{j}}^{b_{j}} \left( s - a_{j} \right)^{n} \left( b_{j} -s \right)^{n} \textrm{d}s \\ & \lesssim_{n,k_{0}} k^{2} \left( \log k \right)^{2n}, \end{align*} where we used that 0 < bj − aj < L. Thus, the result. We are now ready to prove Theorem 3.2. Proof. (of Theorem 3.2): Céa’s lemma (cf. inequality (2.5)) implies \begin{align} \Vert \eta -\hat{\eta} \Vert_{L^{2} \left( \partial K \right)} \leq \dfrac{C_{k}}{c_{k}} \inf_{\hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}} \Vert \eta- \hat{\mu} \Vert_{L^{2} \left( \partial K \right)} \end{align} (4.6) for the unique solution $$\hat{\eta }$$ of the Galerkin formulation (3.4). As we identify $$L^{2} \left ( \partial K \right )$$ with $$L^{2} ( \cup _{j=1}^{J} \mathscr{I}_{j} )$$ through the L−periodic arc length parameterization γ, there holds \begin{align*} \Vert \eta- \hat{\mu} \Vert_{L^{2}(\partial K)} = \Vert \eta- \hat{\mu} \Vert_{L^{2} \left( \cup_{j=1}^{J} \mathscr{I}_{j} \right)} \le \sum_{j=1}^{J} \Vert \eta- \hat{\mu} \Vert_{L^{2} \left( \mathscr{I}_{j} \right)} \end{align*} for any $$\hat{\mu } \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$. Accordingly, the very definition of Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ entails \begin{align} \inf_{\hat{\mu} \in \mathscr{A}_{\mathbf{d}}^{\mathscr{C}}} \Vert \eta- \hat{\mu} \Vert_{L^{2}(\partial K)} \le \sum_{j=1}^{4} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} - p \circ \phi^{-1} \Vert_{L^{2}(\mathscr{I}_{j})} + \sum_{j=5}^{J} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} - p \Vert_{L^{2}(\mathscr{I}_{j})}. \end{align} (4.7) On the other hand, utilizing the change of variables ϕ on the transition intervals $$\mathscr{I}_{j}$$ (j = 1, 2, 3, 4), for any $$p \in \mathbb{P}_{d_{j}}$$, we have \begin{align*} \Vert \eta^{\textrm{slow}} - p \circ \phi^{-1} \Vert^{2}_{L^{2}(\mathscr{I}_{j})} & = \int_{a_{j}}^{b_{j}} \left| \left( \eta^{\textrm{slow}} - p \circ \phi^{-1}\right)(s) \right|{}^{2} \textrm{d}s \\ & = \int_{a_{j}}^{b_{j}} \left| \left( \eta^{\textrm{slow}} \circ \phi - p \right)(s) \right|{}^{2} {D^{1}_{s}}\phi(s) \, \textrm{d}s \\ & \lesssim_{k_{0}} \log k \ \Vert \eta^{\textrm{slow}} \circ \phi - p \Vert^{2}_{L^{2}(\mathscr{I}_{j})}, \end{align*} where we used Proposition 4.3 in conjunction with the fact that kψ < 1. Combining this last estimate with (4.6) and (4.7), we deduce \begin{align*} \Vert \eta -\hat{\eta} \Vert_{L^{2} \left( \partial K \right)} \lesssim_{k_{0}} \dfrac{C_{k}}{c_{k}} \left\{ \sum_{j=1}^{4} \left( \log k \right)^{1/2} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} \circ \phi - p \Vert_{L^{2}(\mathscr{I}_{j})} + \sum_{j=5}^{J} \inf_{p \in \mathbb{P}_{d_{j}}} \Vert \eta^{\textrm{slow}} - p \Vert_{L^{2}(\mathscr{I}_{j})} \right\} \end{align*} and this, on account of Theorem 4.1, implies \begin{align*} \Vert \eta -\hat{\eta} \Vert_{L^{2} \left( \partial K \right)} \lesssim_{n_{1},\ldots,n_{J},k_{0}} \dfrac{C_{k}}{c_{k}} \left\{ \sum_{j=1}^{4} \left( \log k \right)^{1/2} \left| \eta^{\textrm{slow}} \circ \phi \right|{}_{n_{j},\mathscr{I}_{j}} d_{j}^{-n_{j}} + \sum_{j=5}^{J} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}} d_{j}^{-n_{j}} \right\}\!\!. \end{align*} Therefore, to complete the proof, it suffices to show that \begin{align} \left| \eta^{\textrm{slow}} \circ \phi \right|{}_{n_{j},\mathscr{I}_{j}} \lesssim_{n_{j},k_{0}} k \left( \log k \right)^{n_{j}}, \qquad j = 1,2,3,4, \end{align} (4.8) and \begin{align} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}} \lesssim_{n_{j},k_{0}} k \, k^{-1/6}, \qquad j = 5,6, \end{align} (4.9) and (if J = 8) \begin{align} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}} \lesssim_{n_{j},k_{0}} k, \qquad j = 7,8. \end{align} (4.10) While the estimates in (4.8) are given by Corollary 4.1, for the shadow boundary intervals $$\mathscr{I}_{j}$$ (j = 5, 6) we use Theorem 4.2 to deduce \begin{align*} \left| D_{s}^{n_{j}} \eta^{\textrm{slow}}(s,k) \right| \lesssim_{n_{j},k_{0}} k + \sum_{m=4}^{n_{j}+2} \left( k^{-1/3} + \left| w(s) \right| \right)^{-m} \lesssim_{n_{j},k_{0}} k + k^{(n_{j}+2)/3}; \end{align*} this implies \begin{align*} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}}^{2} & = \int_{a_{j}}^{b_{j}} \left| D^{n_{j}}_{s} \eta^{\textrm{slow}}(s) \right|{}^{2} \left( s -a_{j} \right)^{n_{j}} \left( b_{j}-s \right)^{n_{j}} \textrm{d}s \\ & \lesssim_{n_{j},k_{0}} \left( k + k^{(n_{j}+2)/3} \right)^{2} \left( b_{j} -a_{j} \right)^{2n_{j}+1} \\ & \lesssim_{n_{j},k_{0}} \left( k + k^{(n_{j}+2)/3} \right)^{2} \left( k^{-1/3} \right)^{2n_{j}+1} \\ & \lesssim_{n_{j},k_{0}} \left( k \, k^{-1/6} \right)^{2}, \end{align*} which justifies the estimates in (4.9). This completes the proof when J = 6. When J = 8, for the illuminated and shadow region intervals $$\mathscr{I}_{j}$$ (j = 7, 8), we use Theorem 4.2 to estimate \begin{align*} \left| D_{s}^{n_{j}} \eta^{\textrm{slow}}(s,k) \right| \lesssim_{n_{j},k_{0}} k + \sum_{m=4}^{n_{j}+2} \left( k^{-1/3} + \left| w(s) \right| \right)^{-m} \lesssim_{n_{j},k_{0}} k + \sum_{m=4}^{n_{j}+2} \left| w(s) \right|{}^{-m} \lesssim_{n_{j},k_{0}} k \end{align*} so that \begin{align*} \left| \eta^{\textrm{slow}} \right|{}_{n_{j},\mathscr{I}_{j}}^{2} & = \int_{a_{j}}^{b_{j}} \left| D^{n_{j}}_{s} \eta^{\textrm{slow}}(s) \right|{}^{2} \left( s -a_{j} \right)^{n_{j}} \left( b_{j}-s \right)^{n_{j}} \textrm{d}s \\ & \lesssim_{n_{j},k_{0}} k^{2} \left( b_{j} -a_{j} \right)^{2n_{j}+1} \\ & \lesssim_{n_{j},k_{0}} k^{2}, \end{align*} which verifies (4.10). This finishes the proof. 5. Numerical tests In this section, we first discuss the implementation details related with the Galerkin BEMs based on frequency-dependent changes of variables. Next we present numerical tests exhibiting their performance in comparison with the frequency-adapted Galerkin BEMs. Finally, we discuss the accuracy of numerical solutions in the shadow region and compare the results with those in Asheim & Huybrechs (2014). 5.1. Implementation details As we mentioned, the developments central to this paper are independent of the integral equation used. In order to allow a simple performance comparison with the aforementioned algorithms, we base our numerical implementations on the CFIE wherein the integral operator and the right-hand side are given by \begin{align*} \mathscr{R}_{k} = \dfrac{1}{2} \, I + \mathscr{D} - ik \mathscr{S} \quad \textrm{and} \quad f_{k} = \dfrac{\partial u^{\textrm{inc}}}{\partial \nu} - ik u^{\textrm{inc}}. \end{align*} Here, $$\mathscr{S}$$ is the acoustic single-layer integral operator and $$\mathscr{D}$$ is its normal derivative, and they are defined by \begin{align*} \mathscr{S} \eta(x) =&\, \int_{\partial K} \varPhi(x,y) \, \eta(y) \, \textrm{d}s(y), \qquad x \in \partial K,\\ \mathscr{D} \eta(x) =&\, \int_{\partial K} \dfrac{ \partial \varPhi(x,y)}{\partial \nu(x)} \, \eta(y) \, \textrm{d}s(y), \qquad x \in \partial K, \end{align*} where ν(x) is the outward unit normal to ∂K. In this case, the conjugates of the entries of the Galerkin matrices arising in connection with the approximation spaces $$\mathscr{A}_{d}$$ or $$\mathscr{A}_{d}^{\mathscr{C}}$$ are of the form \begin{align} \overline{B_{k}(\hat{\mu}_{j},\hat{\mu}_{j^{\prime}})} = \dfrac{1}{2} \int_{\mathscr{I}_{j} \cap \mathscr{I}_{j^{\prime}}} \overline{\hat{\mu}_{j}(s)} \, \hat{\mu}_{j^{\prime}}(s) \, \textrm{d}s + \iint_{\mathscr{I}_{j} \times \mathscr{I}_{j^{\prime}}} F(s,t,k) \, \overline{\hat{\mu}_{j}(s)} \, \hat{\mu}_{j^{\prime}}(t) \, \textrm{d}t \, \textrm{d}s, \end{align} (5.1) where $$\hat{\mu }_{j}$$ and $$\hat{\mu }_{j^{\prime }}$$ are elements of $$\mathscr{A}_{d}$$ or $$\mathscr{A}_{d}^{\mathscr{C}}$$ supported on $$\mathscr{I}_{j}$$ and $$\mathscr{I}_{j^{\prime }}$$, and F(s, t, k) = ∂ν(s)Φ(γ(s), γ(t)) − ikΦ(γ(s), γ(t)). Note that $$d H^{(1)}_{0}(z)/dz = - H^{(1)}_{1}(z)$$ and \begin{align*} H_{r}^{(1)}(z) = e^{iz} \left\{ \sum_{\ell=0}^{N} \dfrac{c_{r,\ell}}{z^{\ell+1/2}} + \theta_{r,N}(z) \, \dfrac{c_{r,N+1}}{z^{N+3/2}}\right\}\!\!,\qquad r,N \ge 0 \end{align*} with $$\left | \theta _{r,N}(z) \right | < 1$$ for all z ≫ 1 for some constants cr, ℓ (Gradshteyn & Ryzhik, 2000). Therefore, setting \begin{align*} G(s,t,k) = e^{-ik|\gamma(s)-\gamma(t)|} F(s,t,k) \qquad \textrm{and} \qquad \varPsi(s,t) = \alpha \cdot \left(\gamma(t)-\gamma(s)\right)+\left|\gamma(s)-\gamma(t)\right| \end{align*} and using p and q to denote polynomials, the right-hand side of (5.1) takes on the generic forms \begin{align} \dfrac{1}{2} \int_{\mathscr{I}_{j} \cap \mathscr{I}_{j^{\prime}}} p(s) \, q(s) \, \textrm{d}s + \iint_{\mathscr{I}_{j} \times \mathscr{I}_{j^{\prime}}} e^{ik \, \varPsi(s,t)} \, G(s,t,k) \, p(s) \, q(t) \, \textrm{d}t \, \textrm{d}s \end{align} (5.2) for $$\mathscr{A}_{d}$$, and \begin{align} \dfrac{1}{2} \int_{\mathscr{I}_{j} \cap \mathscr{I}_{j^{\prime}}} p\left(\phi_{j}^{-1}(s)\right) \, q\left(\phi_{j^{\prime}}^{-1}(s)\right) \, \textrm{d}s + \iint_{\mathscr{I}_{j} \times \mathscr{I}_{j^{\prime}}} e^{ik \, \varPsi(s,t)} \, G(s,t,k) \, p\left(\phi_{j}^{-1}(s)\right) \, q\left(\phi_{j^{\prime}}^{-1}(t)\right) \, \textrm{d}t \, \textrm{d}s \end{align} (5.3) for $$\mathscr{A}^{\mathscr{C}}_{d}$$. Here ϕj = ϕ if $$\mathscr{I}_{j}$$ is an interval of change of variables, and it is the identity map otherwise (the same comment applies to $$\phi _{j^{\prime }}$$). While the one-dimensional integrals in (5.2) and (5.3) can be evaluated through any suitable quadrature rule, the highly-oscillatory double integrals therein require special treatment. Indeed, Kim (2012) provides a detailed theoretical discussion on the efficient evaluation of the double integrals in (5.2) by analyzing the regularity of the function G and the stationary points of the phase function Ψ for compact, smooth, and strictly convex obstacles. The discussions there are confined to the method in Domínguez et al. (2007) and thus limited to the use of three regions only, but provide an excellent source for the extension and numerical analysis of the methodologies discussed there to the more general case here involving 4m regions. Concerning the double integrals in (5.3) one would be inclined to substitute $$u = \phi _{j}^{-1}(s)$$ and $$v = \phi _{j^{\prime }}^{-1}(t)$$, but this transforms the phase functions into $$\varPsi (\phi _{j}(u),\phi _{j^{\prime \prime }}(v))$$. In this case the associated double integrals cannot be evaluated efficiently through standard quadrature rules for highly oscillatory integrals since they are based on the assumption that the phase functions do not depend on k (here ϕj and $$\phi _{j^{\prime }}$$ depend on k if they correspond to a transition interval). It is therefore reasonable to invert these functions so as to preserve the form of the phase functions. In this connection, observe that if $$\mathscr{I}_{j} = I_{\textrm{IT}_{1}}$$, then ϕj can be expressed as \begin{align} \phi_{j}(s) = t_{1} + \left(\xi_{1} + \left(\xi^{\prime}_{1}-\xi_{1}\right) s^{\prime}\right) \, k^{(s^{\prime}-1)/3}, \qquad s \in [a_{1},b_{1}] = \left[t_{1} + \xi_{1} k^{-1/3}, t_{1} + \xi_{1}^{\prime}\right] \end{align} (5.4) with $$0 \le s^{\prime } = \frac{s-a_{1}}{b_{1}-a_{1}} \le 1$$, and ϕj admits similar representations in the other transition intervals. It follows from (5.4) that ϕj does not oscillate with increasing k and can be inverted numerically. Here we suggest a more direct approach observing first that when $$\xi _{1} = \xi ^{\prime }_{1}$$ \begin{align*} \phi_{j}^{-1}(s) = a_{1} + (b_{1}-a_{1}) \left( 3\dfrac{\log(s-t_{1})-\log \xi_{1}}{\log k} + 1 \right)\!. \end{align*} Otherwise, $$\xi _{1} < \xi ^{\prime }_{1}$$ and equation (5.4) can be rewritten as \begin{align} \dfrac{1}{3} \dfrac{\phi_{j}(s)-t_{1}}{\xi^{\prime}_{1}-\xi_{1}} = \left( \dfrac{1}{3} \dfrac{\xi^{\prime}_{1}}{\xi^{\prime}_{1}-\xi_{1}} + u \right) k^{u}, \qquad s \in [a_{1},b_{1}] \end{align} (5.5) with $$-\frac{1}{3}\le u = \frac{s^{\prime }-1}{3} \le 0$$. In this case, since equations (5.4) and (5.5) are related by affine transformations, the problem of numerically inverting ϕj is equivalent to that of f(u) = (c + u)ku on the interval [−1/3, 0] where c > 1/3 is a constant. Using the affine transformations v = c + u and $$w = v \log k$$, this reduces to inverting g(w) = wew on $$[(c-1/3) \log k, c \log k]$$. The inverse of g is known as the Lambert or omega function and can be evaluated efficiently (Corless et al., 1996). With this information, the highly oscillatory double integrals in (5.3) can be efficiently evaluated based on extensions of the methodologies in Kim (2012) to J = 6 and J = 8 regions. Further, the explicit forms of the functions ϕj (as in (5.4)) also allow for the extension of the regularity analysis and numerical quadrature discussed therein to cover the highly oscillatory double integrals in (5.3). Leaving these details for future work, here we evaluate the double integrals appearing as the entries of Galerkin matrices using a combination of the Nyström and trapezoidal rule discretizations (Colton & Kress, 1992) for the inner integrals and the trapezoidal rule for the outer integrals. This latter rule is also used in the evaluation of one-dimensional integrals. In order to preserve the high-order approximation properties of these numerical integration rules for smooth and periodic integrands, as in Ecevit & Özen (2017), we additionally utilize a smooth partition of unity confined to the regions on the boundary of the scatterer described in Section 3. In each case, based on our experience in Ecevit & Özen (2017), the number of discretization points is chosen approximately as 10 to 12 points per wave length. An important component of our algorithms relates to the choice of bases for the spaces of algebraic polynomials and, in order to minimize the numerical instabilities arising from the use of high-degree polynomials (Ecevit & Özen, 2017), on any generic interval $$\mathscr{I} = [a,b]$$, we use the bases {ρr : r = 0, …, d} and {(ρ∘ϕ−1)r : r = 0, …, d} for $$\mathbb{P}_{d}$$ and $$\mathbb{P}_{d} \circ \phi ^{-1}$$, respectively, where ρ is the affine function that maps the interval $$\mathscr{I}$$ onto [−1, 1]. In the same vein, the choice of the parameters ξ, ξ′, ζ, ζ′ appearing in the definitions of the intervals $$I_{IT_{j}},I_{ST_{j}},I_{sb:{j}}$$ (j = 1, 2) and IIT, IDS related with the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces is of great importance, since a random choice may result in a loss of accuracy due to poor resolution of the boundary layers in the solution. We therefore optimize these parameters for a small wavenumber (k = 10 for all the implementations here) through a simple iterative procedure, and use these values for all larger wavenumbers. For instance, when J = 6, we take $$\xi _{1}^{\prime } = \xi _{2}^{\prime }$$, $$\zeta _{1}^{\prime }=\zeta _{2}^{\prime }$$ and initially require that $$\alpha \cdot \gamma (\xi _{1}^{\prime }) = -1$$ and $$\alpha \cdot \gamma (\zeta _{1}^{\prime }) = 1$$. We then take ξ1 to be the mid-point between t1 and $$\xi _{1}^{\prime }$$, and change it in small increments until the local error in IT1 ∪ SB1 is minimized. We treat the triplet $$(t_{2},\xi _{2},\xi _{2}^{\prime })$$ similarly. Finally, we fix ξ1 and ξ2, and change $$\xi _{1}^{\prime }=\xi _{2}^{\prime }$$ in small increments until the error in IL ∪ IT1 ∪ IT2 ∪ SB1 ∪ SB2 is minimized. The optimization of ζ parameters is realized similarly. The computed values are taken as the initial guess for the next iterate, and the procedure is repeated with smaller increments, with step size half that of the previous one, until the variations in the global error stabilizes. Following the prescriptions in Ecevit & Özen (2017, §4.2), then we introduce a smooth partition of unity confined to the regions$$I_{IT_{j}},I_{ST_{j}}, I_{sb:{j}}$$(j = 1, 2) and IIT, IDS, and optimize the shapes of hat functions therein using a similar iterative procedure. We follow a similar procedure for the optimizaton of the parameters appearing in connection with the $$\mathscr{A}_{\mathbf{d}}$$ spaces. For further details, we refer to Eruslu (2015). Fig. 3. View largeDownload slide Single-scattering geometries and associated incidence directions α. Fig. 3. View largeDownload slide Single-scattering geometries and associated incidence directions α. 5.2. Performance comparison of $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces Here we present numerical tests comparing the performances of frequency-adapted Galerkin BEMs and Galerkin BEMs based on frequency-dependent changes of variables. For the former, the Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}$$ are generated by taking m = 2 and assigning the same polynomial degree d = 3, 6, 9, 12, 15 to each one of the eight direct summands which give rise to the number of DoF 8(d + 1) = 32, 56, 80, 104, 128. Let us mention in this connection that the numerical tests we presented in Ecevit & Özen (2017) are restricted to the choices m = 1 and m = 2 (with some variations) for the very reason that utilization of larger values of m is simply impractical from an implementation point of view. On the other hand, these tests have displayed that the condition numbers of the Galerkin linear systems corresponding to m = 2 are significantly smaller when compared to m = 1, and this enhances numerical stability with increasing polynomial degrees and thereby allows for more accurate approximations. For the latter, the Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ are generated by taking J = 6 and assigning the same polynomial degree d = 4, 8, 12, 16, 20 to each one of the six direct summands resulting in the number of DoF 6(d + 1) = 30, 54, 78, 102, 126. We consider three different single-scattering geometries (see Fig. 3) consisting of the unit circle, the ellipse with major/minor axes (aligned with the x/y-axes) of 2/1, and the kite shaped obstacle given parametrically as $$\{(\cos t + 0.65 \cos 2t - 0.65, 1.5 \sin t) \, : \, t \in [0,2\pi ]\}$$. The unit circle is the standard test example for single-scattering solvers since circles are the only two-dimensional obstacles for which explicit solutions are available through a straightforward Fourier analysis (see equation (5.7) below). As for the ellipse and kite shaped obstacles, we compare the outcome of our numerical implementations with highly accurate reference solutions obtained by a combination of the Nyström and trapezoidal rule discretizations applied to the CFIE (Colton & Kress, 1992). Let us note that the kite shaped obstacle and the associated direction of incidence is chosen to display that the methodologies discussed in this paper can be applied to more general single-scattering problems involving nonconvex smooth scatterers, provided that there are exactly two shadow boundaries and no multiple reflections. The numerical results in Table 1 display the DoF vs. the relative L2 error \begin{align} \dfrac{\Vert \eta - \hat{\eta} \Vert_{L^{2}(\partial K)}}{\Vert \eta \Vert_{L^{2}(\partial K)}} \end{align} (5.6) related with the approximation spaces $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ for the wavenumbers k = 50, 100, 200, 400, 800. For any fixed value of k, the table depicts for both $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces that the error decreases with increasing DoF. The data corresponding to the unit circle shows that, for fixed DoF, the approximation errors associated with the $$\mathscr{A}_{\mathbf{d}}$$ spaces vary depending on the value of k, and those related with the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces tend to decrease with increasing k. Further, $$\mathscr{A}_{\mathbf{d}}$$ spaces provide better approximations for smaller values of k, and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces perform better for larger values of k (when DoF ≥ 54). In the case of the ellipse, for fixed DoF, the approximation errors corresponding to the $$\mathscr{A}_{\mathbf{d}}$$ spaces increase with increasing k, and those related with the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces decrease with increasing k; and when DOF ≥ 54, the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces perform better for all values of k. In light of these observations, we can therefore infer that, for sufficiently large DoF, the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces are better adopted to the high-frequency structure of the scattering problem for convex obstacles. Finally, in the case of the nonconvex kite shaped scatterer, no reasonable performance comparison can be made except that for large DoF the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces seem to provide better approximations. Table 1 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces for the single-scattering configurations in Fig. 3 and for various different wavenumbers k Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Table 1 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces for the single-scattering configurations in Fig. 3 and for various different wavenumbers k Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Circle Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 5.94e − 3 5.49e − 3 4.76e − 3 3.91e − 3 3.20e − 3 56 6.71e − 4 7.32e − 4 9.52e − 4 1.09e − 3 1.08e − 3 80 1.15e − 4 1.69e − 4 1.64e − 4 2.07e − 4 2.54e − 4 104 1.21e − 5 3.41e − 5 4.46e − 5 3.95e − 5 4.56e − 5 128 7.37e − 6 1.20e − 5 3.35e − 5 2.11e − 5 1.73e − 5 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.75e − 3 4.20e − 3 5.10e − 3 4.21e − 3 4.74e − 3 54 8.59e − 4 6.82e − 4 5.63e − 4 4.34e − 4 3.43e − 4 78 2.93e − 4 2.30e − 4 1.95e − 4 1.45e − 4 9.13e − 5 102 6.62e − 5 5.60e − 5 4.16e − 5 2.89e − 5 2.06e − 5 126 3.82e − 5 2.74e − 5 2.08e − 5 1.59e − 5 1.20e − 5 Ellipse Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 7.19e − 3 7.45e − 3 9.01e − 3 1.14e − 2 1.44e − 2 56 1.39e − 3 1.53e − 3 1.68e − 3 1.73e − 3 1.82e − 3 80 4.21e − 4 5.84e − 4 6.79e − 4 7.78e − 4 8.17e − 4 104 1.14e − 4 2.28e − 4 2.87e − 4 3.17e − 4 3.47e − 4 128 3.77e − 5 5.86e − 5 9.46e − 5 1.09e − 4 1.23e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 1.08e − 2 9.36e − 3 1.25e − 2 1.05e − 2 1.17e − 2 54 7.71e − 4 7.35e − 4 6.38e − 4 5.32e − 4 4.87e − 4 78 2.25e − 4 1.94e − 4 1.59e − 4 1.66e − 4 1.23e − 4 102 4.84e − 5 4.46e − 5 4.15e − 5 3.85e − 5 2.99e − 5 126 3.14e − 5 3.09e − 5 2.07e − 5 1.72e − 5 1.47e − 5 Kite Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 32 3.12e − 2 1.92e − 2 2.19e − 2 2.56e − 2 2.91e − 2 56 1.93e − 2 3.59e − 3 2.21e − 3 2.69e − 3 4.98e − 3 80 1.01e − 2 3.08e − 3 4.18e − 4 8.68e − 4 1.63e − 3 104 7.04e − 3 2.88e − 3 1.45e − 4 1.43e − 4 3.84e − 4 128 5.75e − 3 2.46e − 3 1.17e − 4 7.61e − 5 1.86e − 4 DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) k = 50 k = 100 k = 200 k = 400 k = 800 30 5.05e − 2 2.53e − 2 2.16e − 2 1.85e − 2 2.54e − 2 54 1.71e − 2 5.04e − 3 2.91e − 3 4.54e − 3 1.92e − 3 78 6.49e − 3 2.63e − 3 5.77e − 4 5.87e − 4 4.23e − 4 102 1.96e − 3 1.55e − 3 1.63e − 4 1.46e − 4 7.41e − 5 126 7.82e − 4 6.35e − 4 6.51e − 5 5.99e − 5 2.77e − 5 Let us finally note that the error term on the right-hand side of equation (3.6) in Corollary 3.1 suggests (excluding the fact that $$C_{k}/c_{k} = \mathscr{O}(k^{1/3})$$ as $$k \to \infty$$ for CFIE) that DoF should increase in proportion to $$\log k$$ with increasing k in order to maintain an error independent on k. Although the data in Table 1 essentially show for the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces that the error decreases with increasing k for sufficiently large DoF, in Table 2 we display DoF vs. the relative L2 error for the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces when $$\textrm{DoF} \approx 80 \times \log _{10}(k/20)$$. Table 2 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces when DoF increases logarithmically with increasing k: DoF $$\approx 80 \times \log _{10}(k/20)$$ Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 View Large Table 2 DoF vs. relative L2 error related with the $$\mathscr{A}_{\mathbf{d}}$$ and $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ spaces when DoF increases logarithmically with increasing k: DoF $$\approx 80 \times \log _{10}(k/20)$$ Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 Circle k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 5.94e − 3 30 5.75e − 3 100 56 7.32e − 4 54 6.82e − 4 200 80 1.64e − 4 78 1.95e − 4 400 104 3.95e − 5 102 2.89e − 5 800 128 1.73e − 5 126 1.20e − 5 Ellipse k DoF ($$\mathscr{A}_{\mathbf{d}}$$ ) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 7.19e − 3 30 1.08e − 2 100 56 1.53e − 3 54 7.35e − 4 200 80 6.79e − 4 78 1.59e − 4 400 104 3.17e − 4 102 3.85e − 5 800 128 1.23e − 4 126 1.47e − 5 Kite k DoF ($$\mathscr{A}_{\mathbf{d}}$$) Relative L2 error DoF ($$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$) Relative L2 error 50 32 3.12e − 2 30 5.05e − 2 100 56 3.59e − 3 54 5.04e − 3 200 80 4.18e − 4 78 5.77e − 4 400 104 1.43e − 4 102 1.46e − 5 800 128 1.86e − 4 126 2.77e − 5 View Large 5.3. Approximations in the shadow region Here we focus on approximations in the shadow region, and we present numerical tests exhibiting the performance of Galerkin BEMs based on frequency dependent changes of variables. To this end, in Fig. 4, we consider the unit circle parametrized as $$(\cos \theta ,\sin \theta )$$, 0 ≤ θ ≤ 2π, and illuminated by a plane wave incidence with direction α = (1, 0), and we display the pointwise relative error as a function of the radial angle θ for k = 10, 20, 40, 80, 120. In these experiments, we have constructed the Galerkin approximation spaces $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ defined in equation (3.3) by using J = 6 regions. We have assigned the polynomial degrees d = 3, 7, 11 to each of the six regions and this resulted in the total numbers of DoF = 6(d + 1) = 24, 48, 72. Fig. 4. View largeDownload slide Pointwise relative error on the boundary of the unit circle for increasing DoF and increasing wavenumbers k. Fig. 4. View largeDownload slide Pointwise relative error on the boundary of the unit circle for increasing DoF and increasing wavenumbers k. The pointwise relative errors displayed in Fig. 4 show that the quality of approximations in the illuminated region, namely the interval (π/2, 3π/2), improves with increasing DoF and it does not essentially depend on the wavenumber k. The approximations in the shadow region also improve with increasing DoF while degrading for increasing values of the wavenumber k. Except for k = 40, there is essentially no difference between the approximations in the intervals [0, θ0] and [θ1, 2π] for DoF = 48 and 72 (here θ0 = 0.86 and θ1 = 2π − θ0 and they are depicted by the vertical dotted lines in these plots). As we discuss below, this behavior is related with the asymptotic behavior of the amplitude ηslow in the shadow region. On the other hand, let us emphasize that the change of variables used to construct the $$\mathscr{A}_{\mathbf{d}}^{\mathscr{C}}$$ Galerkin approximation spaces are designed to resolve the boundary layers around the shadow boundaries (i.e. the regions where the derivatives of ηslow actually blow up with increasing k—as implied by the estimates in Theorem 4.2) by assigning more DoF in these regions when compared to the deep shadow regions with increasing k. Moreover, from a numerical perspective, it is also important to note that the condition numbers of the Galerkin linear systems increase with increasing DoF, e.g. they are about 1e + 3, 3.2e + 5 and 3.2e + 8 (almost independently of the value of k) for DoF = 24, 48 and 72, respectively. Pointwise relative error on the boundary of the unit circle for increasing DoF and increasing wavenumbers k. As for the asymptotic behavior of the amplitude ηslow in the shadow region, its phase in this region is modulated depending on the wavenumber, the curvature, and the geodesic distance to the shadow boundaries. Moreover, while its size increases linearly in the illuminated region with increasing wavenumber k, it decays exponentially in the shadow region. Furthermore, this decay depends additionally on the curvature and the geodesic distance to the shadow boundary points. A detailed discussion on these issues can be found in Asheim & Huybrechs (2014) where, based on suitable combinations of the GTD and Melrose–Taylor asymptotics, the authors design numerical schemes more appropriate for approximations in the shadow region. To illustrate this decaying behavior, let us consider a circular scatterer of radius r parametrized as $$(r\cos \theta ,r\sin \theta )$$, 0 ≤ θ ≤ 2π, and illuminated by a plane wave incidence with direction α = (1, 0). In this case, simple calculations based on Fourier analysis (Bowman et al., 1970) shows that η can be expressed as a function of the radial angle θ in the form \begin{align} \eta(\theta) = - \dfrac{2i}{\pi r} \sum_{n = -\infty}^{\infty} \dfrac{e^{in(\frac{\pi}{2}+\theta)}}{H^{(1)}_{n}(kr)}, \qquad 0 \le \theta \le 2\pi, \end{align} (5.7) where $$H^{(1)}_{n}$$ is the Hankel function of the first kind and order n. In Fig. 5, we display the absolute value of η as a function of the radial angle θ for increasing values of k for circles of radii r = 1 and r = 1/2. As this figure depicts, the region where |η(θ)|≥ 1 shrinks towards the illuminated region (i.e. the interval (π/2, 3π/2)) with increasing k, and its size is smaller for r = 1. Moreover, $$\max |\eta |$$ is approximately 2k for both r = 1 and r = 1/2, but as k increases from 10 to 1280 the values of $$\min |\eta |$$ decreases from 3.2e − 2 to 1.1e − 10 for r = 1 and from 1.7e − 1 to 6.3e − 8 for r = 1/2. Therefore, it can be inferred that accurately approximating the solution in the (deep) shadow region requires the use of significantly higher numerical precision, not only with increasing wavenumber k, but also increasing radius of curvature. Fig. 5. View largeDownload slide |η(θ)| plotted against the radial angle θ for circles of radii r = 1 (left) and r = 1/2 (right) and increasing values of k. Fig. 5. View largeDownload slide |η(θ)| plotted against the radial angle θ for circles of radii r = 1 (left) and r = 1/2 (right) and increasing values of k. With these observations, let us finally compare the numerical results in Fig. 4 above with those presented by Asheim & Huybrechs (2014). Just as in this figure, their implementations are confined to k = 10, 20, 40, 80, 120. Figure 11 therein concerns scattering by a circle of radius r = 1/2 and displays the pointwise relative errors for DoF = 20 (note that, since the number of wave lengths λ = 2π/k on the surface of the unit circle is twice that of the circle of radius r = 1/2, this corresponds to DoF = 40 for the unit circle). There they approximate the solution in $$[0,\theta ^{\prime }_{0}] \cup [\theta ^{\prime }_{1},2\pi ]$$ using the values of the approximate solution at the radial angles $$\theta ^{\prime }_{0} = \pi /5$$ and $$\theta ^{\prime }_{1} = 2\pi -\theta ^{\prime }_{0}$$ (these angles are exactly the vertical dotted lines on the right pane in Fig. 5 above). Evaluating η(θ) through equation (5.7) to machine accuracy for r = 1/2 reveals that $$0.14 \le |\eta (\theta ^{\prime }_{j})| \le 0.74$$ (j = 0, 1) for k = 10, 20, 40, 80, 120. For r = 1, the vertical dotted lines in the left pane in Fig. 5 as well as those in Fig. 4 correspond exactly to θ0 = 0.86 and θ1 = 2π − θ0, and evaluating η(θ) through formula (5.7) we see that, in this case, 0.14 ≤|η(θj)|≤ 0.67 (j = 0, 1) for k = 10, 20, 40, 80, 120. Table 3 $$\min |\eta (\theta )|$$ for circles of radii r = 1 and r = 1/2 and the ellipse for the wavenumbers k = 10, 20, 40, 80, 120 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 Table 3 $$\min |\eta (\theta )|$$ for circles of radii r = 1 and r = 1/2 and the ellipse for the wavenumbers k = 10, 20, 40, 80, 120 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 k = 10 k = 20 k = 40 k = 80 k = 120 Circle (r = 1) 3.2e − 2 8.4e − 3 1.6e − 3 1.6e − 4 3.3e − 5 Circle (r = 1/2) 1.7e − 1 6.4e − 2 1.7e − 2 3.2e − 3 8.9e − 4 Ellipse 4.1e − 2 1.2e − 2 2.6e − 3 3.7e − 4 1.2e − 4 In light of the discussions in Asheim & Huybrechs (2014), this explains the lack of accuracy in the pointwise relative errors on [0, θ0] ∪ [θ1, 2π] in Fig. 4 above. Let us mention in this connection that the quality of numerical approximations in Fig. 11 therein are virtually independent on frequency. This is, however, not the case in Fig. 14 therein which concerns scattering by an ellipse parametrized as $$(\frac{1}{2}\cos \theta , \sin \theta )$$ that is also rotated by π/6 in the counterclockwise direction and illuminated by a plane wave incidence with direction α = (1, 0). The circumference of this scatterer is approximately 4.97 and they use DoF = 30 (note that, with the same reasoning as above, DoF = 20 for the circle of radius r = 1/2 would correspond to about DoF = 32 for this ellipse). In the illuminated region the error decreases as k increases from 10 to 80, but there are virtually no differences between the errors for k = 80 and 120. In the shadow region, the error degrades and blows up in the deep shadow region just as in Fig. 4 above. Let us finally note that the differences between the sizes of blow up in the pointwise relative errors in the deep shadow regions are related with the minimum values of |η(θ)| depicted in Table 3. 6. Conclusions In this paper we proposed a new class of Galerkin BEMs based on novel changes of variables for the solution of two-dimensional single-scattering problems. The Galerkin approximation spaces, generated in the form of a direct sum of algebraic polynomial spaces weighted by the oscillations in the incident field of radiation, are additionally coupled with novel frequency dependent changes of variables in the transition regions. As we have shown, this construction ensures that the global approximation spaces are perfectly matched with the changes in the boundary layers of the solutions with increasing wavenumber k. Furthermore, for large frequencies, they provide significant savings over their counterparts in Ecevit & Özen (2017) in regards to the total number of DoF necessary to obtain a prescribed accuracy. Acknowledgments We thank the anonymous referees for their constructive comments which significantly improved the presentation in the paper. References Abboud , T. , Nédélec , J.-C. & Zhou , B. ( 1994 ) Méthode des équations intégrales pour les hautes fréquences . C. R. Acad. Sci. Paris Sér. I Math. , 318 , 165 -- 170 . Abboud , T. , Nédélec , J.-C. & Zhou , B. ( 1995 ) Improvements of the integral equation method for high frequency problems . In Proc. of 3rd Int. Conf. on Mathematical Aspects of Wave Propagation Phenomena, SIAM, Philadelphia, PA, USA , 178 -- 187 . Anand , A. , Boubendir , Y. , Ecevit , F. & Reitich , F. ( 2010 ) Analysis of multiple scattering iterations for high-frequency scattering problems. II. The three-dimensional scalar case . Numer. Math. , 114 , 373 -- 427 . Google Scholar CrossRef Search ADS Asheim , A. & Huybrechs , D. ( 2014 ) Extraction of uniformly accurate phase functions across smooth shadow boundaries in high frequency scattering problems . SIAM J. Appl. Math. , 74 , 454 -- 476 . Google Scholar CrossRef Search ADS Banjai , L. & Hackbusch , W. ( 2008 ) Hierarchical matrix techniques for low- and high-frequency Helmholtz problems . IMA J. Numer. Anal. , 28 , 46 -- 79 . Google Scholar CrossRef Search ADS Boffi , D . ( 2010 ) Finite element approximation of eigenvalue problems . Acta Numer ., 19 , 1 -- 120 . Google Scholar CrossRef Search ADS Boubendir , Y. , Ecevit , F. & Reitich , F. ( 2016 ) Acceleration of an iterative method for the evaluation of high-frequency multiple scattering effects. ArXiv e-prints . Bowman , J. J. , Senior , T. B. A. & Uslenghi , P. L. E. ( 1987 ) Electromagnetic and Acoustic Scattering by Simple Shapes . New York : Hemisphere Publishing Corp ., xix + 747 . Bruno , O. P. & Geuzaine , C. A. ( 2007 ) An O(1) integration scheme for three-dimensional surface scattering problems . J. Comput. Appl. Math. , 204 , 463 -- 476 . Google Scholar CrossRef Search ADS Bruno , O. P. & Kunyansky , L. A. ( 2001 ) A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications . J. Comput. Phys. , 169 , 80 -- 110 . Google Scholar CrossRef Search ADS Bruno , O. , Geuzaine , C. & Reitich , F. ( 2005 ) On the O(1) solution of multiple-scattering problems . IEEE Trans. Magn. , 41 , 1488 -- 1491 . Google Scholar CrossRef Search ADS Bruno , O. P. , Geuzaine , C. A. , Monro , J. A. Jr. & Reitich , F. ( 2004 ) Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case . Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. , 362 , 629 -- 645 . Google Scholar CrossRef Search ADS Bruno , O. P. , Domínguez , V. & Sayas , F.-J. ( 2013 ) Convergence analysis of a high-order Nyström integral-equation method for surface scattering problems . Numer. Math. , 124 , 603 -- 645 . Google Scholar CrossRef Search ADS Chandler-Wilde , S. N. , Graham , I. G. , Langdon , S. & Lindner , M. ( 2009 ) Condition number estimates for combined potential boundary integral operators in acoustic scattering . J. Integral Equations Appl. , 21 , 229 -- 279 . Google Scholar CrossRef Search ADS Chandler-Wilde , S. N. , Graham , I. G. , Langdon , S. & Spence , E. A. ( 2012 ) Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering . Acta Numer ., 21 , 89 -- 305 . Google Scholar CrossRef Search ADS Chandler-Wilde , S. N. , Hewett , D. P. , Langdon , S. & Twigger , A. ( 2015 ) A high frequency boundary element method for scattering by a class of nonconvex obstacles . Numer. Math. , 129 , 647 -- 689 . Google Scholar CrossRef Search ADS Colton , D. & Kress , R. ( 1992 ) Inverse Acoustic and Electromagnetic Scattering Theory . Applied Mathematical Sciences, vol. 93. Berlin: Springer, pp. x+305 . Corless , R. M. , Gonnet , G. H. , Hare , D. E. G. , Jeffrey , D. J. & Knuth , D. E. ( 1996 ) On the LambertW function . Adv. Comput. Math. , 5 , 329 -- 359 . Google Scholar CrossRef Search ADS Domínguez , V. , Graham , I. G. & Smyshlyaev , V. P. ( 2007 ) A hybrid numerical-asymptotic boundary integral method for high-frequency acoustic scattering . Numer. Math. , 106 , 471 -- 510 . Google Scholar CrossRef Search ADS Darrigrand , E . ( 2002 ) Coupling of fast multipole method and microlocal discretization for the 3-D Helmholtz equation . J. Comput. Phys. , 181 , 126 -- 154 . Google Scholar CrossRef Search ADS Davies , R. W. , Morgan , K. & Hassan , O. ( 2009 ) A high order hybrid finite element method applied to the solution of electromagnetic wave scattering problems in the time domain . Comput. Mech. , 44 , 321 -- 331 . Google Scholar CrossRef Search ADS Ecevit , F. & Özen , H. Ç. ( 2017 ) Frequency-adapted galerkin boundary element methods for convex scattering problems . Numer. Math. , 135 , 27 -- 71 . Google Scholar CrossRef Search ADS Ecevit , F. & Reitich , F. ( 2009 ) Analysis of multiple scattering iterations for high-frequency scattering problems. I. The two-dimensional case . Numer. Math. , 114 , 271 -- 354 . Google Scholar CrossRef Search ADS Engquist , B. & Majda , A. ( 1977 ) Absorbing boundary conditions for the numerical simulation of waves . Math. Comp. , 31 , 629 -- 651 . Google Scholar CrossRef Search ADS Eruslu , H. H. ( 2015 ) An optimal change of variables scheme for single scattering problems . Master’s Thesis . Istanbul : Boğaziçi University . Galkowski , J. ( 2014 ) Distribution of resonances in scattering by thin barriers . ArXiv e-prints . Galkowski , J. , Müller , E. H. & Spence , E. A. ( 2016 ) Wavenumber-explicit analysis for the Helmholtz h-BEM: error estimates and iteration counts for the Dirichlet problem . ArXiv e-prints . Ganesh , M. & Hawkins , S. C. ( 2011 ) A fully discrete Galerkin method for high frequency exterior acoustic scattering in three dimensions . J. Comput. Phys. , 230 , 104 -- 125 . Google Scholar CrossRef Search ADS Giladi , E . ( 2007 ) Asymptotically derived boundary elements for the Helmholtz equation in high frequencies . J. Comput. Appl. Math. , 198 , 52 -- 74 . Google Scholar CrossRef Search ADS Givoli , D . ( 2004 ) High-order local non-reflecting boundary conditions: a review . Wave Motion , 39 , 319 -- 326 . Google Scholar CrossRef Search ADS Gradshteyn , I. S. & Ryzhik , I. M. ( 2000 ) Table of Integrals, Series, and Products, 6th edn. San Diego, CA: Academic Press Inc., pp. xlvii+1163 . Grote , M. J. & Sim , I. ( 2011 ) Local nonreflecting boundary condition for time-dependent multiple scattering . J. Comput. Phys. , 230 , 3135 -- 3154 . Google Scholar CrossRef Search ADS Han , X. & Tacy , M. ( 2015 ) Sharp norm estimates of layer potentials and operators at high frequency . J. Funct. Anal. , 269 , 2890 -- 2926 . Google Scholar CrossRef Search ADS Hesthaven , J. & Warburton , T. ( 2004 ) High-order accurate methods for time-domain electromagnetics. CMES Comput. Model. Eng. Sci. , 5 , 395 -- 407 . Hewett , D. P. , Langdon , S. & Melenk , J. M. ( 2013 ) A high frequency $hp$ boundary element method for scattering by convex polygons . SIAM J. Numer. Anal. , 51 , 629 -- 653 . Google Scholar CrossRef Search ADS Hewett , D. P. , Langdon , S. & Chandler-Wilde , S. N. ( 2015 ) A frequency-independent boundary element method for scattering by two-dimensional screens and apertures . IMA J. Numer. Anal. , 35 , 1698 -- 1728 . Google Scholar CrossRef Search ADS Huybrechs , D. & Vandewalle , S. ( 2007 ) A sparse discretization for integral equation formulations of high frequency scattering problems . SIAM J. Sci. Comput. , 29 , 2305 -- 2328 . Google Scholar CrossRef Search ADS Kim , T. ( 2012 ) Asymptotic and numerical methods for high-frequency scattering problems. Ph.D. Thesis , University of Bath, UK . Lazergui , S. & Boubendir , Y. ( 2016 ) Asymptotic expansions of the Helmholtz equation solutions using approximations of the Dirichlet to Neumann operator . ArXiv e-prints . Melrose , R. B. & Taylor , M. E. ( 1985 ) Near peak scattering and the corrected Kirchhoff approximation for a convex obstacle . Adv. Math. , 55 , 242 -- 315 . Google Scholar CrossRef Search ADS Schwab , C. ( 1998 ) p- and hp-finite Element Methods: Theory and Applications in Solid and Fluid Mechanics. Numerical Mathematics and Scientific Computation. The Clarendon Press, New York: Oxford University Press, pp. xii+374 . Spence , E. A. , Chandler-Wilde , S. N. , Graham , I. G. & Smyshlyaev , V. P. ( 2011 ) A new frequency-uniform coercive boundary integral equation for acoustic scattering . Comm. Pure Appl. Math. , 64 , 1384 -- 1415 . Google Scholar CrossRef Search ADS Spence , E. A. , Kamotski , I. V. & Smyshlyaev , V. P. ( 2015 ) Coercivity of combined boundary integral equations in high-frequency scattering . Comm. Pure Appl. Math. , 68 , 1587 -- 1639 . Google Scholar CrossRef Search ADS Tong , M. S. & Chew , W. C. ( 2010 ) Multilevel fast multipole acceleration in the Nyström discretization of surface electromagnetic integral equations for composite objects . IEEE Trans. Antennas Propagation , 58 , 3411 -- 3416 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Feb 5, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations