# Time-splitting Chebyshev-Spectral method for the Schrodinger equation in the semiclassical regime with zero far-field boundary conditions.

Byline: Hongsheng Wanga, Yushan Nia and Junwan LibAbstract: Semiclassical limit of Schrodinger equation with zero far-field boundary conditions is investigated by the time- splitting Chebyshev-spectral method. The numerical results of the position density and current density are presented. The time-splitting Chebyshev-spectral method is based on Strang splitting method in time coupled with Chebyshev- spectral approximation in space. Compared with the time-splitting Fourier-spectral method, the time-splitting Chebyshev- spectral method is unnecessary to treat the wave function as periodic and holds the smoothness of the wave function.

For different initial conditions and potential (e.g. constant potential and harmonic potential), extensive numerical test examples in one-dimension are studied. The numerical results are in good agreement with the weak limit solutions. It shows that the time-splitting Chebyshev-spectral method is effective in capturing s-oscillatory solutions of the Schrodinger equation with zero far-field boundary conditions. In addition, the time-splitting Chebyshev-spectral method surpasses the traditional finite difference method in the meshing strategy due to the exponentially high-order accuracy of Chebyshev-spectral method.

Keywords: Schrodinger equation, time-splitting Chebyshev-spectral method, zero far-field boundary conditions, semiclassical limit.

1. INTRODUCTION

Many problems of solid state physics require the (numerical) solution of the following Schrodinger equation [1-3] in the case of a small Planck constant s (0less thansless thanless than1):

where V(x) is a given real-valued electrostatic potential, and us = us(x, t) is the complex-valued wave function which is used to compute

observables (including the primary physical quantities) in classical quantum physics. The introduction of the small Planck constant s in Schrodinger equation with oscillatory initial data offers great help for the equation resolution [4].

The work about the behavior of solutions of the nonlinear Schrodinger equation in the semiclassical limit (s*0) was entirely motivated by a natural mathematicalquestion [5], however, in the recent years many researches show that the semiclassical limit of Schrodinger equation is also of direct importance to basic physics and technology in nonlinear optics [5-7]. As mentioned in [5], the semiclassical limit of nonlinear Schrodinger equation provides an idealized description of optical shocks and wave breaking in the nonlinearpropagation of laser pulses in optical fibers.

In [6, 7], the semiclassical limit of nonlinear Schrodinger equation has also been applied to describe "nonreturn to zero" optical pulses in nonlinear fibers, a problem of central importance of the current technology of long distance (transoceanic) telephone communication. Investigations on the semiclassical limit of Schrodinger equation including theoretical analysis, numericalcalculation and application have become an important subject, particularly by the introduction of tools frommicrolocal analysis, such as defect measures [8], H- measures [9] and Wigner measures [10-12]. These studies have provided technical tools for exploiting properties of Schrodinger equation in the semiclassical regime.

It is well known that Eq. (1) propagates oscillations of wave length O(s) in space and time, which inhibit usfrom converging strongly. This implies that the analysis of the so-called semiclassical limit (s*0) is a mathematically rather complex issue [1, 2, 13], and up to now there has not been a rigorous analytical theory of the semicalssical limit of Schrodinger equation [13], thus it is natural these days to study this problem numerically. Unfortunately, the highly oscillatory nature of the solutions of Schrodinger equation with small s provides severe numerical burdens [2, 13, 14].

The semiclassical initial value problem, namely Eq. (1) and (2), is a classic example of a "stiff" problem, as it contains two vastly different spatial and temporal scales and involves essential competition between effects on these two scales [13].

Even for stable discretization schemes (or under mesh size restrictions which guarantee stability), the oscillations may very well disturb the solution in such a way that the quadratic macroscopic quantities and other physical observables come out completely wrong unless the spatial-temporal oscillations are fully resolved numerically, i.e. using many grid points per wavelength of O(s) [1, 2]. Markowich [1, 15] et al. applied Wigner- transform techniques to the analysis of finite difference methods for Schrodinger equation with small s andzero far-field boundary conditions ( us ( x, t ) decays to zero sufficiently fast as x c oe [1, 2, 5]), they obtainedsharp conditions on the spatial-temporal grid whichguaranteed convergence for values of all observables as Plank constant s tended to zero.

