High order surface radiation conditions for time-harmonic waves in exterior domains

TL;DR

This paper introduces a new family of high-order on-surface radiation conditions for accurately approximating outgoing solutions to the Helmholtz equation in exterior domains, enhancing computational methods for wave problems.

Contribution

It develops a systematic procedure to implement high-order pseudo-differential symbols for surface radiation conditions, improving solution accuracy in exterior Helmholtz problems.

Findings

Numerical results demonstrate improved accuracy of the proposed method.

The approach effectively handles Dirichlet and Neumann boundary conditions.

Potential extensions suggest broader applicability in wave simulations.

Abstract

We formulate a new family of high order on-surface radiation conditions to approximate the outgoing solution to the Helmholtz equation in exterior domains. Motivated by the pseudo-differential expansion of the Dirichlet-to-Neumann operator developed by Antoine et al. (J. Math. Anal. Appl. 229:184-211, 1999), we design a systematic procedure to apply pseudo-differential symbols of arbitrarily high order. Numerical results are presented to illustrate the performance of the proposed method for solving both the Dirichlet and the Neumann boundary value problems. Possible improvements and extensions are also discussed.

Equations

$$\Delta u+k^{2}u=0\text{in Ω},$$

$$u=f_{\mathrm{Dir}}\text{or}{\partial }_{\nu }u=f_{\mathrm{Neu}}\text{on Γ},$$

$$\underset{r\to \infty }{\lim }r\left({\partial }_{r}u-iku\right)=0,$$

$$u(\text{x})={\int }_{\Gamma }\left[u(\text{y}){\partial }_{\nu (\text{y})}\Phi (\text{x},\text{y})-{\partial }_{\nu }u(\text{y})\Phi (\text{x},\text{y})\right]dS(\text{y}),\text{x}\in \Omega ,$$

$$Lu=\Delta u-c^{-2}{\partial }_t^{2}u=0,$$

$$\mathcal{L}u={\partial }_r^{2}u+2\mathcal{H}{\partial }_{r}u+{\Delta }_{\Gamma }u+k^{2}u.$$

$$\mathcal{L}=({\partial }_r-{\Lambda }^{-})({\partial }_r-{\Lambda }^{+})u.$$

$${\lambda }^{\pm }\sim \sum _{n=-1}^{+\infty }{\lambda }_{-n}^{\pm }$$

$$\mathcal{L}u={\partial }_r^{2}u-\left({\Lambda }^{+}+{\Lambda }^{-}\right){\partial }_{r}u+\left({\Lambda }^{-}{\Lambda }^{+}-\{{\partial }_r{\lambda }^{+}\}\right)u.$$

$${\Lambda }^{+}+{\Lambda }^{-}=-2\mathcal{H},$$

$${\Lambda }^{-}{\Lambda }^{+}-\{{\partial }_r{\lambda }^{+}\}=k^2+{\Delta }_{\Gamma }.$$

$$\sum _{n=-1}^{+\infty }\left({\lambda }_{-n}^{+}+{\lambda }_{-n}^{-}\right)=-2\mathcal{H},$$

$$\sum _{n=-2}^{+\infty }\sum _{j=-1}^{n+1}{\lambda }_{-j}^{-}{\lambda }_{j-n}^{+}-\sum _{m=-1}^{+\infty }{\partial }_r{\lambda }_{-m}^{+}=k^2+{\Delta }_{\Gamma }.$$

$${\lambda }_{-n}^{+}+{\lambda }_{-n}^{-}=0n\ne 0,$$

$${\lambda }_0^{+}+{\lambda }_0^{-}=-2\mathcal{H}.$$

$${\lambda }_1^{-}{\lambda }_1^{+}=k^2,$$

$${\lambda }_1^{-}{\lambda }_0^{+}+{\lambda }_0^{-}{\lambda }_1^{+}-{\partial }_r{\lambda }_1^{+}=0.$$

$${\lambda }_1^{\pm }=\pm ik,$$

$${\lambda }_0^{\pm }=-\mathcal{H}.$$

$${\lambda }_1^{-}{\lambda }_{-1}^{+}+{\lambda }_0^{-}{\lambda }_0^{+}+{\lambda }_{-1}^{-}{\lambda }_1^{+}-{\partial }_r{\lambda }_0^{+}={\Delta }_{\Gamma }$$

