Transparent boundary conditions for iterative high-order parabolic equations

Journal of Computational Physics 313 (2016) 144–158

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Transparent boundary conditions for iterative high-order parabolic equations P.S. Petrov a,b,∗ , M. Ehrhardt c a b c

Il’ichev Pacific Oceanological Institute, 43 Baltiyskaya str., Vladivostok, 690041, Russia Far Eastern Federal University, 8 Sukhanova str., 690950, Vladivostok, Russia Bergische Universität Wuppertal, Gaußstrasse 20, D-42119 Wuppertal, Germany

a r t i c l e

i n f o

Article history: Received 7 July 2015 Received in revised form 13 February 2016 Accepted 15 February 2016 Available online 18 February 2016 Keywords: Parabolic equation method Multiple-scale method Transparent boundary conditions

a b s t r a c t Recently a new approach to the construction of high-order parabolic approximations for the Helmholtz equation was developed. These approximations have the form of the system of iterative parabolic equations, where the solution of the n-th equation is used as an input term for the (n + 1)-th equation. In this study the transparent boundary conditions for such systems of coupled parabolic equations are derived. The existence and uniqueness of the solution of the initial boundary value problem for the system of iterative parabolic equations with the derived boundary conditions are proved. The well-posedness of this problem is also established and an unconditionally stable finite difference scheme for its solution is proposed. © 2016 Elsevier Inc. All rights reserved.

1. Introduction

The wide-angle parabolic equations (WAPEs) presently are considered a main computational tool for many problems of wave propagation [1,2]. The most important applications probably include geophysics [3], radiowave propagation problems [1], and acoustics [2], to mention only a few. The WAPEs are traditionally derived by means of the operator square root approximation with a Padé series (hereafter they are referred to as Padé WAPE). Recently another approach to the wideangle parabolic approximations was proposed in [4]. This derivation relies upon the systematic use of the multiple-scale expansion method, and the resulting high-order parabolic approximations have the form of the system of parabolic equations (PEs), where the input term of the n-th PE is obtained from the solution of (n − 1)-th PE [4] (hence we call them iterative PEs, and this reflects the nature of the numerical algorithm required to solve such system). These new parabolic approximations have some advantages over the classical Padé-style approach to the WAPE derivation, which are discussed in great detail in [4]. The most important feature of iterative PE approach is the clear and consistent derivation of interface and boundary conditions for the case when the interfaces and boundaries have complicated shape [4]. By contrast, in classical Padé WAPE theory it is still quite unclear which conditions should be used for the boundary conditions at the curved or sloping interfaces, excepting for a number of very special cases [5,6]. This shortcoming restricts one to the use of the staircase approximation of a complicated interface or boundary in various applications. This issue is especially important in the acoustical applications, where the interfaces between the media are ubiquitous [2].

*

Corresponding author at: Il’ichev Pacific Oceanological Institute, 43 Baltiyskaya str., Vladivostok, 690041, Russia. E-mail addresses: [email protected] (P.S. Petrov), [email protected] (M. Ehrhardt).

http://dx.doi.org/10.1016/j.jcp.2016.02.037 0021-9991/© 2016 Elsevier Inc. All rights reserved.

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

145

In order to solve numerically the practical problems of wave propagation on the unbounded domains using these new parabolic approximations from [4] one has to truncate the domain with artificial boundaries. For more than two decades many research efforts were made to develop the methods of the artificial domain truncation for the Schrödinger-type equations (including the acoustical parabolic equation and optical paraxial equation). This domain truncation may be accomplished either by imposing the transparent boundary conditions (TBC) or by extending the computational domain with so-called perfectly matching layers (PMLs). The theory of the TBCs for the Schrödinger-like equations originated from the papers of Baskakov and Popov [7] and Papadakis [8]. For a review of different approaches to TBCs and PMLs see [9] and the numerous references therein. There are also some works concerning the TBCs for the conventional Taylor and Padé WAPEs, e.g. TBCs for the high-order Taylor WAPEs were proposed in [10], the TBCs for the rational–linear WAPEs were derived in [11,12] while the case of general Padé WAPE is dealt with in [13–15]. In this paper we derive the TBCs for the parabolic approximations proposed in [4]. These conditions are basically a natural but non-trivial generalization of the classical Baskakov–Popov TBCs [7]. We also show that the initial-boundary value problems (IBVPs) for the coupled PEs constituting the system from [4] with the derived TBCs are well posed and that its solution coincides with the solution of the same system on the unbounded domain. The derived TBCs may be also used for the solution of the wide angle mode parabolic equations. We also propose a finite-difference scheme for the solution of the PEs system from [4] supplied with the derived TBCs. It is again a generalization of the numerical scheme of Baskakov and Popov [7]. In the interior of the computational domain the PEs are discretized using a second-order implicit Crank–Nicholson finite-difference method which is unconditionally stable for the unbounded domain or homogeneous Dirichlet conditions [9]. The incorporation of the TBCs into a numerical scheme may however render it only conditionally stable [16]. Sun and Wu [17] proved however that the Baskakov and Popov discretization of the TBCs leads to an unconditionally stable scheme. We adapted their proof for our new numerical scheme which also turns out to be unconditionally stable. 2. Wide angle parabolic approximations Let us consider the problem of sound propagation in a 2D acoustical waveguide  = {(x, z)| z ≥ 0} consisting of the water layer and one or more layers of bottom (which is assumed to be liquid), where z is the depth and x is the horizontal variable (here we use the acoustical notation following [4], although the same results may be reproduced for the open waveguides in optics and radiowave propagation theory). The acoustical pressure p (x, z) due to a point source located at x = 0, z = z s then satisfies the Helmholtz equation

∂ ∂x



1 ∂p



ρ ∂x

+

∂ ∂z



1 ∂p

ρ ∂z



+

1

ρ

κ2 p =

1

ρ

where the medium parameters are the density Dirichlet-type boundary condition

δ(x)δ( z − zs ) ,

(2.1)

ρ = ρ (x, z) and the wavenumber κ = κ (x, z) = ω2 /c 2 . A pressure-release

p (x, 0) = 0

(2.2)

is usually imposed at the ocean surface z = 0, while for sufficiently large values of depth (say for z ≥ L) the medium is assumed to be homogeneous, i.e. ρ (x, z) = ρb and κ (x, z) = κb for all z ≥ L. In order √ to correctly set a BVP for the equation (2.1) one also requires certain radiation conditions to be satisfied at infinity R = x2 + z2 → ∞ (see e.g. [2]). The propagation problems in the shallow-water acoustics usually feature additional complication associated with the presence of interfaces, i.e. surfaces z = H (x) where media parameters have finite jump discontinuities (e.g. water-bottom interface). The following coupling conditions are imposed at the interface z = H (x) (see e.g. [2]):

