Discrete transparent boundary conditions for Schrödinger-type equations for non-compactly supported initial data

Applied Numerical Mathematics 58 (2008) 660–673 www.elsevier.com/locate/apnum

Discrete transparent boundary conditions for Schrödinger-type equations for non-compactly supported initial data Matthias Ehrhardt 1 Institut für Mathematik, Technische Universität Berlin, Straße des 17, Juni 136, D-10623 Berlin, Germany Available online 16 February 2007

Abstract Transparent boundary conditions (TBCs) are an important tool for the truncation of the computational domain in order to compute solutions on an unbounded domain. In this work we want to show how the standard assumption of ‘compactly supported data’ could be relaxed and derive TBCs for a generalized Schrödinger equation directly for the numerical scheme on the discrete level. With this inhomogeneous TBCs it is not necessary that the initial data lies completely inside the computational region. However, an increased computational effort must be accepted. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. PACS: 02.70.Bf; 31.15.Fx Keywords: Schrödinger-type equation; Transparent boundary condition; Unbounded domain; Non-compactly supported initial data; Finite differences

1. Introduction Transparent boundary conditions (TBCs) are an important tool for the truncation of the computational domain in order to compute solutions on an unbounded domain (see the reviews in [12,13,28]). In this work we want to show how the standard assumption of ‘compactly supported initial data’ could be relaxed and derive inhomogeneous TBCs for a finite difference discretization of so-called standard and wide angle “parabolic” equations [16]. These models appear as one-way approximations to the Helmholtz equation in cylindrical coordinates with azimuthal symmetry and include as a special case the Schrödinger equation. With this TBCs it is not necessary that the starting field is completely inside the computational region. This is the case in underwater acoustics if the source is close to the bottom or in radiowave propagation problems when computing coverage diagrams of airborne antennas [17]. In the past two decades “parabolic” equation (PE) models have been widely used for wave propagation problems in various application areas, e.g. seismology [7], optics and plasma physics (cf. the references in [6]). Further applications to wave propagation problems can be found in radio frequency technology [24]. Here we will be mainly interested in their application to underwater acoustics, where PEs have been introduced by Tappert [26]. An account on the vast recent literature is given in the survey article [15] and in the book [14]. Thus in the sequel we will use a notation E-mail address: [email protected]. URL: http://www.math.tu-berlin.de/~ehrhardt/. 1 Supported by the DFG Research Center M ATHEON “Mathematics for key technologies” in Berlin. 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.02.002

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

661

common to the application in underwater acoustics. Nevertheless our approach is generally applicable to all one-way wave propagation problems in 2D. In oceanography one wants to calculate the underwater acoustic pressure p(z, r) emerging from a time-harmonic point source located in the water at (zs , 0). Here, r > 0 denotes the radial range variable and 0 < z < zb the depth variable. The water surface is at z = 0, and the (horizontal) sea bottom at z = zb . We denote the local sound speed by c(z, r), the density by ρ(z, r), and the attenuation by α(z, r)  0. n(z, r) = c0 /c(z, r) is the refractive index, with a reference sound speed c0 . The reference wave number is k0 = 2πf/c0 , where f denotes the (usually low) frequency of the emitted sound. In the far field approximation (k0 r 1) the outgoing acoustic field  (1) ψ(z, r) = k0 r p(z, r) e−ik0 r satisfies the one-way Helmholtz equation: √ ψr = ik0 ( 1 − L − 1)ψ, r > 0. √ Here, 1 − L is a pseudo-differential operator, and L the Schrödinger operator   L = −k0−2 ρ ∂z ρ −1 ∂z + 1 − N 2 (z, r),

(2)

(3)

where N(z, r) = n(z, r) + iα(z, r)/k0 denotes the complex refractive index. √ “Parabolic” approximations of (2) consist in formally√approximating the pseudo-differential operator 1 − L by rational functions of L. The linear approximation of 1 − λ by 1 − λ/2 gives the narrow angle or standard “parabolic” equation (SPE) of Tappert: ik0 Lψ, r > 0. (4) 2 ◦ This Schrödinger equation is a reasonable description √ of waves with a propagation direction within about 15 of the horizontal. Rational approximations of the form 1 − λ ≈ f (λ) = (p0 − p1 λ)/(1 − q1 λ) with real p0 , p1 , q1 yield the wide angle “parabolic” equations (WAPE)   p0 − p 1 L − 1 ψ, r > 0. (5) ψr = ik0 1 − q1 L ψr = −