$${\lambda }_{-1}^{\pm }=\pm \frac{i}{2k}\left({\Delta }_{\Gamma }-{\mathcal{H}}^2-{\partial }_r\mathcal{H}\right)=\pm \frac{i}{2k}\left({\Delta }_{\Gamma }+{\mathcal{H}}^2-\mathcal{K}\right)$$

$${\lambda }_{-n-1}^{+}=\frac{i}{2k}\left[{\partial }_r{\lambda }_{-n}^{+}+\sum _{j=1}^{n-1}{\lambda }_{-j}^{+}{\lambda }_{j-n}^{+}\right],n\ge 1.$$

$${\lambda }_{-n}^{+}(r)={\left(\frac{r_{\mathrm{o}}}{r}\right)}^{n+1}{\lambda }_{-n}^{+}(r_{\mathrm{o}})n=-1,0,1,2,...$$

$$\mathcal{H}(r)=\frac{r_{\mathrm{o}}}{r}\mathcal{H}(r_{\mathrm{o}})\text{and}\mathcal{K}(r)={\left(\frac{r_{\mathrm{o}}}{r}\right)}^2\mathcal{K}(r_{\mathrm{o}})$$

$${\lambda }_{-(n+1)}^{+}(r)=\frac{i}{2k}{\left(\frac{r_{\mathrm{o}}}{r}\right)}^{n+2}\left[-\frac{(n+1)}{r_{\mathrm{o}}}{\lambda }_{-n}^{+}(r_{\mathrm{o}})+\sum _{j=1}^{n-1}{\lambda }_{-j}^{+}(r_{\mathrm{o}}){\lambda }_{j-n}^{+}(r_{\mathrm{o}})\right],n=1,2,...$$

$${\lambda }_1^{+}$$

$${\lambda }_0^{+}$$

$${\lambda }_{-1}^{+}$$

$${\lambda }_{-(n+1)}^{+}$$

$$[{\lambda }^{+}{]}_N=\sum _{n=-1}^{N-1}{\lambda }_{-n}^{+}N\ge 0,$$

$$[{\lambda }^{+}{]}_0$$

$$[{\lambda }^{+}{]}_1$$

$$[{\lambda }^{+}{]}_2$$

$$[{\lambda }^{+}{]}_3$$

$$[{\lambda }^{+}{]}_4$$

$$({\Delta }_{\Gamma }u)({\text{x}}_j)=\frac{2}{{\pi }^2}\sum _i\frac{{\left(\sin {\alpha }_i+\sin {\beta }_i\right)}^2}{{\left(A_j{A}_{j(i)}\right)}^{1/2}}\left(u({\text{x}}_{j(i)})-u({\text{x}}_j)\right),$$

$$\text{H}({\text{x}}_j)=-\frac{1}{{\pi }^2}\sum _i\frac{{\left(\sin {\alpha }_i+\sin {\beta }_i\right)}^2}{{\left(A_j{A}_{j(i)}\right)}^{1/2}}\left({\text{x}}_{j(i)}-{\text{x}}_j\right).$$

Full text

††journal:  

url]https://sites.google.com/site/acostasebastian01

1 Introduction

This paper is concerned with the approximation of the outgoing solution to the Helmholtz equation in a domain $\Omega$ exterior to a closed surface $\Gamma$. We consider both the Dirichlet and the Neumann boundary value problems,

[TABLE]

where $r=|\text{\bf x}|$ and $\text{\bf x}\in\mathbb{R}^{3}$. The main underlying assumptions are that $k>0$ is constant in $\Omega$, that $\Gamma$ is sufficiently smooth, and that the domain enclosed by $\Gamma$ is strictly convex. These assumptions imply that the wave field in $\Omega$ is purely outgoing, ruling out the possibility of waves reflecting back and forth between disjoint subsets of $\Gamma$.

When exact analytical solutions to the boundary value problem (1) are not available, boundary integral equations and their discretizations offer some of the most robust choices to approximate the solution numerically [1, 2, 3, 4]. Another approach is to truncate and discretize the domain surrounding $\Gamma$ and impose an absorbing condition on an artificial boundary to approximate the outgoing nature of the solution. Some reviews of these techniques are found in [5, 6, 7, 8, 9]. A third line of work is the development of on–surface radiation conditions (OSRC). These conditions can be seen as the application of absorbing boundary conditions directly on the surface $\Gamma$ in order to approximate the Neumann data when the Dirichlet data is prescribed, or vice versa. Once the Cauchy data $(u,\partial_{\nu}u)$ on $\Gamma$ is suitably approximated, then the Green’s representation