p |z= H (x)+0 = p |z= H (x)−0 ,

  1 ∂p   = . ρ ∂ n z= H (x)+0 ρ ∂ n z= H (x)−0 1 ∂p

(2.3)

The bottom relief described by the function z = H (x) in practical problems is often very complicated. Usually the BVP (2.1)–(2.2)–(2.3) is too complicated to be solved directly, and the high-order parabolic equations are used [2] to approximate the solution of the Helmholtz equation (2.1). While sacrificing some relatively unimportant propagation features, thus we obtain mathematical formulation of the problem which is much more efficient and easier for the numerical implementation. Among the shortcomings of the traditional approach to the PEs derivation (the one based on the operator square root approximation) is its inability to systematically account for the sloping and variable bottoms. Strictly speaking the very Helmholtz operator factorization which leads to the equations for the forward- and backward-propagating waves (containing the operator square root) relies upon the assumption of waveguide range-independence [2]. It is also well-known that even the simplest boundary and interface conditions for the case of the sloping bottom may lead to the ill-posedness of the IBVP for the parabolic equation [5]. Although some efforts we made to derive the proper interface and boundary conditions for the PEs in the case of the sloping bottoms [6,5], to our knowledge they were only partially successful. For example, in the paper [6] authors consider the case of the rational–linear wide-angle PE and the sloping pressure-release bottom, while for the general n-th order Padé WAPE and the arbitrarily sloping penetrable bottom no interface conditions were proposed

146

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

so far. Usually the piecewise-linear approximation of the bottom is used for the computations involving these high-order WAPEs [2]. At the same time, the consistent interface and boundary conditions may be successfully obtained within another approach to the WAPE derivation [4]. This approach is based on the multiple-scale asymptotic expansion of the acoustic pressure p (x, z). Here we outline the main results of this theory following closely [4]. We assume that the medium properties variation in x is much slower than in z and that κ = κ0 + ν (x, z) (here |ν (x, z)|  κ0 ). Under these assumptions the complex acoustical pressure p (x, z) may be approximated by the formula (see [4])

pn (x, z) = exp(iκ0 x)

j =n 

A j (x, z) ,

j =0

where

2i

1

ρ

κ0

∂Aj ∂ + ∂x ∂z



1 ∂Aj

ρ ∂z



      ∂ 1 ∂ A j −1 1 1 κ0 + ν A j + + i = 0, ρ ρ ∂x ρ ∂x x

(2.4)

and A −1 (x, z) = 0. System (2.4) consists of Schrödinger-type equations, and the input term of the j-th equation is determined from the solution of the ( j − 1)-th equation. Hence the whole system may be solved iteratively. The quantity pn (x, z) is called hereafter the wide-angle parabolic approximation of the order n to the solution p (x, z) of the equation (2.1). It requires the solution of n + 1 parabolic equations (2.4) to be computed. It is interesting to note that the system of PEs (2.4) for the case of the Gaussian beams propagation in the free space was derived earlier in the paper of Grikurov and Kiselev [18], similar approach was also outlined in the early (unpublished) works of Malyuzhinets. Note that the solutions to (2.4) are sought as usual [19] in the space C ([0, xmax ], L 2 ([0, ∞))). To approximate the point-source solution of the Helmholtz equation, i.e. the solution of (2.1) with the right-hand side replaced by δ(x, z − z s ), we use the following Cauchy initial conditions A j (0, z) for the system (2.4):

A 0 (0, z) = S ( z) , A j (0, z) = 0,

j = 1, 2, . . . ,

(2.5)

where S ( z) is any starter from standard PE theory (a nice collection of starters can be found in [2], see also examples in Section 6). We also note that the same asymptotic expansion as was used to derive (2.4) allows to obtain the consistent interface conditions for the parabolic equations (2.4). This feature of our approach allows to establish some important results on the convergence of our approximations and the energy flux conservation (see extensive discussion in [4]). 3. Transparent boundary conditions for the system of PE The modeling of sound propagation in the typical problems of underwater acoustics requires the solution of the PEs system (2.4) in the domain  = {(x, z)| z ≥ 0, 0 ≤ x ≤ xmax }, where the field vanishes at the ocean surface z = 0 and the medium is assumed to have constant acoustical properties for sufficiently large values of depth z (i.e. there exists L such that ρ (x, z) = ρb , κ (x, z) = κb and ν (x, z) = νb for z ≥ L). For computing the numerical solution of (2.4) it is therefore natural to truncate the computational domain at z = L introducing a TBC at this fictitious boundary. Note that we also require that all interfaces (2.3) where media discontinuities occur are located above this boundary (i.e. we require H (x) < L for all x). Without any loss of generality we consider in this paper the following problem. Let A = ( A 0 (x, z), A 1 (x, z), . . . , A n (x, z)) be a solution to the reference IBVP for the system

2iκ0 A j ,x + A j ,zz + ν A j + A j −1,xx = 0 , A 0 (0, z) = S ( z) ,

A j (0, z) = 0 ,

j = 1, 2, . . . ,

A j (x, 0) = 0 , lim | A j (x, z)| = 0

(3.1)

z→∞

(here the subscript indices x, z denote partial derivatives) in the domain  with the initial conditions S ( z) compactly supported on [0, L ]. We seek to construct the artificial boundary conditions of the form

B( A j ) = 0 ,

(3.2)

for (3.1) at z = L such that the solution A = to the IBVP for the system (3.1) on the truncated domain t = {(x, z)|0 ≤ z ≤ L , 0 ≤ x ≤ xmax } with initial conditions (2.5) and boundary conditions (3.2) at z = L coincides t

( A t0 (x, z),

A t1 (x, z), . . . ,

A nt (x, z))

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

147

with the solution of the reference IBVP. Recall that conditions (3.2) ensuring At = A for all (x, z) ∈ t are called transparent boundary conditions (TBCs) for (3.1). Note that we use the simplified form (3.1) of the equation (2.4) since in the halfspace z ≥ 0 ρ and κ are assumed to be independent of z, x. From the TBC derivation in the next section it will be clear that it can be used without any changes for general case of the system (2.4). We assume also that the initial condition S ( z) vanishes for all z ≥ L. 3.1. Construction of the TBCs The reference IBVP in the halfspace z ≥ 0 is obviously equivalent to the two coupled systems of IBVPs:

