Hybrid Spectral Difference Methods for Elliptic Equations on Exterior Domains with the Discrete Radial Absorbing Boundary Condition

The hybrid spectral difference methods (HSD) for the Laplace and Helmholtz equations in exterior domains are proposed. We consider the fictitious domain method with the absorbing boundary conditions (ABCs). The HSD method is a finite difference version of the hybridized Galerkin method, and it consists of two types of finite difference approximations; the cell finite difference and the interface finite difference. The fictitious domain is composed of two subregions; the Cartesian grid region and the boundary layer region in which the radial grid is imposed. The boundary layer region with the radial grid makes it easy to implement the discrete radial ABC. The discrete radial ABC is a discrete version of the Bayliss–Gunzburger–Turkel ABC without pertaining any radial derivatives. Numerical experiments confirming efficiency of our numerical scheme are provided.

As shown in Fig. 1 we consider a fictitious domain 1 ∪ 2 with the inner boundary 1 and outer boundary 2 for the exterior domain problem. The outer boundary 2 is a circle with radius r 0 > 1, and an absorbing boundary condition (ABC) is imposed on it, which replaces the decay condition at infinity. In 1 the Cartesian coordinate system is used, and the rectangular or trapezoidal partition is used for implementation of the hybrid difference method (see Fig. 3). The subdomain 2 is a thin layer boundary region with single or double layers, on which the polar coordinate system is used. The width of the boundary layer region is of size h, the maximum diameter of a partition of 1 . Since our ABC is constructed based on the multipole expansion of a solution the layer boundary region with the polar coordinate and radial grid can be more efficient for implementing high order ABCs.
In 1 ∪ 2 the equation (1.1) can be reformulated as follows with a proper ABC: (1.2e) Here, the unit outward normal vector ν i on 12 is emanating from i , and (r, t) is the polar coordinate for 2 . The interface 12 is a polygon, and it is designed to change according to the mesh refinement in 1 . As a mesh becomes finer 12 converges to a circle. Figure 1 represents a simple grid generation of 1 and 2 . For 1 the mesh is composed of rectangles (in the interior of 1 ) and trapezoids (adjacent to 12 ). The radial mesh for 2 should be constructed subordinately to the mesh in 1 . The advantage of considering two separate subdomains is that the inner boundary of a general shape can be treated somehow well with the Cartesian grid and the high order ABC can be more easily applied on a radial grid. In this paper a new absorbing boundary condition (that is named the discrete radial ABC) is introduced for the Laplace and Helmholtz equations. To derive the discrete ABC we follow the idea of the Bayliss-Gunzburger-Turkel (BGT) [3] absorbing boundary condition. To the author's knowledge there has not been much attention to the ABCs for frequency domain problems. However, the frequency domain approach on a parallel computing environment can provide very efficient numerical solvers. On the other hand there have been many efforts for construction of the ABCs for the time dependent wave equations. Engquist and Majda [5] made a pioneering work, in which the ABC is obtained by considering waves that enter or leave the domain and annihilating those scattered waves entering the domain from the outside. Bayliss and Turkel [2] constructed a sequence of absorbing boundary conditions for the wave equation, based on annihilation of first several negative radial power terms in the multipole expansion of a solution, and those idea were exploited later to derive the BGT condition with Gunzburger for the static elliptic and wave equations. The high order derivative term of the ABC in [2] can be impractical for implementation, and Hagstrom and Hariharan [7] developed the ABCs that contain only first order derivatives at a cost of auxiliary functions. Huan and Thompson [9] developed these ideas further to obtain the ABCs applicable in the finite element method in a symmetric way. Higdon [8] obtained absorbing boundary conditions by invoking to the finite time difference rather than the differential equation, which improved convergence quality of numerical solutions. Our discrete radial ABC is also obtained by treating the annihilation condition of the BGT condition directly in the radial finite difference equation. For absorbing boundary conditions with elliptic or general shaped obstacles we refer to [1,13,14]. For a brief review on the various ABCs we refer to [6] and references therein.
For a numerical method of partial differential equations we consider the hybrid spectral difference method (HSD). The HSD can be understood as the finite difference version of the hybridized Galerkin method. The HSD is simple to implement and it has several advantages over the conventional finite difference methods. The method can be applied to nonuniform grids, retaining the optimal order of convergence, and numerical methods with an arbitrarily high-order convergence can be obtained easily. Problems on a complicated geometry can be treated reasonably well, and the boundary condition can be imposed exactly on the exact boundary. The mass conservation property holds in each cell and flux continuity holds across inter-cell boundaries. Embedded static condensation property considerably reduces global degrees of freedom [10][11][12].
The aim of this paper is to introduce the discrete radial ABC and the hybrid spectral difference method for elliptic equations. The hybrid spectral difference method can manage reasonably well obstacles of a general shape in scattering problems. The fictitious boundary can be taken as a circle, independently of the shape of obstacles, which makes easy the implementation of the discrete radial ABC of a high order.
The paper is organized as follows. In Sect. 2 we review the Bayliss-Gunzburger-Turkel absorbing condition, and the discrete radial ABC is introduced by following the annihilation idea of them. The high order radial derivatives in the BGT condition may not be approximated efficiently in the finite difference setting when the wave number is large. With the discrete radial ABC an exact annihilation is accomplished. In Sect. 3 unique solvability and convergence analysis on an annulus region are presented, and the proofs follow closely the ideas in [3]. In Sect. 4 the HSD on the rectangular grid with the Cartesian coordinate is introduced. The HSD in the polar grid can be derived in a similar manner and we delete details. Section 5 is devoted to numerical experiments for the exterior harmonic and Helmholtz equations. Some comparisons between the BGT condition and discrete radial ABC are made. It is shown that the discrete radial ABC produces quite accurate numerical solutions even for the wave equation of a large wave number. A scattering problem with a scatterer of a general shape is tested, which justify that our new method can be an effective numerical solver for a class of scattering problems.