[TABLE]

provides an approximate solution in $\Omega$. Here $\Phi$ is the fundamental solution to the Helmholtz equation. The OSRC is a semi–analytical approach meant to sacrifice some accuracy in favor of tremendous computational speed. The OSRC may also render ultra fast methods to explore parameter spaces for optimization problems, for the construction of reduced bases, or inexpensive preconditioners to solve boundary integral equations numerically via Krylov iterative techniques [10, 11, 12]. In general, the OSRC provides a crude approximation of the solution. The state–of–the–art lacks error estimates, and the only notion of convergence is obtained in a microlocal sense. Therefore, the OSRC should not be employed if a high degree of accuracy is needed.

The original conception of the OSRC is due to Kriegsmann et al. [13]. Over the years, improvements and extensions have been devised, including the works of Jones [14, 15, 16], Ammari [17, 18], Calvo et al. [19, 20], Antoine, Barucq et al. [21, 22, 23, 24, 25, 26, 27, 28] and Darbas et al. [29, 10, 24, 11, 30, 12, 31, 32]. See also [33, 34, 35, 36, 37, 38, 39, 40, 28]. An overview of the OSRC theory, with an extensive list of references, was published in [22]. One of the open problems discussed by Antoine therein was the lack of a systematic procedure to increase the order of the OSRC beyond second order. Our present work is a novel attempt to partially solve this problem. Other efforts have been carried out in the past. Teymur extended the formulation of Jones, which assumes that the wave field has a well defined phase front conforming to the surface $\Gamma$, to furnish OSRC of high order [37]. Explicit expressions were provided up to fourth order. Calvo [20, 19] started from the eigenfunction representation of the Dirichlet–to–Neumann map for a circle and used rational approximations for the logarithmic derivative of the Hankel functions to obtained applicable OSRC of high order.

As reviewed in Section 2, the approach developed in the present paper is primarily influenced by [21] which is based on rigorous microlocal analysis to express the outgoing Dirichlet–to–Neumann (DtN) map in a classical pseudo–differential sense. This pseudo–differential calculus renders a recursive formula to obtain a differential expansion of the DtN map to an arbitrarily high order. In Section 3, we show that the recursive formula becomes algorithmically tractable for the spherical geometry which we formally generalize to other arbitrary convex surfaces. In Section 4, we propose a numerical implementation of the high order OSRC based on a triangulation for the surface $\Gamma$. A few numerical results and comparisons with exact solutions are shown in Sections 5 and 6. Limitations, possible improvements and conclusions are discussed in Section 7.

2 Approximation of the outgoing Dirichlet–to–Neumann map

We follow closely the approach described in [21], but we only consider the pseudo–differential analysis with respect to the time variable $t$ and its Fourier dual $k$. We start with the wave equation,

[TABLE]

where the space derivatives can be written in terms of a local, tangential system of coordinates on the surface $\Gamma$ to obtain the symbol (with respect to time) of the wave operator,

[TABLE]

Here $\mathcal{H}$ is the mean curvature of $\Gamma$, and $\partial_{r}$ represents the derivative in the outward normal direction on $\Gamma$. Multiplication by $k^{2}$ is an operator of order $+2$. Here $\Delta_{\Gamma}$ represents the Laplace–Beltrami operator on the surface $\Gamma$ which we consider to be of order [math] with respect to $k$. A concise review of the differential geometry of surfaces, including the definition of curvatures and the Laplace–Beltrami operator, is found in Section 2 of [21]. We seek to decompose the wave field propagating across a surface into incoming and outgoing components. This can be done for convex surfaces. Based on Nirenberg’s factorization theorem [41, 21], there are two pseudo–differential operators $\Lambda^{-}$ and $\Lambda^{+}$ of order $+1$, such that

[TABLE]

We refer to $\Lambda^{+}$ and $\Lambda^{-}$ as the outgoing and incoming Dirichlet–to–Neumann (DtN) operators for the radiating boundary value problem defined in the exterior of $\Gamma$. The operators $\Lambda^{\pm}$ have respective symbols $\lambda^{\pm}$ that admit the following expansion