and

⎧ 2iκ0 A tj ,x + A tj ,zz + ν A tj + A tj −1,xx = 0 , (x, z) ∈ t , ⎪ ⎪ ⎪ t ⎨ A 0 (0, z) = S ( z) , A tj (0, z) = 0, j = 1, 2, . . . , ⎪ A tj (x, 0) = 0, j = 0, 1, . . . , ⎪ ⎪ ⎩ t A j ,z (x, L ) = A j ,z (x, L ) ,

(3.3)

⎧ 2iκ0 A rj ,x + A rj ,zz + νb A rj + A rj −1,xx = 0 , ⎪ ⎪ ⎪ ⎨ A r (0, z) = 0, j = 0, 1, . . . , j t r ⎪ A ⎪ j (x, L ) = A j (x, L ) , ⎪ ⎩ limz→∞ | A rj (x, z)| = 0 .

(3.4)

(x, z) ∈ r ,

Here (3.3) is the problem on the truncated domain t and (3.4) is the bottom exterior problem on the halfspace r = [0, xmax ] × [ L , ∞) (in the works related to the Schrödinger equations (3.4) is usually called right exterior problem [9]). For the given boundary input A tj (x, L ) the bottom exterior problem may be solved explicitly and A tj ,z (x, L ) is then computed from the solution. Thus we consider the solution of (3.4) as an operator relating A tj (x, L ) to A tj ,z (x, L ). This relation is called Dirichlet-to-Neumann map, it could be used as a boundary condition for (3.3). Note that according to (3.3), (3.4) there are no density discontinuities at z = L. The technique allowing to account for the density jumps is described in [11,12]. For a given function A tj (x, L ) the problem (3.4) may be solved explicitly. We now apply the Laplace transform L : f (x) →

ˆf (ξ ) to all equations and BCs in (3.4) and arrive at the following BVP:

⎧ ˆ r + (2iκ0 ξ + νb ) Aˆ r = −ξ 2 Aˆ r , ⎪ A ⎪ j j −1 ⎨ j ,zz ˆ r (ξ, L ) = Aˆ t (ξ, L ) , A j j ⎪ ⎪ ⎩ lim | Aˆ r (ξ, z)| = 0 . z→∞

z ∈ [ L , ∞) , (3.5)

j

ˆ r (ξ, z) = u j (t ), 2iκ0 ξ + νb = − w 2 , ξ 2 = v, Aˆ t (ξ, L ) = a j in (3.5), we rewrite Introducing a new variable t = z − L and setting A j j the BVP (3.5) as

⎧  2 ⎪ ⎨ u j − w u j = − vu j −1 , u j (0) = a j , ⎪ ⎩ limt →∞ |u j | = 0 .

t ∈ [0, ∞) , (3.6)



Note that v , w are independent of t, and w = + −2iκ0 ξ − νb denotes the branch of the square root with the positive real part. In (3.6) w is considered a parameter. This system of BVPs may be easily solved using the variation of parameters, and its solution has the general form



u j (t ) = e− wt a j + a j −1 v P 1 (t , w ) + a j −2 v 2 P 2 (t , w ) + · · · + a0 P j (t , w ) , where P k (t , w ) for k ≥ 1 are polynomials in t such that the function qk (t ) = homogeneous BCs, i.e.:

e− wt P

(3.7) k (t , w )

is a solution to the BVP with

⎧  2 − wt P ⎪ k−1 (t , w ) , ⎨ qk − w qk = −e qk (0) = 0 , ⎪ ⎩ limt →∞ |q| = 0 .

Substituting qk (t ) into the latter BVP, we obtain a BVP for P k (t , w )

⎧   ⎪ ⎨ P k − 2w P k = − P k−1 , P k (0, w ) = 0 , ⎪ ⎩ limt →∞ | P k (t , w )e− wt | = 0 .

(3.8)

148

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

Note that w is considered a parameter in P k (t , w ), and primes in (3.8) denote the derivatives of the polynomials with respect to t. It is easy to verify explicitly that the BVPs (3.8) have solutions of the form

P k (t , w ) =

k 

αk, j

j =1

tj w 2k− j

(3.9)

,

¯ k+1 = {αk+1, j } may be computed from the set α¯ k by solving the following system of linear where the set of coefficients α algebraic equations:

⎧ ⎪ ⎨ −2(k + 1)αk+1,k+1 = −αk,k , −2( j + 1)( j + 2)αk+1, j +2 − 2( j + 1)αk+1, j +1 = −αk, j , ⎪ ⎩ 2αk+1,2 − 2αk+1,1 = 0 .

j = 1, . . . , k − 1 ,

The solution of the system (3.10) by Gaussian elimination is an easy task; it requires the recursive evaluation of αk+1, j+1 by the formulae

⎧ 1 ⎪ ⎨ αk+1,k+1 = 2(k+1) αk,k , αk+1, j = j+2 1 αk+1, j+1 + ⎪ ⎩ αk+1,1 = αk+1,2 .

1 2j

j = 2, . . . , k ,

αk−1, j ,

(3.10)

αk+1, j from

(3.11)

Now we have all necessary ingredients to solve the BVP (3.5). In the Laplace domain (ξ, z) we have the following ˆ r in the halfspace z ≥ 0: expression for A j

ˆ r (ξ, z) A j

=e

⎛ j 

− w (ξ )( z− L ) ⎝

⎞

2k

ξ P k (z −

k =0

ˆ t (ξ, L )⎠ L , w (ξ )) A j −k

,

where we write P k ( z − L , w (ξ )) instead of P k ( z − L , w ) to stress the dependence of the polynomial coefficients on ξ . To obtain a Dirichlet-to-Neumann (DtN) condition, we differentiate the last equation with respect to z

∂ Aˆ rj (ξ, z)

⎛ j  = e− w ( z − L ) ⎝ ξ 2k ( P  ( z − L , w ) − w P k ( z − L , w )) Aˆ t

⎠.

j −k (ξ, L )

k

∂z

⎞

k =0

(3.12)

Let us recall that the coupling condition in the IBVP (3.3) reads as A tj ,z (x, L ) = A rj ,z (x, L ). Next we substitute the expression for A rj ,z (ξ, z) from (3.12) into this condition. Observing that

⎛ ⎞  j   w (ξ )e− w (ξ )(z− L ) ⎝ ξ 2k P k ( z − L , w (ξ )) Aˆ tj −k (ξ, L )⎠  k =0

