Efficient schemes for the damped nonlinear Schrödinger equation in high dimensions

Applied Mathematics Letters 102 (2020) 106158

Contents lists available at ScienceDirect

Applied Mathematics Letters www.elsevier.com/locate/aml

Efficient schemes for the damped nonlinear Schrödinger equation in high dimensions Jiaxiang Cai ∗, Haihui Zhang School of Mathematical Science, Huaiyin Normal University, Huaian, Jiangsu 223300, China

article

info

Article history: Received 9 October 2019 Received in revised form 21 November 2019 Accepted 21 November 2019 Available online 27 November 2019 Keywords: Schrödinger equation Integrating factor method Conservative system Dissipative equation Fast Fourier transform

abstract Temporal second- and high-order efficient one-step schemes are proposed for the damped nonlinear Schrödinger equation in high dimensions by combining the integrating factor method with the discretization matrices in diagonalized forms. Our schemes decouple the nonlinear part at tn+1 from the linear part and also have a compact representation, so they greatly save computational cost and storage. Numerical results are reported to exhibit the effectiveness and accuracy of the proposed schemes. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The nonlinear Schr¨ odinger-type equations are the very important models of mathematical physics, and they have been also wildly appeared in material science, plasma physics, quantum physics, nonlinear optics, water bimolecular dynamics and so on. The damped nonlinear Schr¨odinger (DNLS) equation is 2

iψt + ∆ψ + β|ψ| ψ + iγψ = 0,

(1)

where ψ(x, t) is a complex function with x = (x1 , . . . , xd )⊤ ∈ Rd and t ∈ (0, T ], the parameter β ̸= 0 is a real constant (β > 0 and β < 0 represent the focusing and defocusing cases, respectively) and the real constant γ > 0 represents the dissipation rate of system. For γ = 0, the equation is conservative. In practical computation, the equation is usually solved on a bounded Ω which is an interval in one dimension (d = 1), a rectangle in two dimensions (d = 2) or a box in three dimensions (d = 3). In this work, we focus on developing highly efficient schemes for the equation in d dimensions with initial condition ψ0 (x) = ψ(x, 0), ∗ Corresponding author. E-mail addresses: [email protected], [email protected] (J. Cai).

https://doi.org/10.1016/j.aml.2019.106158 0893-9659/© 2019 Elsevier Ltd. All rights reserved.

2

J. Cai and H. Zhang / Applied Mathematics Letters 102 (2020) 106158

x ∈ Ω and (l1 , . . . , ld )⊤ -periodic boundary conditions. The above initial-periodic value problem has a mass dissipation law ∫ 2 −2γt M (t) = e M (0), M (t) = |ψ(x, t)| dx. (2) Ω

There have been some mathematical and numerical studies on the DNLS equation in the literature. On the mathematical analysis, Temam [1] studied the regularity of solution and existence of attractors. On the numerical studies, sine spectral methods [2–4] and finite difference methods [5–9] have been adopted to solve the equation. Most of the existing methods are fully implicit and their linear and nonlinear terms are coupled at tn+1 so that they are time-consuming and require amount of memory for the problem in three or higher dimensions. Integrating factor (IF) or exponential integrator [10,11] methods and compact implicit integrating factor (cIIF) method [12] are popular for advection–diffusion equations. The idea of these methods is to evaluate the linear part exactly and then to integrate the nonlinear part with a quadrature in time. Both IF and cIIF methods improve the stability constraint associated with the linear part, so large time step can be used in computation. Besides, cIIF method has a significant improvement on memory and CUP savings. The recent works on IF and cIIF methods and their applications can be referred to [13–16]. In this work, inspired by the IF and cIIF methods together with the fast discrete Fourier transform (DFFT), we firstly develop a one-step second-order efficient stabilized compact exponential scheme for the DNLS equation. Secondly, the classical high-order IF and cIIF schemes are inconvenient to use since they are multi-step. To this end, we propose a fourth-order one-step scheme by employing the composition method. The remarkable advantages of our proposed schemes are as follows: • The nonlinear term at tn+1 is decoupled from the linear term. The schemes are highly efficient for the problem in high dimensions since only a local nonlinear system needs to be solved at each spatial grid point. • The computational cost can be further reduced through the use of DFFT by taking advantage of the circular structure of discretized matrices. • The one-step feature of high-order scheme makes it easy to use. This paper is organized as follows. Section 2 shows how to design the temporal second-order cIIF2 scheme for the DNLS equation in three dimensions. The composition method is employed to develop one-step highorder schemes in Section 3. Numerical results are reported in Section 4, and finally a conclusion is drawn in Section 5. 2. cIIF scheme for the DNLS equation in three dimensions In this section, we only show the cIIF scheme for the DNLS equation (1) in three dimensions, i.e., d = 3, x = (x, y, z) and Ω = {xa < x < xb , ya < y < yb , za < z < zb }. Our method can be applied to the equation in one, two or higher dimensions, analogously. Let us discretize the spatial domain by a rectangular mesh as follows: (xj , yk , zl ) = (xa + (j − 1)hx , ya + (k − 1)hy , za + (l − 1)hz ) where hu = (ub − ua )/Nu , u = x, y, z, 1 ≤ j ≤ Nx , 1 ≤ k ≤ Ny and 1 ≤ l ≤ Nz , and tn = nτ with τ = T /N . Denote the numerical approximation n n of ψ(xj , yk , zl , tn ) by ψj,k,l and the unknowns by Ψ n = (ψj,k,l ) ∈ CNx ×Ny ×Nz . x Nx ×Nx y Ny ×Ny z Nz ×Nz Let D ∈ R , D ∈ R and D ∈ R be the differential matrices associated with derivatives ∂xx , ∂yy and ∂zz , respectively. Due to the periodic boundary conditions, these matrices are usually circulant skew-symmetric for the second-order central differential discretization, Fourier pseudospectral discretization, high-order compact finite difference discretization, wavelet discretization and so on. We further define the following operators: x (Dx ⃝Ψ )j,k,l =

