Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations

Journal of Computational and Applied Mathematics (

)

–

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations Tiffany Jones *, Qin Sheng Department of Mathematics and Center for Astrophysics, Space Physics and Engineering Research, Baylor University, One Bear Place, Waco, TX 76798-7328, USA

article

info

Article history: Received 23 November 2017 Received in revised form 4 March 2018 MSC: 65M06 65M12 65M50 65Z05 78A15 78M20 Keywords: Wave equation High oscillations Compact schemes Decomposition Asymptotic stability Self-focusing phenomenon

a b s t r a c t This paper concerns numerical stabilities of a decomposed compact finite difference method for solving Helmholtz partial differential equation problems with large wave numbers. Radially symmetric transverse fields and standard polar coordinates are considered. A decomposition is implemented to remove the anticipated singularity in the transverse direction. Compact scheme structures are introduced in the transverse direction to raise the accuracy for laser and subwavelength optical computations. It is proven that, while the highly accurate compact algorithm shies away from the conventional stability in the von Neumann sense, it is asymptotically stable with index one. Additionally, a discussion on theoretical and numerical stabilities for this method is presented. Numerical experiments further demonstrate the high reliability of the scheme when implemented in highly oscillatory optical self-focusing wave propagation simulations. © 2018 Elsevier B.V. All rights reserved.

1. Introduction The field of computational partial differential equations is an ever-evolving area of numerical mathematics. This growth is fueled by the physical significance of many such equations. And yet, even with this driving force behind the various advances of recent years, the theory of numerical PDEs is still incomplete [1]. A crucial property required for the study of numerical PDEs is stability. In this paper we will center our discussion on the study of two types of stability for a decomposed compact finite difference scheme used for solving a highly oscillatory Helmholtz equation. These two types of stability will emphasize different physical quantities. Before continuing this discussion, we will first introduce our partial differential equation. Maxwell’s field equations play a major role in electrodynamic theory, electrical circuit designs, and subwavelength metal optics [2–6]. These fundamental fields provide the ideal theoretical framework for many innovative technologies. Maxwell’s equations describe the propagation and interaction of electric and magnetic fields. However, these first-order partial differential equations are unsuitable for immediate use in conventional initial–boundary value problem computations. On the other hand, when light beam propagations are considered, these equations can be decoupled to yield the following time-dependent Helmholtz equation [2–4,7,8]:

( ) wtt = c 2 wxx + wyy + wzz , (x, y) ∈ D2 , z > 0, t > t0 ,

*

(1.1)

Corresponding author. E-mail address: [email protected] (T. Jones).

https://doi.org/10.1016/j.cam.2018.03.031 0377-0427/© 2018 Elsevier B.V. All rights reserved.

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

2

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

where w = w (x, y, z , t) is the electric field intensity function, z is the beam propagation direction of the optical beam, D2 is the transverse domain, and c is the speed of light in vacuum. When considering monochromatic light, we may 2π iκ0 t denote , where Ψ is the complex wavefunction, κ0 is the wave frequency in free space, and √ w(x, y, z , t) = Ψ (x, y, z)e i = −1 [2,9]. Utilizing this substitution, we observe that

Ψxx + Ψyy + Ψzz = −κ 2 Ψ , (x, y) ∈ D2 , z > 0,

(1.2)

where κ = 2π κ0 /c > 102 is the wave number of the optical beam. Finally, let U be the complex envelope of Ψ , meaning U(x, y, z) = Ψ (x, y, z)eiκ z . Hence, from (1.2) we acquire 2iκ Uz = Uxx + Uyy + Uzz , (x, y) ∈ D2 , z > 0.

(1.3)

A multitude of recent publications investigate fast and effective numerical methods for computing oscillatory solutions of (1.1)–(1.3) with relatively small wave numbers. However, the exploration of highly accurate and efficient algorithms for solving highly oscillatory beam propagation problems beyond finite-difference time-domain methods is still developing [7,10]. Fourier integral formalism and diffraction continue to prevail as solution procedures for paraxial optical systems [11]. We consider cases where radially symmetric electric fields are expected in the transverse direction. The singularity resulting from the polar coordinate transform is successfully removed through a domain decomposition strategy in the transverse direction [8,12,13]. We are primarily interested in schemes that are highly accurate in the transverse direction for subwavelength metal optical applications, e.g., those in nanophotonics, metamaterials, and plasmonics. We must also ensure computational efficiency and effectiveness. Considerable attention is given to the investigation of the numerical stabilities for the constructed decomposed numerical method. In the next section, we first transform (1.3) to a polarized form. The polarized problem is decomposed via a separation of the transverse domain. The singularity emerging from polarization is removed, and a global accuracy of the approximation is ensured through asymptotic expansions and auxiliary comparisons. A tridiagonal linear system is generated once a full discretization is implemented. We demonstrate in Section 3 that our underlying scheme is not linearly stable according to a conventional von Neumann stability analysis. Leading to Section 4, where a more appropriate asymptotic stability definition is discussed for highly oscillatory wave phenomena. It is shown that the constructed numerical method is unconditionally stable in the asymptotic sense. It should be noted that throughout this paper ∥·∥ represents any appropriate matrix norm. In Section 5, the numerical method is alternatively validated in computational experiments, especially when circularly polarized light is concerned. Computer simulations for the highly oscillatory optical beam produce meaningful propagation results, correctly reproducing the anticipated self-focusing phenomenon. Consideration is also given to the numerical order of accuracy for the method, and the accompanying numerical convergence figures are examined. Finally, concluding remarks are given in Section 6. 2. Compact decomposition algorithm For many monochromatic light beam propagation applications, electric fields in transverse direction, or optical tunnel gaps, may be taken to be radially symmetric [4,9,14]. When D2 is radially symmetric, then (1.3) can be conveniently reformulated as 2iκ Uz = Urr +

1 r

Ur + Uzz , 0 < r ≤ R, z > 0,

via a standard polar transformation. Letting β = −i/(2κ ) results in Uz = β

( Urr +

1 r

) Ur + Uzz

, 0 < r ≤ R, z > 0.

(2.1)

Next, adopt a transparent, or reflective, boundary condition [15,16]: Ur |r =R = 0, z ≥ 0,