Their results show that, for the best combination of the time and space discretizations, one needs the following constraint: mesh size h = o(s) and time step k = o(s). Failure to satisfy this constraint leads to wrong numerical observables. It appears clearly that the finite difference methods for this problem have very high requirement in the meshing strategy when s inclines to zero. Since the solution of Schrodinger equation in the semiclassical regime is highly oscillatory but smooth, spectral methodis perfect to compute the spatial derivatives of us tohigh accuracy on a simple domain [16], and is helpful to solve the Schrodinger equation in the semiclassical regime. Feit [17] et al. used spectral method to determine the eigenvalues and eigenfunctions of Schrodinger equation.

Pathria [18] et al. studied the time-splitting spectral discretizations of Schrodinger equation for s = O(1) which did not give any clue about the semiclassical limit. Recently, Bao et al. systemically investigated the time-splitting Fourier -spectral method [2, 19] and time-splitting sine-spectral method [20] for the semiclassical limit of Schrodinger equation with zero far-field boundary conditions and applied themsuccessfully to simulate Bose-Einstein condensation21]. Bao [22] also proposed the time-splittingChebyshev-spectral method to preliminarily solve Schrodinger equation with nonzero far-field conditions. Moreover, Zhang [14, 23] et al. presented a time- splitting and space-time adaptive wavelet scheme to solve Schrodinger equation with zero far-field boundary conditions.

As a rule, Fourier-spectral method is used to compute the spatial derivatives of the periodic function [16]. For the mathematically nonperiodic function, if it is very close to zero (for example, exponentially decaying to zero) at the ends of the computing interval, one approach would be to regard the function as periodic and to use Fourier-spectral method approximately after choosing an appropriately large computing interval [16]. Based on this approach, although the solution of Schrodinger equation in the semiclassical regime is nonperiodic, Bao [2, 19] presented the time-splitting Fourier-spectral method to solve Schrodinger equation with zero far-field boundary conditions.

However, this kind of approximation on periodicity usually damagesthe smoothness of the wave function, and sacrifices the accuracy advantages of the spectral methods because the accuracy of spectral methods depends on the smoothness of the functions being solved. Sometimes Gibbs phenomenon even occurs [16]. Now that the solution of Schrodinger equation with zero far-field boundary conditions is nonperiodic, it is a better choice to use the time-splitting Chebyshev-spectral method which is more applicable to the nonperiodic function [16].

In this paper, we derive the algorithmic formula of the time-splitting Chebyshev-spectral method and use it to solve Schrodinger equation with zero far-field boundary conditions which is previously solved by the time-splitting Fourier-spectral method [2, 19]. Extensive numerical test examples in one-dimension are studied for different initial conditions and potential (e.g. constant potential and harmonic potential). Comparisons between the numerical results and the weak limits verify the availability of the time-splitting Chebyshev-spectral method for Schrodinger equation with zero far-field boundary conditions. We also present the meshing strategy of the time-splitting Chebyshev-spectral method, and it surpasses the traditional finite difference method due to the exponentially high-order accuracy of Chebyshev- spectral method.

2. TIME-SPLITTING CHEBYSHEV-SPECTRAL METHOD

In this section, we derive the algorithmic formula of the time-splitting Chebyshev -spectral method. We consider Eq. (1) and (2) with zero far-field boundary conditions in the case of one-dimension (d = 1):

It must be pointed out that the solution of Schrodinger equation (3) with zero far-field boundary conditions (5) is nonperiodic. Compared with the time- splitting Fourier-spectral method, we use time-splitting Chebyshev-spectral method to solve this nonperiodic problem.

By classical quantum physics [24] the wave function us is an auxiliary quantity used to compute the primary physical quantities, which are quadratic functions of us,e.g. the position densityqs ( x, t ) =| us ( x, t ) |2 (6) and the current density where "---" denotes complex conjugation2.1. Grid PartitioningFor simplicity, we choose the spatial grid points asChebyshev-Gauss-Lobatto interpolation points [25, 26], then let the grid points beWe choose the time step size k and let time stepsbe t := nk n = 0,1,2, ......2.2. Strang Splitting MethodAs preparatory steps, we begin by introducing the splitting method [27] for a general evolution equation where f(u) is a operator and the splitting f(u) = Au+ Bucan be quite arbitrary. For a given time step k, let tn =nk, n = 0,1,2, ... and un be the approximation of u(t ).