In the sequel we will √ repeatedly require f  (0) = p0 q1 − p1 < 0. With the choice p0 = 1, p1 = 3/4, q1 = 1/4 ((1,1)Padé approximant of 1 − λ ) one obtains the WAPE of Claerbout. If we assume ρ = ρ(z) one can apply the operator 1 − q1 L to (5):       1 − q1 V + q1 k0−2 ρ∂z ρ −1 ∂z ψr = ik0 p0 − 1 − (p1 − q1 )V + (p1 − q1 )k0−2 ρ ∂z ρ −1 ∂z ψ. (6) In this paper we shall focus on boundary conditions (BCs) for the SPE (4) and the WAPE (5). At the water surface one usually employs a Dirichlet (“pressure release”) BC: ψ(z = 0, r) = 0. In the z-direction one wishes to restrict the computational domain by introducing an artificial boundary at or below the sea bottom. Papadakis [22] derived impedance BCs (or TBCs) for the SPE and the WAPE: complementing the WAPE (5) with a TBC at zb allows to recover—on the finite computational domain (0, zb )—the exact half-space solution on 0 < z < ∞. For an overview paper we refer the reader to [23]. As the SPE is a Schrödinger equation, similar strategies have been developed independently for quantum mechanical applications. While TBCs fully solve the problem of cutting off the z-domain for the analytical equation, all available numerical discretizations suffer from reduced accuracy (in comparison to the discretized half-space problem) and render the overall numerical scheme only conditionally stable [27]. In [4] discrete transparent boundary conditions (DTBCs) for a Crank–Nicolson finite difference discretization of the WAPE were constructed such that the overall scheme is unconditionally stable and as accurate as the discretized half-space problem. In this work we show how to remove one essential restriction of the (discrete) TBCs: we derive (discrete) TBCs for the case that the initial data which models a point source located at (zs , 0), is not completely contained in the computational domain. However, an increased computational effort must be accepted. While most authors [18,20] (cf. also Section 8.6 in [16]) use an exterior solution to ‘homogenize’ the problem, we directly solve the inhomogeneous problem. The usage of this DTBC is especially beneficial when one needs several computations with the same initial

662

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

data. Then the calculation of the additional term has only to be done once. The same applies when the initial field is concentrated far outside the computational domain. This is the case e.g. in radiowave propagation when computing coverage diagrams of airborne antennas; the (high) source need not be included in the domain and the computation needs far less time. The paper is organized as follows: In Section 2 we review the derivation of the continuous TBC and in Section 3 we derive the discrete TBCs for the SPE and the WAPE. Finally, in Section 4 we conclude with a numerical example. 2. Transparent boundary conditions In this section we first will derive the continuous transparent boundary condition at the sea-bottom interface. For simplicity of the presentation we will restrict ourselves to the TBC for the SPE. The basic idea of the derivation is to explicitly solve the equation in the bottom region, which is the exterior of the computational domain (0, zb ). Therefore we assume that the bottom region is homogeneous, i.e. all physical parameters are constant for z > zb (denoted with a subscript ‘b’). As the density is typically discontinuous at the water-bottom interface (z = zb ), one requires continuity of the pressure and the normal particle velocity (matching conditions): ψz (zb− , r) ψz (zb+ , r) = , (7) ψ(zb− , r) = ψ(zb+ , r), ρw ρb where ρw = ρ(zb− , r) and ρb denotes the constant density of the bottom. In order to derive the TBC we assume that the starting field ψ I (z) is continuous and consider the SPE (4) in the bottom region:   (8) ψzz + 2ik0 ψr − k02 1 − Nb2 ψ = 0, z > zb . The Laplace transformation in range of ψ is given by ∞ ˆ ψ(z, s) = ψ(z, r) e−sr dr,

(9)

0

where we set s = η + iξ , ξ ∈ R, and η > 0 is fixed, with the idea to later perform the limit η → 0. Now the exterior problem (8) is transformed to    k0  1 − Nb2 ψˆ = c2 ψ I (z), z > zb , ψˆ zz + c2 s + i (10) 2 √ where we set c = (1 + i) k0 . The basic idea is to solve this inhomogeneous second order differential equation (10) explicitly. The homogeneous solution is (11) ψˆ hom (z, s) = C1 (s) eicq(s)(z−zb ) + C2 (s) e−icq(s)(z−zb ) , z > zb ,

with the abbreviation q(s) = + s + ik0 (1 − Nb2 )/2. A particular solution of (10) is given by  z z c par icq(s)(z−ζ ) I icq(s)(ζ −z) I ψˆ (z, s) = e ψ (ζ ) dζ − e ψ (ζ ) dζ , (12) 2iq(s) zb

zb