[TABLE]

where $\lambda^{\pm}_{-n}$ are symbols of degree $-n$ with respect to $k$. The expansion is understood in the sense of classical pseudo–differential operators. See [21, 42, 43, 44] for details. We develop (5) to obtain

[TABLE]

By comparing (4) and (7), and matching terms with same degree of radial derivatives, we find that,

[TABLE]

Now we plug the expansion (6) to write the system (8)–(9) in terms of pseudo–differential symbols. We obtain,

[TABLE]

Now in (10) we match symbols of the same order with respect to $k$ to find that

[TABLE]

Similarly, matching the symbols of order $+2$ and $+1$ in (11), we find that

[TABLE]

Since $\lambda_{1}^{-}=-\lambda_{1}^{+}$ and $\lambda_{0}^{-}=-2\mathcal{H}-\lambda_{0}^{+}$ then we arrive at

[TABLE]

Now we proceed to match the symbols of order [math] from (11) to obtain

[TABLE]

Using (14)–(15) and the fact that $\partial_{r}\mathcal{H}=\mathcal{K}-2\mathcal{H}^{2}$ (see [45, Probl. 11 §3.5]) we arrive at

[TABLE]

where $\mathcal{K}$ is the Gauss curvature of $\Gamma$. Using (11) we leave the rest of the symbols expressed recursively as follows,

[TABLE]

3 Formal Generalization from Spherical Case

The only impediment to algorithmically apply (17) is the lack of knowledge of the dependence of $\lambda_{-n}^{+}=\lambda_{-n}^{+}(r)$ in order to evaluate $\partial_{r}\lambda_{-n}^{+}$. Here we seek to provide an explicit expression for the symbols $\lambda^{+}_{-n}$ as a function of $r$ for the case where $\Gamma$ is a sphere of radius $r_{\rm o}$. We take a family of parallel spherical surfaces $\Gamma(r)$ parametrized by their radius $r$. Our practical application of the recursive formula (17) is based on the following lemma.

Lemma 1**.**

For the family of spherical surfaces $\Gamma(r)$ of radius $r$, the symbols $\lambda_{-n}^{+}$ have the following dependence on the parameter $r$,

[TABLE]

where $\lambda_{-n}^{+}(r_{\rm o})$ is independent of $r$.

Proof.

The claim is clearly true for $n=-1$ even if $\Gamma$ is not a sphere. Now for the spherical case, the mean and Gauss curvatures have the following properties,

[TABLE]

where $\mathcal{H}(r_{\rm o})=1/r_{\rm o}$ (mean curvature) and $\mathcal{K}(r_{\rm o})=1/r_{\rm o}^{2}$ (Gauss curvature) for a sphere. Therefore, the claim is true for $n=0$. For a sphere, the Laplace–Beltrami operator satisfies $\Delta_{\Gamma(r)}=(\frac{r_{\rm o}}{r})^{2}\Delta_{\Gamma}$. Hence, the claim is also valid for $n=1$. Now we proceed by induction using the recursive formula (17). Directly by inductive hypothesis, we have

[TABLE]

∎

Using Lemma 1, we can express $\partial_{r}\lambda_{-n}^{+}$ explicitly in order to apply DtN operator $\lambda^{+}$ using the series (6). We obtain the following representation for the symbols on the surface $\Gamma$,

[TABLE]

where we have identified $r=\mathcal{H}^{-1}$. Although the recursive formula (19d) is only valid for the spherical case, we formally generalize it to other geometries. In order to conform to the literature, we will use the following terminology in the sequel. Our approximation $[\lambda^{+}]_{N}$, of order $N$ with respect to $k$, for the symbol $\lambda^{+}$ of the outgoing DtN operator $\Lambda^{+}$ is given by

[TABLE]

where $\lambda_{-n}^{+}$ are defined by (19). At least for the spherical case, we get $|\lambda^{+}-[\lambda^{+}]_{N}|\sim\mathcal{O}(k^{-N})$ as $k\to\infty$, which reveals the order of the approximation in a microlocal sense. Notice that $[\lambda^{+}]_{N}$ is not expected to converge as $N\to\infty$ for fixed $k$. However, as illustrated in Section 5, we do see more accurate results as $N$ increases within a reasonable range. We refer to $[\lambda^{+}]_{N}$ as the Dirichet–to–Neumann On–Surface Radiation Condition of order $N$ (DtN–OSRCN). Analogously, we refer to $[\lambda^{+}]_{N}^{-1}$ as the Neumann–to–Dirichlet On–Surface Radiation Condition of order $N$ (NtD–OSRCN).