and j 

2k



ξ P k (z −

k =1

=

j 

 ξ

2k

k =1

=

j 

 ξ

2k

=

j 

ξ 2k

k =1

z= L

  

ˆ t (ξ, L ) L, w ) A j −k  

z= L

k ( z − L )m ∂  αk,m 2k−m ∂z w m =1

αk,1 w 2k−1

k =1

= w (ξ ) Aˆ tj (ξ, L ) ,

+

k  m =2

    

ˆ t (ξ, L ) A j −k z= L

 ( z − L )m  mαk,m  w 2k−m 

ˆ t (ξ, L ) A j −k z= L

αk,1 ˆ t A j −k (ξ, L ) ,

w 2k−1

we obtain the following DtN TBC in the Laplace domain:

 ∂ Aˆ tj (ξ, z)    ∂z

= − w (ξ ) Aˆ tj (ξ, L ) + z= L

j  k =1

ξ 2k

αk,1 ˆ t A j −k (ξ, L ) .

w 2k−1

(3.13)

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

149

For the sake of completeness we note that an analogous condition for the left artificial boundary at z = 0 reads

 ∂ Aˆ tj (ξ, z)    ∂z

= w (ξ ) Aˆ tj (ξ, 0) −

j 

ξ 2k

k =1

z =0

αk,1 w (ξ )2k−1

ˆ t (ξ, 0) . A j −k

(3.14)

In order to obtain the TBCs in the physical domain we now apply the inverse Laplace transform L−1 : ˆf (ξ ) → f (x) to the conditions (3.13)–(3.14) using the well-known properties of L−1 :

d2k f

L−1 (ξ 2k ˆf (ξ )) =

,

dx2k

L−1 ( gˆ (ξ − a) ˆf (ξ )) = eax

x

g (x − y )e−ay f ( y )dy = eax L−1 ( gˆ (ξ )L(e−ax f (x))) ,

0

L −1

√ +

 x  y yk y 2 ξˆ d dy f (ξ ) = . . . f ( y 1 )dy 1 dy 2 . . . dyk , √ dx x− y ξk 0

0

0

0

and observing that



π

1 w (ξ )2k−1

=

ei 4 (2k−1) 1

(2κ0 )k− 2

ξ − 2iνκb0

k . ξ − 2iνκb0 +

Thus we straightforwardly arrive at the following TBCs

∂ A tj (x, z) ∂n

 =−

2κ0 −i π i 2νκb x d e 4e 0 π dx

x √ 0

 j

−

−k

y 2

0

ν

e

... 0

A tj ( y , z)e

x− y

 y yk

αk,1 (−2iκ0 )

k =1



dy

−i 2κb y 1 0

ν

−i 2κb y 0

∂ 2k A tj −k ( y 1 , z) ∂ y 2k 1

0

 dy 1 dy 2 . . . dyk

(3.15)

at z = L and z = 0 (n denotes the outward unit normal vector at z = L, z = 0 respectively). It is reasonable to tabulate the first coefficients αk,1 :

α1,1 =

1 2

,

α2,1 =

1 8

,

α3,1 =

1 16

,

α4,1 =

5 128

,

α5,1 =

7 256

,...

since in practice it is enough to keep only the first few terms of the sum in (3.15). The TBCs (3.15) simplify significantly if we choose the reference wavenumber κ0 in such a way that νb = 0 (this condition is fulfilled if we set κ0 = κb in (2.4)). Under this assumption the multiple integral on the right-hand side of (3.15) vanishes, and the TBC becomes

∂ A tj (x, z) ∂n

 =−

 ∂ k A t ( y , z) j 2 κ0 − i π  d dy j −k e 4 αk,1 (−2iκ0 )−k , √ k π dx x− y ∂y x

k =0

0

at z = 0, L ,

(3.16)

where α0,1 = −1. Note that the TBCs (3.15) are a natural generalization of the TBC for the single parabolic equation (see e.g. [9]), while its simplified form (3.16) is a generalization of the TBC for the PE with the potential vanishing at the artificial boundary [7,17]. 4. Well-posedness of the IBVP and the uniqueness of its solution In this section we first prove the uniqueness of the solution of the IBVP supplied with the TBC (3.15) and then show that the resulting problems for A tj (x, z) are well-posed. For the sake of simplicity we restrict our attention to the problem with the right transparent boundary condition:

150

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

2iκ0 A tj ,x + A tj ,zz + ν A tj + A tj −1,xx = 0 , (x, z) ∈ t , A t0 (0, z) = S ( z) ,

A tj (0, z) = 0, j = 1, 2, . . . ,

A tj (x, 0) = 0, j = 0, 1, . . . ,

∂ A tj (x, z) ∂z



=−

j 2 κ0 − i π  e 4

π

αk,1

d

(−2i)k 0k k =0

κ dx

 x ∂ k A t ( y , z) j −k ∂ yk 0

√

dy x− y

at z = L .

(4.1)

4.1. Existence and uniqueness The uniqueness proof relies on the following standard result for the homogeneous Schrödinger-type equation [12] Lemma 4.1. (See [12].) An IBVP

2iκ0 B x + B zz + ν B = 0 , (x, z) ∈ t , B (0, z) = S ( z) ,

S ( z) = 0 for z ≥ L ,

B (x, 0) = 0 ,

 x ∂ B (x, z) 2 κ0 − i π d B ( y , z)dy =− e 4 α0,1 at z = L √ ∂z π dx x− y

(4.2)

0

for the standard PE with the TBC at z = L and compactly supported initial data S ( z) ∈ H 2 [0, L ] has unique solution B (x, z) ∈ L 2 [0, L ]. This lemma immediately implies the following result: Proposition 4.1. If there exist two solutions A 1j (x, z) and A 2j (x, z) of the IBVP (4.1) then A 1j (x, z) = A 2j (x, z) for all (x, y ) ∈ t . Proof. If A 1j (x, z) and A 2j (x, z) satisfy (4.1) then A 1j (x, z) − A 2j (x, z) satisfies IBVP (4.2) with S ( z) = 0 for all z. From the uniqueness of the solution of IBVP (4.2) it follows that A 1j (x, z) − A 2j (x, z) ≡ 0.

2

To prove the existence of the solution of IBVP (4.1) and its well-posedness we first consider a general inhomogeneous IBVP of the form

2iκ0 A x + A zz + ν A + ζ = 0 , A (0, z) = 0 , A (x, 0) = 0 , lim | A (x, z)| = 0 .

z→∞