Absorbing Boundary Conditions on the Boundary Layer
In this section we review the Bayliss-Gunzbuger-Turkel (BGT) absorbing boundary condition, and a new absorbing boundary condition (the discrete radial ABC) is proposed. The discrete radial ABC bears the same idea of the BGT condition, and it can be more easily implemented for a high order method. Moreover, convergence looks better than that of the BGT condition in finite difference settings.

The BGT Absorbing Boundary Condition for the Laplace Equation
For the Laplace equation on an exterior domain, it has the multipole expansion, where φ k (t) ∈ span{cos kt, sin kt}. Simple calculation yields that Let us define B m := m−1 j=1 r ∂ ∂r + j . Then, it is easy to see that Therefore, Hence, we have the error estimate; where r 0 is the radius of 2 .

The BGT Absorbing Boundary Condition for the Helmholtz Eqution
For the Helmholtz equation on an exterior domain the solution has the multipole expansion such that Here, H 0 and H 1 are the Hankel functions of the first kind of orders 0 and 1, respectively. The asymptotic expansion is due to [4]. Let us introduce the operator, By the same way as the above, it can be shown that Therefore, Trivially, we have The error estimate holds;

The Discrete Radial Absorbing Boundary Condition
We consider a multipole expansion of a solution in the form: are the grid points on the ray with an angle t emanating from the origin with (τ m , t) ∈ 2 . Note that r 0 = τ m is the radius of the circle 2 . We consider the discrete radial ABC operator as follows: such that the coefficients, {c j } m j=1 satisfy  Suppose |r 0 − τ 1 | 1, then R h m u satisfy the following estimate: The boundary condition in (1.2) on 2 will be replaced by R h m u = 0. Tables 1, 2, 3 show numerical approximation properties of the BGT and discrete radial ABCs. We consider the radial grid on annulus 1 < r < r 0 (= 3). The finite difference approximation As shown in Tables 1, 2 the finite difference (FD) approximation of the BGT condition induces a lot of error for a high wave number ω since the FD approximation of high order derivatives in the BGT condition is not accurate enough. For this reason it is recommended in [3] to use the ABC condition where the constants k 1 , k 2 , k 3 depend only on ω and r 0 when u satisfies the Helmholtz equation. It is obtained by substituting u rr in On the other hand the discrete ABC (2.2) approximate the absorbing boundary condition exactly as shown in Table 3.