In order to easily compare with other results in the literature of OSRC, we display explicitly the DtN-OSRCN up to the fourth order, that is, for $N=0,1,2,3,4$:

[TABLE]

4 Discrete Implementation

In this section, we describe a numerical implementation of the proposed on–surface radiation conditions defined by (20). We numerically solve both the Dirichlet and the Neumann boundary value problems (1), using a discrete version of the DtN–OSRCN and of NtD–OSRCN, respectively. The approach is based on having a triangulation that approximates the surface $\Gamma$. In the following subsections, we describe the discrete approximations of the Laplace–Beltrami operator, of the mean and Gauss curvatures of $\Gamma$, and of the OSRC.

4.1 Discrete Laplace–Beltrami operator and curvatures

The subject of approximating geometrical properties and operators on triangulated surfaces has received considerable attention in the past decades. See [46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58] for some recent advances in this area which we use as our guiding references.

The discrete Laplace–Beltrami operator is described as follows. Let $\{\text{\bf x}_{j}\}_{j=1}^{J}$ be the collection of vertices of the triangulation $\mathcal{T}$ of the surface $\Gamma$. For a fixed vertex $\text{\bf x}_{j}$, let $\{\text{\bf x}_{j(i)}\}_{i=1,2,...}$ be the neighboring vertices of $\text{\bf x}_{j}$, and let $\{e_{i}\}_{i=1,2,...}$ be the edges of the triangulation that connect $\text{\bf x}_{j}$ and $\text{\bf x}_{j(i)}$. Now for each edge $e_{i}$, let $\alpha_{i}$ and $\beta_{i}$ be angles at the vertex $\text{\bf x}_{j}$ of the triangles on either side of the edge $e_{i}$. For a smooth function $u$ defined on the surface $\Gamma$, we use the following discrete Laplace–Beltrami operator,

[TABLE]

where $A_{j}$ is the area associated with the vertex $\text{\bf x}_{j}$ defined as $1/3$ of the total area of the triangular elements sharing the vertex $\text{\bf x}_{j}$.

For the discrete mean curvature, we first define the discrete mean curvature vector as

[TABLE]

Then the mean curvature is given by

[TABLE]

where $\textbf{n}(\text{\bf x}_{j})$ is the normal vector of the vertex $\text{\bf x}_{j}$.

For the discrete Gauss curvature at each vertex $\text{\bf x}_{j}$ of the triangulation, we follow the approach described in [52, 54, 48] based on the Gauss–Bonnet theorem. Let $\alpha_{i}$ be the angle between two successive edges sharing the vertex $\text{\bf x}_{j}$. Then, the approximated Gauss curvature is given by

[TABLE]

where $A_{j}$ is the area associated with the vertex $\text{\bf x}_{j}$ as defined above.

4.2 Discrete on–surface radiation conditions

Given the triangulation $\mathcal{T}$ of the surface $\Gamma$, we have all the necessary definitions in Subsection 4.1 to implement a discrete version of the formulas (19) to apply the discrete DtN–OSRCN given by (20). On the triangulation $\mathcal{T}$, the operator $[\lambda^{+}]_{N}$ can be represented by a matrix of size $M\times M$, where $M$ is the number of vertices of the triangulation $\mathcal{T}$. This is equivalent to working on a piecewise linear function space supported by a mesh of triangular elements. We are not concerned with increasing the polynomial order of approximation because the error from the OSRC is known to dominate the calculations once the mesh is fine enough to resolve the wavelength $2\pi/k$.

The recursive formula (19d) and the fact that (19c) contains the Laplace–Beltrami operator $\Delta_{\Gamma}$ imply that the symbols $\lambda_{-n}^{+}$ apply differential operators of increasingly higher order as $n$ increases. A direct application of the operators $\lambda_{-n}^{+}$ quickly becomes ill–conditioned (unstable) at the numerical level. Therefore, we propose a regularized application of the discrete DtN–OSRCN and NtD–OSRCN operators to obtain $u\approx\lambda^{+}f_{\rm Dir}$ or $v\approx(\lambda^{+})^{-1}f_{\rm Neu}$. This application is defined in precise term by the Algorithm 1.

