- Email: [email protected]
Numerical methods for the simulation of trapped nonlinear Schrödinger systems
Applied Mathematics and Computation 144 (2003) 215–235 www.elsevier.com/locate/amc
Numerical methods for the simulation €dinger systems of trapped nonlinear Schro Vıctor M. Perez-Garcıa a
a,*
, Xiao-yan Liu
b
Departamento de Matem aticas, Escuela T ecnica Superior de Ingenieros Industriales, Universidad de Castilla-La Mancha, Avenida de Camilo Jose Cela, s/n, 13071 Ciudad Real, Spain b Department of Mathematics, Northeast Normal University, Changchun 130024, China
Abstract We propose, analyze and compare the efficiency and accuracy of different numerical schemes for the solution of the nonlinear Schr€ odinger equation with a trapping potential. In particular we study schemes of finite difference, pseudospectral and spectral types for the space discretization together with explicit symplectic, multistep, split-step and standard variable-step integrators to solve the time evolution. All of these schemes are evaluated comparatively and some recommendations based on their accuracy and computational efficiency are made. Ó 2002 Elsevier Inc. All rights reserved. Keywords: Pseudospectral schemes; Nonlinear Schr€ odinger equations; Finite difference schemes; Symplectic schemes
1. Introduction The nonlinear Schr€ odinger equation (NLSE) is one of the most important models of mathematical physics, with applications to different fields such as plasma physics, nonlinear optics, water waves, biomolecule dynamics and many other fields (see e.g. Refs. [1–5]). In many of those fields the equation appears as an asymptotic limit for a slowly varying dispersive wave envelope propagating in a nonlinear medium. A new burst of interest on this equation has been motivated by the recent achievement of Bose–Einstein condensation *
Corresponding author. E-mail address: [email protected] (V.M. Perez-Garcıa).
0096-3003/02/$ - see front matter Ó 2002 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00402-2
216
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
using ultracold neutral bosonic gases. In this case the NLSE is obtained as a mean field model for the dynamics of the quantum gas [6], which, in reduced units reads 1 1 2 iWt ¼ DW þ x2 W þ U jWj W; 2 2
ð1Þ
where W : X R2 ! C, D ¼ ðo2 =ox21 Þ þ þ ðo2 =ox2n Þ, x2 ¼ x21 þ þ x2n , n ¼ 1; 2; 3 being the spatial dimensionality of the problem and U 2 R. In BEC problems one is interested on the solution of the Cauchy problem for initial data in L2 in unbounded domains. From now on we will restrict our attention to the case n ¼ 1, i.e. the problem to be solved is i
ow 1 o2 w 1 2 2 ¼ þ x w þ U jwj w; ot 2 ox2 2
wðx; 0Þ ¼ w0 ðxÞ; lim wðx; tÞ ¼ 0:
x! 1
ð2aÞ ð2bÞ ð2cÞ
The effect of the potential term V ðxÞ ¼ x2 =2 is to localize the solutions near x ¼ 0. This is why this term is called a ‘‘trap’’. In what follows we will refer to this problem as harmonically trapped nonlinear Schr€ odinger equation (HTNLSE). This one-dimensional situation is relevant at least for applications in the fields of Bose–Einstein condensation [7–10] and nonlinear optics [11] and has received attention in the last years. In addition to its specific interest, the analysis of the one-dimensional problem is the first step for the analysis of more complicated multidimensional configurations. Concerning the HTNLSE, it is known that this equation has at least two conserved quantities which are Z 2 N ðwÞ ¼ jwðx; tÞj dx; ð3aÞ H EðwÞ ¼
Z 2 1 ~ 1 1 2 4 rw þ x2 jwj þ U jwj dx: 2 2 2
ð3bÞ
In our paper we will consider solutions with N ðwÞ ¼ 1 since other values may be scaled into the nonlinear coefficient U. Despite the many applied works which rely on numerical simulations of the HTNLSE, this equation has not received specific attention in the mathematical literature. This is in contrast with the fact that the NLSE without trapping potential is a very well studied problem for which different numerical schemes have been analyzed. For instance, the usual NLS equation has been solved using finite differences (FDs) in space plus time variables [12–15] (in some of
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
217
those cases the time-discretization was conservative), FD in space plus symplectic integration in time [16,17], finite elements [18–21] and pseudospectral methods for the discretization of the spatial variable [22–24]. Of remarkable success are the methods based on split-step decompositions and FFT [25,26], which have been adopted as standard integrators by a relevant part of the scientific community (see e.g. [27]). In this paper we consider several numerical schemes for Eqs. (2a)–(2c) and compare their accuracy and performance by using numerical simulations. To our knowledge and despite the practical importance of these problems, this analysis has not been considered previously in the scientific computing literature. In particular for the space discretization we will consider second order FD schemes, spectral and pseudospectral methods with different basis functions. For the time evolution we use several methods such as split-step, centered (conservative) FDs, modified Euler method and explicit symplectic integrators.
2. Exact solutions of the HTNLSE The simplest solutions of Eqs. (2a)–(2c) are the so-called solitary waves or stationary solutions, which are of the form wl ðx; tÞ ¼ /l ðxÞeilt ;
ð4Þ
where l 2 R and /l ðxÞ : R ! R are to be determined from the equations 1 o2 1 2 2 x l/l ðxÞ ¼ þ þ U j/ j ð5aÞ /l ðxÞ; l 2 ox2 2 lim /l ðxÞ ¼ 0:
x! 1
ð5bÞ
In the linear case (U ¼ 0) this problem has a numerable set of stationary real 2 solutions /k ðxÞ ¼ Hk ðxÞex =2 , with k 2 N and Hk ðxÞ being the Hermite polynomials [29] and lk ¼ k þ 1=2. For the nonlinear repulsive case (U > 0) it is has been proven in limited situations [11] and conjectured for the general case that the situation is the same, i.e. that for fixed N ¼ 1 there is a numerable set of solutions fln ; /ln g. In this paper we will be interested on the ground state solution of Eqs. (5a) and (5b) which is the stationary solution with smallest value of the eigenvalue l for a fixed value of the L2 -norm. It is known that the ground state of Eqs. (5a) and (5b) is p a ffiffiunique positive radial solution [30] which for r 0 behaves as ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /l ðrÞ ð1= uÞ ðl ð1=2Þr2 Þþ and for r ! 1 is bounded by /l ðrÞ 6 2 ½C0 þ CðlÞrer =2 [29]. General properties of the solutions of Eqs. (5a) and (5b) have been studied in a more general context [31].
218
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Although many properties of the ground state are known and there are several asymptotic formulae available for its approximation, a closed analytical form giving its shape is not available, thus the most accurate way of obtaining it is by means of numerical methods. This turns out to be a difficult problem since the approximation of /0 ðxÞ requires the solution of a nonlinear eigenvalue problem (l0 is also unknown) for which simple shooting methods are not useful. This is because the asymptotic behavior for large x values is a delicate balance between the increasing x2 term and the other ones present in the equation. Thus, to obtain accurate approximations of the ground states, which will be later used as initial data we have used the so called Sobolev-gradients method as described in Ref. [33]. Following this scheme we have obtained the ground state for three values of the nonlinearity U ¼ 2; 10; 100 and N ð/0 Þ ¼ 1 after a relaxation time of T ¼ 100:000 with s ¼ 0:2, on the uniformly spaced grid Xh ¼ ½Xmax ; Xmax þ h; . . . ; Xmax with Xmax ¼ 40 and Nx ¼ 8192 points. Working in this way we get /0 ðxÞ with approximate maximum errors of the order of 104 –105 (this depends on the specific value of U), which is enough for our present study. Through this work when approximations of /0 ðxÞ for points xj 62 Xh have been necessary we have used Lagrange interpolation with cubic polynomials to get those values from their four nearest neighbors. The dynamics of the ground state although could reveal some simple properties of the numerical schemes to be later used is in a certain sense trivial. To build more interesting solutions of Eqs. (2a)–(2c) we will use a property of stationary solutions of HTNLSE which allows to build time dependent solutions from stationary ones [32]. For our present purposes we present a limited formulation which is that given a stationary solution of (5a) and (5b) then one may build infinitely many new time dependent solutions of the form wðx; t; l; X0 ; V0 Þ ¼ /l ðx X ðtÞ; tÞeiðltþxX ðtÞÞ ; _
ð6Þ
where X ðtÞ satisfies X€ þ X ¼ 0 with arbitrary initial data X ðt0 Þ ¼ X0 ; X_ ðt0 Þ ¼ V0 . The existence of these nontrivial solutions provides a way to study quantitatively the errors of the different numerical methods.
3. Explicit symplectic schemes 3.1. The scheme Although the simplest explicit schemes are usually shown to be unstable, there are explicit symplectic schemes which do not display such unwanted
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
219
behavior. This is why these schemes, which are simple to implement have been used in several applications to linear problems [34–36]. Let us first define the discrete grid Xh ¼ fx1 ; . . . ; xnx g ¼ fX ; X þ h; . . . ; X h; X g with Nx elements. We rewrite our equations by defining wj ðtÞ ¼ wðxj ; tÞ and decomposing wj ðtÞ ¼ aj ðtÞ þ ibj ðtÞ. The nonlinear equation (2a) can be written on the discrete grid in the following way oa ¼ Hb; ot
ð7aÞ
ob ¼ Ha; ot
ð7bÞ >
>
where a ¼ ða1 ; a2 ; . . . ; an Þ ; b ¼ ðb1 ; b2 ; . . . ; bn Þ ; H ¼ HK þ HV where HK stands for the effect of the operator ðo2 =ox2 Þ on the particular grid point and HV is the diagonal matrix 01 2 1 x þ U ða21 þ b21 Þ 2 1 1 2 2 2 B C x þ U ða2 þ b2 Þ 2 2 B C HV ¼ B C: . .. @ A 1 2 2 2 x þ U ðan þ bn Þ 2 n ð8Þ Eqs. (7a) and (7b) form a hamiltonian system which may be accurately integrated by using symplectic schemes. For instance, a second order explicit scheme is given by A ¼ an þ p1 H sbn ;
ð9aÞ
n
B ¼ b þ q1 H sA;
ð9bÞ
anþ1 ¼ A þ p2 H sB;
ð9cÞ
bnþ1 ¼ B þ q2 H sanþ1 ; ð9dÞ pffiffiffi pffiffiffi pffiffiffi pffiffiffi where p1 ¼ 1 ð1= 2Þ; p2 ¼ ð1= 2Þ; q1 ¼ ð1= 2Þ; q2 ¼ 1 þ ð1= 2Þ are the coefficients defining the method. The matrix form of this scheme is
nþ1 n M11 M12 a a ¼ ; ð10Þ M21 M22 bn bnþ1 where
M¼
¼
1 q2 H s
0 1
1 0
p2 H s 1
1 þ q1 p2 H 2 s2 ðq1 þ q2 ÞH s þ q1 q2 p2 H 3 s3
1 0 q1 H s 1
1 0
p1 H s 1
ðp1 þ p2 ÞH s þ p1 p2 q1 H 3 s3 : 1 þ ðp1 q1 þ p1 q2 þ p2 q2 ÞH 2 s2 þ p1 p2 q1 q2 H 4 s4
220
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
To define a fully discrete problem we must discretize the spatial derivatives. Using FDs and zero Dirichlet boundary conditions on the boundaries, the operator HV may be written into the standard form 1 0 2 1 0 0 0 B 1 2 1 0 0 C B .. C .. C B 1 B 0 . C . ð11Þ HK ¼ 2 B . C: .. 2h B .. . 0 C C B @ 0 0 1 2 1 A 0
0
0
1
2
3.2. Stability of symplectic schemes Eq. (10) is equivalent to a linear mapping Z kþ1 ¼ MZ k where k ¼ 0; 1; 2; . . ., > corresponds to t ¼ t0 ; t0 þ s; t0 þ 2s; . . ., and Z ¼ ða; bÞ . The stability of such a mappings depends on the eigenvalues of M. In particular, stable behavior occurs if these eigenvalues have modulus 6 1. The eigenvalues of M, defining x0 ¼ sh, are given by ! M11 ðx0 Þ k M12 ðx0 Þ det ¼ k2 ðM11 ðx0 Þ þ M22 ðx0 ÞÞk M21 ðx0 Þ M22 ðx0 Þ k þ M11 ðx0 ÞM22 ðx0 Þ M12 ðx0 ÞM21 ðx0 Þ ¼ 0: Here we have used the definitions M11 ðx0 Þ ¼ 1 þ q1 p2 x02 ; M12 ðx0 Þ ¼ p1 x0 þ p2 x0 þ p1 p2 q1 x03 ; M21 ðx0 Þ ¼ q1 x0 þq2 x0 þq1 q2 p2 x03 ; M22 ðx0 Þ ¼ 1þp1 q1 x02 þp1 q2 x02 þ 02 04 p2 q2 x þ p1 p2 q1 q2 x . The matrix M is symplectic and then jMj ¼ 1, i.e.
M11 ðx0 Þ M12 ðx0 Þ det ¼ M11 ðx0 ÞM22 ðx0 Þ M12 ðx0 ÞM21 ðx0 Þ ¼ 1: ð12Þ M21 ðx0 Þ M22 ðx0 Þ Thus, we get k2 ðM11 ðx0 Þ þ M22 ðx0 ÞÞk þ 1 ¼ 0:
ð13Þ
In order to guarantee jkj 6 1, we need to satisfy jk1 þ k2 j ¼ jM11 ðx0 Þ þ M22 ðx0 Þj 6 jk1 j þ jk2 j 6 2;
ð14Þ
thus we find that the stability condition is kx0 k 6 2:26 [37]. That is to say if kx0 k 6 xs ¼ 2:26 using an appropriate matrix norm then any eigenvalue jkj 6 1. Then xs ¼ 2:26 is the stability limit for the 2-order explicit symplectic scheme. Now we can estimate a maximum stable time step ss from xs by using the fact that for any natural norm the spectral radius of the matrix qðH Þ 6 kH k and thus
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
ss ¼
xs : qðH Þ
221
ð15Þ
For instance, if we consider the FD discretization matrix for HK and using Chebyshev approximation we get qðH Þ <
2 1 U 2 þ max jxj j þ : h2 2 j¼1;...;nx h
ð16Þ
This equation provides a bound for the maximum time step which may be used for the stable integration of Eq. (2a) by the explicit second order explicit scheme ss 6
2 h2
xs : þ 12 X 2 þ Uh
ð17Þ
It is remarkable that not only the usual term 2=h2 due to the FD discretization of the problem is present but also two terms due to the nonlinear effects and boundary effects which also restrict the time step. In fact, for large integration regions (X) and values of the nonlinear parameter (U) these terms would be responsible for the stability constraints of the scheme, a result which is very different for that of the free (nontrapped) LSE. The fact that the nonlinear term contributes with a term of type U =h is a first indication of the relevant role of the nonlinear parameter and its possible connection with unstable behavior, a connection which will be more clear later. This type of stability restrictions appear also when higher order schemes are used. The only difference is on the particular value of the parameter xs , which should be replaced by the corresponding value of the stability condition of the scheme for the free problem. 3.3. Results We have computed the solutions of Eq. (2a) using the explicit symplectic scheme (9a)–(9d) (similar simulations where done with a fourth order explicit scheme with very similar qualitative results). We have used as initial state the displaced ground state /0 ðx 1Þ obtained above on a grid with Xmax ¼ 40 and Nx ¼ 4096 points. All simulations were repeated for U ¼ 2, 10 and 100. The corresponding stability limits for these values of the nonlinear parameter obtained from Eq. (15) are ss ðU ¼ 2Þ ¼ 3:68 104 , ss ðU ¼ 10Þ ¼ 3:45 104 and ss ðU ¼ 100Þ ¼ 2:03 104 , respectively. It can be seen in Fig. 1 that when U ¼ 2 taking s ¼ 3 104 gives reasonable results (Fig. 1(a)–(c)). However when U ¼ 10 and s ¼ 3 104 the error in the solution becomes already very large for t ’ 20 (see Fig. 1(d)–(f)). For these values of the nonlinear parameter, despite they lay above the stability limit the errors accumulate in such a way that for moderate times the approximate solution is not good. Taking smaller time steps, e.g. s ¼ 5 105 as shown in
222
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Fig. 1. Results from typical simulations using the second order explicit symplectic schemes for U ¼ 2 and U ¼ 10, respectively, up to T ¼ 100. Shown are (a)–(c) simulations with Nx ¼ 4096, s ¼ 0:0003: (a) solution jwðx; tÞj2 for t ¼ 10, (b) error jeðx; tÞj2 for t ¼ 100, (c) evolution of the maximum error EðtÞ ¼ kwðx; tÞ wnj k1 . (d)–(f) Same as (a)–(c) but for U ¼ 10 with s ¼ 0:0003. (g)– (i) Same as (d)–(f) but for U ¼ 10 and s ¼ 0:00005.
Fig. 1(g)–(i) solves this problem at the cost of a higher computational cost. We have verified that for U ¼ 100 very small time steps must be used to get reasonable results, so small in fact that this scheme becomes unusable for the numerical approximation of the solution for long time simulations.
4. A conservative finite difference scheme for the HTNLS 4.1. Introduction The lack of stability of the symplectic scheme used in Section 3.2 was partially due to the explicit nature of the method which is known to impose strict stability constraints on these type of problems. Although there is the possibility
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
223
of using implicit symplectic FD schemes in the space and time variables such as those used in [16,17] the computational implementation of those schemes is computationally expensive mainly due to the need for solving a system of nonlinear equations at each time step. This is why our next step will be to consider an extension of a conservative FD scheme which has been used successfully in several contexts [15] and which is linearly implicit and leads to very efficient computational implementations. 4.2. The finite difference conservative scheme A natural extension of the linearly implicit, conservative numerical scheme of Ref. [15] is the following " nþ1 # nþ1 nþ1 n1 n1 n1 wjnþ1 wjn1 1 wjþ1 2wj þ wj1 wjþ1 2wj þ wj1 i ¼ þ 2 2s 2h2 2h2 nþ1
1 wj þ xj 2
þ wjn1 wnþ1 þ wjn1 n 2 j þ U jwj j : 2 2
ð18Þ
If, as in the previous case, zero Dirichlet boundary conditions are used and introducing aj ¼ 1=h2 þ i=s x2j =2 and the matrix 1 0 1 0 0 a1 þ U jw1 j2 2h2 C B 1 a2 þ U jw2 j2 2h12 0 C B 2h2 C B .. .. C B C B 0 . . A¼B C; . . C B . . . . 0 C B C B 2 n 1 1 a þ U jw j 0 A @ Nx 1 Nx 1 2h2 2h2 2 n 1 0 0 aNx þ U jwNx j 2h2 ð19Þ
the numerical scheme may be written in matrix form as Awnþ1 ¼ A wn1 ;
ð20Þ
where w ¼ ðw1 ; . . . ; wM ÞT . Following the methods developed in Ref. [15], it is straightforward to prove that this scheme is stable and convergent and that no artificial blow-up is possible. The main advantages of this scheme are that it is linearly implicit and conservative. This scheme has been successfully applied when there is no potential to different nonlinear Schr€ odinger problems. By means of numerical simulations we have found that in spite of its stability properties the scheme does not perform well for long-time simulations of Eq. (2a). Although the accuracy of the scheme is not bad for short time simulations in agreement with the Oðs2 þ h2 Þ truncation error of the scheme, it
224
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
does not work properly for long simulations unless a very small time step is chosen, with a behavior similar to that of the symplectic scheme. An example of this behavior will be presented in a later section.
5. Split-step schemes for the HTNLS 5.1. Introduction The fact that simple numerical schemes such as those presented up to now require very small time steps to provide accurate results in long-time integrations of the HTNLS equation has motivated the search for other type of schemes of different nature such as those which will be presented in this section. Split-step schemes for integrating partial differential equations are based on the observation that many problems may be decomposed into exactly solvable parts and the exact solution of the full problem can be approximated as a composition of the solutions of the individual problems. For instance the solution of a partial differential equation of the type ot uðx; tÞ ¼ N ðt; x; u; ox ; . . . ; Þu ¼ ðA þ BÞu;
ð21Þ
may be approximated from the exact solutions of the subproblems ot u ¼ Au;
ð22aÞ
ot u ¼ Bu:
ð22bÞ
This way of obtaining solutions by partitioning of the original problem has been studied both in the field of the approximation of ordinary [28] and partial differential equations. In the case of the ordinary linear or NLSE this type of schemes have been successfully applied previously [26]. 5.2. Fourier pseudospectral scheme Let us first split the operator in Eq. (2a) following the standard approach, which consists in taking A¼
1 o2 ; 2 ox2
1 B ¼ x2 þ U jwj2 : 2
ð23aÞ ð23bÞ
Let us first note that wðx; t þ sÞ ¼ esA=2 esB esA=2 wðx; tÞ þ Oðs2 Þ:
ð24Þ
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
225
To obtain the exact flows we must solve the problems 1 o2 w ot w ¼ ; 2 ox2
1 2 2 x þ U jwj w: ot w ¼ 2
ð25aÞ ð25bÞ
Eq. (25a) is easy R þ1to solve in Fourier space. Defining the Fourier transform as usual w^ðkÞ ¼ 1 eikx wðxÞ dx, Eq. (25a) reads 1 ot w^ ¼ k 2 w^; 2
ð26Þ
2 whose solution is w^ðtÞ ¼ w^ðt0 Þek ðtt0 Þ . Solution of Eq. (25b) is also straight2 1 2 forward wðt; xÞ ¼ wðt0 ; xÞeð2x þU jwj Þðtt0 Þ . Thus, we have the explicit form of the operators
esA ¼ F1 expðsk 2 ÞF;
ð27aÞ
esB ¼ esð2x
ð27bÞ
1 2 þU jwj2
Þ:
This scheme has several advantages. First, it is norm preserving and secondly its computational implementation using the fast Fourier transform (FFT) leads to very fast scheme. This is mainly due to the fact that the computational cost to compute one FFT on a grid of Nx points is OðNx log Nx Þ. 5.3. A pseudospectral method with Hermite basis functions In the case of the harmonically trapped NLS equation it is possible to design a pseudospectral method based on a different separation. If we split now the operators A¼
1 o2 1 2 þ x; 2 ox2 2
ð28aÞ
2
B ¼ U jwj ;
ð28bÞ
we may again use scheme (24). In this case to obtain the exact flows we must solve the problems 1 o2 w 1 2 x ot w ¼ þ w; ð29aÞ 2 ox2 2 2
ot w ¼ U jwj w:
ð29bÞ 2
The second one is solved easily to give wðx; tÞ ¼ eU jwj ðtt0 Þ wðx; t0 Þ. To solve the first one we just notice that the eigenfunctions and eigenvalues of the linear operator A are just
226
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235 2 =2
/n ðxÞ ¼ hn ðxÞex kn ¼ n þ
1 2
;
ð30aÞ ð30bÞ
with n 2 N and hn ðxÞ being the Hermite polynomials (normalized to one) [29]. These functions form a complete basis for L2 ðRÞ, basis any L2 P1thus usingxthis 2 =2 solution of (29a) may be written as wðx; tÞ ¼ k¼1 cn /n ðxÞe . The pseudospectral Pmethod proceeds by taking the truncated expression 2 M wðx; tÞ ’ k¼1 cn /n ðxÞex =2 so that our approximate solution for esA is given by 1 0 es=2 0 0 0 0 0 C B 0 e3s=2 0 0 0 0 C B C 5s=2 1 B 0 sA 0 e 0 0 0 e ¼ HM B ð31Þ C HM ; C B .. A @ . 0 0 0 0 0 eðMþ1=2Þs where H is the operator which obtains the truncated spectral decomposition of the function w on the orthonormal basis (30a), i.e. 0 1 0 1 ð/0 ; wÞ c0 B ð/1 ; wÞ C B c1 C B C B C HM w ¼ B ð32Þ C ¼ B .. C: .. @ A @ . A . cM ð/M ; wÞ When choosing a particular spatial grid the linear functional HM may be expressed in matrix form. In our case we have used the trapezoidal rule for computing numerically the integrals in Eq. (32) (i.e. the scalar products ð/j ; wÞ). Incidentally, we must say that for the integration of functions with fast asymptotic decay this formula is of very high order. When this is done the discrete expression of the operator is 0 1 2 2 2 h0 ðx1 Þex1 =2 h0 ðx2 Þex2 =2 h0 ðxNx ÞexNx =2 B h ðx Þex21 =2 h ðx Þex22 =2 h ðx Þex2Nx =2 C C; ð33Þ 1 1 1 2 1 Nx HM ¼ hB @ A 2 2 2 hM ðx1 Þex1 =2 hM ðx2 Þex2 =2 hM ðxNx ÞexNx =2 which operates over the vector wj ; j ¼ 1; . . . ; Nx . The inverse operator is given t by H1 M ¼ HM . Although these expressions look simple, its practical evaluation in the computer is not easy. It is well known that the evaluation of very high order polynomials (keep in mind that hM ðxÞ is a polynomial of degree M), leads to numerical instabilities due to roundoff errors and cancellations. This is also the case with Hermite polynomials. Although it is usually suggested that evaluation of polynomial by the Horner method aids in supressing this type of
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
227
instabilities, we have found that the most robust strategy to evaluate the polynomial on the lattice points xj consists of applying the recursion relation Hkþ1 ðxj Þ ¼ 2xj Hk ðxj Þ 2kHk1 ðxj Þ;
k 2 N;
ð34Þ
to obtain the expression of the nonnormalized Hermite polynomial. To use Eq. (34) one first evaluates the two lowest order Hermite polynomials H0 ðxÞ ¼ 1 and H1 ðxÞ ¼ x and then generates the values for all the Hk ðxÞ. Finally, to get an orthonormal basis ðhn ðxÞÞ, one mustpnormalize the resulting ffiffiffi 2 2 basis to one using the fact that ðHk ex =2 ; Hk ex =2 Þ ¼ p2n n!. 5.4. Computational complexity of the pseudospectral methods In principle both of the psudospectral methods presented here have similar computational costs per time step since they involve the computation of an integral transformation (OðNx2 Þ) plus the exponentials for advancing on time. In practice, the availability of the fast Fourier transformation for the evaluation of the Fourier transform, whose cost per step is OðNx log Nx Þ could seem to imply a priori that this method should be more efficient. However, the fact that both methods could require a different number of grid points Nx to obtain good results defers the question of which of them is superior to a later section.
6. Results We have computed the solutions of Eq. (2a) using the pseudospectral schemes here described and compared them against the exact time-dependent solutions based on the ground state presented in Section 2. We have studied the dependence of the error of the solutions on the relevant parameters of the method, i.e. the number of grid points Nx , the time-step s and the number of Hermite modes used in the expansion M (in the case of the Fourier scheme the number of Fourier modes used is the same as the grid size). Some results are summarized in the following figures. In all the cases studied we choose as computational domain x 2 ½40; 40 because in this region the solution is zero to machine precision on the boundaries for all the values of U which we have considered. In this way the effect of the zero boundary condition used in our simulation is, if not completely removed, at least minimized. All the simulations where performed with optimized (vectorized) MATLAB codes in personal computers and workstations. All times given are referred to the test machine, which was an IBM Intellistation 6849U equipped with a Pentium 4 1.4 GHz processor and 512 Mb RAM running MATLAB 6.1 under Windows XP. In all the numerical simulations to be presented in this section we have used as initial data /0 ðx 1Þ (i.e. the ground state displaced initially to X0 ¼ 1). We
228
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
have verified that the results are not essentially affected by this choice, e.g. for taking X0 ¼ 5 all the results are qualitatively unaffected. To evaluate the accuracy of the different schemes we have followed the evolution of the error during the simulation. We have considered long time simulations up to T ¼ 100, which corresponds to 16 full oscillation periods of the exact solution and thus allows to get an idea of the real stability properties of the numerical scheme. Typical behaviors are presented in Fig. 2. It is seen in Fig. 2 how, for this specific choice of parameters, the pseudospectral method with Hermite basis (to be abbreviated hereafter as ‘‘PH2 scheme’’) performs better than the pseudospectral method with trigonometric basis (PFFT2) and FD schemes (FD2). As discussed before, for the parameters used here, the FD scheme despite its computation time is the largest of all three schemes, fails to approximate the solution, being necessary much smaller time
Fig. 2. Results from typical simulations for the most representative schemes for U ¼ 2 up to T ¼ 100. Shown are (a)–(c) simulations of the PH2 scheme with m ¼ 20, Nx ¼ 512, s ¼ 0:04: (a) solution jwðx; tÞj2 for t ¼ 100, (b) error as a function of x for t ¼ 100, (c) evolution of the maximum error EðtÞ ¼ kwðx; tÞ wnj k1 . (d)–(f) Same as (a)–(c) but for the PFFT2 scheme with s ¼ 0:04 and Nx ¼ 4096. (g)–(i) Same as (a)–(c) but for the FD2 scheme with s ¼ 0:02; Nx ¼ 1024.
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
229
Fig. 3. Graph of computation time vs. error for the simulations presented in Fig. 1.
steps to get reasonable results. Moreover, as it can be seen from Fig. 3, not only the precision of the PH2 is good but also the computational performance of the scheme is the best of all the schemes analyzed in this paper. Because of the failure of FD schemes to approximate properly the solutions with reasonable computational cost (as it happened with the symplectic scheme) in what follows we will concentrate on the accuracy and performance of the pseudospectral schemes. We have analyzed the accuracy (Fig. 4) and computation time (Fig. 5) for both pseudospectral schemes and their dependence on the relevant parameters s, Nx and M (the later only for the PH2 method). First, in Fig. 4 we study the dependence of the error of the PH2 scheme as a function of the number of
Fig. 4. Dependence of the error of the PH2 scheme as a function of M with s ¼ 0:01; Nx ¼ 1024 for (a) U ¼ 2, (b) U ¼ 10, (c) U ¼ 100.
230
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Fig. 5. Computation time for the PH2 scheme as a function of M with s ¼ 0:01; Nx ¼ 1024. Since Nx is fixed, the dependence is approximately linear and the ‘‘zero’’ value of the computation time corresponds to the decomposition of initial data on Hermite basis and the creation of the appropriate matrices which are later used in the computation.
modes used in the expansion for our test problem. It is seen how the larger is U the more the modes which are necessary to approximate the error with maximum precision (the saturation of the curve) which is dictated by the minimization procedure (which in the case U ¼ 10 is more efficient). It is clearly seen how the U ¼ 100 is the most difficult situation for the numerical scheme, in fact for this combination of Nx ; s and for all the range of M values the error is only moderately small (Figs. 4(c) and 5). The analogous of the previous simulations for the case of the PFFT2 scheme is the dependence of the error with the number of grid points (which is equal to the number of Fourier modes). This analysis is presented in Fig. 6. In this case we do not plot the error for U ¼ 100 because it is very large for low Nx values and this choice of s. It is clear from Fig. 6 that the scheme is sensitive on the specific choice of Nx and that the highest precisions are achieved only for Nx as large as possible. To end this section, we present in Figs. 7 and 8, some selected simulations of the dependence of the error with s which is the only remaining relevant parameter. It can be seen how the larger the value of U the more sensitive is the numerical scheme to the choice of this parameter, specially the U ¼ 100 case requires a very small choice of s in order to get reliable results. As discussed above the pseudospectral method based on the Hermite expansion is more accurate than the FFT-based one, in fact the later scheme failed to obtain accurate results when U ¼ 100 for time steps down to Dt ¼ 5 105 .
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
231
Fig. 6. (a) and (b) Dependence of the error of the PFFT2 scheme as a function of Nx with s ¼ 0:01 for (a) U ¼ 2, (b) U ¼ 10. (c) Computation time for the previous cases as a function of Nx . The curve shows the expected superlinear (OðNx log Nx Þ) behavior.
Fig. 7. Dependence of the error of the PH2 scheme as a function of the time step s for m ¼ 100; X ¼ 40; X0 ¼ 1; Nx ¼ 4096. (a) U ¼ 2, (b) U ¼ 10, (c) U ¼ 100.
7. A Galerkin spectral scheme 7.1. The scheme From the previous subsections it is more or less clear that the basis hermite functions is specially suitable to approximate the solutions of this problem. One might wonder if this could be used as the basis for a nodeless spectral discretization of Galerkin type. The idea is to semidiscretize the problem by
232
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Fig. 8. Dependence of the error of the FFT2 scheme as a function of the time step s for X ¼ 40; X0 ¼ 1; Nx ¼ 4096. (a) U ¼ 2, (b) U ¼ 10.
P using the truncated expansion wðx; tÞ ’ M k¼1 ck ðtÞ/k ðxÞ. Inserting this expression into Eq. (2a) one finds M M X X dck ðtÞ 1 /k ðxÞ ¼ kþ ck ðtÞ/k ðxÞ i dt 2 k¼1 k¼1 þU
M X
ck ðtÞ/k ðxÞcl ðtÞ/l ðxÞcm ðtÞ/m ðxÞ:
ð35Þ
k;l;m¼1
Computing the scalar product by /j ðxÞ we get
M X dcj ðtÞ 1 ¼ jþ Ajklm ck ðtÞcl ðtÞcm ðtÞ; cj ðtÞ þ U i dt 2 k;l;m¼1
ð36Þ
R where Ajklm ¼ /j ðxÞ/k ðxÞ/l ðxÞ/m ðxÞ dx. Due to the fact that the basis functions are real this is a completely symmetric fourth order tensor. Eq. (36) is also a Hamiltonian system which may be obtained from M U X 1X 1 2 H ða; bÞ ¼ nþ Aklmn ðam an ak al an ðtÞ þ b2n ðtÞ þ 2 n¼1 2 4 k;l;m;n þ bk bl bm bn þ 2al am bk bn Þ;
ð37Þ
where c ¼ a þ ib. 7.2. Computational considerations The main a priori advantage of scheme (36) is that it is a fully Galerkin method and thus no spatial nodes are needed to obtain the solution (thus there is no need for any grid). In this sense, the boundary conditions are also automatically satisfied by the numerical solutions and there are no undesired effects arising from artificial boundary conditions.
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
233
However, there are several difficulties in applying directly this scheme. The first one is the computation of the tensor Ak;l;m;n . Although this is in principle an straightforward problem––it involves only evaluation of polynomial integrals multiplied by Gaussian functions––we must keep in mind that its computation involves integrals of products of four polynomials of high order. As discussed before this is a difficult problem to solve numerically (even harder than before because now we have a product of four polynomials that must be evaluated) and thus we have opted by computing the integrals using symbolic computation. Specifically we have used MAPLE. The computational cost of this operation increases exponentially with M. In fact, even in our test machine the computation of the tensor for M ¼ 20 took several hours of CPU time and the case M ¼ 30 took a couple of days to complete. As we have seen before, there are situations where M ¼ 50 or more modes must be kept. In those cases even the memory requirements for such a tensor are of the order of 50 Mb and each evaluation of the right hand side of Eq. (36) costs of the order of 5M 4 operations (i.e. for M ¼ 50, the number of operations to be done to evaluate once the RHS is about 3 107 ). Thus we see that the direct use of this scheme is computationally very expensive. But, even if one disregards the previous considerations with the hope that at least it could be possible to get a better than before approximation of the solutions of the problem, this is not true. This is so because the time discretization of Eq. (36) poses important problems. First of all explicit schemes do not preserve the hamiltonian structure of the problem and lead to computational instabilities. We have tried to solve these equations using both simple and advanced schemes such as the explicit midpoint rule and the Dormand– Prince variable step integrators with the result that there are divergences for times below our typical integration time T ¼ 100. Other possibilities such as implicit symplectic integrators could be used but at the cost of many evaluations of the right hand side of Eq. (36) per time step, which will definitely lead to nonpractical methods.
8. Conclusions In this work we have analyzed several schemes for the numerical solution of the HTNLSE. By using as test data exact solutions of the HTNLSE we have shown that several strategies which are appropriate for the integration of the usual NLSEs are not useful for these type of equations. We have proposed a pseudospectral scheme which is computationally very efficient and provides very good approximation of the solutions with long-time stability of the results. The philosophy here used may be applied to the construction of numerical schemes of higher order in time by just using the appropriate compositions of operators [38].
234
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Acknowledgements VMPG is partially supported by the Ministerio de Ciencia y Tecnologıa under grant BFM2000-0521. XL acknowledges support from the ‘‘Programa de profesores invitados’’ of Universidad de Castilla-La Mancha. We would like to acknowledge Luis V azquez (U. Complutense) for a critical reading of the manuscript.
References [1] A. Hasegawa, Optical Solitons in Fibers, Springer-Verlag, Berlin, 1989. [2] R.K. Dodd, J.C. Eilbeck, J.D. Gibbon, H.C. Morris, Solitons and Nonlinear Wave Equations, Academic Press, New York, 1982. [3] A.S. Davydov, Solitons in Molecular Systems, Reidel, Dordrecht, 1985. [4] C. Sulem, P. Sulem, The Nonlinear Schr€ odinger Equation, Springer, Berlin, 2000. [5] L. V azquez, L. Streit, V.M. Perez-Garcıa (Eds.), Nonlinear Schr€ odinger and Klein-Gordon systems: Theory and Applications, World Scientific, Singapore, 1996. [6] F. Dalfovo, S. Giorgini, L.P. Pitaevskii, S. Stringari, Theory of Bose–Einstein condensation in trapped gases, Rev. Mod. Phys. 71 (1999) 463–512. [7] V.M. Perez-Garcıa, H. Michinel, H. Herrero, Bose–Einstein solitons in highly asymmetric traps, Phys. Rev. A 57 (1998) 3837–3845. [8] T. Busch, J.R. Anglin, Dark-bright solitons in inhomogeneous Bose–Einstein condensates, Phys. Rev. Lett. 87 (2001) 010401:1–4. [9] T. Busch, J.R. Anglin, Motion of dark solitons in trapped Bose–Einstein condensates, Phys. Rev. Lett. 84 (2000) 2298–2301. [10] S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G.V. Shlyapnikov, M. Lewenstein, Dark solitons in Bose–Einstein condensates, Phys. Rev. Lett. 83 (1999) 5198–5201. [11] M. Kunze, T. K€ upper, V.K. Mezentsev, E.G. Shapiro, S. Turitsyn, Nonlinear solitary waves with gaussian tails, Physica D 128 (1999) 273–295. [12] M. Delfour, M. Fortin, G. Payre, Finite difference solutions of a nonlinear Schr€ odinger equation, J. Comp. Phys. 44 (1981) 277–288. [13] J.M. Sanz-Serna, J.G. Verwer, Conservative and nonconservative schemes for the solution of the nonlinear Schr€ odinger equation, IMA J. Numer. Anal. 6 (1986) 25–42. [14] L.S. Peranich, A finite difference scheme for solving a nonlinear Schr€ odinger equation with a linear damping term, J. Comp. Phys. 68 (1987) 501–505. [15] F. Zhang, V.M. Perez-Garcıa, L. Vazquez, Numerical simulation of nonlinear Schr€ odinger systems: a new conservative scheme, Appl. Math. Comput. 71 (1995) 164–177. [16] Y. Tang, L. Vazquez, F. Zhang, V.M. Perez-Garcıa, Symplectic methods for the nonlinear Schr€ odinger equation, Comp. Math. Appl. 32 (1996) 73–83. [17] Y. Tang, V.M. Perez-Garcıa, L. Vazquez, Symplectic methods for the Ablowitz–Ladik model, Appl. Math. Comput. 82 (1997) 17–38. [18] O. Karakashian, C. Makridakis, A space-time finite element method for the nonlinear Schr€ odinger equation: the discontinuous Galerkin method, Math. Comput. 67 (1998) 479–499. [19] Y. Tourigny, J. Morris, An investigation into the effect of product approximation in the numerical solution of the cubic nonlinear Schr€ odinger equation, J. Comp. Phys. 76 (1988) 103– 130. [20] B.M. Herbst, F. Varadi, M.J. Ablowitz, Symplectic methods for the nonlinear Schr€ odinger equation, Math. Comput. Simul. 37 (1994) 353–369.
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
235
[21] R. MacLahan, Symplectic integration of Hamiltonian wave equations, Numer. Math. 66 (1994) 465–492. [22] D. Pathria, J.Ll. Morris, Pseudo-spectral solution of nonlinear Schr€ odinger equations, J. Comp. Phys. 87 (1990) 108–125. [23] J.A.C. Weideman, B.M. Herbst, Split-step methods for the solution of the nonlinear Schr€ odinger equation, SIAM J. Numer. Anal. 23 (1986) 485–507. [24] A.D. Bandrauk, H. Shen, High-order split-step exponential methods for solving coupled nonlinear Schr€ odinger equations, J. Phys. A Math. Gen. 27 (1994) 7147–7155. [25] B.M. Herbst, J.Ll. Morris, A.R. Mitchell, Numerical experience with the nonlinear Schr€ odinger equation, J. Comput. Phys. 60 (1985) 282–305. [26] T.R. Taha, M.J. Ablowitz, Analytical and numerical aspects of certain nonlinear evolution equations II. Numerical nonlinear Schr€ odinger equation, J. Comp. Phys. 55 (1984) 203–230. [27] G. Agrawall, Nonlinear Fiber Optics, Academic Press, New York, 1995. [28] S. Hairer, Nordsett, Wanner, in: Solving Ordinary Differential Equations, vol. 1, Springer, Berlin, 1996. [29] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions, Dover, New York, 1964. [30] L. Erbe, M. Tang, Uniqueness of positive radial solutions of Du þ f ðjxj; uÞ ¼ 0, Diff. Int. Eq. 11 (1995) 725–743. [31] Z.-Q. Wang, Existence and symmetry of multi-bump solutions for nonlinear Schr€ odinger equations, J. Diff. Eq. 159 (1999) 102–137. [32] J.J. Garcıa-Ripoll, V.M. Perez-Garcıa, V. Vekslerchik, Construction of exact solutions by spatial translations in inhomogeneous nonlinear Schr€ odinger equations, Phys. Rev. E 64 (2001) 056602-1:6. [33] J.J. Garcıa-Ripoll, V.M. Perez-Garcıa, Minimization of Schr€ odinger functionals using sobolev gradients: Applications to quantum mechanics and nonlinear optics, SIAM J. Sci. Comp. 23 (2001) 1315–1324. [34] P.Z. Ding, C.X. Wu, Y.K. Mu, Y.X. Li, M.X. Jin, Square-preserving and symplectic structure and scheme for quantum system, Chin. Phys. Lett. 13 (1996) 245–248. [35] Z.Y. Zhou, P.Z. Ding, S.P. Pan, Study of a symplectic scheme for the time evolution of an atom in an external field, J. Kor. Phys. Soc. 32 (1998) 417–424. [36] X. Liu, X. Liu, Z. Zhou, P. Ding, S. Pan, Numerical solution of one-dimensional timeindependent Schr€ odinger equation by using symplectic schemes, Int. J. Quant. Chem. 79 (2000) 343–349. [37] S.K. Gray, D.E. Manolopoulos, Symplectic integrators tailored to the time-dependent Schr€ odinger equation, J. Chem. Phys. 104 (1996) 7099–7112. [38] S. Blanes, P.C. Moan, Practical symplectic partitioned Runge-Kutta and Runge–Kutta– Nystrom methods, J. Comput. Appl. Math. 142 (2002) 313–330.
Numerical methods for the simulation €dinger systems of trapped nonlinear Schro Vıctor M. Perez-Garcıa a
a,*
, Xiao-yan Liu
b
Departamento de Matem aticas, Escuela T ecnica Superior de Ingenieros Industriales, Universidad de Castilla-La Mancha, Avenida de Camilo Jose Cela, s/n, 13071 Ciudad Real, Spain b Department of Mathematics, Northeast Normal University, Changchun 130024, China
Abstract We propose, analyze and compare the efficiency and accuracy of different numerical schemes for the solution of the nonlinear Schr€ odinger equation with a trapping potential. In particular we study schemes of finite difference, pseudospectral and spectral types for the space discretization together with explicit symplectic, multistep, split-step and standard variable-step integrators to solve the time evolution. All of these schemes are evaluated comparatively and some recommendations based on their accuracy and computational efficiency are made. Ó 2002 Elsevier Inc. All rights reserved. Keywords: Pseudospectral schemes; Nonlinear Schr€ odinger equations; Finite difference schemes; Symplectic schemes
1. Introduction The nonlinear Schr€ odinger equation (NLSE) is one of the most important models of mathematical physics, with applications to different fields such as plasma physics, nonlinear optics, water waves, biomolecule dynamics and many other fields (see e.g. Refs. [1–5]). In many of those fields the equation appears as an asymptotic limit for a slowly varying dispersive wave envelope propagating in a nonlinear medium. A new burst of interest on this equation has been motivated by the recent achievement of Bose–Einstein condensation *
Corresponding author. E-mail address: [email protected] (V.M. Perez-Garcıa).
0096-3003/02/$ - see front matter Ó 2002 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(02)00402-2
216
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
using ultracold neutral bosonic gases. In this case the NLSE is obtained as a mean field model for the dynamics of the quantum gas [6], which, in reduced units reads 1 1 2 iWt ¼ DW þ x2 W þ U jWj W; 2 2
ð1Þ
where W : X R2 ! C, D ¼ ðo2 =ox21 Þ þ þ ðo2 =ox2n Þ, x2 ¼ x21 þ þ x2n , n ¼ 1; 2; 3 being the spatial dimensionality of the problem and U 2 R. In BEC problems one is interested on the solution of the Cauchy problem for initial data in L2 in unbounded domains. From now on we will restrict our attention to the case n ¼ 1, i.e. the problem to be solved is i
ow 1 o2 w 1 2 2 ¼ þ x w þ U jwj w; ot 2 ox2 2
wðx; 0Þ ¼ w0 ðxÞ; lim wðx; tÞ ¼ 0:
x! 1
ð2aÞ ð2bÞ ð2cÞ
The effect of the potential term V ðxÞ ¼ x2 =2 is to localize the solutions near x ¼ 0. This is why this term is called a ‘‘trap’’. In what follows we will refer to this problem as harmonically trapped nonlinear Schr€ odinger equation (HTNLSE). This one-dimensional situation is relevant at least for applications in the fields of Bose–Einstein condensation [7–10] and nonlinear optics [11] and has received attention in the last years. In addition to its specific interest, the analysis of the one-dimensional problem is the first step for the analysis of more complicated multidimensional configurations. Concerning the HTNLSE, it is known that this equation has at least two conserved quantities which are Z 2 N ðwÞ ¼ jwðx; tÞj dx; ð3aÞ H EðwÞ ¼
Z 2 1 ~ 1 1 2 4 rw þ x2 jwj þ U jwj dx: 2 2 2
ð3bÞ
In our paper we will consider solutions with N ðwÞ ¼ 1 since other values may be scaled into the nonlinear coefficient U. Despite the many applied works which rely on numerical simulations of the HTNLSE, this equation has not received specific attention in the mathematical literature. This is in contrast with the fact that the NLSE without trapping potential is a very well studied problem for which different numerical schemes have been analyzed. For instance, the usual NLS equation has been solved using finite differences (FDs) in space plus time variables [12–15] (in some of
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
217
those cases the time-discretization was conservative), FD in space plus symplectic integration in time [16,17], finite elements [18–21] and pseudospectral methods for the discretization of the spatial variable [22–24]. Of remarkable success are the methods based on split-step decompositions and FFT [25,26], which have been adopted as standard integrators by a relevant part of the scientific community (see e.g. [27]). In this paper we consider several numerical schemes for Eqs. (2a)–(2c) and compare their accuracy and performance by using numerical simulations. To our knowledge and despite the practical importance of these problems, this analysis has not been considered previously in the scientific computing literature. In particular for the space discretization we will consider second order FD schemes, spectral and pseudospectral methods with different basis functions. For the time evolution we use several methods such as split-step, centered (conservative) FDs, modified Euler method and explicit symplectic integrators.
2. Exact solutions of the HTNLSE The simplest solutions of Eqs. (2a)–(2c) are the so-called solitary waves or stationary solutions, which are of the form wl ðx; tÞ ¼ /l ðxÞeilt ;
ð4Þ
where l 2 R and /l ðxÞ : R ! R are to be determined from the equations 1 o2 1 2 2 x l/l ðxÞ ¼ þ þ U j/ j ð5aÞ /l ðxÞ; l 2 ox2 2 lim /l ðxÞ ¼ 0:
x! 1
ð5bÞ
In the linear case (U ¼ 0) this problem has a numerable set of stationary real 2 solutions /k ðxÞ ¼ Hk ðxÞex =2 , with k 2 N and Hk ðxÞ being the Hermite polynomials [29] and lk ¼ k þ 1=2. For the nonlinear repulsive case (U > 0) it is has been proven in limited situations [11] and conjectured for the general case that the situation is the same, i.e. that for fixed N ¼ 1 there is a numerable set of solutions fln ; /ln g. In this paper we will be interested on the ground state solution of Eqs. (5a) and (5b) which is the stationary solution with smallest value of the eigenvalue l for a fixed value of the L2 -norm. It is known that the ground state of Eqs. (5a) and (5b) is p a ffiffiunique positive radial solution [30] which for r 0 behaves as ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /l ðrÞ ð1= uÞ ðl ð1=2Þr2 Þþ and for r ! 1 is bounded by /l ðrÞ 6 2 ½C0 þ CðlÞrer =2 [29]. General properties of the solutions of Eqs. (5a) and (5b) have been studied in a more general context [31].
218
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Although many properties of the ground state are known and there are several asymptotic formulae available for its approximation, a closed analytical form giving its shape is not available, thus the most accurate way of obtaining it is by means of numerical methods. This turns out to be a difficult problem since the approximation of /0 ðxÞ requires the solution of a nonlinear eigenvalue problem (l0 is also unknown) for which simple shooting methods are not useful. This is because the asymptotic behavior for large x values is a delicate balance between the increasing x2 term and the other ones present in the equation. Thus, to obtain accurate approximations of the ground states, which will be later used as initial data we have used the so called Sobolev-gradients method as described in Ref. [33]. Following this scheme we have obtained the ground state for three values of the nonlinearity U ¼ 2; 10; 100 and N ð/0 Þ ¼ 1 after a relaxation time of T ¼ 100:000 with s ¼ 0:2, on the uniformly spaced grid Xh ¼ ½Xmax ; Xmax þ h; . . . ; Xmax with Xmax ¼ 40 and Nx ¼ 8192 points. Working in this way we get /0 ðxÞ with approximate maximum errors of the order of 104 –105 (this depends on the specific value of U), which is enough for our present study. Through this work when approximations of /0 ðxÞ for points xj 62 Xh have been necessary we have used Lagrange interpolation with cubic polynomials to get those values from their four nearest neighbors. The dynamics of the ground state although could reveal some simple properties of the numerical schemes to be later used is in a certain sense trivial. To build more interesting solutions of Eqs. (2a)–(2c) we will use a property of stationary solutions of HTNLSE which allows to build time dependent solutions from stationary ones [32]. For our present purposes we present a limited formulation which is that given a stationary solution of (5a) and (5b) then one may build infinitely many new time dependent solutions of the form wðx; t; l; X0 ; V0 Þ ¼ /l ðx X ðtÞ; tÞeiðltþxX ðtÞÞ ; _
ð6Þ
where X ðtÞ satisfies X€ þ X ¼ 0 with arbitrary initial data X ðt0 Þ ¼ X0 ; X_ ðt0 Þ ¼ V0 . The existence of these nontrivial solutions provides a way to study quantitatively the errors of the different numerical methods.
3. Explicit symplectic schemes 3.1. The scheme Although the simplest explicit schemes are usually shown to be unstable, there are explicit symplectic schemes which do not display such unwanted
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
219
behavior. This is why these schemes, which are simple to implement have been used in several applications to linear problems [34–36]. Let us first define the discrete grid Xh ¼ fx1 ; . . . ; xnx g ¼ fX ; X þ h; . . . ; X h; X g with Nx elements. We rewrite our equations by defining wj ðtÞ ¼ wðxj ; tÞ and decomposing wj ðtÞ ¼ aj ðtÞ þ ibj ðtÞ. The nonlinear equation (2a) can be written on the discrete grid in the following way oa ¼ Hb; ot
ð7aÞ
ob ¼ Ha; ot
ð7bÞ >
>
where a ¼ ða1 ; a2 ; . . . ; an Þ ; b ¼ ðb1 ; b2 ; . . . ; bn Þ ; H ¼ HK þ HV where HK stands for the effect of the operator ðo2 =ox2 Þ on the particular grid point and HV is the diagonal matrix 01 2 1 x þ U ða21 þ b21 Þ 2 1 1 2 2 2 B C x þ U ða2 þ b2 Þ 2 2 B C HV ¼ B C: . .. @ A 1 2 2 2 x þ U ðan þ bn Þ 2 n ð8Þ Eqs. (7a) and (7b) form a hamiltonian system which may be accurately integrated by using symplectic schemes. For instance, a second order explicit scheme is given by A ¼ an þ p1 H sbn ;
ð9aÞ
n
B ¼ b þ q1 H sA;
ð9bÞ
anþ1 ¼ A þ p2 H sB;
ð9cÞ
bnþ1 ¼ B þ q2 H sanþ1 ; ð9dÞ pffiffiffi pffiffiffi pffiffiffi pffiffiffi where p1 ¼ 1 ð1= 2Þ; p2 ¼ ð1= 2Þ; q1 ¼ ð1= 2Þ; q2 ¼ 1 þ ð1= 2Þ are the coefficients defining the method. The matrix form of this scheme is
nþ1 n M11 M12 a a ¼ ; ð10Þ M21 M22 bn bnþ1 where
M¼
¼
1 q2 H s
0 1
1 0
p2 H s 1
1 þ q1 p2 H 2 s2 ðq1 þ q2 ÞH s þ q1 q2 p2 H 3 s3
1 0 q1 H s 1
1 0
p1 H s 1
ðp1 þ p2 ÞH s þ p1 p2 q1 H 3 s3 : 1 þ ðp1 q1 þ p1 q2 þ p2 q2 ÞH 2 s2 þ p1 p2 q1 q2 H 4 s4
220
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
To define a fully discrete problem we must discretize the spatial derivatives. Using FDs and zero Dirichlet boundary conditions on the boundaries, the operator HV may be written into the standard form 1 0 2 1 0 0 0 B 1 2 1 0 0 C B .. C .. C B 1 B 0 . C . ð11Þ HK ¼ 2 B . C: .. 2h B .. . 0 C C B @ 0 0 1 2 1 A 0
0
0
1
2
3.2. Stability of symplectic schemes Eq. (10) is equivalent to a linear mapping Z kþ1 ¼ MZ k where k ¼ 0; 1; 2; . . ., > corresponds to t ¼ t0 ; t0 þ s; t0 þ 2s; . . ., and Z ¼ ða; bÞ . The stability of such a mappings depends on the eigenvalues of M. In particular, stable behavior occurs if these eigenvalues have modulus 6 1. The eigenvalues of M, defining x0 ¼ sh, are given by ! M11 ðx0 Þ k M12 ðx0 Þ det ¼ k2 ðM11 ðx0 Þ þ M22 ðx0 ÞÞk M21 ðx0 Þ M22 ðx0 Þ k þ M11 ðx0 ÞM22 ðx0 Þ M12 ðx0 ÞM21 ðx0 Þ ¼ 0: Here we have used the definitions M11 ðx0 Þ ¼ 1 þ q1 p2 x02 ; M12 ðx0 Þ ¼ p1 x0 þ p2 x0 þ p1 p2 q1 x03 ; M21 ðx0 Þ ¼ q1 x0 þq2 x0 þq1 q2 p2 x03 ; M22 ðx0 Þ ¼ 1þp1 q1 x02 þp1 q2 x02 þ 02 04 p2 q2 x þ p1 p2 q1 q2 x . The matrix M is symplectic and then jMj ¼ 1, i.e.
M11 ðx0 Þ M12 ðx0 Þ det ¼ M11 ðx0 ÞM22 ðx0 Þ M12 ðx0 ÞM21 ðx0 Þ ¼ 1: ð12Þ M21 ðx0 Þ M22 ðx0 Þ Thus, we get k2 ðM11 ðx0 Þ þ M22 ðx0 ÞÞk þ 1 ¼ 0:
ð13Þ
In order to guarantee jkj 6 1, we need to satisfy jk1 þ k2 j ¼ jM11 ðx0 Þ þ M22 ðx0 Þj 6 jk1 j þ jk2 j 6 2;
ð14Þ
thus we find that the stability condition is kx0 k 6 2:26 [37]. That is to say if kx0 k 6 xs ¼ 2:26 using an appropriate matrix norm then any eigenvalue jkj 6 1. Then xs ¼ 2:26 is the stability limit for the 2-order explicit symplectic scheme. Now we can estimate a maximum stable time step ss from xs by using the fact that for any natural norm the spectral radius of the matrix qðH Þ 6 kH k and thus
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
ss ¼
xs : qðH Þ
221
ð15Þ
For instance, if we consider the FD discretization matrix for HK and using Chebyshev approximation we get qðH Þ <
2 1 U 2 þ max jxj j þ : h2 2 j¼1;...;nx h
ð16Þ
This equation provides a bound for the maximum time step which may be used for the stable integration of Eq. (2a) by the explicit second order explicit scheme ss 6
2 h2
xs : þ 12 X 2 þ Uh
ð17Þ
It is remarkable that not only the usual term 2=h2 due to the FD discretization of the problem is present but also two terms due to the nonlinear effects and boundary effects which also restrict the time step. In fact, for large integration regions (X) and values of the nonlinear parameter (U) these terms would be responsible for the stability constraints of the scheme, a result which is very different for that of the free (nontrapped) LSE. The fact that the nonlinear term contributes with a term of type U =h is a first indication of the relevant role of the nonlinear parameter and its possible connection with unstable behavior, a connection which will be more clear later. This type of stability restrictions appear also when higher order schemes are used. The only difference is on the particular value of the parameter xs , which should be replaced by the corresponding value of the stability condition of the scheme for the free problem. 3.3. Results We have computed the solutions of Eq. (2a) using the explicit symplectic scheme (9a)–(9d) (similar simulations where done with a fourth order explicit scheme with very similar qualitative results). We have used as initial state the displaced ground state /0 ðx 1Þ obtained above on a grid with Xmax ¼ 40 and Nx ¼ 4096 points. All simulations were repeated for U ¼ 2, 10 and 100. The corresponding stability limits for these values of the nonlinear parameter obtained from Eq. (15) are ss ðU ¼ 2Þ ¼ 3:68 104 , ss ðU ¼ 10Þ ¼ 3:45 104 and ss ðU ¼ 100Þ ¼ 2:03 104 , respectively. It can be seen in Fig. 1 that when U ¼ 2 taking s ¼ 3 104 gives reasonable results (Fig. 1(a)–(c)). However when U ¼ 10 and s ¼ 3 104 the error in the solution becomes already very large for t ’ 20 (see Fig. 1(d)–(f)). For these values of the nonlinear parameter, despite they lay above the stability limit the errors accumulate in such a way that for moderate times the approximate solution is not good. Taking smaller time steps, e.g. s ¼ 5 105 as shown in
222
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Fig. 1. Results from typical simulations using the second order explicit symplectic schemes for U ¼ 2 and U ¼ 10, respectively, up to T ¼ 100. Shown are (a)–(c) simulations with Nx ¼ 4096, s ¼ 0:0003: (a) solution jwðx; tÞj2 for t ¼ 10, (b) error jeðx; tÞj2 for t ¼ 100, (c) evolution of the maximum error EðtÞ ¼ kwðx; tÞ wnj k1 . (d)–(f) Same as (a)–(c) but for U ¼ 10 with s ¼ 0:0003. (g)– (i) Same as (d)–(f) but for U ¼ 10 and s ¼ 0:00005.
Fig. 1(g)–(i) solves this problem at the cost of a higher computational cost. We have verified that for U ¼ 100 very small time steps must be used to get reasonable results, so small in fact that this scheme becomes unusable for the numerical approximation of the solution for long time simulations.
4. A conservative finite difference scheme for the HTNLS 4.1. Introduction The lack of stability of the symplectic scheme used in Section 3.2 was partially due to the explicit nature of the method which is known to impose strict stability constraints on these type of problems. Although there is the possibility
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
223
of using implicit symplectic FD schemes in the space and time variables such as those used in [16,17] the computational implementation of those schemes is computationally expensive mainly due to the need for solving a system of nonlinear equations at each time step. This is why our next step will be to consider an extension of a conservative FD scheme which has been used successfully in several contexts [15] and which is linearly implicit and leads to very efficient computational implementations. 4.2. The finite difference conservative scheme A natural extension of the linearly implicit, conservative numerical scheme of Ref. [15] is the following " nþ1 # nþ1 nþ1 n1 n1 n1 wjnþ1 wjn1 1 wjþ1 2wj þ wj1 wjþ1 2wj þ wj1 i ¼ þ 2 2s 2h2 2h2 nþ1
1 wj þ xj 2
þ wjn1 wnþ1 þ wjn1 n 2 j þ U jwj j : 2 2
ð18Þ
If, as in the previous case, zero Dirichlet boundary conditions are used and introducing aj ¼ 1=h2 þ i=s x2j =2 and the matrix 1 0 1 0 0 a1 þ U jw1 j2 2h2 C B 1 a2 þ U jw2 j2 2h12 0 C B 2h2 C B .. .. C B C B 0 . . A¼B C; . . C B . . . . 0 C B C B 2 n 1 1 a þ U jw j 0 A @ Nx 1 Nx 1 2h2 2h2 2 n 1 0 0 aNx þ U jwNx j 2h2 ð19Þ
the numerical scheme may be written in matrix form as Awnþ1 ¼ A wn1 ;
ð20Þ
where w ¼ ðw1 ; . . . ; wM ÞT . Following the methods developed in Ref. [15], it is straightforward to prove that this scheme is stable and convergent and that no artificial blow-up is possible. The main advantages of this scheme are that it is linearly implicit and conservative. This scheme has been successfully applied when there is no potential to different nonlinear Schr€ odinger problems. By means of numerical simulations we have found that in spite of its stability properties the scheme does not perform well for long-time simulations of Eq. (2a). Although the accuracy of the scheme is not bad for short time simulations in agreement with the Oðs2 þ h2 Þ truncation error of the scheme, it
224
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
does not work properly for long simulations unless a very small time step is chosen, with a behavior similar to that of the symplectic scheme. An example of this behavior will be presented in a later section.
5. Split-step schemes for the HTNLS 5.1. Introduction The fact that simple numerical schemes such as those presented up to now require very small time steps to provide accurate results in long-time integrations of the HTNLS equation has motivated the search for other type of schemes of different nature such as those which will be presented in this section. Split-step schemes for integrating partial differential equations are based on the observation that many problems may be decomposed into exactly solvable parts and the exact solution of the full problem can be approximated as a composition of the solutions of the individual problems. For instance the solution of a partial differential equation of the type ot uðx; tÞ ¼ N ðt; x; u; ox ; . . . ; Þu ¼ ðA þ BÞu;
ð21Þ
may be approximated from the exact solutions of the subproblems ot u ¼ Au;
ð22aÞ
ot u ¼ Bu:
ð22bÞ
This way of obtaining solutions by partitioning of the original problem has been studied both in the field of the approximation of ordinary [28] and partial differential equations. In the case of the ordinary linear or NLSE this type of schemes have been successfully applied previously [26]. 5.2. Fourier pseudospectral scheme Let us first split the operator in Eq. (2a) following the standard approach, which consists in taking A¼
1 o2 ; 2 ox2
1 B ¼ x2 þ U jwj2 : 2
ð23aÞ ð23bÞ
Let us first note that wðx; t þ sÞ ¼ esA=2 esB esA=2 wðx; tÞ þ Oðs2 Þ:
ð24Þ
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
225
To obtain the exact flows we must solve the problems 1 o2 w ot w ¼ ; 2 ox2
1 2 2 x þ U jwj w: ot w ¼ 2
ð25aÞ ð25bÞ
Eq. (25a) is easy R þ1to solve in Fourier space. Defining the Fourier transform as usual w^ðkÞ ¼ 1 eikx wðxÞ dx, Eq. (25a) reads 1 ot w^ ¼ k 2 w^; 2
ð26Þ
2 whose solution is w^ðtÞ ¼ w^ðt0 Þek ðtt0 Þ . Solution of Eq. (25b) is also straight2 1 2 forward wðt; xÞ ¼ wðt0 ; xÞeð2x þU jwj Þðtt0 Þ . Thus, we have the explicit form of the operators
esA ¼ F1 expðsk 2 ÞF;
ð27aÞ
esB ¼ esð2x
ð27bÞ
1 2 þU jwj2
Þ:
This scheme has several advantages. First, it is norm preserving and secondly its computational implementation using the fast Fourier transform (FFT) leads to very fast scheme. This is mainly due to the fact that the computational cost to compute one FFT on a grid of Nx points is OðNx log Nx Þ. 5.3. A pseudospectral method with Hermite basis functions In the case of the harmonically trapped NLS equation it is possible to design a pseudospectral method based on a different separation. If we split now the operators A¼
1 o2 1 2 þ x; 2 ox2 2
ð28aÞ
2
B ¼ U jwj ;
ð28bÞ
we may again use scheme (24). In this case to obtain the exact flows we must solve the problems 1 o2 w 1 2 x ot w ¼ þ w; ð29aÞ 2 ox2 2 2
ot w ¼ U jwj w:
ð29bÞ 2
The second one is solved easily to give wðx; tÞ ¼ eU jwj ðtt0 Þ wðx; t0 Þ. To solve the first one we just notice that the eigenfunctions and eigenvalues of the linear operator A are just
226
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235 2 =2
/n ðxÞ ¼ hn ðxÞex kn ¼ n þ
1 2
;
ð30aÞ ð30bÞ
with n 2 N and hn ðxÞ being the Hermite polynomials (normalized to one) [29]. These functions form a complete basis for L2 ðRÞ, basis any L2 P1thus usingxthis 2 =2 solution of (29a) may be written as wðx; tÞ ¼ k¼1 cn /n ðxÞe . The pseudospectral Pmethod proceeds by taking the truncated expression 2 M wðx; tÞ ’ k¼1 cn /n ðxÞex =2 so that our approximate solution for esA is given by 1 0 es=2 0 0 0 0 0 C B 0 e3s=2 0 0 0 0 C B C 5s=2 1 B 0 sA 0 e 0 0 0 e ¼ HM B ð31Þ C HM ; C B .. A @ . 0 0 0 0 0 eðMþ1=2Þs where H is the operator which obtains the truncated spectral decomposition of the function w on the orthonormal basis (30a), i.e. 0 1 0 1 ð/0 ; wÞ c0 B ð/1 ; wÞ C B c1 C B C B C HM w ¼ B ð32Þ C ¼ B .. C: .. @ A @ . A . cM ð/M ; wÞ When choosing a particular spatial grid the linear functional HM may be expressed in matrix form. In our case we have used the trapezoidal rule for computing numerically the integrals in Eq. (32) (i.e. the scalar products ð/j ; wÞ). Incidentally, we must say that for the integration of functions with fast asymptotic decay this formula is of very high order. When this is done the discrete expression of the operator is 0 1 2 2 2 h0 ðx1 Þex1 =2 h0 ðx2 Þex2 =2 h0 ðxNx ÞexNx =2 B h ðx Þex21 =2 h ðx Þex22 =2 h ðx Þex2Nx =2 C C; ð33Þ 1 1 1 2 1 Nx HM ¼ hB @ A 2 2 2 hM ðx1 Þex1 =2 hM ðx2 Þex2 =2 hM ðxNx ÞexNx =2 which operates over the vector wj ; j ¼ 1; . . . ; Nx . The inverse operator is given t by H1 M ¼ HM . Although these expressions look simple, its practical evaluation in the computer is not easy. It is well known that the evaluation of very high order polynomials (keep in mind that hM ðxÞ is a polynomial of degree M), leads to numerical instabilities due to roundoff errors and cancellations. This is also the case with Hermite polynomials. Although it is usually suggested that evaluation of polynomial by the Horner method aids in supressing this type of
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
227
instabilities, we have found that the most robust strategy to evaluate the polynomial on the lattice points xj consists of applying the recursion relation Hkþ1 ðxj Þ ¼ 2xj Hk ðxj Þ 2kHk1 ðxj Þ;
k 2 N;
ð34Þ
to obtain the expression of the nonnormalized Hermite polynomial. To use Eq. (34) one first evaluates the two lowest order Hermite polynomials H0 ðxÞ ¼ 1 and H1 ðxÞ ¼ x and then generates the values for all the Hk ðxÞ. Finally, to get an orthonormal basis ðhn ðxÞÞ, one mustpnormalize the resulting ffiffiffi 2 2 basis to one using the fact that ðHk ex =2 ; Hk ex =2 Þ ¼ p2n n!. 5.4. Computational complexity of the pseudospectral methods In principle both of the psudospectral methods presented here have similar computational costs per time step since they involve the computation of an integral transformation (OðNx2 Þ) plus the exponentials for advancing on time. In practice, the availability of the fast Fourier transformation for the evaluation of the Fourier transform, whose cost per step is OðNx log Nx Þ could seem to imply a priori that this method should be more efficient. However, the fact that both methods could require a different number of grid points Nx to obtain good results defers the question of which of them is superior to a later section.
6. Results We have computed the solutions of Eq. (2a) using the pseudospectral schemes here described and compared them against the exact time-dependent solutions based on the ground state presented in Section 2. We have studied the dependence of the error of the solutions on the relevant parameters of the method, i.e. the number of grid points Nx , the time-step s and the number of Hermite modes used in the expansion M (in the case of the Fourier scheme the number of Fourier modes used is the same as the grid size). Some results are summarized in the following figures. In all the cases studied we choose as computational domain x 2 ½40; 40 because in this region the solution is zero to machine precision on the boundaries for all the values of U which we have considered. In this way the effect of the zero boundary condition used in our simulation is, if not completely removed, at least minimized. All the simulations where performed with optimized (vectorized) MATLAB codes in personal computers and workstations. All times given are referred to the test machine, which was an IBM Intellistation 6849U equipped with a Pentium 4 1.4 GHz processor and 512 Mb RAM running MATLAB 6.1 under Windows XP. In all the numerical simulations to be presented in this section we have used as initial data /0 ðx 1Þ (i.e. the ground state displaced initially to X0 ¼ 1). We
228
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
have verified that the results are not essentially affected by this choice, e.g. for taking X0 ¼ 5 all the results are qualitatively unaffected. To evaluate the accuracy of the different schemes we have followed the evolution of the error during the simulation. We have considered long time simulations up to T ¼ 100, which corresponds to 16 full oscillation periods of the exact solution and thus allows to get an idea of the real stability properties of the numerical scheme. Typical behaviors are presented in Fig. 2. It is seen in Fig. 2 how, for this specific choice of parameters, the pseudospectral method with Hermite basis (to be abbreviated hereafter as ‘‘PH2 scheme’’) performs better than the pseudospectral method with trigonometric basis (PFFT2) and FD schemes (FD2). As discussed before, for the parameters used here, the FD scheme despite its computation time is the largest of all three schemes, fails to approximate the solution, being necessary much smaller time
Fig. 2. Results from typical simulations for the most representative schemes for U ¼ 2 up to T ¼ 100. Shown are (a)–(c) simulations of the PH2 scheme with m ¼ 20, Nx ¼ 512, s ¼ 0:04: (a) solution jwðx; tÞj2 for t ¼ 100, (b) error as a function of x for t ¼ 100, (c) evolution of the maximum error EðtÞ ¼ kwðx; tÞ wnj k1 . (d)–(f) Same as (a)–(c) but for the PFFT2 scheme with s ¼ 0:04 and Nx ¼ 4096. (g)–(i) Same as (a)–(c) but for the FD2 scheme with s ¼ 0:02; Nx ¼ 1024.
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
229
Fig. 3. Graph of computation time vs. error for the simulations presented in Fig. 1.
steps to get reasonable results. Moreover, as it can be seen from Fig. 3, not only the precision of the PH2 is good but also the computational performance of the scheme is the best of all the schemes analyzed in this paper. Because of the failure of FD schemes to approximate properly the solutions with reasonable computational cost (as it happened with the symplectic scheme) in what follows we will concentrate on the accuracy and performance of the pseudospectral schemes. We have analyzed the accuracy (Fig. 4) and computation time (Fig. 5) for both pseudospectral schemes and their dependence on the relevant parameters s, Nx and M (the later only for the PH2 method). First, in Fig. 4 we study the dependence of the error of the PH2 scheme as a function of the number of
Fig. 4. Dependence of the error of the PH2 scheme as a function of M with s ¼ 0:01; Nx ¼ 1024 for (a) U ¼ 2, (b) U ¼ 10, (c) U ¼ 100.
230
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Fig. 5. Computation time for the PH2 scheme as a function of M with s ¼ 0:01; Nx ¼ 1024. Since Nx is fixed, the dependence is approximately linear and the ‘‘zero’’ value of the computation time corresponds to the decomposition of initial data on Hermite basis and the creation of the appropriate matrices which are later used in the computation.
modes used in the expansion for our test problem. It is seen how the larger is U the more the modes which are necessary to approximate the error with maximum precision (the saturation of the curve) which is dictated by the minimization procedure (which in the case U ¼ 10 is more efficient). It is clearly seen how the U ¼ 100 is the most difficult situation for the numerical scheme, in fact for this combination of Nx ; s and for all the range of M values the error is only moderately small (Figs. 4(c) and 5). The analogous of the previous simulations for the case of the PFFT2 scheme is the dependence of the error with the number of grid points (which is equal to the number of Fourier modes). This analysis is presented in Fig. 6. In this case we do not plot the error for U ¼ 100 because it is very large for low Nx values and this choice of s. It is clear from Fig. 6 that the scheme is sensitive on the specific choice of Nx and that the highest precisions are achieved only for Nx as large as possible. To end this section, we present in Figs. 7 and 8, some selected simulations of the dependence of the error with s which is the only remaining relevant parameter. It can be seen how the larger the value of U the more sensitive is the numerical scheme to the choice of this parameter, specially the U ¼ 100 case requires a very small choice of s in order to get reliable results. As discussed above the pseudospectral method based on the Hermite expansion is more accurate than the FFT-based one, in fact the later scheme failed to obtain accurate results when U ¼ 100 for time steps down to Dt ¼ 5 105 .
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
231
Fig. 6. (a) and (b) Dependence of the error of the PFFT2 scheme as a function of Nx with s ¼ 0:01 for (a) U ¼ 2, (b) U ¼ 10. (c) Computation time for the previous cases as a function of Nx . The curve shows the expected superlinear (OðNx log Nx Þ) behavior.
Fig. 7. Dependence of the error of the PH2 scheme as a function of the time step s for m ¼ 100; X ¼ 40; X0 ¼ 1; Nx ¼ 4096. (a) U ¼ 2, (b) U ¼ 10, (c) U ¼ 100.
7. A Galerkin spectral scheme 7.1. The scheme From the previous subsections it is more or less clear that the basis hermite functions is specially suitable to approximate the solutions of this problem. One might wonder if this could be used as the basis for a nodeless spectral discretization of Galerkin type. The idea is to semidiscretize the problem by
232
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Fig. 8. Dependence of the error of the FFT2 scheme as a function of the time step s for X ¼ 40; X0 ¼ 1; Nx ¼ 4096. (a) U ¼ 2, (b) U ¼ 10.
P using the truncated expansion wðx; tÞ ’ M k¼1 ck ðtÞ/k ðxÞ. Inserting this expression into Eq. (2a) one finds M M X X dck ðtÞ 1 /k ðxÞ ¼ kþ ck ðtÞ/k ðxÞ i dt 2 k¼1 k¼1 þU
M X
ck ðtÞ/k ðxÞcl ðtÞ/l ðxÞcm ðtÞ/m ðxÞ:
ð35Þ
k;l;m¼1
Computing the scalar product by /j ðxÞ we get
M X dcj ðtÞ 1 ¼ jþ Ajklm ck ðtÞcl ðtÞcm ðtÞ; cj ðtÞ þ U i dt 2 k;l;m¼1
ð36Þ
R where Ajklm ¼ /j ðxÞ/k ðxÞ/l ðxÞ/m ðxÞ dx. Due to the fact that the basis functions are real this is a completely symmetric fourth order tensor. Eq. (36) is also a Hamiltonian system which may be obtained from M U X 1X 1 2 H ða; bÞ ¼ nþ Aklmn ðam an ak al an ðtÞ þ b2n ðtÞ þ 2 n¼1 2 4 k;l;m;n þ bk bl bm bn þ 2al am bk bn Þ;
ð37Þ
where c ¼ a þ ib. 7.2. Computational considerations The main a priori advantage of scheme (36) is that it is a fully Galerkin method and thus no spatial nodes are needed to obtain the solution (thus there is no need for any grid). In this sense, the boundary conditions are also automatically satisfied by the numerical solutions and there are no undesired effects arising from artificial boundary conditions.
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
233
However, there are several difficulties in applying directly this scheme. The first one is the computation of the tensor Ak;l;m;n . Although this is in principle an straightforward problem––it involves only evaluation of polynomial integrals multiplied by Gaussian functions––we must keep in mind that its computation involves integrals of products of four polynomials of high order. As discussed before this is a difficult problem to solve numerically (even harder than before because now we have a product of four polynomials that must be evaluated) and thus we have opted by computing the integrals using symbolic computation. Specifically we have used MAPLE. The computational cost of this operation increases exponentially with M. In fact, even in our test machine the computation of the tensor for M ¼ 20 took several hours of CPU time and the case M ¼ 30 took a couple of days to complete. As we have seen before, there are situations where M ¼ 50 or more modes must be kept. In those cases even the memory requirements for such a tensor are of the order of 50 Mb and each evaluation of the right hand side of Eq. (36) costs of the order of 5M 4 operations (i.e. for M ¼ 50, the number of operations to be done to evaluate once the RHS is about 3 107 ). Thus we see that the direct use of this scheme is computationally very expensive. But, even if one disregards the previous considerations with the hope that at least it could be possible to get a better than before approximation of the solutions of the problem, this is not true. This is so because the time discretization of Eq. (36) poses important problems. First of all explicit schemes do not preserve the hamiltonian structure of the problem and lead to computational instabilities. We have tried to solve these equations using both simple and advanced schemes such as the explicit midpoint rule and the Dormand– Prince variable step integrators with the result that there are divergences for times below our typical integration time T ¼ 100. Other possibilities such as implicit symplectic integrators could be used but at the cost of many evaluations of the right hand side of Eq. (36) per time step, which will definitely lead to nonpractical methods.
8. Conclusions In this work we have analyzed several schemes for the numerical solution of the HTNLSE. By using as test data exact solutions of the HTNLSE we have shown that several strategies which are appropriate for the integration of the usual NLSEs are not useful for these type of equations. We have proposed a pseudospectral scheme which is computationally very efficient and provides very good approximation of the solutions with long-time stability of the results. The philosophy here used may be applied to the construction of numerical schemes of higher order in time by just using the appropriate compositions of operators [38].
234
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
Acknowledgements VMPG is partially supported by the Ministerio de Ciencia y Tecnologıa under grant BFM2000-0521. XL acknowledges support from the ‘‘Programa de profesores invitados’’ of Universidad de Castilla-La Mancha. We would like to acknowledge Luis V azquez (U. Complutense) for a critical reading of the manuscript.
References [1] A. Hasegawa, Optical Solitons in Fibers, Springer-Verlag, Berlin, 1989. [2] R.K. Dodd, J.C. Eilbeck, J.D. Gibbon, H.C. Morris, Solitons and Nonlinear Wave Equations, Academic Press, New York, 1982. [3] A.S. Davydov, Solitons in Molecular Systems, Reidel, Dordrecht, 1985. [4] C. Sulem, P. Sulem, The Nonlinear Schr€ odinger Equation, Springer, Berlin, 2000. [5] L. V azquez, L. Streit, V.M. Perez-Garcıa (Eds.), Nonlinear Schr€ odinger and Klein-Gordon systems: Theory and Applications, World Scientific, Singapore, 1996. [6] F. Dalfovo, S. Giorgini, L.P. Pitaevskii, S. Stringari, Theory of Bose–Einstein condensation in trapped gases, Rev. Mod. Phys. 71 (1999) 463–512. [7] V.M. Perez-Garcıa, H. Michinel, H. Herrero, Bose–Einstein solitons in highly asymmetric traps, Phys. Rev. A 57 (1998) 3837–3845. [8] T. Busch, J.R. Anglin, Dark-bright solitons in inhomogeneous Bose–Einstein condensates, Phys. Rev. Lett. 87 (2001) 010401:1–4. [9] T. Busch, J.R. Anglin, Motion of dark solitons in trapped Bose–Einstein condensates, Phys. Rev. Lett. 84 (2000) 2298–2301. [10] S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, A. Sanpera, G.V. Shlyapnikov, M. Lewenstein, Dark solitons in Bose–Einstein condensates, Phys. Rev. Lett. 83 (1999) 5198–5201. [11] M. Kunze, T. K€ upper, V.K. Mezentsev, E.G. Shapiro, S. Turitsyn, Nonlinear solitary waves with gaussian tails, Physica D 128 (1999) 273–295. [12] M. Delfour, M. Fortin, G. Payre, Finite difference solutions of a nonlinear Schr€ odinger equation, J. Comp. Phys. 44 (1981) 277–288. [13] J.M. Sanz-Serna, J.G. Verwer, Conservative and nonconservative schemes for the solution of the nonlinear Schr€ odinger equation, IMA J. Numer. Anal. 6 (1986) 25–42. [14] L.S. Peranich, A finite difference scheme for solving a nonlinear Schr€ odinger equation with a linear damping term, J. Comp. Phys. 68 (1987) 501–505. [15] F. Zhang, V.M. Perez-Garcıa, L. Vazquez, Numerical simulation of nonlinear Schr€ odinger systems: a new conservative scheme, Appl. Math. Comput. 71 (1995) 164–177. [16] Y. Tang, L. Vazquez, F. Zhang, V.M. Perez-Garcıa, Symplectic methods for the nonlinear Schr€ odinger equation, Comp. Math. Appl. 32 (1996) 73–83. [17] Y. Tang, V.M. Perez-Garcıa, L. Vazquez, Symplectic methods for the Ablowitz–Ladik model, Appl. Math. Comput. 82 (1997) 17–38. [18] O. Karakashian, C. Makridakis, A space-time finite element method for the nonlinear Schr€ odinger equation: the discontinuous Galerkin method, Math. Comput. 67 (1998) 479–499. [19] Y. Tourigny, J. Morris, An investigation into the effect of product approximation in the numerical solution of the cubic nonlinear Schr€ odinger equation, J. Comp. Phys. 76 (1988) 103– 130. [20] B.M. Herbst, F. Varadi, M.J. Ablowitz, Symplectic methods for the nonlinear Schr€ odinger equation, Math. Comput. Simul. 37 (1994) 353–369.
V.M. Perez-Garcıa, X.-y. Liu / Appl. Math. Comput. 144 (2003) 215–235
235
[21] R. MacLahan, Symplectic integration of Hamiltonian wave equations, Numer. Math. 66 (1994) 465–492. [22] D. Pathria, J.Ll. Morris, Pseudo-spectral solution of nonlinear Schr€ odinger equations, J. Comp. Phys. 87 (1990) 108–125. [23] J.A.C. Weideman, B.M. Herbst, Split-step methods for the solution of the nonlinear Schr€ odinger equation, SIAM J. Numer. Anal. 23 (1986) 485–507. [24] A.D. Bandrauk, H. Shen, High-order split-step exponential methods for solving coupled nonlinear Schr€ odinger equations, J. Phys. A Math. Gen. 27 (1994) 7147–7155. [25] B.M. Herbst, J.Ll. Morris, A.R. Mitchell, Numerical experience with the nonlinear Schr€ odinger equation, J. Comput. Phys. 60 (1985) 282–305. [26] T.R. Taha, M.J. Ablowitz, Analytical and numerical aspects of certain nonlinear evolution equations II. Numerical nonlinear Schr€ odinger equation, J. Comp. Phys. 55 (1984) 203–230. [27] G. Agrawall, Nonlinear Fiber Optics, Academic Press, New York, 1995. [28] S. Hairer, Nordsett, Wanner, in: Solving Ordinary Differential Equations, vol. 1, Springer, Berlin, 1996. [29] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions, Dover, New York, 1964. [30] L. Erbe, M. Tang, Uniqueness of positive radial solutions of Du þ f ðjxj; uÞ ¼ 0, Diff. Int. Eq. 11 (1995) 725–743. [31] Z.-Q. Wang, Existence and symmetry of multi-bump solutions for nonlinear Schr€ odinger equations, J. Diff. Eq. 159 (1999) 102–137. [32] J.J. Garcıa-Ripoll, V.M. Perez-Garcıa, V. Vekslerchik, Construction of exact solutions by spatial translations in inhomogeneous nonlinear Schr€ odinger equations, Phys. Rev. E 64 (2001) 056602-1:6. [33] J.J. Garcıa-Ripoll, V.M. Perez-Garcıa, Minimization of Schr€ odinger functionals using sobolev gradients: Applications to quantum mechanics and nonlinear optics, SIAM J. Sci. Comp. 23 (2001) 1315–1324. [34] P.Z. Ding, C.X. Wu, Y.K. Mu, Y.X. Li, M.X. Jin, Square-preserving and symplectic structure and scheme for quantum system, Chin. Phys. Lett. 13 (1996) 245–248. [35] Z.Y. Zhou, P.Z. Ding, S.P. Pan, Study of a symplectic scheme for the time evolution of an atom in an external field, J. Kor. Phys. Soc. 32 (1998) 417–424. [36] X. Liu, X. Liu, Z. Zhou, P. Ding, S. Pan, Numerical solution of one-dimensional timeindependent Schr€ odinger equation by using symplectic schemes, Int. J. Quant. Chem. 79 (2000) 343–349. [37] S.K. Gray, D.E. Manolopoulos, Symplectic integrators tailored to the time-dependent Schr€ odinger equation, J. Chem. Phys. 104 (1996) 7099–7112. [38] S. Blanes, P.C. Moan, Practical symplectic partitioned Runge-Kutta and Runge–Kutta– Nystrom methods, J. Comput. Appl. Math. 142 (2002) 313–330.