(4.3)

The following result may be then established [19,20]: Lemma 4.2. If the input term ζ (x, z) ∈ C ([0, xmax ], H 2 (R)) and the potential ν (x, z) ∈ C ([0, xmax ], H 2 (R)) and ν (x, z) = 0 for z ≥ L then the IBVP (4.3) has a unique classical solution

A (x, z) ∈ C ([0, xmax ], H 2 (R)) ∩ C 1 ([0, xmax ], L 2 (R)) . Note that S.-i. Doi [20] proved that the solution also belongs to C ([0, xmax ], H k (R)) for sufficiently smooth initial data and potential. Hereafter we assume that these conditions are fulfilled, i.e. that ζ (x, z) ∈ C ([0, xmax ], C M (R)) and ν (x, z) ∈ C ([0, xmax ], C M (R)) for sufficiently large M such that A jxx (x, z) in (3.1) belongs to C ([0, xmax ], H 2 (R)) for all j = 1, 2, . . . , n. Although this requirement is somewhat too restrictive, the more tedious formulation of necessary conditions for S ( z) and ν (x, z) would lead us too far from the subject of this study. The function S (z) in most PE starters is in C ∞ (R), but the potential may actually be piecewise smooth in z (usually piecewise linear), and the well-posedness proof given here does not cover this case. Summarizing the arguments above we formulate Lemma 4.3. For ν (x, z) ∈ C ([0, xmax ], C M (R)) and S ( z) ∈ C M (R) where M is sufficiently large the IBVP (3.1) has a unique classical solution A = ( A 0 (x, z), A 1 (x, z), . . . , A n (x, z)) where A j (x, z) ∈ C ([0, xmax ], H 2 (R)).

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

151

This lemma immediately implies Proposition 4.2. If the conditions of Lemma 4.3 are fulfilled and additionally ν (x, z) = 0 and S ( z) = 0 for all z ≥ L then there exist a classical solution of the IBVP (4.1) (which is unique according to Proposition 4.1). This solution coincides with the solution of (3.1) for all (x, z) ∈ t . Proof. By construction of the conditions (3.16), any solution A j (x, z) of the halfspace problem (3.1) satisfies the TBC (3.16). Thus the solution from the Lemma 4.3 also solves (4.1). 2 4.2. Well-posedness Since the solution of the IBVP with the TBC (4.1) is merely a restriction of the halfspace problem (3.1) solution to the truncated domain t , the well-posedness of (4.1) also follows from the well-posedness of the IBVP (3.1). The latter fact is established by the following proposition Proposition 4.3. For the solution of the IBVP (3.1) under the assumptions of Lemma 4.3 the following inequality holds for j = 1, 2, . . .

x 2κ0 ( N (x) − N (0)) ≤

η(ξ )dξ ,

(4.4)

0

where

N (x) =  A j (x, z) L 2 z

⎞1/2 ⎛∞  2 = ⎝ | A j (x, z)| dz⎠ , 0

η(x) = ∂ 2 A j−1 (x, z)/∂ x2 L 2z

⎞1/2 ⎛∞  = ⎝ |∂ 2 A j −1 (x, z)/∂ x2 |2 dz⎠ . 0

Proof. For convenience we use the notation from (4.3): A (x, z) = A j (x, z) and ζ (x, z) = ∂ 2 A j −1 (x, z)/∂ x2 . Let us recast the parabolic equation from (4.3) in the form

2κ0 A x = i A zz + iν A + iζ , multiply it by A ∗ (the star denotes complex conjugation) and take the real part of the resulting expression. After obvious transformations we obtain

κ0

∂ | A |2 = (i A ∗ A zz ) + (iν | A |2 ) + (i A ∗ ζ ) . ∂x

Now let us integrate the last equation with respect to z:

κ

d

 A 2L 2 0 z dx

⎞ ⎞ ⎛∞ ⎛ ∞  ∞  ∗ 2 ∗ =  ⎝ i A A zz dz⎠ + (ν )| A | dz +  ⎝i A ζ dz⎠ . 0

0

0

Firstly using integration by parts we find that

⎞ ⎛ ⎞ ⎛∞  ∞  ∞  ⎝ i A ∗ A zz dz⎠ =  ⎝ i A ∗ A z  − i A ∗z A z dz⎠ 0

0

⎛ =  ⎝i

∞

⎞

0

| A z |2 dz⎠ = 0 .

0

Secondly we observe that (ν ) = 0 for the medium with no attenuation, while for the lossy one (ν ) < 0. We therefore conclude that

⎛

κ0

d dx

 A 2L 2 ≤  ⎝i

∞

⎞

A ∗ ζ dz⎠ ≤  A  L 2 ζ  L 2 . z

z

0

z

152

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

This inequality may be rewritten using our notation as

2κ0 N  (x) N (x) ≤ N (x)η(x) , and this concludes the proof of the proposition

x 2κ0 ( N (x) − N (0)) = 2κ0

x



N (ξ )dξ ≤ 0

2

η(ξ )dξ . 0

5. Numerical scheme for the solution of the IBVP and its stability In this section we construct a finite difference numerical scheme for the solution of the IBVP (4.1). It is a simple generalization of the scheme devised by Sun and Wu [17] which in turn reduces to the numerical method of Baskakov and Popov [7]. This scheme enjoys a remarkable property of the unconditional stability, and this fact was first proven in [17]. Our scheme also inherits this property as we will show here. We introduce a uniform grid xn = nx, zm = m z, n = 0, 1, . . . , N, m = 0, 1, . . . , M, where  zM = L, xN = xmax . Hereafter upper indices are used to denote values of any given function U (x, z) in the grid points:

1

U n+1/2,m ≡

U n,m ≡ U (xn , zm ) ,

2



U (xn+1 , zm ) + U (xn , zm ) .

For the parabolic equations (2.4) in the problem (4.1) we use the Crank–Nicholson finite difference scheme modified to include the input term: n+1,m

2iκ0

Aj

− Anj ,m

x

n+1/2,m+1

Aj

+

n+1/2,m

− 2A j 

n+1/2,m−1

+ Aj

n+1/2,m

+ ν n+1/2,m A j

 z2

+ ( A j −1,xx )n+1/2,m = 0 ,

(5.1)



where the input term ( A j −1,xx )n+1/2,m = ( A j −1,xx )n+1,m + ( A j −1,xx )n,m /2 is computed from the solution of the equation for A j −1 using the standard 3-point finite difference approximation n+1,m

( A j −1,xx )n,m ≈

A j −1