In this paper, we use the following regularizer,

[TABLE]

where $\text{\bf x}_{j(i)}$, $\alpha_{i}$, $\beta_{j}$, $A_{i}$ and $A_{j(i)}$ are defined in Subsection 4.1. The constant $C_{j}$ normalizes $\mathcal{A}$ so that the entries in each row sum to $1$. In other words, $\mathcal{A}$ is a local averaging operator mimicking the off–diagonal entries of the Laplace–Beltrami operator (22). Other approaches, such as global smoothing operators, may work too. For instance, good results are obtained by using,

[TABLE]

where $c_{n}\sim\left(n/k\right)^{2}$. In practice, the application of (27) requires the inversion of an $M\times M$ matrix representing the discretization of a coercive operator on $\Gamma$. In this paper we show results using the application of (26), which is computationally lighter than (27). The regularizers (26) and (27) damp out the spurious oscillations that may appear due to non–smooth imperfections in approximating the surface $\Gamma$, the curvatures $\mathcal{H}$ and $\mathcal{K}$, as well as the operator $\Delta_{\Gamma}$. It turns out that approximation of geometrical properties on triangulations is slowly convergent and quite sensitive to the valence of vertices and to the quality of triangles [48, 50, 49, 53, 54, 55].

5 Numerical Results

In this section we present a few numerical results obtained from the implementation of the proposed high order OSRC described in Section 4. We discuss three realizations for the surface $\Gamma$: the unit sphere, a marshmallow–like surface and a couple of more challenging non–convex surfaces in the likeness of a squash and a red blood cell. These four shapes are shown in Figure 1. The approximate mean curvature, defined by (24), is shown in Figure 2. Although the pseudo–differential approximation of the DtN operator described in Section 2 is valid for convex surfaces, we still apply it to non–convex surfaces to explore its performance.

For comparison, we want to have an exact solution to the boundary values problems (1), regardless of the geometry of the surface $\Gamma$. Hence, we define very specific boundary data in order to engineer an exact solution. First consider the field,

[TABLE]

where $\Phi(\text{\bf x})=e^{ik|\text{\bf x}|}/(4\pi|\text{\bf x}|)$ is the outgoing fundamental solution to the Helmholtz equation. Hence, $F$ represents the superposition of $L$ point–sources with respective locations at $\text{\bf c}_{l}$. If the points $\text{\bf c}_{l}$ for $l=1,2,...,L$ are located inside the surface $\Gamma$, then the exact solution to the boundary values problems (1) with Dirichlet data $f_{\rm Dir}=F|_{\Gamma}$ or Neumann data $f_{\rm Neu}=\partial_{\nu}F$ is given by $F|_{\Omega}$. Figure 3 displays the $L^{2}$–norm relative errors for both DtN–OSRCN and NtD–OSRCN as $N$ increases, for $k=10\pi$, and using (26) as the regularizer. In order to create an interesting (non–trivial) problem, we used $L=5$ point sources to define (28). These 5 point sources are located on the corners of a regular pentagon of circumradius equal to $1/2$ and lying on the $\{z=0\}$ plane. The errors are defined as follows,

[TABLE]

We observe in Figure 3 that as the order $N$ of the OSRCN increases, the error decreases. However, at some point the error becomes stagnant. This is not the result of the finite–dimensional approximation obtained by the discrete triangulation. A refinement of the mesh does not diminish the error any further. Nevertheless, we are pleased to see that the higher order OSRC diminishes the error in a systematic manner. In the spherical case, we are able to approximate the solution with an error below $1\%$. For the other three surfaces, the error drops below $2\%$, $4\%$ and $10\%$, respectively.