Solvability and Convergence Analysis
In this section we provide the unique solvability and convergence analysis for the discrete radial ABC of the exterior Laplace equation on an annulus region. The analysis closely follows those ideas in [3]. It must be possible to extend those analysis to the case of the Helmholtz equation. For simplicity of analysis we consider the fictitious domain as an annulus, 1 < r < r 0 . Proof The solution of the Laplace equation on the annulus region has the representation in the polar coordinate system as follows: where (r, t) is the polar coordinate and φ k (t) = c k cos kt + d k sin kt. We will show that the trivial boundary condition induces the trivial solution. On the boundary we have that Then, we have φ k = μ k φ k , which implies that φ k = 0 for m ≤ k < ∞ if μ k = 1. It follows immediately that φ −k = 0 for m ≤ k < ∞. As a conclusion φ k = 0, −∞ < k < ∞.

Remark 3.2 In the proof of the above theorem we have
Then, there exists k 0 > 0 such that |μ k | < 1 for all k ≥ k 0 . Therefore, the assumption in Theorem 3.1 is reasonable. Proof The exact solution u and its approximation w by the discrete radial ABC have the multipole expansions, From the Eqs. (3.1) and (3.1b) we have The above can be summarized as follows: Therefore, Simple calculation yields that Then, we have The Q * 2 grid: Hence, We arrive at the conclusion that

Hybrid Finite Difference Methods
As mentioned in Sect. 1 the domain = 1 ∪ 2 is a multiply connected domain with the inner boundary 1 , the outer boundary 2 and the interface 12 . To define the hybrid difference method (HDM) we need a quasi-rectangular partition of the subdomain 1 . Here, the quasi-rectangle includes a rectangle, a trapezoid (Fig. 3) or a trapezoid-like mesh with one rounded side (Fig. 4). For the subdomain 2 we consider a radial subdivision by polar rectangles. In this section we describe the HDM in the Cartesian coordinate system because the hybrid difference in the polar coordinate can be obtained analogously. The hybrid spectral difference method (HSD) is a high order version of the hybrid difference method.
Let K h denote the skeleton of a mesh generation T h of , and let N ( ) and N (K h ) denote the set of grid points in the closure of a domain and that on its skeleton, respectively. It is worth to note that the mesh is composed mostly of rectangles. However, we need trapezoidal meshes to match the interface 12 and we may need rounded trapezoidal meshes to match 1 if 1 is curved. We call the grid configuration in Fig. 2 as the Q * 2 grid and the grid configuration in Fig. 3 as the Q * 4 grid. The Q * m grid denotes the grid configuration obtained from the Q m grid by deleting four vertices, where the Q m grid represents the grid on which the binomials of degree up to m (Q m ) can be uniquely defined. In the HDM the nodes are the Gauss-Legendre points in a cell and its edges.
Let us begin with the simplest hybrid difference method (HDM) in the Cartesian coordinate system. The derivation is based on the grid configuration in Fig. 2 (4.1) on the cell R 1 , and the interface finite difference approximates the flux continuity at the intercell point so that Here, ∂ h ν u = ∇ h u · ν, and ν is the outward unit normal vector from each cell. The cell finite difference (4.1) can be solved for u(η 4 ) in terms of the cell boundary values, {u(η 1 ), u(η 3 ), u(η 5 ), u(η 8 )} for each cell. Then, the interface difference (4.2) yields the global stiffness system with unknowns {u(η) : η ∈ N (K h )}. Therefore, the static condensation property is naturally embedded in the HDM.
To derive a higher order method we consider a more complicated cell configuration as in Fig. 3. Then, the cell finite difference solves the cell problem: and the interface finite difference solves the interface condition For a detailed derivation of the high order finite difference formulas for h and ∂ h ν we refer to [11,12]. Let us define the Gaussian quadratures on a reference cell R (with |R| = h × k) Table 4 Reduction in degrees of freedom by static condensation and its boundary ∂ R on the Q * m grid as follows: Here, {σ j } m−1 j=1 is the Gaussian weights on the reference interval [−1, 1]. There exists a natural composite quadrature defined on the whole domain as follows.
Then, the Eqs. (4.3) and (4.4) can be unified in a variational form: Some numerical analysis of (4.5) can be found in [11], where analysis is performed for the Poisson equation. The hybrid difference method is understood as the finite difference version of the hybridized finite element method. The main advantage of the HDM is that the finite differences are of one dimensional nature (it is not related to two or three dimensional polynomial interpolation). Therefore, we can handle the boundary data exactly even on a domain with a curved boundary by extending the line of derivative evaluation up to the given curved boundary and taking the intersection as a nodal point (see Fig. 4). Secondly, the method possesses the embedded static condensation property, which considerably reduces global degrees of freedom in the global stiffness system.