,m −1,m − 2 Anj − + Anj − 1 1

.

x2

n ,0

n ,1

n, M

The scheme (5.1) is marching in x, i.e. it computes the vector Anj = ( A j , A j , . . . , A j condition (4.1) at x = 0 yields the starter 0,m

A0

= S ( zm ) ,

0,m

Aj

= 0,

) from the vector Anj −1 . The initial

A0j :

j = 1, 2, . . . .

(5.2)

The Dirichlet initial condition A j (x, 0) = 0 at z = 0 turns into n,0

Aj

= 0,

j = 0, 1, 2, . . . .

(5.3)

We now describe the approximation of the TBC (3.16) in (4.1). Let us first recast the convolution terms of the TBC in the following form

 x ∂ k A t ( y , z) j −k

d

∂ yk

dx 0

√

dy x− y

=

dk

 x ∂ A t ( y , z) j −k

dxk

∂y

√

0

dy x− y

.

For the latter integral we use the approximation developed by Baskakov and Popov [7]: (n+ 1)x

0

where

 ∂ A j ( y , L ) dy ≈ γk Anj +1−k,M , √ ∂y x− y n

(5.4)

k =0

√ √

γ0 = 2/ x, while for k = 1, 2, . . .

γk = √

x

√

k+1+

√  √ k

−2

k−1+

√  √

To further simplify the notation, we define βk = in the BC at z = L:

k



k+1+

2κ0 −i π4

π e

√

.

k−1

γk and introduce the following approximation for the input terms

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158



x(n+1) k t ∂ A

j 2 κ0 − i π  d 4 e αk,1 (−2iκ0 )−k π dx

j −k

∂ yk

k =1



0

j

≈

( y , z)

αk,1 (−2iκ0 )−k ( D x )k

n 

√

153

dy x− y

βk Anj +1−s, M ,

(5.5)

s =0

k =1

where ( D x )k is the k-th degree of the central difference operator D x defined by

D x U n,m =

U n+1,m − U n−1,m 2 x

for any function U n,m = U (xn , zm ) on the grid xn , zm . The right-hand side of (5.5) is denoted hereafter as W nj . The discretized form of the TBC (3.16) from the IBVP (4.1) can be written as

1 2 z

n+1/2, M +1

Aj

n+1/2, M −1

− Aj



=−

n 

n+1/2−s, M

βs A j

n+1/2

−Wj

n+1/2, M +1

Expressing the quantity A j n +1 , M

(5.6)

Aj

− Anj , M

from (5.6) and substituting it into the formula (5.1) (for m = M), we finally arrive at

n+1/2, M −1

Aj

n+1/2, M

− Aj

n+1/2, M

+ ν n+1/2, M A j  x  z2  n  2  n+1/2−s, M n+1/2 − βs A j +Wj + ( A j −1,xx )n+1/2, M = 0 , z

2iκ0

.

k =0

+2

(5.7)

k =0

Equations (5.1), (5.7) and (5.3) constitute together a linear system which allows to obtain Anj +1 from Anj . Sun and Wu [17] have proved the following result for the numerical scheme (5.1)–(5.7)–(5.3) with no input (i.e. for the case of A j −1,xx = 0 and W j = 0): Proposition 5.1. In the case of the IBVP (4.2) the numerical scheme (5.1)–(5.7)–(5.3) is unconditionally stable, i.e. the following inequality holds for all n:

Bn 2 ≤ B0 2 , where Bn = ( B n,1 , . . . , B n, M ), and  · 2 denotes the discrete 2 -norm:

Bn 22 =

M 

| B n,m |2 .

m =0

It is easy to see that in the case of the equation for A j (xn , zm ) the input terms may be considered as additional initial conditions (introduced at x = x1 , x = x2 , . . . ). The resulting solutions for all of these additional initial conditions are combined in the inhomogeneous equation for A j and by Proposition 5.1 we conclude the following: Proposition 5.2. The numerical solution of the IBVP (4.1) obtained using the numerical scheme (5.1)–(5.7)–(5.3) is unconditionally stable, and the following inequality is satisfied for any N = 1, 2, . . . :

BN 2 ≤ B0 2 + C 1

N  n =0

| W nj |2 + C 2

N  M 

|( A j −1,xx )n, M |2 ,

n =0 m =0

where the constants C 1 and C 2 are independent of W j , A j −1,xx and N. The proof of this result obviously follows from Proposition 5.1. 6. Test examples In this section we conduct some numerical experiments in order to test the derived TBCs for the equations (2.4). Since the performance of the TBCs does not depend on the inhomogeneities in the interior of the computational domain, we start with the simplest case of sound propagation in the homogeneous medium. In this case an analytical solution is available, and the accuracy of the derived TBCs can be directly estimated. The in the second example we deal with more practical 2D wedge problem which is usually used for benchmarking various PE models.

154

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

Fig. 1. 2 -errors E Nj computed for different mesh sizes x for the solutions A j (x, y ), j = 0, 1, 2.

6.1. Homogeneous halfspace Let us consider a halfspace z ≥ 0 with the parameters c (x, z) = c and ρ (x, z) = 1. The source of the frequency f is located at x = 0, z = z s (hence κ0 = 2π f /c). We set up an artificial boundary at z = L = 2z s and solve the problem of sound propagation using the system of the IBVPs (4.1) and the finite difference scheme (5.1)–(5.7)–(5.3). In this example we assigned the following values to the parameters: c = 1500 m/s, z s = 100 m, f = 100 Hz. We use the Gaussian initial condition [2] of the form

¯ A 0 (0, z) = S ( z) = Ae

− (z−z2s )

2

(6.1)

σ

in the IBVP (4.1) for A 0 in order to simulate an acoustical field produced by the point source. It easy to see that for such initial conditions all the equations of the system (3.1) in the homogeneous halfspace z ≥ 0 may be solved analytically using the Fourier transform in z (cf. [2]). The solution has a very simple form, e.g. A 0A and A 1A (subscript A for “analytical”) writes as

A 0A (x, z) = U 0 (x, z − z s ) − U 0 (x, z + z s ) , A 1A (x, z) = U 1 (x, z − z s ) − U 1 (x, z + z s ) ,

(6.2)

where U 0 (x, z), U 1 (x, z) are the solutions for the unbounded media (i.e. for −∞ ≤ z ≤ ∞) with the source located at z = 0:



¯ U 0 (x, z) = A U 1 (x, z) = ix





σ 2 κ0 κ0 z 2 − , exp σ 2 κ0 + 2ix σ 2 κ0 + 2ix