for z > zb , i.e. the general solution in the bottom region is ˆ s) = ψˆ hom (z, s) + ψˆ par (z, s) ψ(z,  z c −icq(s)zb −icq(s)ζ I = C1 (s) e + e ψ (ζ ) dζ eicq(s)z 2iq(s) zb

 + C2 (s) e

icq(s)zb

c − 2iq(s)

∞ e zb

icq(s)ζ

c ψ (ζ ) dζ + 2iq(s)



∞

I

e z

icq(s)ζ

ψ (ζ ) dζ e−icq(s)z . I

(13)

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

663

√ Here, + denotes the branch of the square root with nonnegative real part. We note that the last term in (13) is bounded for fixed s and z √ → ∞. Since the solutions have to decrease as z → ∞, the idea is to eliminate the growing factor −icq(s)z (1−i) k0 q(s)z by simply choosing e =e c C2 (s) = 2iq(s)

∞

eicq(s)(ζ −zb ) ψ I (ζ ) dζ.

(14)

zb

Consequently, we obtain C1 (s) as: ˆ b+ , s) − C1 (s) = ψ(z

c 2iq(s)

∞

eicq(s)(ζ −zb ) ψ I (ζ ) dζ.

(15)

zb

From this we get with the matching conditions (7) the following representation of the transformed TBC : ψˆ z (zb− , s) = ic

= ic

ρw c 2 ρw q(s) C1 (s) − ρb 2 ρb

∞

eicq(s)(ζ −zb ) ψ I (ζ ) dζ

zb

ρ ρw ˆ b− , s) − c2 w q(s) ψ(z ρb ρb

∞

eicq(s)(ζ −zb ) ψ I (ζ ) dζ.

(16)

zb

It remains to inverse transform (16). If we further assume that ψ I is continuously differentiable, then integration by parts yields: ψˆ z (zb , s) =

∞

 ic ρw  2 ic ρw ˆ b , s) − ψ I (zb ) − q (s)ψ(z q(s) ρb q(s) ρb

eicq(s)(ζ −zb ) ψzI (ζ ) dζ.

zb

The inverse Laplace transformation using the convolution theorem gives ic ρw ψz (zb , r) = √ π ρb

r 0



∞ ψr (zb , τ ) ρw −1 1 L eicq(s)(ζ −zb ) ψzI (ζ ) dζ . dτ − ic √ ρb q(s) r −τ

(17)

zb

Finally, if ψzI is integrable for z > zb , Levy proved [19] that the integration and the inverse Laplace transform can be interchanged in (17) to obtain the TBC ice−ibr ρw ψz (zb , r) = √ π ρb

r 0

with b

= k0 (Nb2

ψr (zb , τ ) ice−ibr ρw dτ − √ √ πr ρb r −τ

∞ ψzI (z) e

ik0 (z−zb )2 2r

dz,

(18)

zb

− 1)/2. This TBC was derived by Levy [19] for the case Nb = 1.

Remark 1. Clearly, if ψ I (z) = 0 for z > zb then (18) reduces to the well-known TBC derived by Papadakis [21,22]:  r 2k0 − iπ ibr ρw d ψ(zb , τ ) e−ibτ e 4 e ψz (zb , r) = − dτ. (19) √ π ρb dr r −τ 0

Equivalently, it can be written as a Neumann-to-Dirichlet map [3,9] π

e 4 i ρb ψ(zb , r) = − √ 2πk0 ρw

r 0

eibτ ψz (zb , r − τ ) √ dτ. τ

(20)

664

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

Let us note that the asymptotic behaviour (for r → ∞) of the convolution kernel in the TBC (19) is O(r −3/2 ), which can be seen after an integration by parts. This decay rate and the one of the corresponding discrete kernel will be discussed in Section 3.1. Remark 2. The results of this section can be extended to the case that the refractive index in the bottom Nb depends on the range r. Using the dephasing function approach [2] (also known as gauge change in quantum mechanics) we introduce in (8) the new variable

 r  k0  ϕ(z, r) = ψ(z, r) exp i 1 − Nb2 (τ ) dτ , z > zb , (21) 2 0

to get the simple equation ϕzz + 2ik0 ϕr = 0,

(22)

z > zb .

Remark 3 (TBCs for the Standard Schrödinger equation). Since the standard Schrödinger equation is of great importance in applications we state for completeness the TBC. It is a special case of (19) using t instead of r as the temporal variable. We consider the Schrödinger equation in the form h¯ 2 ψxx + V (x, t)ψ, 2 ψ(x, 0) = ψ I (x),

i hψ ¯ t =−

x ∈ R, t > 0,

(23)

where h¯ denotes the (scaled) Planck constant, ψ I ∈ L2 (R), V (·, t) ∈ L∞ (R) and V (x, ·) is piecewise continuous. Again we assume that the potential is constant for x < xl or x > xr and the denote these values with Vl , Vr , respectively. Proceeding as before we obtain the TBCs at the left and right artificial boundaries x = xl,r :  Vl,r t 2 −i π −i Vl,r t d ψ(xl,r , τ ) ei h¯ τ h e 4e ¯ ψx (xl,r , t) = ± dτ √ dt h¯ π t −τ 0

−i

Vl,r h¯ t

e ∓ √ h¯ πt

∞ ψxI (x) e

i(x−xl,r )2 2h¯ t

(24)

dx.

zb

All these BCs (19), (20) are nonlocal in the range variable r; in range marching algorithms they thus require storing the bottom boundary data of all previous range levels. Also they involve a mildly singular convolution kernel. As motivated in the introduction we will not discretize the continuous TBCs. Instead we will show now how to derive the TBCs on a fully discrete level by mimicking the derivation of the continuous TBC. This approach avoids any unphysical reflections and conserves the stability properties of the underlying finite difference scheme. 3. Discrete transparent boundary conditions In this section we will present how to derive the discrete TBC for a Crank–Nicolson finite difference scheme for the SPE and the WAPE. Here we shall only consider uniform grids in z and r. While a uniform range discretization is crucial for our construction of discrete TBCs, this construction is independent of the (possibly nonuniform) zdiscretization on the interior domain. For simplicity we consider the uniform grid zj = j h, rn = nk (h = z, k = r) and the approximation ψj(n) ∼ ψ(zj , rn ). The discretized WAPE (6) then reads:    (n+ 1 ) (n) 1 − q1 Vj 2 + q1 k0−2 ρj D 0h ρj−1 D 0h Dk+ ψj 2

2

(n)

(n+1)

   ψj + ψj (n+ 1 ) = ik0 p0 − 1 − (p1 − q1 )Vj 2 + (p1 − q1 )k0−2 ρj D 0h ρj−1 D 0h 2 2 2

,

(25)

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673 (n+ 12 )

with Vj

665

:= V (zj , rn+ 1 ) and the usual difference operators 2

(n+1)

Dk+ ψj(n)

=

ψj

(n)

− ψj

D 0h ψj(n)

ψ

(n) j + 12

−ψ

(n) j − 12

q=

2 q1 k −1 , k p1 − q 1 0

, . = k h 2 It is well known [1] that this scheme is second order in h and k and unconditionally stable. To derive the DTBC we will now mimic the derivation of the analytic TBCs from the previous section on a discrete level. We solve the discrete exterior problem in the bottom region, i.e. the Crank–Nicolson finite difference scheme (25) for j  J − 1 :   (n+1)   (n+1) (n)  (n)  Rδb + q 2h ψj (26) − ψj = i Rκb + 2h ψj + ψj , with   δb = 1 − q1 1 − Nb2 ,

2k0 h2 , p1 − q 1 k   k  κb = k0 p0 − 1 − (p1 − q1 ) 1 − Nb2 , 2 (n)

(n)

R=

(n)

(n)

where 2h ψj = ψj +1 − 2 ψj + ψj −1 denotes the second order difference operator, and R is proportional to the parabolic mesh ratio. To solve this difference scheme (26) we use the Z-transform: ∞   (n)  (n) Z ψj = ψˆ j (ζ ) := ψj ζ −n ,

ζ ∈ C, |ζ | > Rψ ,

(27)

n=0

where Rψ denotes the convergence radius of this Laurent series. Note that we denoted in (27) the transformation variable with ζ in order to assign z for the depth variable. Thus, the difference equation (26) is transformed to

2h ψˆ j (ζ ) + iR

δb (ζ − 1) − iκb (ζ + 1) ψˆ j (ζ ) = γj (ζ ), ζ + 1 + iq(ζ − 1)

j  J − 1,

(28)

with the inhomogeneity γj (ζ ) =

ζ iRζ [δb − iκb ]

2 ψ (0) + ψ (0) . ζ +1 h j ζ + 1 + iq(ζ − 1) j

(29)

Eq. (28) is an inhomogeneous second order difference equation of the form Uj +1 + aUj + b Uj −1 = γj ,

j  J − 1.

(30)

The two linearly independent homogeneous solutions take the form α j , β j , j  J with αβ = b and a particular solution (30) can be found with the ansatz of “variation of constants”. According to Ehrhardt and Arnold [9] the general solution to the inhomogeneous equation (30) is given by  j j   1 j j j −m j −m Uj = cα + dβ + α γm − β γm , j  J − 1, (31) α−β m=J

m=J

which is the discrete analogue to the solution formula (12) in the continuous case. Now we use (31) to design a boundary condition at j = J . For that purpose we assume |α| < 1, |β| > 1 (recall that b = αβ = 1 for the Crank–Nicolson scheme for solving the wide angle parabolic equation). Proceeding analogously to the continuous case we have to eliminate the growing factor β j by choosing d appropriately as d=

∞ 1  −m β γm . α−β

(32)

m=J

We obtain from (31)  j ∞  1  −m 1 j α γm α + β j −m γm , Uj = c + α−β α−β m=J

m=j +1

j  J − 1.

(33)

666

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

The value of c can be expressed with UJ −1 :  J −1 ∞ 1  −m β UJ −1 β γm , c = J −1 − α α−β α m=J

and inserting this into (33) with j = J : ∞ 1  J −m UJ = cα + β γm α−β J

m=J

yields   ∞ ∞  1  J −m α β γm = α UJ −1 − β −1 β −m γJ +m , UJ = αUJ −1 − 1 − β α−β m=J

m=0

or equivalently bUJ −1 = βUJ +

∞ 

b−m α m γJ +m .

(34)

m=0

Finally, we want to apply these results to the discretized WAPE (26) and derive the DTBC at j = J in the situation, (0) when the initial data ψj does not vanish for j  J − 1. In this case the Z-transformed exterior Crank–Nicolson scheme reads:   δb (ζ − 1) − iκb (ζ + 1) ψˆ j +1 (ζ ) − 2 − iR ψˆ j (ζ ) + ψˆ j −1 (ζ ) = γj (ζ ), (35) ζ + 1 + iq(ζ − 1) j  J − 1, where the inhomogeneity γj is given by (29). An inverse Z-transformation of γj (ζ ) yields:     iR[δb − iκb ] −1 + iq n (0) (n) −1 n 2 (0) ψj . γj (ζ ) = (−1) h ψj + γj = Z 1 + iq 1 + iq

(36)

Thus, we can use the general solution formula (34) to obtain the transformed discrete TBC: ψˆ J −1 (ζ ) = ν2 (ζ )ψˆ J (ζ ) +

∞ 

ν1m (ζ )γJ +m (ζ ),

(37)

m=0

where ν1 , ν2 = ν1−1 are the two solutions of the quadratic equation   iR δb (ζ − 1) − iκb (ζ + 1) 2 ν + 1 = 0. ν −2 1− 2 ζ + 1 + iq(ζ − 1)

(38)

Remark 4. Setting γj ≡ 0 (37) reduces to the transformed discrete TBC derived by Arnold and Ehrhardt [4]. In order to formulate the discrete TBC we define (pm ) := (1 + iq)Z −1 {ν1m (ζ )}, m ∈ N0 , and set ((n) ) := (1 + iq)Z −1 {ν2 (ζ )}. We obtain by inverse transforming (37) (n)

(0) (n) (1 + iq)ψJ(n) −1 −  ψJ =

n−1 

(n−k) ψJ(k) + (1 + iq)γJ(n) +

k=0

∞  n 

(k) (n−k) pm γJ +m ,

(39)

m=1 k=0

n  1, with the convolution coefficients (n) given by   i i (n) −iξ  = 1 + iq + (γ − iσ )e δn0 − H (−1)n einξ 2 2 

n−1   iξ n−m −n −iξ −1 −iϕ −e λ − ελ Pm (μ) , Pn (μ) + e λ Pn−1 (μ) + ωe m=0

(40)

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

667



γ = Rδb ,

F E H2 , μ= √ , , ω= + G |E| EG ϕ 1 i ϕ = arg E, ε = |E| 2 ei 2 , 2  − 4q + i(σ + 4) , F = γ (γ − 4q) + σ (σ + 4),  − 4q − i(σ + 4) , H = γ + iσ + (γ − iσ )e−iξ .

σ = −Rκb ,

1 − iq , 1 + iq  E = (γ + iσ ) γ  G = (γ − iσ ) γ

ξ = arg

λ=

+

In (40) δn0 denotes the Kronecker symbol and Pn (μ) the Legendre polynomials (P0 ≡ 1, P−1 ≡ 0). In the nondissipative case (αb = 0) we have |λ| = 1, μ ∈ [−1, 1], and hence |Pn (μ)|  1. In the dissipative case αb > 0 we have |λ| > 1, μ becomes complex and |Pn (μ)| typically grows with n. In order to evaluate (n) in a numerically stable fashion it is therefore necessary to use the damped Legendre polynomials λ−n Pn (μ) for a numerical evaluation in (40). The convolution coefficients (n) behave asymptotically as ∼ const(−1)n einξ , n → ∞, (n) = (41) (n)

(n+1)

which may lead to subtractive cancellation in (39) (note that ψJ ≈ ψJ we introduce the summed coefficients s (n) := (n) + eiξ (n−1) ,

in a reasonable discretization). Therefore

n  1, s (0) := (0) ,

(42)

which are calculated as   Pn (μ) − Pn−2 (μ) i (n) iξ s = (1 + iq)e + (γ − iσ ) δn1 + ελ−n , n  1. 2 2n − 1 Alternatively, they can be calculated directly with the recurrence formula 2n − 3 −1 (n−1) n − 3 −2 (n−2) s (n) = − , n  3, μλ s λ s n n once s (1) , s (2) are computed from (43). Finally, we obtain the DTBC for non-compactly supported initial data (n  1): (n) (1 + iq) ψJ −1

−s

(0)

(n) ψJ

=

n−1 

(k)

(n−1)

(43)

(44)

(0)

s (n−k) ψJ − (1 − iq) ψJ −1 + 2iq(−1)n 2h ψJ

k=0

+

∞ 

(n)

(0) pm γJ +m +

m=1

∞  n−1   (k+1)  (k) (n−1−k) pm γJ +m . + eiξ pm

(45)

m=1 k=0

3.1. Decay rate of the convolution kernel Now we want to discuss the decay rate of the convolution kernel of the discrete TBC (45 in non-dissipative case of the SPE. Since |λj | = 1 the following relation was proven using asymptotic properties of the Legendre polynomials [9] √ √ 2 sin θ −n sin[(n − 12 )θ − π4 ] (n) , n → ∞, (46) λ s ∼ −ε √ √ π (n − 12 ) n where ε, μ = cos θ , 0 < θ < π are given in (40). It shows that s (n) = O(n−3/2 ), which agrees with the decay of the convolution kernel in the differential TBC (19). To show the decay property of the continuous TBC we consider the TBC (19) with N 2 ≡ 1 (k0 = 1) and obtain after an integration by parts ψz (zb , r) = c

d dr 

r 0

d =c dr

ψ(zb , τ ) dτ √ r −τ r

r− r

ψ(zb , τ ) 1 dτ − √ 2 r −τ

r− r

0

ψ(zb , τ ) ψ(zb , r − r) dτ + √ 3/2 (r − τ )

r



668

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

with

 c=−

2k0 −i π ρw (i − 1)k0 ρw e 4 = √ . π ρb ρb π

To compare the discrete convolution in the discrete TBC with the continuous convolution in the differential TBC (19) we consider the following discretization c − 2

r− r

0

−3/2 ψ(zb , τ ) c  (k)  dτ ≈ − ψJ (n − k) r

r. 3/2 2 (r − τ ) n−1

(47)

k=1

If we compare this with (45) (for compactly supported initial data) written in the form  n−1 (n) (n) ψJ − ψJ −1 1  (n−k) (k) (n−1) (n) (0) (n) =− s ψJ − ψJ −1 + s ψJ − ψJ , n  1,

z

z

(48)

k=1

we would roughly expect c z −3/2 (1 − i)k0 ρw z −3/2 = , n → ∞. (49) n s (n) ∼ − √ √ √ n 2 π ρb r 2 r This asymptotic decay can also be shown in the general case of the WAPE [4]. We will illustrate this finding in the numerical example at the end of Section 4. 3.2. Implementation In practical situations the sum (over m) in (45) of course has to be finite (e.g. up to an index m = M). This means that the initial condition is still compactly supported, but possibly outside of the computational interval. The (n) coefficients pm , m = 1, 2, . . . , M, can be calculated recursively by “continued convolution”, i.e.   (n) p1 = (1 + iq)Z −1 ν1 (ζ ) ,

(n)

p2 =

n 

(n−k) (k) p1 ,

(n)

p3 =

p1

k=0

n 

(n−k) (k) p1 ,

p2

etc.

(50)

k=0

Since this computation is rather costly (even when using fast convolution algorithms with FFTs we seek for another  n−1 (k+1) (k) (n−1−k) way to calculate ∞ + eiξ pm )γJ +m , n  1. The key idea is to use (38) for ν1 (ζ ) in order to m=1 k=0 (pm (n)

construct a recurrence relation for the pm (w.r.t. m). (38) for ν1 (ζ ) gives   δb (ζ − 1) − iκb (ζ + 1) m m+1 ν1 (ζ ) = 2 − iR ν1 (ζ ) − ν1m−1 (ζ ), m  1. ζ + 1 + iq(ζ − 1)

(51)

An inverse Z-transformation of the term in the brackets in (51) is (1 + iq)[δb + iκb ]δn0 − 2[δb + iκb ](   Z −1 [. . .] = 2δn0 − iR −1 − q 2

−1+iq n 1+iq )

,

(52)

and thus we obtain for m  1: (n)

(n) pm+1 = c1 pm − c2

n 

(n)

(n−k) (−1)k eikξ pm − pm−1 ,

(53)

k=0

with c1 = 2(1 + iq) + iRe−iξ [δb + iκb ],

c2 = 2iR[δb + iκb ],

(n) (n) and the starting sequences p0 = δn0 and p1 = (1 + iq)Z −1 {ν1 (ζ )}, (n) (n) (n−1) (−1) we consider qm := pm + eiξ pm , pm = 0 and obtain (n)

(n)

(n) (n) − c 2 pm − qm−1 , qm+1 = c1 qm

m  1,

(54) n  0. To circumvent the convolution in (53) (55)

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

669

to use in the DTBC (45) of the form (n)

(n)

(1 + iq) ψJ −1 − s (0) ψJ =

n−1 

(k)

(n−1)

(0)

s (n−k) ψJ − (1 − iq) ψJ −1 + 2iq(−1)n 2h ψJ

k=0

+

M 

(n)

(n)

(0) pm γJ +m + SM ,

n  1,

(56)

m=1

where M  n 

(n)

SM :=

(n−k)

(k) qm γJ +m ,

n  1.

(57)

m=1 k=1

The calculation of (57) with the aid of the recursion formula (55) is done by the following algorithm (1) (2) (3) (4)

(n)

q0 = (1 + iq)δn0 + (1 − iq)δn1 (n) (n) (n−1) q1 = p1 + eiξ p1 = s˜ (n)  (n) (k) (n−k) n S1 = k=1 qm γJ +m for m = 1, . . . , M − 1 do (n) (n) (n) (n) qm+1 = c1 qm − c2 pm − qm−1  (n) (n) (k) (n−k) Sm+1 = Sm + nk=1 qm+1 γJ +m+1 (0)

(0)

pm+1 = qm+1 for n = 1, . . . , N do (n) (n) (n−1) pm+1 = qm+1 − eiξ pm+1 end

end (n)

Here N denotes the maximum time index and s˜J the summed coefficients given by (43) but with the other sign in front of ‘ε’. The computational effort of the above implementation of the discrete TBC is O(M · N ) for the whole simulation, i.e. the same effort as when enlarging the computational domain sufficiently. Alternatively, a second possible implementation is to consider the transformed discrete TBC (37) and to calculate numerically the inverse Z-transform of the finite sum once: 

M    ν m (ζ ) γJ +m (ζ ) . (58) F (n) = Z −1 Fˆ (ζ ) = Z −1 (1 + iq) 1

m=0

The discrete TBC then reads (n)

(n)

(1 + iq)ψJ −1 − (0) ψJ =

n−1 

(k)

(n−k) ψJ + F (n) ,

n  1,

(59)

k=0

with the coefficient F (n) given by F

(n)

τn = 2π

2π

  Fˆ τ eiϕ einϕ dϕ,

n ∈ N0 , τ > 0.

(60)

0

Since this inverse Z-transformation cannot be done explicitly, we use a numerical inversion technique based on FFT (details can be found in [9]). 3.3. Discrete treatment of the density jump So far we did not consider the (typical) density jump at the sea bottom in the discrete TBC (56). In [4] we used an offset grid, i.e. z˜ j = (j + 12 )h, ψ˜ jn ∼ ψ(˜zj , rn ), j = −1(1)J , where the water-bottom interface with the density

670

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

jump lies between the grid points j = J − 1 and J . For discretizing the matching conditions in this case one wants to find suitable approximations for ψ and ρ at the interface zb , Ψ ∼ ψ(zb ) and ρeff = ρ(zb ), such that both sides of the discretized second matching condition (7) (n) (n) 1 ψ˜ J − Ψ 1 Ψ − ψ˜ J −1 = ρw h/2 ρb h/2

1 ψ˜ J − ψ˜ J −1 . ρeff h (n)

are equal to

(n)

(61)

This approach results in an effective density ρeff = (ρw + ρb )/2 (based on a different derivation this was also used by (n) (n) (n) Collins [8]). At the surface we use instead of ψ0 = 0 the offset BC ψ˜ 0 = −ψ˜ −1 . Finally it remains to reformulate the discrete TBC (56) such that the density jump is taken into account. We rewrite the discretization of the second depth derivative at j = J from (25):     ρb  (n) (n)  (n) (n)  h2 ρJ D 0h ρJ−1 D 0h ψ˜ J = 2h ψ˜ J + 1 − (62) ψ˜ J − ψ˜ J −1 . ρeff 2 2 Comparing the r.h.s. of (62) to (26) we observe that only one additional term appears, and instead of (28) we get    δb (ζ − 1) − iκb (ζ + 1) ˆ ρb  ˆ ψ˜ J (ζ ) = ψ˜ˆ J +1 (ζ ) − 1 − iR (63) ψ˜ J (ζ ) − ψ˜ˆ J −1 (ζ ) + γ˜J (ζ ) ζ + 1 + iq(ζ − 1) ρeff with the inhomogeneity   iRζ [δb − iκb ] ζ ρb  (0) (0) (0) (0)  (0) ˜ ˜ ˜ ˜ γ˜J (ζ ) = (64) ψJ +1 − ψJ − ψ˜ . ψJ − ψJ −1 + ζ +1 ρeff ζ + 1 + iq(ζ − 1) J  m Using ψ˜ˆ J +1 (ζ ) = ν1 (ζ ) ψ˜ˆ J (ζ ) − ∞ m=1 ν1 γ˜J +m (ζ ), where ν1 (z) denotes the solution of (38), and considering the −1 fact that ν1 (ζ ) + ν1 (ζ ) − 1 is equal to the term in the squared brackets in (63) we obtain the Z-transformed discrete TBC: ∞   ρb  ˆ ν1m γ˜J +m (ζ ). (65) ψ˜ J (ζ ) − ψ˜ˆ J −1 (ζ ) = ψ˜ˆ J (ζ ) − ν2 (ζ )ψ˜ˆ J (ζ ) − ρeff m=0

Hence, the discrete TBC including the density jump reads     ρb (n) ρb (n) − (0) ψ˜ J (1 + iq) ψ˜ J −1 + (1 + iq) 1 − ρeff ρeff =

n−1 

(n−k) ψ˜ J + (1 + iq) γ˜J + (k)

(n)

k=0

∞  n 

(n−k)

(k) pm γ˜J +m ,

n  1,

(66)

m=1 k=0

and formulated for the summed coefficients:     ρb (n) ρb (0) ˜ (n) ˜ (1 + iq) −s ψJ ψ + (1 + iq) 1 − ρeff J −1 ρeff   n−1  ρb (n−1) ρb (n−1) (n−k) ˜ (k) ˜ ψ˜ J ψJ − (1 − iq) ψ = s − (1 − iq) 1 − ρeff J −1 ρeff k=0

+ 2iq(−1)

n

(0)

2h ψ˜ J

+

∞  m=1

(0) (n) pm γ˜J +m

+

∞  n 

(n−k)

(k) qm γ˜J +m

(67)

m=1 k=1

with the coefficients s (n) given by (43). Setting ρb = ρeff this discrete TBC (67) reduces to (45). In the derivation of the discrete TBC we needed a uniform range discretization and a uniform z-discretization for the exterior problem on z > zb , i.e. j  J − 1. For the interior problem, however, a nonuniform discretization (even adaptive in r) may be used, and this would not change our discrete TBC (66). For any given interior z-discretization and a uniform grid spacing h in the exterior domain, the discrete TBC will always yield, on the interior domain, the same solution as the corresponding discrete half-space solution (up to round-off errors). We refer the reader to [4] for a concise discussion of strategies for the ‘best exterior discretization’ (i.e. in the sea bottom).

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

671

4. Numerical example In this section we present a simple model example to illustrate the numerical results when using our new discrete TBC for the SPE (4). We use an initial field, that is partially outside the computational domain. Due to its construction, our DTBC yields exactly (up to round-off errors) the numerical whole-space solution restricted to the computational interval [0, zb ]. Example 1. This example shows a simulation of a right travelling Gaussian beam [ψ I (z) = exp(i100z−30(z−0.8)2 )] at three consecutive times evolving under the SPE with N 2 ≡ 1 (k0 = 1) and with the rather coarse depth discretization of 161 grid points for the interval 0  z  1 (i.e. z = 1/160) and the range step r = 2 × 10−5 . For the exterior (computational) domain we choose the same depth step z and use 60 grid points which results in the exterior interval 1 < z  1.38125. In the following Fig. 1 we plotted the absolute value of the initial data and the solution obtained with the discrete TBCs (45) at the range steps r = 0.002, r = 0.004, r = 0.006. One clearly sees in Fig. 1 that the solution is solely propagated to the right and no artificial reflections are caused. In this example the computation using the inhomogeneous DTBCs (56) needs approximately the same CPU-time than just enlarging the domain to the interval 0  z  1.8 using a simple Neumann boundary condition at z = 1.8. From Fig. 1 one can guess that the solution at r = 0.004 has already reached the boundary at z = 1.8. Hence it is worthwhile in this example to use the inhomogeneous DTBCs whenever the solution for r > 0.004 is needed.

Fig. 1. Solution |ψ(z, r)| at range r = 0, r = 0.002, r = 0.004, r = 0.006: the solution with the new discrete TBCs (45) coincides with the whole-space solution and does not introduce any numerical reflections.

672

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

Fig. 2. Discrete TBC: convolution coefficients s (n) compared to the r.h.s. of (46) and the decaying rate (49).

Finally, we compare in Fig. 2 the summed convolution coefficients s (n) for increasing range levels n with the r.h.s. of (46). Fig. 2 displays a quite good agreement of the asymptotic behaviour of the s (n) with the predicted one in (46). The values of s (n) oscillate around the rate (49). 5. Conclusion We have derived discrete transparent boundary conditions for starting fields that are (partially) supported outside of the computational domain. While the discrete TBC solves the problem of cutting off the computational domain, the resulting numerical effort of this approach is not completely settled yet and subject to further investigations. In (m) (m) particular one has to compare an “optimal” computation algorithm for the coefficients pn or qn with simulations on a sufficiently enlarged computational domain. In a subsequent paper we will introduce a suitable approximation to the discrete convolutions appearing in the inhomogeneous DTBC (67) in the spirit of the ideas in [5,10]. In [11,20,25] (semi)-discrete TBCs for split-step methods for high-order WAPE like a general (p, q)-Padé approximation to the one-way Helmholtz equation (2) are developed. It seems to be possible to extend all the technical tools to this general case. This topic will also be regarded in a forthcoming paper. References [1] G.D. Akrivis, V.A. Dougalis, G.E. Zouraris, Error estimates for finite difference methods for a wide-angle “parabolic” equation, SIAM J. Numer. Anal. 33 (1996) 2488–2509. [2] 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. [3] A. Arnold, Numerically absorbing boundary conditions for quantum evolution equations, VLSI Design 6 (1998) 313–319. [4] A. Arnold, M. Ehrhardt, Discrete transparent boundary conditions for wide angle parabolic equations in underwater acoustics, J. Comput. Phys. 145 (1998) 611–638. [5] A. Arnold, M. Ehrhardt, I. Sofronov, Discrete transparent boundary conditions for the Schrödinger equation: Fast calculation, approximation, and stability, Comm. Math. Sci. 1 (2003) 501–556. [6] A. Bamberger, B. Engquist, L. Halpern, P. Joly, Parabolic wave equation approximations in heterogeneous media, SIAM J. Appl. Math. 48 (1988) 99–128. [7] J.F. Claerbout, Fundamentals of Geophysical Data Processing, McGraw-Hill, New York, 1976. [8] M.D. Collins, A higher-order parabolic equation for wave propagation in an ocean overlying an elastic bottom, J. Acous. Soc. Am. 86 (1989) 1459–1464.

M. Ehrhardt / Applied Numerical Mathematics 58 (2008) 660–673

673

[9] M. Ehrhardt, A. Arnold, Discrete transparent boundary conditions for the Schrödinger equation, Riv. Mat. Univ. Parma 6 (2001) 57–108. [10] M. Ehrhardt, A. Arnold, Discrete transparent boundary conditions for wide angle parabolic equations: Fast calculation and approximation, in: D.G. Simons (Ed.), Proceedings of the Seventh European Conference on Underwater Acoustics, July 3–8, 2004, Netherlands Organisation for Applied Scientific Research (TNO) and Delft University of Technology, Delft, The Netherlands pp. 9–14. [11] 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. [12] D. Givoli, Non-reflecting boundary conditions, J. Comput. Phys. 94 (1991) 1–29. [13] D. Givoli, Numerical methods for problems in infinite domains, Studies in Applied Mechanics, vol. 33, Elsevier, Amsterdam, 1992. [14] F.B. Jensen, W.A. Kuperman, M.B. Porter, H. Schmidt, Computational Ocean Acoustics, AIP Press, New York, 1994. [15] D. Lee, S.T. McDaniel, Ocean acoustic propagation by finite difference methods, Comput. Math. Appl. 14 (1987) 305–423. [16] M.F. Levy, Parabolic Equation Models for Electromagnetic Wave Propagation, IEE Electromagnetic Waves Series, vol. 45, IEE, London, 2000. [17] M.F. Levy, Fast PE models for mixed environments, in: AGARD Conference on Propagation Assessment in Coastal Environments, Bremerhaven, Germany, 19–22 September, 1994, AGARD—Advisory Group for Aerospace Research & Development, Neuilly-sur-Seine, France, North Atlantic Treaty Organization, pp. 8-1–8-6. [18] M.F. Levy, Transparent boundary conditions for parabolic equation: Solutions of radiowave propagation problems, IEEE Trans. Antennas Propag. 45 (1997) 66–72. [19] M.F. Levy, Non-local boundary conditions for radiowave propagation, in: G. Cohen (Ed.), Third International Conference on Mathematical and Numerical Aspects of Wave Propagation Phenomena, Juan-les-Pins, France, 24–28 April 1995, SIAM, Philadelphia, PA, 1995, pp. 499– 505. [20] D. Mikhin, Exact discrete nonlocal boundary conditions for high-order Padé parabolic equations, J. Acous. Soc. Am. 116 (2004) 2864–2875. [21] J.S. Papadakis, Impedance formulation of the bottom boundary condition for the parabolic equation model in underwater acoustics, NORDA Parabolic Equation Workshop, NORDA Tech. Note 143, 1982. [22] J.S. Papadakis, Impedance bottom boundary conditions for the parabolic-type approximations in underwater acoustics, in: R. Vichnevetsky, D. Knight, G. Richter (Eds.), Advances in Computer Methods for Partial Differential Equations VII, IMACS, New Brunswick, NJ, 1992, pp. 585–590. [23] J.S. Papadakis, Exact nonreflecting boundary conditions for parabolic-type approximations in underwater acoustics, J. Comput. Acous. 2 (1994) 83–98. [24] A.V. Popov, N.Y. Zhu, Low-order modes in smoothly curved, oversized and lossy waveguides of arbitrary cross section: Parabolic equation approach, AEÜ Int. J. Electr. Commun. 54 (2000) 175–182. [25] F. Schmidt, T. Friese, D. Yevick, Transparent boundary conditions for split-step Padé approximations of the one-way Helmholtz equation, J. Comput. Phys. 170 (2001) 696–719. [26] F.D. Tappert, The parabolic approximation method, in: J.B. Keller, J.S. Papadakis (Eds.), Wave Propagation and Underwater Acoustics, Lecture Notes in Physics, vol. 70, Springer, New York, 1977, pp. 224–287. [27] D.J. Thomson, M.E. Mayfield, An exact radiation condition for use with the a posteriori PE method, J. Comput. Acous. 2 (1994) 113–132. [28] S.V. Tsynkov, Numerical solution of problems on unbounded domains. A review, Appl. Numer. Math. 27 (1998) 465–532.