Nx ∑

x y Dj,p Ψp,k,l , (Dy ⃝Ψ )j,k,l =

p=1

Obviously, all these operators are commutative.

Ny ∑ p=1

y z Dk,p Ψj,p,l , (Dz ⃝Ψ )j,k,l =

Nz ∑ p=1

z Dl,p Ψj,k,p .

(3)

J. Cai and H. Zhang / Applied Mathematics Letters 102 (2020) 106158

3

In light of these definitions, one obtains the following semi-discretization form of (1) i 2

dΨ 2 y z x +D ⃝ y + D ⃝) z Ψ + iγΨ + β|Ψ | Ψ = 0, + (Dx ⃝ dt

(4)

2

where |Ψ | Ψ = (|ψj,k,l | ψj,k,l ) ∈ CNx ×Ny ×Nz . Note that the above circulant skew-symmetric matrices can u be diagonalized, i.e., Du = (F u )−1 ΛD F u , u = x, y, z, where F u represents the matrix of discrete Fourier u transform coefficients and the diagonal matrix ΛD = diag(F u du ) with du being the first column of Du [17]. y y x z to Eq. (4) reads Applying the operator F x ⃝F ⃝F z ⃝ ( ) 2 ˜ Dy Dz ˜ /dt + ΛDx ⃝ ˜ + iγ Ψ ˜ + β |Ψ x +Λ y +Λ z idΨ ⃝ ⃝ Ψ | Ψ = 0,

(5)

2 2 ˜ y y y y ˜ = F x ⃝F x z x z where Ψ ⃝F z ⃝Ψ and |Ψ | Ψ = F x ⃝F ⃝F z ⃝(|Ψ | Ψ ). It follows from (5) that ) ( ( x ) d ˜ 2 ˜ Dy Dz ˜ . (6) |Ψ | Ψ Ψj,k,l = i ΛD + Λ + Λ + iγ Ψ + iβ j,k,l j,j k,k l,l dt j,k,l ) ( x Dy Dz Let Aj,k,l = i ΛD j,j + Λk,k + Λl,l + iγ . To apply the integrating factor method to (6), one multiplies the

equation by e−Aj,k,l t to obtain ( ) ˜ j,k,l ) d(e−Aj,k,l t Ψ 2 ˜ = iβe−Aj,k,l t |Ψ | Ψ . dt j,k,l

(7)

Integrating on both sides of (7) from tn to tn+1 leads to ( ) ∫ τ 2 ˜ n+1 Aj,k,l τ ˜ n Aj,k,l τ −Aj,k,l s ˜ Ψj,k,l = e Ψj,k,l + iβe e |Ψ (tn + s)| Ψ (tn + s) 0

ds

(8)

j,k,l

By utilizing the trapezoidal quadrature to approximate the integrand in (8), one has a temporal second-order scheme ( ) ( ) ( ) iβτ iβτ 2 2 n+1 ˜ ˜ n+1 n Aj,k,l τ n n n+1 ˜ ˜ Ψj,k,l + Ψj,k,l = e |Ψ | Ψ + |Ψ | Ψ . (9) 2 2 j,k,l j,k,l In practical computation, the following equivalent matrix form (cIIF2) is more convenient to use, that is, ( ) iβτ n+1 iβτ ˜ 2 n 2 n+1 ∗ Aτ n n ˜ ˜ |Ψ | Ψ |Ψ ˜ | Ψ n+1 , (10) Ψ = (e ) ⊙ Ψ + + 2 2 where A = (Aj,k,l )Nx ×Ny ×Nz , (e∗ )Aτ = (eτ Aj,k,l )Nx ×Ny ×Nz and the operator ‘⊙’ means element-by-element multiplication of the two matrices. Remark 2.1. 1. The nonlinear term at tn+1 is decoupled from the linear term in the scheme (10). One only needs to solve a local nonlinear system at each grid point when the fix-point iteration is used. Therefore, the scheme is advantageous in both CPU time and memory savings. 2. τ < 2/β can ensure the convergence of the fix-point iteration. 3. High-order cIIF schemes In general, a temporal rth-order (r > 2) scheme can be obtained by evaluating the integrand in (8) with a multi-step approximation. For example, one may let 2 F (s) = e−Aj,k,l s (|Ψ (tn + ˜ s)| Ψ (tn + s))j,k,l ,

J. Cai and H. Zhang / Applied Mathematics Letters 102 (2020) 106158

4

Table 1 The solution errors of various schemes (the results of the existing methods are quoted from [21]). Method

cIIF2

[8]

[22]

[3]

[21]

L∞ L2

2.23e−10 2.70e−10

1.98e−2 2.82e−2

1.84e−2 2.66e−2

2.00e−9 2.61e−9

1.34e−9 2.21e−9

and then use a (r−1)th-order Lagrange polynomial at interpolation point tn+1 , tn , . . . , tn+2−r to approximate F (s). Clearly, the obtained rth-order multi-step scheme needs the initial values Ψ 0 , Ψ 1 , . . . , Ψ r−1 . However, one-step high-order scheme is much easier to use than the multi-step one in computations. Next, we will compose the second-order cIIF2 scheme into a temporal fourth-order scheme. Let ϕτ : yn → yn+1 be a basic one-step integrator for the ODEs y˙ = f (y), y ∈ Rd . There exist the following results: Definition 3.1 ([18]). A one-step integrator ϕτ is called self-adjoint, if it satisfies ϕ−τ ϕτ = ϕτ ϕ−τ = I. Theorem 3.1 ([18,19]). (Yoshida’s Triple Jump) Let ϕτ be a one-step self-adjoint integrator of even order m. Then the composition integrator Φτ = ϕγ3 τ ◦ ϕγ2 τ ◦ ϕγ1 τ with γ1 = γ3 = 1/(2 − 21/(m+1) ) and γ2 = −21/(m+1) /(2 − 21/(m+1) ) is of order m + 2. The result in Theorem 3.1 is helpful for composing the cIIF2 scheme into a fourth-order scheme since cIIF2 is self-adjoint/symmetric. Here, the self-adjoint/symmetric feature of cIIF2 can be verified by Definition 3.1 directly or that the scheme is unchanged by exchanging the indices n and n+1 together with the replacement 2 of τ by −τ . Denote cIIF2 by Φτ,h : Ψ n → Ψ n+1 where h = (hx , hy , hz ). It follows from Theorem 3.1 that 4 Φτ,h = Φγ21 τ,h ◦ Φγ22 τ,h ◦ Φγ21 τ,h

(11)

is a fourth-order scheme (cIIF4). Remark 3.1. 1. Higher order integrators can be obtained by repeating the above procedure. 2. Suzuki’s Fractals [18] can also be employed to compose the cIIF2 into a fourth-order cIIF scheme. 4. Numerical results In this section, we will test effectiveness and accuracy of cIIF2 and cIIF4 through some numerical experiments. In our simulations, the differences between the numerical solution and the exact solution at nth time level is measured by the maximal error norm L∞ and the average error norm L2 , and the change in ∑ 2 n mass at nth time level is monitored by |M n − M (tn )| where M n = hx hy hz j,k,l |ψj,k,l | . In the following simulations, let Du , u = x, y, z be the Fourier pseudo-spectral matrices [20]. DNLS equation in one dimension: Consider the equation with parameters β = 2 and γ =5e-3 and initial condition ψ0 (x) = sech(x)e2ix , x ∈ [−25, 25]. The simulation are done with τ =1e-5, h = 50/1024 and T = 1. Here, we use the numerical ‘exact’ solution obtained by cIIF4 with a very fine mesh to replace the exact solution of the equation to obtain the solution error. Table 1 displays the solution errors of various temporal second-order schemes. It is clear that the solution of cIIF2 is much more precise than those of others. The exact integration of the linear part may the reason that cIIF2 provides a better solution than the spectral schemes in [3,21]. The total mass and mass error are shown in Table 2. The total mass decreases as time evolves which coincides with the exact one. The change in mass at a given time indicates the mass is preserved well.

J. Cai and H. Zhang / Applied Mathematics Letters 102 (2020) 106158

5

Table 2 The results on total mass for cIIF2. Time

t=0

t = 0.2

t = 0.4

t = 0.8

t=1

Mn |M n − M (tn )|

2 0

1.99920015997937 7.0388e−13

1.99840063983082 1.4502e−12

1.99680255863813 2.9170e−12

1.99600403725843 3.6828e−12

Table 3 Numerical comparison: γ = 0.1, T = 1. (τ, h)

Method

L∞

L2

(5e−3, π/16)

cIIF2 cIIF4 [4] [9]

5.6315e−6 2.2078e−9 5.3916e−5 8.7130e−3

8.8694e−5 3.4772e−8 8.4488e−4 1.3723e−1

(2.5e−3, π/32)

cIIF2 cIIF4 [4] [9]

1.4078e−6 1.3815e−10 1.3462e−5 2.1840e−3

2.2173e−5 2.1757e−9 2.1122e−4 3.4340e−2

Table 4 The total mass and the mass error at different times: γ = 0.1, τ =5e−3, h = 2π/64, T = 1. Time

t=0

t = 0.4

t = 0.8

t=1

Exact cIIF2 Error

248.0502134424 248.0502134424 0

228.9792067534 228.9794183500 2.1160e−4

211.3744487367 211.3748105136 3.6178e−4

203.0863380528 203.0867565069 4.1845e−4

cIIF4 Error

248.0502134424 0

228.9792067728 1.9415e−8

211.3744487676 3.0932e−8

203.0863380875 3.4619e−8

DNLS equation in three dimensions: The DNLS equation in three dimensions admits an analytical solution ψ(x, t) = Kei(k·x−δ(t))−γt , x ∈ Ω = (0, 2π)3 where k = (k1 , k2 , k3 ), K is an arbitrary constant and { 2 −2γt (k · k)t + β|K| − 1), γ ̸= 0 2γ (e δ(t) = (12) 2 (k · k)t − β|K| t, γ ̸= 0. In the following simulation, we choose parameters β = 2, K = 1, k = (1, 1, 1) and hx = hy = hz = h = 2π/64. Firstly, we compare the numerical results of our schemes with those of the schemes in [4,9]. The results for γ = 0.1 at T = 1 are listed in Table 3. It is clear that the second-order cIIF2 scheme provides much more accurate solution than the existing schemes, and the solution of the fourth-order cIIF4 scheme is the best. Secondly, Fig. 1 displays the maximal error vs. time step and CPU time, respectively. The results verify that the solutions of cIIF2 and cIIF4 converge to the exact one with second and fourth order in the maximal error norm, respectively. The right graph shows, although cIIF4 repeats cIIF2 three times with time steps {γ1 τ, γ2 τ, γ1 τ } at each time step, the former are much more efficient than the latter for a given small solution error. Finally, we exhibit the total mass and mass error in Table 4. One can see that the values of mass decrease as time evolves and they are very close to the exact one at a given time. 5. Conclusion The DNLS equation is an important model in physics and it has attracted much attention on numerical studies recently. The linear part and nonlinear part at tn+1 of most of the existing methods are coupled, so that they are inefficient for the problem in high dimensions. In this paper, by using the integrating factor method, we integrate the linear part of the DNLS equation exactly, and then propose a second-order and a fourth-order cIIF schemes. The schemes are highly efficient for the multi-dimensional DNLS equation. Numerical experiments are carried out to exhibit the good solution accuracy and excellent preservation of

J. Cai and H. Zhang / Applied Mathematics Letters 102 (2020) 106158

6

Fig. 1. Numerical results for cIIF2 (blue line) and cIIF4 (red dashed line). Left: time step vs. maximal solution error (star: γ = 0; plus: γ = 0.1; pentagram: γ = 1; Blue and red lines are second- and fourth-order reference lines, respectively). Right: maximal solution error vs. CPU time.

mass dissipation. In this work, we study the equation with periodic boundary conditions. For the equation imposed on other boundary conditions such as Dirichlet boundary condition, one can move the values of ψ at nodes on the boundary to the right-hand side of (4), and then use the diagonalizable representation u of differential matrices Du , u = x, y, z, i.e., Du = (P u )−1 ΛD P u . Again, the efficient cIIF schemes can be proposed. CRediT authorship contribution statement Jiaxiang Cai: Formal analysis, Investigation, Methodology, Resources, Validation, Writing - original draft, writing - review & editing. Haihui Zhang: Formal analysis. Acknowledgments The work was supported by the Natural Science Foundation of Jiangsu Province of China (BK20181482), Qing Lan Project of Jiangsu Province of China, and Jiangsu Overseas Visiting Scholar Program for University Prominent Young & Middle-aged Teachers and Presidents, China. The authors would like to express sincere gratitude to the referee for the insightful comments and suggestions which help to improve the paper. References [1] R. Temam, Infinite-Dimensional Dynamical Systems in Mechanics and Physics, second ed., Springer Science & Business Media, New York, 2012. [2] X. Antoine, W. Bao, C. Besse, Computational methods for the dynamics of the nonlinear Schr¨ odinger/Gross–Pitaevskii equations, Comput. Phys. Comm. 184 (2013) 2621–2633. [3] W. Bao, D. Jaksch, An explicit unconditionally stable numerical method for solving damped nonlinear Schr¨ odinger equations with a focusing nonlinearity, SIAM J. Numer. Anal. 41 (2003) 1406–1426. [4] C. Jiang, Y. Song, Y. Wang, A linearized and conservative Fourier pseudo-spectral method for the damped nonlinear Schr¨ odinger equation in three dimensions, 2018, arXiv:1807.00091v1. [5] R. Ciegis, V. Pakalnyte, The finite difference scheme for the solution of weakly damped nonlinear Schr¨ odinger equation, Int. J. Appl. Sci. Comput. 8 (2001) 127–134.

J. Cai and H. Zhang / Applied Mathematics Letters 102 (2020) 106158

7

[6] S.R.K. Iyengar, G. Jayaraman, V. Balasubramanian, Variable mesh difference schemes for solving a nonlinear Schr¨ odinger equation with a linear damping term, Comput. Math. Appl. 40 (2000) 1375–1385. [7] L.S. Peranich, A finite difference scheme for solving a non-linear Schr¨ odinger equation with a linear damping term, J. Comput. Phys. 68 (1987) 501–505. [8] B.E. Moore, L. Nore˜ na, C.M. Schober, Conformal conservation laws and geometric integration for damped Hamiltonian PDES, J. Comput. Phys. 232 (2013) 214–233. [9] F. Zhang, Long-time behavior of finite difference solutions of three-dimensional nonlinear Schr¨ odinger equation with weakly damped, J. Comput. Math. 22 (2004) 593–604. [10] M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numer. 19 (2010) 209–286. [11] S.M. Cox, P.C. Matthews, Exponential time differencing for stiff systems, J. Comput. Phys. 176 (2002) 430–455. [12] Q. Nie, F.Y.M. Wan, Y. Zhang, X. Liu, Compact integration factor methods in high spatial dimensions, J. Comput. Phys. 227 (2008) 5238–5255. [13] D. Wang, W. Chen, Q. Nie, Semi-implicit integration factor methods on sparse grids for high-dimensional systems, J. Comput. Phys. 292 (2015) 43–55. [14] L. Ju, J. Zhang, Q. Du, Fast and accurate algorithms for simulating coarsening dynamics of Cahn-Hilliard equations, Comput. Mater. Sci. 108 (2015) 272–282. [15] S. Ahmed, X. Liu, High order integration factor methods for systems with inhomogeneous boundary conditions, J. Comput. Appl. Math. 348 (2019) 89–102. [16] T. Jiang, Y. Zhang, Krylov single-step implicit integration factor WENO methods for advection-diffusion-reaction equations, J. Comput. Phys. 311 (2016) 22–44. [17] J. Cai, C. Bai, H. Zhang, Decoupled local/global energy-preserving schemes for the N -coupled nonlinear Schr¨ odinger equations, J. Comput. Phys. 374 (2018) 281–299. [18] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, Springer-Verlag, Berlin, 2006. [19] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150 (1990) 262–268. [20] J. Cai, C. Bai, H. Zhang, Efficient schemes for the coupled Schr¨ odinger-KdV equations: Decoupled and conserving three invariants, Appl. Math. Lett. 86 (2018) 200–207. [21] C. Jiang, W. Cai, Y. Wang, Optimal error estimate of a conformal Fourier pseudo-spectral method for the damped nonlinear Schr¨ odinger equation, Numer. Math. PDE 34 (2018) 1422–1454. [22] H. Fu, W. Zhou, X. Qian, S. Song, L. Zhang, Conformal structure-preserving method for damped nonlinear Schr¨ odinger equation, Chin. Phys. B. 25 (2016) 110201.