6x2 − 6iσ 2 κ0 x + 12iκ0 z2 x − 2κ02 z4 + 6σ 2 κ02 z2 −

κ0 (σ 2 κ0 + 2ix)4

3 2

σ 4 κ02

U 0 (x, z) .

(6.3)

The solutions of higher-order equations have similar form and may be obtained in the same straightforward way. In the book [2], it is shown that the parameters in the initial condition (6.1) must be chosen as

¯= A

√

√

κ0 , σ =

2

κ0

.

We may now compare the numerical solutions A 0 , A 1 , . . . of the IBVP (4.1) against the analytical solutions A 0A , A 1A , . . . (6.2) of the halfspace problem (3.1). The numerical solutions were computed for different pairs of xq ,  zq , where

xq =

1 100.2q

x0 ,

 zq =

1 100.1q

 z0 ,

and  z0 = 1 m, x0 = 1 m. The ratio ( zq )2 /xq = ( z0 )2 /x0 is therefore kept constant (this is standard practice for finite difference schemes for parabolic equations). In order to estimate the quality of approximation, we compute E Nj , the 2 -error at xmax = N x divided by the 2 -norm of the exact solution in the domain [0, L ]:

E Nj =

ANj − A Aj , N 2 A Aj , N 2

,

(6.4)

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

155

Fig. 2. Contour plots of log10 | A j (x, z)| for j = 0 (top), j = 1 (middle), j = 2 (bottom). The numerical solutions of the IBVP (4.1) are shown on the right, while their analytical counterparts are presented on the left. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

A ,n

where A j

      = A Aj xn , z0 , A Aj xn , z0 , . . . , A Aj xn , z0 . Note that since our numerical scheme is a direct generalization of

the Baskakov–Popov scheme, the error estimate of Sun and Wu [17] can be easily extended to our case. More precisely, we have

E Nj ≤ C | z3/2 + x3/2  z−1/2 | .

(6.5)

For constant  z2 /x we obtain an estimate E Nj = O ( z3/2 ) = O (x3/4 ). Our calculations confirm that all the numerical solutions A j converge steadily to their analytical counterparts, and the rate of convergence is the same as predicted by (6.5) (i.e. the error is approximately proportional to x3/4 ). This convergence is illustrated in the Fig. 1, where errors E Nj are plotted against the mesh sizes  zq for j = 0, 1, 2. The computed solutions A 0 , A 1 , A 2 and their analytical counterparts are shown in Fig. 2 in the form of the log-contour plots (i.e. we plotted the contours 0, −1, −2, . . . for the function log10 | A j (x, y )|). These log-contours allow us to detect very small levels of the reflection. Note that although the level of accuracy decreases as j increases, this fact is related mostly to the numerical differentiation on the right-hand side of (2.4). In the paper [4] we suppressed these errors using the Lanczos low-noise differentiation formulae. Here we did not use this technique, since the proper TBCs (3.15) allowed us to use standard 3-point central difference formula.

156

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