(2.2)

where Ur |r =R is the outgoing normal derivative of U along the boundary of the two-dimensional domain D2 . We also consider a typical Gaussian beam type initial function as well as a paraxial approximation for initialization, U(r , 0; z0 ) =

e−r

2 /[2(1+iz

0 )]

1 + iz0

Uzz (r , 0) = 0, 0 ≤ r ≤ R,

, 0 ≤ r ≤ R,

(2.3) (2.4)

where κ z0 > 0 is the anticipated focusing location of the light beam caused by the physical phenomenon referred to as the Kerr effect [2]. Due to the singularity of the model equation as r → 0+ , we desire to decompose the domain D2 . With this in mind, Taylor series expansions are utilized to produce the following result. Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

3

Theorem 2.1 (Approximation for Small Radii [8]). Let 0 < r˜ ≪ R and U be sufficiently smooth in the transverse direction. Then for 0 < r < r˜ and z > 0, 1 r

Ur (r , z) − Urr (r , z) = g(U) −

r4 30

Ur 6 (ξ1 , z) −

r4 24

Ur 7 (ξ0 , z)(ξ2 − ξ1 ),

where r r2 r3 g(U) = − Urrr (0, z) − Ur 4 (0, z) − Ur 5 (0, z) 2 3 8 and 0 < ξk < r˜ , k = 0, 1, 2. Due to the symmetry of the radial domain, odd derivatives across the origin are zero. Thus, we take g(U) = −

r2 3

Ur 4 (0, z).

Through the use of Theorem 2.1, we reconstitute (2.1) into a system of partial differential equations in connected transverse domains. Uz = β (2Urr + Uzz + g(U)) , 0 ≤ r < r˜ , ) ( 1 Uz = β Urr + Ur + Uzz , r˜ ≤ r ≤ R, r

(2.5) (2.6)

where 0 < r < r˜ ≪ R and the error induced is O(r 4 ) as r → 0+ . Furthermore, since the singularity has been removed, for calculations involving (2.5) we may consider a slightly extended inner transverse domain, [0, r˜ ). We also note that (2.5)–(2.6) provide the potential for a multiscaled modeling orientation which is particularly useful to computations involving circularly polarized nanotunnel waves [14]. Now, we construct a numerical solution of the initial–boundary value problem (2.2)–(2.6). Let h = R/(n − 1) < r˜ be sufficiently small, and denote the superimposed corresponding mesh region on [0, R] as Dh = {rk : rk = (k − 1)h, k = 1, 2, . . . , n}. Due to the symmetry of U in the transverse direction, we may require U(−h, z) = U(h, z), z ≥ 0.

(2.7)

In the following difference approximations for derivatives in the transverse direction, we include traditionally truncated terms for additional accuracy in the transverse domain. Uk+1 − 2Uk + Uk−1 h2 Uk+1 − Uk−1 2h

= (Urr )k + = (Ur )k +

2h2 4! h

(Ur 4 )k +

2

3!

(Ur 3 )k +

h

2h4 6! 4

5!

(Ur 6 )k + · · ·

(Ur 5 )k + · · · .

Utilizing these approximations, we obtain the O(h6 ) semidiscretized system shown below from (2.5)–(2.6)

β h2 4β h4 − (Ur 4 )k − (Ur 6 )k 6! ( 36 )! + β (Uzz )k + β g(Uk ) + O h , k = 1, 2, . . . , k˜ , ( ) Uk+1 − 2Uk + Uk−1 β Uk+1 − Uk−1 (Uz )k = β + + β (Uzz )k h2 rk 2h [ 2( ) ( h 1 1 h4 1 −β (Ur 4 )k + (Ur 3 )k + (Ur 6 )k 3! 2 rk 5! 3 )] ( ) 1 + (Ur 5 )k + O h6 , k = k˜ + 1, k˜ + 2, . . . , n, (Uz )k = 2β

Uk+1 − 2Uk + Uk−1 h2

rk

(2.8)

(2.9)

˜ < r˜ ≤ (k˜ + 1)h [13]. Note that (2.2) and (2.7) can as h → 0+ , where k˜ ≥ 1 is a sufficiently small integer for which 0 < kh be used to determine the values of Un+1 and U0 , respectively. In order to eliminate the higher order derivatives in the above system, we may utilize (2.5)–(2.6) to create a compact algorithm. Considering these underlying differential equations is often advantageous when developing proper compact or essentially compact algorithms [17–19]. Continuing differentiations of (2.5) produce Urrrr =

1 2β

Uzrr −

1 2

Uzzrr +

1 3

Ur 4 (0, z), 0 ≤ r < r˜ , z > 0.

(2.10)

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

4

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

In the same manner, from (2.6) we acquire

Urrrr +

1 r

1

Urrr =

β

Urrr =

β

1

1 Ur − Uzzr , r˜ ≤ r ≤ R, z > 0, r2 2 2 + 2 Urr − 3 Ur − Uzzrr , r˜ ≤ r ≤ R, z > 0. r r

Uzr − Uzrr

1 r

Urr +

Combining these equations yields Urrrr +

2 r

1

Urrr =

β

Uzrr + 1

1

βr

Uzr − Uzzrr −

1 r

Uzzr

1

Ur , r˜ ≤ r ≤ R, z > 0. r3 Now we employ the results above by substituting (2.10) and (2.11) into (2.8) and (2.9), respectively. This leads to

+

(

r2

Urr −

Uk+1 − 2Uk + Uk−1

)

(Uz )k+1 − 2(Uz )k + (Uz )k−1 + β (Uzz )k − h2 12 ( 4) (Uzz )k+1 − 2(Uzz )k + (Uzz )k−1 ˜ k) + O h , +β + β g(U 12 k = 1, 2, . . . , k˜ , z > 0; ) ( Uk+1 − Uk−1 Uk+1 − 2Uk + Uk−1 + + β (Uzz )k (Uz )k = β h2 2rk h

(Uz )k = 2β

(2.11)

(2.12)

[ β h (Uz )k+1 − 2(Uz )k + (Uz )k−1 (Uz )k+1 − (Uz )k−1 + 12 βh 2β rk (Uzz )k+1 − 2(Uzz )k + (Uzz )k−1 (Uzz )k+1 − (Uzz )k−1 − − h 2rk ] ( ) Uk+1 − 2Uk + Uk−1 Uk+1 − Uk−1 + − + O h4 , 2 3

−

rk h

2rk

k = k˜ + 1, k˜ + 2, . . . , n, z > 0,

(2.13)

as h → 0 , where +

˜ k ) = g(Uk ) − g(U

h2 18

Ur 4 (0, z) = −

(

h2 18

+

rk2 3

)

Ur 4 (0, z).

Remark 2.1. Since the application of Theorem 2.1 restricts (2.12) to being at most an order four approximation, the O(h6 ) terms have been truncated. Thus (2.12)–(2.13) are order four approximations to (2.5)–(2.6) as well as (2.1). This order of accuracy can be raised if higher order difference approximations and additional expansion terms for g are utilized. However, this may considerably increase the intricacy of the compact algorithm. To fully discretize the differential equations, adopt standard central difference approximations for the z-derivatives, that is, j+1

( ) − Ukj−1 + O τ 2 , 0 < τ ≪ 1, j > 0, 2τ j+1 ( ) U − 2Ukj + Ukj−1 = k + O τ 2 , 0 < τ ≪ 1, j > 0, 2 τ

(Uz )(k,j) = (Uzz )(k,j)

Uk

where τ → 0+ is the variable z-step which can be determined via suitable monitoring functions [3,13]. Therefore, from (2.12)–(2.13) we acquire the full discretization j+1

Uk

] [ ] ( )[ j j j j+1 j j−1 Uk+1 − 2Uk + Uk−1 − Ukj−1 Uk − 2Uk + Uk h2 = 2β 1 − +β 2τ 12τ 2 h2 τ2 [ j+1 ] j+1 j+1 j−1 j−1 j−1 1 Uk+1 − 2Uk + Uk−1 − Uk+1 + 2Uk − Uk−1 − 12 2τ [ j+1 ] j+1 j+1 j−1 j−1 j−1 β Uk+1 − 2Uk + Uk−1 + Uk+1 − 2Uk + Uk−1 + 12 τ2 ( ) ( ) + β g˜ Ukj + O h4 + τ 2 , k ∈ θ0 ;

(2.14)

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)[

]