Remark 4.1
The subdomain 2 can contain one or two layers of polar (quasi) rectangles. To combine a low order HDM and a high order ABC we need two layers, and just one layer is enough for combination of the Q * m grid HDM and the R h n -ABC if n ≤ m + 1.

Remark 4.2
Suppose a domain is the unit square with a N × N -rectangular partition. Table 4 highlights reduction in degrees of freedom by the static condensation. The reduction rate is computed asymptotically under the assumption that N is a big number.

Numerical Experiments
In this section we provide numerical experiments on three regions: • An annulus with 1 < r < 3 2 , • A Chinese coin (see Fig. 1  on the annulus. Here, H 1 is the Hankel function of the first kind with order one. Only the results for the discrete radial ABC (R h m (u) = 0) are presented. This example is to illuminate the sensitivity of the discrete radial ABC depending on change of the wave number ω.

Example 3
The scattering problem with the sound soft boundary condition is solved on the Chinese coin domain. Here, the grids are give as in Fig. 1. The incident wave u in = e −iω(x cos α+y sin α) with the incident angle α = 0 and ω = 5 is tested, and the scattered wave is plotted by increasing the degrees of freedom. The 4th order discrete radial ABC (R h 4 u = 0) on the Q * 4 grid is used for numerical computation.
Example 4 The scatterer of the shape "D" is tested with the incident wave u in = e −iω(x cos α+y sin α) and the incident angles α = 0, π 4 , π 2 , π. The same computational parameters are chosen as in Example 3.
As shown in Fig. 5 the discrete radial and BGT ABCs perform evenly well for the Laplace equation. The higher order ABC shows better convergence than the lower order ABC. As the mesh becomes finer the error seems to be dominated by the truncation error in the absorbing boundary condition. For the Helmholtz equation numerical examples only with the discrete radial ABC are presented. We observe that the BGT absorbing boundary condition produces spurious numerical solutions because of the poor finite difference approximation property as mentioned in Sect. 2, and we omit it here. As shown in Figs. 6 and 7 the discrete radial ABC performs quite robust even for a larger wave number (in view of the level of discretization). Here, the computation is performed on a N × 2N mesh by doubling N up to N = 64 for the radial Q * 2 grid and up to N = 32 for the radial Q * 4 grid. When the wave number is The real part of the scattered waves with the various incident angles, α = 0, π 4 , π 2 , π. The mesh corresponds to the 128 subdivision of the exterior boundary 2 (Example 4) small the higher order ABC performs clearly better than the lower order ABC. However, convergence property becomes similar regardless of orders of the discrete radial ABC when the wave number increases. Even the 2nd order method R h 2 u = 0 performs the best for the cases, ω = 50, 100. Nonetheless, the high order hybrid difference combined with the high order discrete radial ABC (for example, the Q * 4 grid method and R h 4 = 0) produces more reliable solutions for a wide range of wave numbers. From those observation the scattering problems in Examples 3 and 4 are solved by using the ABC, R h 4 (u) = 0 on the Q * 4 grid. Figure 8 clearly shows how the scattered wave converges as the degrees of freedom (dof) increase. Figure 9 represents the real part of the scattered wave with various incident angles, and this example shows that our method can manage a complicated geometry-problem quite well. The numerical results accord well with those existing numerical results.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.