Fig. 3. Transmisson losses TL(x) for z = 30 m (left) and z = 150 m (right) computed using the PE system (3.1) and TBC (3.16) (solid red line) and COUPLE 2-way model (black dash line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Contour plot of the transmission loss distribution TL(x, z) for the wedge problem computed using the fifth-order PE approximation with TBCs (3.16). Note the source position S and the sloped water-bottom interface. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

6.2. Sound propagation in a penetrable wedge Now we consider the problem of upslope sound propagation in a wedge-shaped ocean. This range-dependent benchmark problem is usually used for the validation of PE models (see discussion in [2]). The wedge is comprised by the water column overlying a fluid bottom. The sound speed in the water is c w = 1500 m/s and the density ρ w = 1 g/cm3 , while the respective parameters in the bottom are cb = 1700 m/s and ρb = 1.5 g/cm3 . The initial water depth at the source location is 200 m, it is decreasing linearly to 0 at a range of 4 km (see also Fig. 4), hence the bottom slope is 2.86◦ . The source depth is 100 m and the source frequency 25 Hz. In this example we use Thomson halfspace source [21] as the initial condition S ( z) in the IBVP (4.1). This PE starter is suitable to handle very wide-angle propagation accurately [2]. The computational domain is truncated at the depth z = L = 240 m where the TBC (3.16) is imposed (the truncation depth is chosen in such a way that S ( z) has a local minimum near z = L). For the solution of this problem a uniform grid with  z = 0.25 m and x = 5 m is used. Here we omit the details of numerical implementation of the interface conditions at the water-bottom interface (they are described in [4]). The solutions of IBVP (4.1) √ are computed for j = 0, 1, 2, 3, 4, and the sound pressure p 4 (x, z) is obtained by the summation of A j (x, y ) divided by x (to account for cylindrical spreading). The resulting solution is shown in Fig. 3 and Fig. 4. Usually the quality of the PE solution is estimated by the comparison of transmission loss plots TL(x) (see [2] for definition) with the COUPLE 2-way normal mode solution at z = 30 m and z = 150 m (which is highly accurate, see [2] for details). From both plots in Fig. 3 we can see that the agreement of the two solutions is very good, and no spurious reflection at the artificial boundary can be observed in Fig. 4. Note that wide-angle Thomson starter is essential to obtain an accurate solution near the source. We also point out that in the paper [4] wedge benchmark was also considered. The lack of the proper TBCs for the system (3.1) however did not allow us to achieve acceptable accuracy of the solution at z = 150 m in the previous work. The accuracy for z = 30 m is also notably improved by using the derived TBCs. We would like to stress that the use of TBCs (3.16) in principle allows to truncate the computational domain just below the last media interface leading to a highly efficient numerical scheme without any unnecessary domain extension.

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

157

7. Conclusion and future work In this work the TBCs for the system of the parabolic equations from [4] were derived. The solutions of the equations (2.4) form the high-order parabolic approximation to the solution of the Helmholtz equation describing the propagation of acoustic waves. Therefore, in many practical situations TBCs are required in order to solve them numerically on the artificially truncated domains. The TBCs were derived by the same method based on the Laplace transform as was used in several classical papers (see [9]). We considered the IBVP (4.1) for the system of PEs (2.4) with the presented TBCs and proved the existence and uniqueness of its solution. Also the well-posedness of the latter IBVP was established. We developed a finite difference scheme for the numerical solution of this IBVP with the derived TBC. This scheme, being a natural generalization of the one proposed by Baskakov and Popov, is shown to be unconditionally stable. The efficiency of the presented scheme and the convergence of the computed solutions to the analytical ones (as  z → ∞) is confirmed by a numerical example. It is important to note that since we must differentiate the numerical solution A j (x, z) with respect to x in order to obtain an input term in the equation for A j +1 (x, z), even small reflections at the fictitious boundary may destroy the method’s order of accuracy (small errors are inflated by differentiation). This is the reason why it is necessary to make these undesirable reflections as small as possible. One possible way to almost completely suppress them is to use the semi-discrete or fully-discrete TBCs consistent with the Crank–Nicolson scheme. Such conditions for the usual parabolic equation (or Schrödinger equation) were derived in the works [22,12], cf. the review article [9]. In our future work we will try to generalize these TBCs to the case of the system (2.4) and develop the corresponding stability proofs. The derivation of the approximate microdifferential TBCs for the system WAPA (2.4) is another development line for the present study. For instance, a direct generalization of the approach of Antoine and Besse [23] seems to be quite feasible here. Next, we will generalize the equations (2.4) to the 3D case, and investigate the TBCs for this case using the methods of Feshchenko and Popov [24] or Schädle [25]. Furthermore, it is very promising to extend our approach (both the system (2.4) and the TBCs (3.16)) to the case of the nonlinear Helmholtz equation in the Kerr medium [26]. It is easy to see that the resulting iterative PEs have both coefficients and input terms depending on the solutions of the previous equations. Thus our iterative PEs technique can handle both the nonlinearity and the wide-angle propagation in the uniform way. Acknowledgements The reported study was accomplished during P.S. Petrov’s visit to the Bergische Universität Wuppertal under the DAAD program “Forschungsaufenthalte für Hochschullehrer und Wissenschaftler”. P.S. Petrov was also supported by the RFFI foundation under the contracts No. 14-05-3148614_mol_a and No. 16-05-01074_a, the POI FEBRAS Program “Nonlinear dynamical processes in the ocean and atmosphere”, and the RF President Grant MK-4323.2015.5. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation, The Institution of Electrical Engineers, 2000. F.B. Jensen, M.B. Porter, W.A. Kuperman, H. Schmidt, Computational Ocean Acoustics, 2nd edition, Springer, 2011. J.F. Claerbout, Fundamentals of Geophysical Data Processing with Application to Petroleum Prospect, McGraw–Hill, 1976. M.Y. Trofimov, P.S. Petrov, A.D. Zakharenko, A direct multiple-scale approach to the parabolic equation method, Wave Motion 50 (2013) 586–595. L. Abrahamsson, H.-O. Kreiss, Boundary conditions for the parabolic equation in a range-dependent duct, J. Acoust. Soc. Am. 87 (1990) 2438–2441. V.A. Dougalis, F. Sturm, G.E. Zouraris, On an initial-boundary value problem for a wide-angle parabolic equation in a waveguide with a variable bottom, Math. Methods Appl. Sci. 32 (2009) 1519–1540. V.A. Baskakov, A.V. Popov, Implementation of transparent boundaries for numerical solution of the Schrödinger equation, Wave Motion 14 (1991) 123–128. J.S. Papadakis, Impedance formulation of the bottom boundary condition for the parabolic equation model in underwater acoustics, in: NORDA Parabolic Equation Workshop, NORDA Tech, 1982, note 143. X. Antoine, A. Arnold, C. Besse, M. Ehrhardt, A. Schädle, A review of transparent and artificial boundary condition techniques for linear and nonlinear Schrödinger equations, Commun. Comput. Phys. 4 (2008) 729–796. A.V. Popov, Accurate modeling of transparent boundaries in quasioptics, Radio Sci. 31 (1996) 1781–1790. A. Arnold, M. Ehrhardt, Discrete transparent boundary conditions for wide angle parabolic equations in underwater acoustics, J. Comput. Phys. 145 (1998) 611–638. M. Ehrhardt, Discrete artificial boundary conditions, Ph.D. thesis, Technische Universität Berlin, Berlin, 2001. D. Mikhin, Exact discrete nonlocal boundary conditions for high-order Padé parabolic equations, J. Acoust. Soc. Am. 116 (2004) 2864–2875. D. Mikhin, Analytic discrete transparent boundary conditions for high-order Padé parabolic equations, Wave Motion 45 (2008) 881–894. M. Ehrhardt, A. Zisowsky, Discrete non-local boundary conditions for split-step Padé approximations of the one-way Helmholtz equation, J. Comput. Appl. Math. 200 (2007) 471–490. B. Mayfield, Non-local boundary conditions for the Schrödinger equation, Ph.D. thesis, University of Rhode Island, Providence, RI, 1989. Z.Z. Sun, X. Wu, The stability and convergence of a difference scheme for the Schrödinger equation on an infinite domain using artificial boundary conditions, J. Comput. Phys. 214 (2006) 209–223. V.E. Grikurov, A.P. Kiselev, Gaussian beams at large distances, Radiophys. Quantum Electron. 29 (1986) 233–237. A. Arnold, Mathematical properties of quantum evolution equations, in: N. Abdallah, G. Frosali (Eds.), Quantum Transport, in: Lect. Notes Math., vol. 1946, Springer, Berlin, Heidelberg, 2008, pp. 45–109. S.-i. Doi, Smoothness of solutions for Schrödinger equations with unbounded potentials, Publ. RIMS, Kyoto Univ. 41 (2005) 175–221. D.J. Thomson, Wide-angle parabolic equation solutions to two range-dependent benchmark problems, J. Acoust. Soc. Am. 87 (1990) 1514–1520.

158

P.S. Petrov, M. Ehrhardt / Journal of Computational Physics 313 (2016) 144–158

[22] X. Antoine, C. Besse, Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation, J. Comput. Phys. 188 (2003) 157–175. [23] X. Antoine, C. Besse, Construction, structure and asymptotic approximations of a microdifferential transparent boundary condition for the linear Schrödinger equation, J. Math. Pures Appl. 80 (2001) 701–738. [24] R.M. Feshchenko, A.V. Popov, Exact transparent boundary condition for the parabolic equation in a rectangular computational domain, J. Opt. Soc. Am. A 28 (2011) 373–380. [25] A. Schädle, Non-reflecting boundary conditions for the two-dimensional Schrödinger equation, Wave Motion 35 (2002) 181–188. [26] G. Fibich, S. Tsynkov, High-order two-way artificial boundary conditions for nonlinear wave propagation with backscattering, J. Comput. Phys. 171 (2001) 632–677.