[

)

–

5

]

( j j j j+1 j j−1 Uk+1 − 2Uk + Uk−1 − Ukj−1 Uk − 2Uk + Uk h2 + β =β 1− 2τ h2 τ2 12rk2 [ [ ] ] ( ) j+1 j+1 j−1 j−1 j j Uk+1 − Uk−1 − Uk+1 + Uk−1 Uk+1 − Uk−1 h2 h β 1+ − + rk 2h 24rk 2τ 12rk2 [ j+1 ] j+1 j+1 j−1 j−1 j−1 1 Uk+1 − 2Uk + Uk−1 − Uk+1 + 2Uk − Uk−1 − 12 2τ [ j+1 ] j+1 j+1 j j j j−1 j−1 j−1 β Uk+1 − 2Uk + Uk−1 − 2Uk+1 + 4Uk − 2Uk−1 + Uk+1 − 2Uk + Uk−1 + 12 τ2 [ j+1 ] j−1 j−1 j j j+1 β h Uk+1 − Uk−1 − 2Uk+1 + 2Uk−1 + Uk+1 − Uk−1 + 24rk τ2 ) ( 4 (2.15) + O h + τ 2 , k ∈ θ1 , { } { } where θ0 = 1, . . . , k˜ , θ1 = k˜ + 1, k˜ + 2, . . . , n , and h, τ → 0+ . Note that the values of U 0 and U 1 can be determined j+1

Uk

via (2.3) and (2.4), respectively. Remark 2.2. The approximations used in the beam propagation direction are O(τ 2 ) as τ → 0+ . Since c = τ 2 /h2 should be observed in computations due to the Courant–Friedrichs–Lewy condition, where c is the Courant number, we may choose τ 2 = ch2 to achieve an overall second order accuracy for the scheme. Note that the scheme is O(h4 ) as h → 0+ , which is preferable in computations involving polarized waves in nanomaterial computations where high oscillations are anticipated in the radial direction. Remark 2.3. Although the radial step-size in the inner and outer domains, θ0 and θ1 , are both represented by h, it is reasonable to generate a multiscale method. This approach allows for a smaller step-size in the inner domain, θ0 , which is especially useful due to the expected focusing of the optical beam. This will be further discussed in Section 4. Denote σ =

τ h

, λ =

12βσ 2

τ

, µ± k = 2 ±

h rk

, νk± = 1 ±

h2 12rk2

, γ± = 1 ±

2β

τ

. Drop all truncation errors from (2.14)–

(2.15) and rearrange terms to organize the scheme in a conducive manner. As a result, we generate the subsequent compact decomposition algorithm for the numerical solution of the optical beam problem (2.2)–(2.6), j+1

j+1

ak Uk−1 + bk Uk

+ ck Ukj++11 = dk Ukj −1 + ek Ukj + fk Ukj +1 + 24βτ gkj + lk Ukj−−11 + mk Ukj−1 + nk Ukj−+11 k = 1, 2, . . . , n, j > 0,

(2.16)

together with (2.3), where ak =

ck

{ − γ , − µ− k γ /2,

{ − γ , = − µ+ k γ /2,

k ∈ θ0 , k ∈ θ1 ; k ∈ θ0 , k ∈ θ1 ;

⎧ ⎨ −8β (5 + 12σ 2 ) , ek = τ ) ⎩ ( − −4 λνk + 10β/τ , { ˜ kj ), k ∈ θ0 , j g(U gk = 0, k ∈ θ1 ; { 10γ + , k ∈ θ0 , mk = 10γ + , k ∈ θ1 ;

{ bk

=

10γ − , 10γ − ,

k ∈ θ0 , k ∈ θ1 ;

⎧ 4β ( ) ⎪ 1 − 12σ 2 , ⎨− τ dk = ⎪ ⎩2λν − − λh ν + − 2β µ− , k rk k τ k ⎧ 4β ( ) ⎪ 1 − 12σ 2 , ⎨− k ∈ θ0 , τ fk = ⎪ ⎩2λν − + λh ν + − 2β µ+ , k ∈ θ1 ; k rk k τ k { + γ , k ∈ θ0 , lk = + µ− k ∈ θ1 ; k γ /2, { + γ , k ∈ θ0 , nk = + µ+ γ / 2 , k ∈ θ1 . k

k ∈ θ0 , k ∈ θ1 ; k ∈ θ0 , k ∈ θ1 ;

(

j

j

j

j

j

We also wish to reformulate the scheme into the following equivalent matrix form. Let U j = U1 , U2 , . . . , Uk , . . . , Un−1 , Un

)⊺

∈ Cn be the solution at the jth z-step. Then, AU j+1 = BU j + CU j−1 + p j , j = 1, 2, 3, . . . ,

(2.17)

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

6

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

where A, B, C) ∈ Cn×n , p ∈ Cn ,(and the entries ( ) are defined by (2.16) as follows. A = tridiag 𝒹 k , ℯ k , 𝒻 k , and C = tridiag 𝓁 k , 𝓂 k , 𝓃 k , where 0, ak , an + cn ,

{

𝒶k =

𝒷 k = bk ∀ k, {

a1 + c1 ,

𝒸 k = ck , 0,

0,

{

𝓁 k = lk ,

ln + nn ,

k = 1, k = 2, . . . , n − 1, k = n; k = 1, k = 2, . . . , n − 1, k = n; k = 1, k = 2, . . . , n − 1, k = n;

0,

{

𝒹 k = dk ,

dn + fn ,

ℯ k = ek ∀ k, {

d1 + f1 ,

𝒻 k = fk , 0,

l1 + n1 ,

{

𝓃 k = nk , 0,

( ) 𝒶 k , 𝒷 k , 𝒸 k , B = tridiag

k = 1, k = 2, . . . , n − 1, k = n; k = 1, k = 2, . . . , n − 1, k = n; k = 1, k = 2, . . . , n − 1, k = n;

and 𝓂 k = mk ∀ k. j

Furthermore, p j = (p1 , p2 , . . . , pn )⊤ where the only nontrivial entries are the first k˜ entries defined through gk in (2.16). 3. Conventional linear stability analysis Conventional linear stability analysis for methods such as (2.16) emphasizes the influence of the spacial step size in the radial direction, h. Since τ = σ h for our method, the step size in the propagation direction, τ , is also accounted for in the analysis. To begin, consider a general definition for conventional linear stability. Definition 3.1 (Conventional Linear Stability). Consider a finite difference method with an amplification matrix Φ for solving a partial differential equation, where h is the spacial step size. We say that the numerical method is stable if there exists a constant p ∈ R+ such that

( ) ∥Φ ∥ < 1 + O hp as h → 0+ , where ∥·∥ is an appropriate matrix norm and the constant p is the conventional stability index of the method. While this form of the definition is often difficult to implement, it clearly shows the relationship between conventional stability and the physical step size. An alternative method for determining conventional stability is von Neumann analysis. Fortunately, it is appropriate to use von Neumann analysis to investigate the stability of (2.16). To that end, we require an additional theorem. Theorem 3.1 (Simple von Neumann Polynomial Criteria [20]). Let φd (z ) = φd∗ (0)φd (z )−φd (0)φd∗ (z )

∑d

n=0 an z

n

, φd∗ (z ) =

∑d

n=0 ad−n z

n

, and φd−1 =

, where d ∈ Z and z ∈ C. φd is a simple von Neumann polynomial if and only if either z ⏐ ∗ ⏐ (a) |φd (0)| < ⏐φd (0)⏐ and φd−1 is a simple von Neumann polynomial, or (b) φd−1 is identically zero and φd′ is a Schur polynomial. +

Remark 3.1. Notice that φd−1 is intended to be an order d − 1 polynomial. Thus, an iterative process can be derived from the above theorem. The process may be repeated up to d − 1 times until the resulting polynomial is either linear or identically equal to zero. This process will be used in Case I. Remark 3.2. If φd−1 is an order d − 2 polynomial which is not identically zero, then the above theorem implies that the scheme is unstable in the von Neumann sense. This fact will be used in Case II. Theorem 3.2. The compact decomposition scheme (2.16) is unconditionally unstable in the von Neumann sense. Proof. Notice that the z-evolution system (2.16) consistently approximates the Helmholtz equation (2.1) together with the boundary condition (2.2). Therefore, von Neumann analysis is fitting in the investigation of stability for our scheme [15,21]. j We start by considering the perturbation function εk = ξ j eα ikh , where ξ is the amplification factor, ε is the perturbed solution, and α ∈ R. Next, we substitute this function into the perturbation equation, j+1

j+1

ak εk−1 + bk εk

+ ck εkj++11 = dk εkj −1 + ek εkj + fk εkj +1 + lk εkj−−11 + mk εkj−1 + nk εkj−+11 , k = 1, 2, . . . , n, j > 0.

(3.1)

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

7

This results in

) ) ( ( ξ 2 ak e−αih + bk + ck eαih = ξ dk e−αih + ek + fk eαih + lk e−αih + mk + nk eαih . Case I (k ∈ θ0 ). From the above equation, we have 0 = ξ

2

2β

[(

)

(

α ih

−α ih

2β

(

)

)]

1− + 10 1 − e +e τ τ [( ) ] ) ) ( ) ( 8β ( 4β −ξ − eα ih + e−α ih − 1 − 12σ 2 5 + 12σ 2 τ τ [( ( ) )] ) 2β ( α ih 2 β −α ih − 1+ + 10 1 + e +e , τ τ ) ] [ ] [ ( ] 4 [ i 0 = ξ2 2 1 + 12σ 2 (1 − cos α h) + (5 + cos α h) (5 + cos α h) − ξ i τκ τκ [ ( ) ] i − 2 1− (5 + cos α h) . τκ Utilizing Theorem 3.1,

[ (

)

]

] 4i [ 12σ 2 (1 − cos α h) + (5 + cos α h) z φ2 (z ) = 2 1 + (5 + cos α h) z 2 − τκ τκ [ ( ) ] i − 2 1− (5 + cos α h) , τκ ) ] [ ( ] 4i [ i ∗ φ2 (z ) = − 2 1 + 12σ 2 (1 − cos α h) + (5 + cos α h) z (5 + cos α h) z 2 + τκ τκ [ ( ) ] i + 2 1− (5 + cos α h) . τκ i

Thus, we may conclude that φ1 (z ) ≡ 0. Hence, we require φ2′ to be a Schur polynomial, meaning all roots rn are such that |rn | < 1. To that end, consider

) ] [ ( ] 4i [ i φ2′ (z ) = 4 1 + 12σ 2 (1 − cos α h) + (5 + cos α h) (5 + cos α h) z − τκ τκ Note the only root r , 4i

r =

=

τκ

[

12σ 2 (1 − cos α h) + (5 + cos α h)

]

[ (

] (5 + cos α h) ( )] 1 − cos α h 2 [1 + iκτ ] . 1 + 12σ 5 + cos α h

4 1+

i

τκ

)

[

1 1 + κ 2τ 2

Observe that 0 ≤ 1 − cos α h ≤ 2 and 4 ≤ 5 + cos α h ≤ 6, so 0 ≤ which is a constant, we have 1

[

1 − cos α h

≤ 12 . Since σ 2 = τ 2 /h2 is the Courant number,

)]

|1 + iκτ | √ [ ] ( ) 1 1 + 12σ 2 + 36σ 4 1 / 2 2 2 2 ≤ 1 + 6 σ 1 + κ τ = . 1 + κ 2τ 2 1 + κ 2 σ 2 h2 ( ) ( ) In order for |r | < 1 we must require κ 2 > 12 3σ 2 + 1 /h2 . But this is contradictory since 1/h2 = O n2 as h → 0+ , or equivalently n → ∞. Next, we consider our outer scheme. Case II (k ∈ θ1 ). Continuing to use Theorem 3.1, we find that [( ) ( ) ( )] ) ) 2β ( α ih h 2β ( α ih 2β 2 −α ih −α ih 0 = ξ 1− e +e + 1− e −e + 10 1 − τ 2rk τ τ [( ( ) ) 2 2 ) 24βσ h 4β ( α ih −ξ 1− − e + e−α ih τ τ 12rk2 ( ( ) ) ( ( ) )] ) 12βσ h2 2β h ( α ih 48βσ h2 40β −α ih + 1+ − e − e + − 1 − − rk τ rk τ τ 12rk2 12rk2 |r | =

1 + κ 2τ 2

1 + 12σ

(

1−cos α h 5+cos α h

2

5 + cos α h

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

8

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

[( −

1+

2β

)

τ

[ (

0 = ξ2 2 1 +

eα ih + e−α ih +

)

(

i

(

h

1+

2rk

)

h

(

2β

)

τ 1

)

(

eα ih − e−α ih + 10 1 +

)

(

2β

)

)]

τ

–

,

]

i− sin α h (5 + cos α h) + rk κτ [ ( ( ) ) −4i h2 +ξ 6σ 2 1 − 1 − cos α h + 5 + cos α h ( ) ( ) κτ 12rk2 ( ( ) ) ] 2 h2 h − 6σ 1 + − sin α h κ rk τ 12rk2 ) ( ) ] [ ( h 1 i i+ sin α h . − 2 1− (cos α h + 5) + κτ rk κτ

κτ

Similar to Case I,

[ ( ) ( ) ] i h 1 φ2 (z ) = 2 1 + i− sin α h z 2 (5 + cos α h) + κτ rk κτ [ ( ( ( ( ) ) ) ) ] 2 h2 2 h −4i h + 6 σ 1 + 6σ 2 1 − − − sin α h z 1 − cos α h + 5 + cos α h ( ) ( ) κτ κ rk τ 12rk2 12rk2 [ ( ) ( ) ] i 1 h − 2 1− i+ sin α h , (cos α h + 5) + κτ rk κτ ( [ ( ) ) ] 1 i h −i + φ2∗ (z ) = − 2 1 + sin α h z 2 (5 + cos α h) + κτ rk κτ ) ) ) [ ( ( ( ( ) ] 2 2 h 4i h h2 1 − cos α h + 5 + cos α h − − + 6σ 2 1 − 6 σ 1 + sin α h z ( ) ( ) κτ κ rk τ 12rk2 12rk2 [ ( ) ( ) ] i h 1 + 2 1− −i − sin α h . (cos α h + 5) + κτ rk κτ After straightforward multiplications and cancellations, we arrive at 8

(

i

1−

)[

6σ

(

h2

)

] h − (5 + cos α h) sin α h κ rk κτ τ 12rk2 )( ( )[ ( ) ] 1 − cos α h 8h i h2 − 1− 6σ 2 1 − + 1 (5 + cos α h) sin α h. κτ rk κτ 5 + cos α h 12rk2

φ1 (z ) = −

1+

Since φ1 (z ) is a nontrivial constant polynomial, we conclude that our scheme is unstable in the von Neumann sense unconditionally. This completes our proof. □ 4. Asymptotic stability analysis In general, relating stability to the physical step size is appropriate. However, since (2.2)-(2.6) is a highly oscillatory problem, there is an alternative physical quantity that may be utilized. We wish to consider a definition for the stability of a scheme which emphasizes the influence of the wave number, κ. Definition 4.1 (Asymptotic Stability). Consider a finite difference method with an amplification matrix Φ for solving an oscillatory problem associated with a high wave number, that is κ ≫ 1. We say that the numerical method is asymptotically stable if there exists a constant p ∈ R+ such that

∥Φ ∥ < 1 + O

(

1

κp

)

as κ → ∞,

where ∥·∥ is an appropriate matrix norm and the constant p is called the asymptotic stability index of the method. Theorem 4.1. The compact decomposition scheme (2.16) is unconditionally asymptotically stable with index p = 1. Proof. Without loss of generality, we begin with a perturbation of our numerical solution AU˜ j+1 = BU˜ j + C U˜ j−1 + p j , j = 1, 2, 3, . . . . By subtracting the above from the numerical solution (2.17), we have Aϵ

j+1

= Bϵ j + C ϵ j−1 ,

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

9

where ϵ j represents the difference in the jth z-step. Thus

ϵ j+1 = A−1 Bϵ j + A−1 C ϵ j−1 . Therefore

 j+1   −1 j      ϵ  = A Bϵ + A−1 C ϵ j−1  ≤ A−1 Bϵ j  + A−1 C ϵ j−1        ≤ A−1 B ϵ j  + A−1 C  ϵ j−1    )     ( ≤ A−1 B + A−1 C  max{ϵ j  , ϵ j−1 }   )   ( = A−1 B + A−1 C  ϵ j  .     This leads to the amplification matrix Φ , where ∥Φ ∥ ≤ A−1 B + A−1 C  . For simplicity in notations, we define the following.

(

)

A˜ = tridiag a˜ k , b˜ k , c˜k ,

(4.1)

where

⎧ 1, ⎪ ⎪ ⎨ h a˜ k = 1 − , ⎪ 2rk ⎪ ⎩ 2,

k = 2, . . . , k˜ k = k˜ + 1, . . . , n − 1, c˜k = k = n;

⎧ 2, ⎪ ⎪ ⎨1,

k = 1, k = 2, . . . , k˜ ,

h ⎪ ⎪ ⎩1 + , 2rk

k = k˜ + 1, . . . , n − 1;

and b˜ k = 10 ∀ k. Notice that A = γ − A˜ and C = γ + A˜ . Since |¯z /z | = 1 for all z ∈ C,

⏐   ⏐ +⏐ ⏐1 −  −1   1 −1  ⏐ ⏐ A C  =  A˜ × γ + A˜  = ⏐ γ ⏐ ∥I∥ = ⏐⏐ γ−  ⏐γ− ⏐ ⏐1 +

i

κτ i

κτ

⏐ ⏐ ⏐ ⏐ = 1. ⏐

By the same token, every entry of the matrix B contains β and τ . so we may rewrite B as B = κ. This leads to

(4.2) β˜ B, τ

where B˜ is independent of

⏐⏐  ⏐ ⏐    ⏐ −1     −1   1 −1 β  ⏐⏐ 1 ⏐⏐ ⏐ ⏐− i ⏐  A B =  A˜ × B˜  ≤ ⏐ A˜  B˜  . ⏐  ⏐ γ− i ⏐ 1 + κτ ⏐ 2κτ ⏐ τ ⏐ ⏐ ⏐ ⏐ 1 Since ⏐ 1 i ⏐ = ⏐⏐ 1 ⏐⏐ = < 1 and A˜ , B˜ are independent of κ, then 1 1+ κτ

⏐1+ κτi ⏐

1+ 2 2 κ τ

( )     −1  1 −1   ˜  A B < 1  ˜ , as κ → ∞. A  B = O 2κτ κ

(4.3)

Lastly, combining (4.2) and (4.3) yields

    ∥Φ ∥ ≤ A−1 B + A−1 C  < 1 + O

( ) 1

κ

, as κ → ∞.

This completes our proof. □ Remark 4.1. Due to its definition, the asymptotic stability is not affected by particular choices of grids following required approximation consistencies. This feature can be utilized to accommodate favorable multiscale settings, in which a microscale is used for the inner domain θ0 , and a macroscale is considered for the outer domain θ1 . The use of fully-adaptive grids is also possible. A natural question may be to consider whether linearly unstable yet asymptotically stable decomposed compact schemes work satisfactorily for simulating highly oscillatory waves in applications, especially those associated with subwavelength computations. To that end, in the next section we will present numerical results from a sequence of typical computational experiments involving the extremely sensitive optical focusing prediction [7,9]. 5. Computational experiments Optical beams are often modeled by radially symmetric initial–boundary value problems, such as (2.1)–(2.4). Cases in which the wave focuses at κ z0 units away from the lens are of particular interest. This phenomenon is often induced by beam-material interactions, specifically exposure of the media to intense electromagnetic radiation resulting in a variance of the refractive index of the material involved. Often referred to as the Kerr effect, this physical phenomenon is frequently used Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

10

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

Fig. 5.1. The real part (left), imaginary part (right), and modules (solid blue in both images) of the initial data, U(r , 0) where −R ≤ r ≤ R. Note that the complex wave function is initially highly oscillatory. n = 16001 mesh points in the transverse, or radial, direction are utilized, and the parameter z0 = 100 is considered. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

in testing new computational algorithms and evaluating their effectiveness, accuracy, and efficiency [2,12]. Our numerical experiments and verifications are constructed to (employ)this self-focusing phenomenon. 1/2 Let κ ≥ 104 . For 1 ≤ z0 ≤ 100, we set R = 4 1 + z02 . These allow a maximal focusing distance up to 100κ which is sufficient for most realistic applications [7,8]. A greater R may also be considered, though it often demands additional computations. On the other hand, smaller R values may not be adequate for simulating details of the focusing phenomenon [9,19]. To summarize, our numerical solutions must quickly reproduce the anticipated self-focusing phenomenon at the correct location, κ z0 , away from the initial point, to be considered satisfactory. For any suitably large integer n, we have defined h = R/(n − 1) as the uniform mesh step in the radial, or transverse, direction. The compact domain decomposition algorithm (2.16) has been shown to be implicit and stable in the asymptotic sense. Thus, we may utilize a wide range of Courant numbers from (5, 2000), and still preserve the accuracy of the scheme. Remark 5.1. While further investigation is needed, we have found requiring n/κ = o(1) produces satisfactory results in computations. This impression is further reaffirmed by considering (4.3), which gives

( )    (n)  −1  1   A B < 1  =O , as κ, n → ∞. A˜ −1  B˜  = O 2κτ κτ κ By requiring n/κ = o(1), we can bound the amplification matrix,

    ∥Φ ∥ ≤ A−1 B + A−1 C  < 1 + O (ϵ) , as ϵ → 0. Therefore, the numerical analysis and computational observations are in agreement. Computations are primarily conducted on a Matlab platform with the parallel computing toolbox on Kodiak, a high performance Cray CS400-AC cluster running CentOS at Baylor University. Kodiak contains 112 computer nodes, each with 256 GB of RAM and dual 18-core Intel E5-2695 V4 processors operating at 2.1 GHz. For networked storage and message passing, an Intel Omni-Path 100 Gb/s network is used. The shared storage capacity is 720 TB, and the majority of computational tasks are finished in minutes. In computations, we have observed that the choice of the domain separation parameter k˜ is often delicate. We frequently let k˜ be between 5 and 20, leading to very favorable numerical solutions in the full transverse domain. To begin our experiments, Fig. 5.1 displays the real and imaginary parts of the electric field initial function U(r , 0), where z0 = 100. The modulus is also shown in each image. Due to the symmetry of the electric field in the radial direction, we also smoothly extend the initial data to the larger interval [−R, R] . Following the implementation of our algorithm (2.16), we arrive at the desired focusing solution. Fig. 5.2 displays the propagation of the solution in the z-direction on the left and the associated contour maps on the right. The solution U(r , κ z0 ) at the focusing location κ z0 = 106 is also emphasized. For clearer viewing of the behavior near the focusing location, only a portion of the radius and propagation are shown. The numerical algorithm (2.16) remains stable past the focusing location. Additionally, it is interesting to view the focused solution in light of the highly oscillatory initial data from Fig. 5.1. In Fig. 5.3, we show the trajectory of the maximal modulus of the resultant electric field function in order to better ⏐ ⏐ see

⏐ j⏐

the more dynamic features of the beam propagation. The maximal modulus of U j is defined by Mj = max1≤k≤n ⏐Uk ⏐ , j =

0, 1, 2, . . ., m, where m is a large integer such that mτ > κ z0 . It is clear to see that the numerical solution correctly reproduces the desired focusing characteristics. The numerical method demonstrates satisfactory numerical stability, both before and after the focusing location. This is particularly favorable since well over 4,000,000 iterations of the algorithm Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

11

Fig. 5.2. The real part (top row), imaginary part (middle row), and modulus (bottom row) of the evolution of the numerical solution, U(r , z) with z0 = 100 and κ = 104 , as well as corresponding contour maps. A large Courant number, τ 2 /h2 ≈ 100, is chosen with a total of 16001 mesh points in the transverse direction. Note that for the imaginary part the valley and peak occur just before and after the focusing location, respectively.

Fig. 5.3. Trajectories of the maximal modulus function are shown. The peak of the curve represents the focusing of the beam when κ z0 = 106 . The focusing distance and magnitude are extremely accurate, with an error of only 0.028% and 0.0022% respectively.

are computed as shown. The method has also been shown to be reliable in target-oriented operations when compared to standard results [4,22,23]. We wish to more closely investigate the behavior of our numerical algorithm at the focusing location. In Fig. 5.4, we provide a more detailed form of the focusing solution. The figures to the left are produced by rotating the vector focusing solution U(r , κ z0 ) about the origin, since radially symmetric fields are utilized. The figures to the right are the contour maps of the left-hand counterparts. Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

12

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

Fig. 5.4. The real part (top row), imaginary part (middle row), and modulus (bottom row) of the rotated focusing solution, U(r , κ z0 ) with z0 = 100 and κ = 104 , are shown to the left. The corresponding contour maps are to the right. A large Courant number, τ 2 /h2 ≈ 100, is chosen with a total of 16001 mesh points in the transverse direction. The figure is constructed using 1/40th of the total radius. Note that although the imaginary part remains oscillatory, its magnitude is relatively small.

Finally, the first 20 radial solution values at the focusing location κ z0 = 106 are shown in Fig. 5.5. With this view of the interface between the split transverse domains, θ0 and θ1 , observe that the solutions are free from nonphysical oscillations across this interface. An index of k˜ = 11 is used in the figure. The scale of Fig. 5.5 has also been restricted once again to allow for a more refined view of the focusing solution. Lastly, we explore the numerical stability of the algorithm further by calculating the order of accuracy at corresponding mesh points. Definition 5.1 (Order of Accuracy [24]). Consider a finite difference method with a sufficiently small spacial step size h and numerical solution u˜ h . Let u˜ h/2 and u˜ h/4 be the numerical solutions corresponding to the spacial step sizes h/2 and h/4, respectively. We say that the numerical scheme has an order of accuracy p such that

( log2

u˜ h − u˜ h/2 u˜ h/2 − u˜ h/4

)

= p + O(h) as h → 0+ .

p is also referred to as the convergence factor of the method. Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

13

Fig. 5.5. The real part (left), imaginary part (middle), and modulus (right) of the first 20 solution values located next to the origin at the focusing location. The index value k˜ = 11 is used in simulations.

Fig. 5.6. The modulus for the order of accuracy in the radial direction (left) and propagation direction (right). The pointwise convergence factors have been smoothed using an order three Savitzky–Golay filter (shown in blue), a linear least squares regression (shown in red), and the mean value (shown in green). The solution at the focusing location (shown in purple) has been vertically scaled to improve visibility. Note τ = 0.25 and h = 0.1, 0.05, 0.025 are used in calculations for the radial direction, while h = 0.025 and τ = 1, 0.5, 0.25 are used for the propagation direction. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

For the images in Fig. 5.6, we first calculate the appropriate pointwise convergence factors. These values are then smoothed using an order three Savitzky–Golay filter, linear least squares regression, and mean value. The affiliated focusing solution values, which have been vertically stretched to improve visibility, are also shown as a reference. Remark 5.2. Recall that the numerical algorithm (2.16) has a theoretical convergence of order four in the radial direction and order two in the propagation direction, since the error incurred is O(h4 + τ 2 ) as h, τ → 0+ . Notice that numerical order of accuracy is slightly less than anticipated in both the radial and propagation directions. Fortunately, this is not entirely unexpected due to the conventional linear instability discussed in Section 3. This implies a deeper relationship between conventional linear stability and asymptotic stability that is beyond the scope of this paper. 6. Conclusions and expectations A domain decomposed compact algorithm (2.16) for the numerical solution of the polarized Helmholtz equation on a radially symmetric electric field has been constructed and examined. Theoretical fourth order accuracy is maintained for the highly oscillatory component of the wave solution in the transverse direction. This particular feature is meaningful in simulations of plasmonic nanomaterials, which involve circularly polarized optical waves. Although it is shown that numerical method is unconditionally unstable in the conventional von Neumann sense, it is unconditionally stable in the asymptotic sense for large wave numbers. Additionally, computational experiments confirm the viability of the algorithm and the partial retention of the desired spacial orders of accuracy. The relationship between these two types of stability demands further investigation. The novel numerical method is simple, straightforward in structure, and proven to be easy to use in sophisticated laser optics and nanomaterials environments. But its computational speed, on the other hand, cannot contend with optimized integral methods, such as those based on Hankel or Fourier transformations. Since our computations are mainly conducted on uniform grids, the computational speed of the algorithm may be improved through the use of a multiscale approach in the radial direction. While it is possible to develop an asymptotically stable multiscaled methodology similar to (2.16), we may also consider implementing adaptive stepping in the future, particularly in the wave propagation direction. Certain Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.

14

T. Jones, Q. Sheng / Journal of Computational and Applied Mathematics (

)

–

crucial issues still need to be explored, including when and where an adaptation is implemented, which physical quantities to monitor, the construction of effective dual monitoring functions, and the stability of new algorithms when fully nonuniform finite differences are employed [3,25]. Immediate endeavors may include adopting a multiscale approach across the decomposed domain in the radial direction, utilizing dual (r , z)-stretching transformations [3], and dual moving meshes via proper arc-length monitoring functions [12,13]. The compact algorithm may also be reinforced through the use of decompositions, such as eikonal or asymptotic decompositions [15,25] as well as cosine splitting [10,13] for nonlinear nanooptical waves. Furthermore, investigations of in-depth stability analysis in connection with the highly oscillatory characteristics of optical waves with both linear and circular polarizations need to be implemented. Acknowledgments The authors would like to thank the editor and anonymous referees for their time spent and extremely valuable remarks given. Their suggestions have significantly improved the quality and presentation of this paper. The second author also appreciates the support of the U.S. Air Force Research Laboratory, as part of research Grant No. FA-8650-11-D-5400. References [1] A. Iserles, The joy and pain of skew symmetry, Found. Comput. Math. 16 (2016) 1607–1630. [2] Y.B. Band, Light and Matter: Electromagnetism, Optics, Spectroscopy and Lasers, John Wiley & Sons, West Sussex, 2006. [3] L.P. Gonzalez, S. Guha, J.W. Rogers, Q. Sheng, An effective z-stretching method for paraxial light beam propagation simulations, J. Comput. Phys. 227 (2008) 7264–7278. [4] M.N.O. Sadiku, Numerical Techniques in Electromagnetics, second ed., CRC Press, London and New York, 2000. [5] P.P. Crooker, W.B. Colson, J. Blau, Representation of a gaussian beam by rays, Amer. J. Phys. 74 (2006) 722–727. [6] R.K. Tyson, Principles and Applications of Fourier Optics, IOP Publishing, Bristol, 2014. [7] G.D. Gillen, K. Gillen, S. Guha, Light Propagation in linear Optical Media, first ed., CRC Press, London and New York, 2007. [8] T.N. Jones, L.P. Gonzalez, S. Guha, Q. Sheng, A continuing exploration of a decomposed compact method for highly oscillatory wave problems, J. Comput. Appl. Math. 299 (2016) 207–220. [9] S. Guha, Validity of the paraxial approximation in the focal region of a small-f -number lens, Opt. Lett. 26 (2001) 1598–1600. [10] B. Engquist, A. Fokas, E. Hairer, A. Iserles, Highly Oscillatory Problems, London Math Society, London, 2009. [11] J.A. Stratton, L.J. Chu, Diffraction theory of electromagnetic waves, Phys. Rev. 56 (1939) 99–107. [12] M.A. Beauregard, Q. Sheng, A fully adaptive method to approximate reaction diffusion equations of the quenching type over circular domains, Numer. Methods Partial Differential Equations 30 (2) (2014) 472–489. [13] Q. Sheng, Adaptive decomposition finite difference methods for solving singular problems, Front. Math. China 4 (2009) 599–626. [14] M.T. Sheldon, J. van de Groep, A.M. Brown, A. Polman, H.A. Atwater, Plasmoelectric potentials in metal nanostructures, Science 346 (2014) 828–831. [15] Q. Sheng, S. Guha, L.P. Gonzalez, A short note on the asymptotic stability of certain oscillation-free eikonal splitting schemes, Appl. Math. Lett. 25 (2012) 1539–1543. [16] D. Ruprecht, A. Schädle, F. Schmidt, L. Zschiedrich, Transparent boundary conditions for time-dependent problems, SIAM J. Sci. Comput. 30 (2008) 2358–2385. [17] W. E, J.G. Liu, Essentially compact schemes for unsteady viscous incompressible flows, J. Comput. Phys. 126 (1996) 122–138. [18] Q. Sheng, Adi, lod and modern decomposition methods for certain multiphysics applications, J. Algorithms Comput. Technol. 9 (2015) 105–120. [19] T.W.H. Sheu, L.W. Hsieh, C.F. Chen, Development of a three-point sixth-order Helmholtz scheme, J. Comput. Acoust. 16 (2008) 343–359. [20] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, second ed., SIAM, Philadelphia, 2004. [21] H. Brunner, A. Iserles, S.P. Nørsett, The spectral problem for a class of highly oscillatory Fredholm differential equations, IMA J. Numer. Anal. 30 (2010) 108–130. [22] G.D. Gillen, C.M. Seck, S. Guha, Analytical beam propagation model for clipped focused-gaussian beams using vector diffraction theory, Opt. Express 18 (2010) 4023–4040. [23] E. Larsson, A domain decomposition method for the helmholtz equation in a multilayer domain, SIAM J. Sci. Comput. 20 (1999) 1713–1731. [24] A. Iserles, A First Course in the Numerical Analysis of Differential Equations, second ed., Cambridge University Press, Cambridge, 2009. [25] Q. Sheng, S. Guha, L.P. Gonzalez, An exponential transformation based splitting method for fast computations of highly oscillatory solutions, J. Comput. Appl. Math. 235 (2011) 4455–4463.

Please cite this article in press as: T. Jones, Q. Sheng, Numerical stabilities study of a decomposed compact method for highly oscillatory Helmholtz equations, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.03.031.