As defined by the Algorithm 1, the application of the NtD–OSRC requires the inversion of the operators $[\lambda_{+}]_{N}$. However, this inversion is not a heavy computational burden because this auxiliary problem is two–dimensional and local in nature. So its discretization leads to medium–size sparse matrices. In our experience, the memory demand is low and Gaussian elimination works very well. Better efficiency can be achieved by using permutation schemes to pre–order the matrices for Gaussian elimination. The sparsity patterns of the matrices representing the discretization of $[\lambda^{+}]_{N}$ for $N=2,4,6$ are displayed in Figure 4. The top three panels display the raw patterns induced by the connectivity of the triangular mesh. Due to the nature of the recursive formula (19), the pattern of $[\lambda^{+}]_{2N+1}$ is the same as the pattern of $[\lambda^{+}]_{2N}$. The appearance of the Laplace–Beltrami operator $\Delta_{\Gamma}$ in (19c) causes the matrix to lose sparsity as $N$ increases. Therefore, we computed the reverse Cuthill–McKee permutation (built in MATLAB) of the matrix representing the Laplace–Beltrami operator $\Delta_{\Gamma}$, and applied the same permutation to $[\lambda^{+}]_{N}$. The sparsity patterns for the re–ordered matrices are shown in the bottom panels of Figure 4. Notice the considerable improvement in reducing the bandwidth of these matrices. Significant improvements in computational efficiency can be expected when using Gaussian elimination to apply the NtD–OSRCN for fine meshes.

6 Far–field Pattern

Now we illustrate the ability of the proposed high order OSRC to approximate the far–field pattern radiated from the surface $\Gamma$. The far–field pattern $u_{\infty}$ of the solution to the boundary value problem (1) can be obtained from the Green’s formula (2) and the following asymptotic behavior of the fundamental solution

[TABLE]

We get

[TABLE]

For the Dirichlet problem, $u$ is prescribed and $\partial_{\nu}u$ is approximated using the DtN–OSRCN operator. For the Neumann problem, $\partial_{\nu}u$ is prescribed and $u$ approximated using the NtD–OSRCN operator. Then the computation of the far–field pattern reduces to numerical quadrature to approximate the integral in (31). The behavior of the error in approximating the far–field pattern with the proposed method is displayed in Figure 5 for increasing order $N$. This relative error is measured in the $L^{2}$–norm on the unit sphere (domain of the far–field pattern) for both the DtN–OSRCN and NtD–OSRCN methods for the same surfaces and parameters discussed in Section 5 and Figure 3.

Figures 6 and 7 display the far–field patterns for $k=10\pi$ and $k=20\pi$, respectively, for the squash–like surface. A comparison is shown between the exact and approximate solution using the DtN–OSRC8. The Dirichlet data is prescribed as in Section 5. The presence of the sources inside $\Gamma$ produces an interesting interference pattern that the OSRC is able to capture quite accurately.

7 Conclusions

We have constructed an on–surface radiation condition (OSRCN) of arbitrarily high order $N$. This condition is applicable on smooth convex surfaces of arbitrary shape. The high order is obtained through a recursive formula for higher order pseudo–differential symbols of the outgoing Dirichlet–to–Neumann (DtN) map. From the numerical experiments, we observe a systematic reduction of the error as the order $N$ increases which is the main achievement of this paper. This is a rather impressive property of the DtN–OSRCN, considering that we are applying local, two–dimensional, explicit operations to solve a three–dimensional elliptic boundary value problem.

The major limitation of our proposed high order OSRC is caused by the formal generalization of the recursive formula (19d) from the spherical case to other geometries. See Section 3 for details. The symbols (19a)–(19c) are valid for arbitrary convex geometry. But the recursion (19d) is not. However, we suspect that some of the correct geometric information from (19c) gets propagated into the higher symbols through (19d). The other limitations are similar to those of previous formulations of OSRC. In general, the OSRC provides a relatively crude approximation of the solution lacking a–priori error estimates. The microlocal analysis employed to define the geometric OSRC does not imply that the DtN–OSRCN converges to the exact DtN operator as $N$ increases for fixed $k$. The only notion of convergence is obtained as $k\to\infty$. And even in that case, the convergence in the microlocal sense allows for a mismatch in the form of a smoothing operator. This may explain the stagnation of the error observed in Figure 3.

Finally, we propose several directions of research that may improve or extend our work: (i) A generalization of Lemma 1 to ellipsoidal geometry may provide more accurate results including the explicit appearance of the principal curvatures in the recursive formula. (ii) We expect the general ideas of this paper to apply to electromagnetic waves; an area where many engineering application reside. (iii) In this paper, we are not particularly concerned with the error associated with the triangulation of the surface $\Gamma$ and the discrete approximation of geometric properties. Variational and finite element formulations compatible with geometric structures may offer great improvements [59]. (iv) Following the approaches from [34, 28], problems with multiple disjoint obstacles may also be addressed.