Also known as the symmetric operator splitting, a standard second-order Strang splitting method [28-32in time is as follows:then we can obtain un+1 from un by solving Eq. (10).Eq. (3) multiplying by 1/s, we getthen from time t = tn to time t = tn+1 Eq. (11) can be solved by the following three steps via Strang splitting method , namely Eq. (10).First, we solvefor half of one time step (of length k), followed by solvingwith the boundary condition (5) for one time step k.Finally, we solveNote that the ODE (13), (15) can be solved exactly:for half of one time step againtherefore we only need to solve Eq. (14) subject to the boundary condition (5).2.3. Chebyshev-Spectral MethodEq. (14) with the boundary condition (5) can be discretized in space by Chebyshev-spectral method [25, 26] and integrated in time exactly by applying a diagonalization technique for the ODE system in phase space.

Using Chebyshev interpolation [25, 26], letwhere Tm(x) is the m-th Chebyshev polynomialPlugging the expansion (17) into the boundary condition (5), we havePlugging the expansion (17) into Eq. (14), we can obtain the following ODE system by using the propertyof the expansion coefficients of Chebyshev interpolation [25]:

Removing all aM-1(t) and aM(t) in Eq. (22) byemploying Eq. (21), letthus the ODE system (22) can now be written aswhere T = (tjk) is an (M-1)x(M-1) constant matrix with entries2.4. Solving via Diagonalization TechniqueThe matrix T has M-1 distinct negative eigenvalues [25], therefore it is diagonalizable, i.e. there exists an invertible matrix P and a diagonal matrix D such thatEq. (24) multiplying by the matrix P-1, we obtainSince the matrix D is diagonal, the equations in ODE system (28) are actually independent of each other. After a simple calculation, we can getCombining Eq. (27) and Eq. (29), we haveChoosing t = tn in Eq. (30), we obtainSubstituting Eq. (32) into Eq. (30), we get

This implies that we can get .(tEq. (33).2.5. Algorithmic Formula n+1 ) from .(t ) viaFinally, we state the result for the algorithmic formula of the time-splitting Chebyshev -spectral method.From time t = tn to t = tn+1, we combine Eq. (16), (17), (21) and (33) via the standard second-order Strang splitting and obtain the time-splitting Chebyshev-spectral method. The steps for obtaining In the above progress, we used the second-order Strang splitting method. For different precision requirement, we could choose other high-order splitting methods [33].

3. EXAMPLES AND NUMERICAL RESULTS

In this section, for different initial conditions and potential, we study numerically the meshing strategy oftime-splitting Chebyshev-spectral method for Schrodinger equation with zero far-field boundary conditions and investigate the semiclassical limit (s*0) of Schrodinger equation in one-dimension.In our computations, the initial condition is always chosen in the classical WKB form:with A0 and S0 independent of s, real valued and withA0(x) decaying to zero sufficiently fast as x c oe .

Weremark that the initial condition (41) is oscillatory, and Schrodinger equation (3) propagates oscillations of wave length O(s) in space and timeWe choose an appropriately long interval [a, b] for the computations such that the zero far-field boundaryconditions (5) do not introduce a significant error relative to the whole space problemAs already pointed out, the main goal when solvingSchrodinger equation is to compute macroscopic quantities (e.g. the position density and the currentdensity) associated to the wave function, in the following computations we therefore give the numericalresults of the position density qs(x, t) and the current density Js(x, t) rather than the results of wave function

Two sets of numerical experiments are designed in the followingExample 1 (Markowich, Pietra and Pohl [1, 15]) The initial condition is taken asWe solve on the x-interval [0, 1], i.e. a = 0 and b =1. Let V(x) = 100 be a constant potential.In order to test the meshing strategy of the time- splitting Chebyshev-spectral method and to investigate the semiclassical limit (s*0), we compute the following three cases with different combinations of s and MThe weak limits q0(x, t), J0(x, t) of qs(x, t), Js(x, t), respectively, as s*0, have been given in [1, 2, 15]. As a reference purposes, we also plot them at t = 0.54 in Figures 1, 2 and 3.

In case (I), as shown in Figure 1, s1 = 0.0256 is too large compared to zero, and the error between numerical results and weak limits is large.In case (II), we perform tests similar to those in case(I). Figure 2 shows the corresponding results. Due to s2= 0.0064 close to zero, the error between numericalresults and weak limits is much smaller than the error in Figure 1.In case (III), s3 = 0.0008 is sufficiently small, thus the numerical solutions capture the correct weak limits in Figure 3.Although the macroscopic quantities qs(x, 0.54) and Js(x, 0.54) for the initial condition (42) have no oscillations in Figure 1, 2 and 3, respectively, the wave function us has oscillations in the phase [1].

A choice of the discretization parameters that does not take care ofthe oscillations can give the wrong results [1]. As far as this example is concerned, in order to guarantee good approximations for s small, the traditional finite difference method needs the following constraint [1]: h= o(s), k = o(s). However, in Figure 1, 2 and 3, we can observe numerical convergence (in the weak sense) tothe limit solution as s*0 under the meshing strategysplitting Chebyshev-spectral method has much better resolution for oscillatory solutions of Schrodinger equation than the finite difference method.Example 2 (Gasser, Markowich [10]) The initial condition is taken asThe weak limits q0(x, t), J0(x, t) of qs(x, t), Js(x, t), respectively, as s*0, have been given in [2, 10]. As areference purposes, we also plot them at t1 = 0.52 inFigure 4 and t2 = 3.6 in Figure 5, respectively. We solve on the x-interval [-2, 2], i.e. a = -2 and b =2. Let V(x) = x2/2, which is a harmonic oscillator.

Once more, in order to test the meshing strategy of the time-splitting Chebyshev-spectral method and to investigate the semiclassical limit (s*0), we compute the following two cases with different combinations of s and M To obtain a better visualization in Figures 4 and 5, we depict the solutions in a subinterval instead of in the whole computational interval [-2, 2]. In Figures 4 and 5, we can see numerical convergence to the weak limits for different t with s = 0.0025 which is very close to zero.

Remark 1. From the numerical results of these two examples, in which different initial conditions and potential are studied, we can see that the time-splitting Chebyshev-spectral method gives very promisingresults under the meshing strategy " b - a = O(s ) , k =MO(s)" for Schrodinger equation with zero far-field boundary conditions.

4. CONCLUSIONS

Semiclassical limit (s*0) of Schrodinger equation with zero far-field boundary conditions is investigated by the time-splitting Chebyshev-spectral method. The time-splitting Chebyshev-spectral method is based on Strang splitting method in time coupled with Chebyshev-spectral approximation in space. Compared with the time-splitting Fourier-spectral method, the time-splitting Chebyshev-spectral method is unnecessary to treat the wave function as periodic and holds the smoothness of the wave function. It must be pointed out that the solution of Schrodinger equation with zero far-field boundary conditions is nonperiodic and the accuracy of spectral methods depends on the smoothness of the functions being approximated.

Therefore, the time-splitting Chebyshev-spectral method can elude the risk of sacrificing the accuracy advantages of the spectral methods and the risk of Gibbs phenomenon. For different initial conditions and potential (e.g. constant potential and harmonic potential), our numerical results capture the correct weak limits as s*0 and show that the time-splitting Chebyshev-spectral method gives promising results for s small when the meshing strategy satisfies" b - a = O(s ) , k = O(s)", while the frequently used finiteMdifference methods require mesh size and time step much smaller than Plank constant s.

It appears clearly that the time-splitting Chebyshev-spectral method has much better resolution for oscillatory solutions of Schrodinger equation than the finite difference methods. In addition, the application of the diagonalization technique greatly simplifies the calculation process, in which the equations in ODE system are transformed to be independent of each other. Hence, the time-splitting Chebyshev-spectral method is very effective in capturing s-oscillatory solutions of Schrodinger equation with zero far-field boundary conditions in the semiclassical regime and it is a better choice compared with the time-splitting Fourier-spectral method and the finite difference methods in the case of zero far-field boundary conditions.

ACKNOWLEDGEMENTS

This work was supported by the National NaturalScience Foundation of China (Grant No. 10576010).

The authors would like to heartfully thank Professor Weizhu Bao of National University of Singapore for his help.

REFERENCES

[1] Markowich PA, Pietra P, Pohl C. Numerical approximation of quadratic observables of Schrodinger-type equations in the semi-classical limit. Numer Math 1999; 81: 595-30. http://dx.doi.org/10.1007/s002110050406

[2] Bao W, Jin S, Markowich PA. On time-splitting spectral approximations for the Schrodinger equation in the semiclassical regime. J Comput Phys 2002; 175: 487-24. http://dx.doi.org/10.1006/jcph.2001.6956

[3] Alonso MA, Forbes GW. New approach to semiclassical analysis in mechanics. J Math Phys 1999; 40: 1699-18. http://dx.doi.org/10.1063/1.532829

[4] Fannjiang A, Jin S, Papanicolaou G. High frequency behavior of the focusing nonlinear Schrodinger equation with random inhomogeneities. SIAM J Appl Math 2003; 63: 1328-58. http://dx.doi.org/10.1137/S003613999935559X

[5] Bronski JC, McLaughlin DW. Semiclassical behavior in the NLS equation: Optical shocks-focusing instabilities. In: Ercolani N M, Gabitov I R, Levermore C D, Serre D, editiors. Singular limits of dispersive waves; 1991: NATO Adv. Sci. Inst. Ser. B Phys., 320. Plenum 1994: pp. 21-38.

[6] Forest MG, McLaughlin KT-R. Onset of oscillations in nonsoliton pulses in nonlinear dispersive fibers. J Nonlinear Sci 1998; 8: 43-62.http://dx.doi.org/10.1007/s003329900043

[7] Kodama Y, Wabnitz S. Analytical theory of guiding-center nonreturn-to-zero and return-to-zero signal transmission in normally dispersive nonlinear optical fibers. Optics Lett 1995;20: 2291-93. http://dx.doi.org/10.1364/OL.20.002291

[8] Gerard P. Microlocal defect measures. Comm PDE 1991; 16:1761-94. http://dx.doi.org/10.1080/03605309108820822

[9] Tartar L. H-measures, a New Approach for Studying Homogenization, Oscillations and Concentration Effects in Partial Differential Equations. Proc Roy Soc Edinburgh Sect A 1990; 115: 193-30. http://dx.doi.org/10.1017/S0308210500020606

[10] Gasser I, Markowich PA. Quantum hydrodynamics, Wigner transform and the classical limit. Asymptotic Anal 1997; 14:97-116.

[11] Gerard P, Markowich PA, Mauser NJ, Poupaud F.Homogeneization limits and Wigner transforms. Comm Pure Appl Math 1997; 50: 323-79. http://dx.doi.org/10.1002/(SICI)1097-0312(199704)50:4less than323::AID-CPA4greater than3.0.CO;2-C

[12] Markowich PA, Mauser NJ, Poupaud F. A Wigner function approach to semiclassical limits: electrons in a periodic potential. J Math Phys 1994; 35: 1066-94. http://dx.doi.org/10.1063/1.530629

[13] Miller PD, Kamvissis S. On the semiclassical limit of the focusing nonlinear Schrodinger equation. Phys Lett A 1998;247: 75-86.http://dx.doi.org/10.1016/S0375-9601(98)00565-9

[14] Zhang R, Zhang K. Application of time-splitting and wavelet based space-time adaptive method to solving Schrodinger equations. J Jilin Univ 2004; 42: 176-78.

[15] Markowich PA, Pietra P, Pohl C, Stimming HP. A Wigner- measure analysis of the Dufort-Frankel scheme for the Schrodinger equation. SIAM J Numer Anal 2002; 40: 1281-10. http://dx.doi.org/10.1137/S0036142900381734

[16] Trefethen Lloyd N. Spectral Methods in MATLAB: Philadelphia 2000.http://dx.doi.org/10.1137/1.9780898719598

[17] Feit MD, Fleck JA, Steiger A. Solution of the Schrodinger equation by a spectral method. J Comput Phys 1982; 47:412-33.http://dx.doi.org/10.1016/0021-9991(82)90091-2

[18] Pathria D, Morris JL. Pseudospectral solution of nonlinear Schrodinger equations. J Comput Phys 1990; 87: 108-25. http://dx.doi.org/10.1016/0021-9991(90)90228-S

[19] Bao W, Jin S, Markowich PA. Numerical study of time- splitting spectral discretizations of nonlinear Schrodinger equations in the semi-classical regimes. SIAM J Sci Comput2003; 25: 27-64. http://dx.doi.org/10.1137/S1064827501393253[20] Bao W, Jaksch D. An explicit unconditionally stable numerical method for solving damped nonlinear Schrodinger equations with a focusing nonlinearity. SIAM J Numer Anal2003; 41: 1406-26. http://dx.doi.org/10.1137/S0036142902413391

[21] Bao W, Jaksch D, Markowich PA. Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation. J Comput Phys 2003; 187: 318-42. http://dx.doi.org/10.1016/S0021-9991(03)00102-5

[22] Bao W. Numerical methods for the nonlinear Schrodinger equation with nonzero far-field conditions. Methods Appl Anal2004; 11: 367-87.

[23] Zhang R, Zhang K, Zhou YS. Numerical study of time- splitting, space-time adaptive wavelet scheme for Schrodinger equations. J Comput Appl Math 2006; 195: 263-73. http://dx.doi.org/10.1016/j.cam.2005.03.086

[24] Landau LD, Lifschitz EM. Lehrbuch der Theoretischen Physik - Quantenmechanik: Akademic-Verlag 1985.

[25] Gottlieb D, Orszag SA. Numerical Analysis of Spectral methods: Philadelphia 1977.

[26] Canuto C, Hussaini MY, Quarteroni A, Zhang TA. SpectralMethods in Fluid Dynamics: Springer-Verlag 1988.

[27] Fornberg B, Driscoll TA. A fast spectral algorithm for nonlinear wave equations with linear dispersion. J Comput Phys 1999; 155: 456-67. http://dx.doi.org/10.1006/jcph.1999.6351

[28] Strang G. On the construction and composition of difference schemes. SIAM J Numer Anal 1968; 5: 506-17. http://dx.doi.org/10.1137/0705041

[29] Gradinaru V. Strang splitting for the time-dependentSchrodinger equation on sparse grids. SIAM J Numer Anal2007; 46: 103-23. http://dx.doi.org/10.1137/050629823

[30] Gradinaru V, Tubingen. Fourier transform on sparse grids: code design and application to the time dependent Schrodinger equation. Computing 2007; 80: 1-22. http://dx.doi.org/10.1007/s00607-007-0225-3

[31] Jahnke T, Lubich C. Error bounds for exponential operator splitting. BIT 2000; 40: 735-44. http://dx.doi.org/10.1023/A:1022396519656

[32] Sportisse B. An analysis of operator splitting techniques in the stiff case. J Comput Phys 2000; 161: 140-68. http://dx.doi.org/10.1006/jcph.2000.6495

[33] Yoshida H. Construction of higher order symplectic integrators. Phys Lett A 1990; 150: 262-68. http://dx.doi.org/10.1016/0375-9601(90)90092-3aDepartment of Mechanics and Engineering Science, Fudan University, Shanghai 200433, ChinabSchool of Materials Science and Engineering, Shanghai University, Shanghai 200072, ChinaAddress corresponding to this author at the Department of Mechanics andEngineering Science, Fudan University, Shanghai 200433, China; Tel: +86-21-65642745; Fax: +86-21-65642745; E-mail: niyushan@fudan.edu.cnPACS numbers: 02.60.Cb, 02.60.Lj, 02.70.Hm, 03.65.Sq.http://dx.doi.org/10.6000/1927-5129.2013.09.11(c) 2013 Wang et al.; Licensee Lifescience Global.

This is an open access article licensed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted, non-commercial use, distribution and reproduction in any medium, provided the work is properly cited.

Printer friendly Cite/link Email Feedback | |

Author: | Wanga, Hongsheng; Nia, Yushan; Lib, Junwan |
---|---|

Publication: | Journal of Basic & Applied Sciences |

Date: | Dec 31, 2013 |

Words: | 3945 |

Previous Article: | Two integral operators defined with Bessel functions on the classN(th ). |

Next Article: | Fundamental problems of internal gravity waves dynamics in ocean. |

Topics: |