Article Contents
Article Contents

# Numerical comparison of mass-conservative schemes for the Gross-Pitaevskii equation

The authors acknowledge the support by the Swedish Research Council (grant 2016-03339). P. Henning also acknowledges the support of the Göran Gustafsson Foundation

• In this paper we present a numerical comparison of various mass-conservative discretizations for the time-dependent Gross-Pitaevskii equation. We have three main objectives. First, we want to clarify how purely mass-conservative methods perform compared to methods that are additionally energy-conservative or symplectic. Second, we shall compare the accuracy of energy-conservative and symplectic methods among each other. Third, we will investigate if a linearized energy-conserving method suffers from a loss of accuracy compared to an approach which requires to solve a full nonlinear problem in each time-step. In order to obtain a representative comparison, our numerical experiments cover different physically relevant test cases, such as traveling solitons, stationary multi-solitons, Bose-Einstein condensates in an optical lattice and vortex pattern in a rapidly rotating superfluid. We shall also consider a computationally severe test case involving a pseudo Mott insulator. Our space discretization is based on finite elements throughout the paper. We will also give special attention to long time behavior and possible coupling conditions between time-step sizes and mesh sizes. The main observation of this paper is that mass conservation alone will not lead to a competitive method in complex settings. Furthermore, energy-conserving and symplectic methods are both reliable and accurate, yet, the energy-conservative schemes achieve a visibly higher accuracy in our test cases. Finally, the scheme that performs best throughout our experiments is an energy-conserving relaxation scheme with linear time-stepping proposed by C. Besse (SINUM, 42(3):934–952, 2004).

Mathematics Subject Classification: Primary: 35Q55, 65Y20, 65P10, 81Q05; Secondary: 65M60, 65Z99.

 Citation:

• Figure 1.  $L^1$-errors of density and $H^1$-errors for all five methods for the single soliton test case

Figure 2.  We observe remarkably small $L^2$-errors for the LCN-FEM when the spatial and temporal discretizations are coupled as $h/\tau = 2$. This feature virtually disappears in the error plot of the density, implying that for this ratio the error in phase is minimal

Figure 3.  Figures (a)-(c) show the solution to the two stationary soliton test case

Figure 4.  Long time behaviour for Model Problem 3.2

Figure 5.  Energy evolution for the Two-Step FEM and the LCN-FEM for the time discretizations $\tau = 2^{-8},2^{-9},\dots,2^{-20}$, note the different time scales. The Two-Step FEM is highly unstable for all ratios of $h$ and $\tau$. The LCN-FEM can be highly unstable depending on the ratio of $h$ and $\tau$

Figure 6.  Solution plots of density $|u_{h,\tau}|^2$ in Model Problem 4.1 for $\tau = 2^{-8}$ and $h = 0.06$

Figure 8.  Solution plots of density $|u_{h,\tau}|^2$ in Model Problem 4.1 for $\tau = 2^{-9}$ and $h = 0.06$

Figure 9.  For Model Problem 4.1, the $L^1$-errors in density and the $H^1$-errors in $u$ for the four methods versus time-step size for the spatial discretization $h = 0.06$

Figure 7.  Energy evolution of the IM-FEM and RE-FEM for two time-step sizes and the spatial discretization $h = 0.06$. The corresponding solutions to Model Problem 4.1 are depicted in Fig.6a-8d

Figure 11.  Numerical approximations at $T = 10$ and mesh resolution $h = 12\cdot2^{-6}$. The left picture depicts the reference solution. The middle figure shows the solution obtained with the Crank-Nicolson FE-method for $\tau = T/256$. We note that the approximations obtained with IM-FEM and RE-FEM for the same mesh size can visually not be distinguished from the CN-FEM result. However, the approximation obtained with LCN-FEM for $\tau = T/256$ as depicted in the right figure clearly suffers from an energy pollution and looks distorted compared to the other solutions

Figure 10.  Both figures depict a comparison between the isolines of the reference solution (black) and the isolines of a RE-FEM approximation (left, red lines) and a IM-FEM approximation (right, red), respectively. The snapshots are at maximum time $T = 10$ and the approximations are computed for $\tau = T/256$ and $h = 12\cdot2^{-6}$

Figure 12.  Density plots, $|u_{h,\tau}|^2$, for $\tau = 0.515\cdot2^{-11}$ and $h = 0.01$

Table 1.  Note that an optimal convergence rate for the error $e_h = u_{h,\tau}^N - u^N$ in $L^2$ is of order $\mathcal{O}(\tau^2+h^2)$ and an optimal convergence rate in $H^1$ is of order $\mathcal{O}(\tau^2+h)$. As everywhere in the paper, this refers to the case of piecewise linear FEM (simplicial Langrange elements) for the space discretization

 IM-FEM CN-FEM RE-FEM LCN-FEM Two-Step FEM Mass-Conservative Yes Yes Yes Yes Yes Energy-Conservative No Yes Yes No No Symplectic Yes No No No No Linear time-steps No No Yes Yes Yes $L^2$-convergence rates Optimal Optimal - Optimal Optimal Conditional Unconditional Unconditional Conditional $H^1$-convergence rates Optimal - - Optimal Optimal Conditional Unconditional Conditional

Table 4.  Average CPU-times over five runs, the number of time-steps is 1024 and $N$ denotes the degrees of freedom. The computational complexity is the same for the linear methods and, in 1D, roughly a tenth of that of the linear methods

 CPU [s] times 1D $N$ IM-FEM CN-FEM RE-FEM LCN-FEM TwoStep FEM 1600 47.92 46.58 5.78 5.40 5.50 3200 96.59 92.57 11.52 10.75 10.91 6400 197.89 192.53 23.38 21.83 22.18

Table 2.  $L^1$-errors of the density for the two stationary soliton test case at time $T = 2$

 $||\rho-\rho_{h,\tau}||_{L^1}$ $\tau$ IM-FEM CN-FEM RE-FEM LCN-FEM TwoStepFEM $\tau$ SP2 $2^{-3 }$ 2.684e14 3.368e12 12.068 21.033 11.757 $2^{-7}$ 11.582 $2^{-4 }$ 23.911 22.546 13.774 18.639 11.0617 $2^{-8}$ 6.2246 $2^{-5 }$ 23.325 19.659 7.544 18.145 19.100 $2^{-9}$ 2.1756 $2^{-6 }$ 15.885 12.738 4.648 18.525 15.223 $2^{-10}$ 0.5953 $2^{-7 }$ 8.498 5.355 1.265 4.868 7.733 $2^{-11}$ 0.1523 $2^{-8 }$ 3.356 1.442 0.319 8.109 2.301 $2^{-12}$ 0.0383 $2^{-9 }$ 0.959 0.360 0.080 3.537 0.577 $2^{-13}$ 0.0096 $2^{-10 }$ 0.252 0.091 0.021 0.723 0.145 $2^{-14}$ 0.0024

Table 3.  $H^1$-errors for the two stationary soliton test case at time $T = 2$

 $||u-u_{h,\tau}||_{H^1}$ $\tau$ IM-FEM CN-FEM RE-FEM LCN-FEM TwoStepFEM $\tau$ SP2 $2^{-3 }$ 3.795e10 4.270e9 11.445 24.614 10.629 $2^{-7}$ 9.8755 $2^{-4 }$ 10.368 11.108 15.194 22.414 11.157 $2^{-8}$ 4.6047 $2^{-5 }$ 9.544 11.474 13.620 29.591 10.877 $2^{-9}$ 1.7293 $2^{-6 }$ 10.400 9.935 3.768 36.725 10.542 $2^{-10}$ 0.4877 $2^{-7 }$ 6.175 4.662 1.022 9.990 6.601 $2^{-11}$ 0.1258 $2^{-8 }$ 2.704 1.261 0.260 12.261 2.043 $2^{-12}$ 0.0317 $2^{-9 }$ 0.806 0.315 0.066 3.277 0.514 $2^{-13}$ 0.0079 $2^{-10 }$ 0.214 0.080 0.018 0.632 0.129 $2^{-14}$ 0.0020

Table 5.  $H^1$-errors and CPU-times for different time-step sizes and mesh size $h = 24/51200$

 $||u-u_{h,\tau}||_{H^1}$, CPU [s] $\tau$ IM-FEM CN-FEM RE-FEM LCN-FEM Two-Step FEM $2^{-9}$ 2.401 1220 3.045 1130 3.045 130 95.175 130 9.338 140 $2^{-10}$ 2.474 1960 0.980 1800 0.979 210 38.324 210 9.148 240 $2^{-11}$ 0.710 3920 0.249 3540 0.249 370 1.948 370 8.942 450 $2^{-12}$ 0.179 7800 0.062 7270 0.062 700 0.295 680 8.899 870 $2^{-13}$ 0.044 10750 0.015 9800 0.015 1240 0.080 1190 8.683 1610

Table 6.  $L^2$ and $H^1$-errors and $L^1$-errors of the density for the rotating condensate experiment (i.e. Model Problem 4.2) at time $T = 10$ and for the mesh size $h = 12\cdot2^{-6}$

 $||u-u_{h,\tau}||_{L^2}$, $||\rho-\rho_{h,\tau}||_{L^1}$, $||u-u_{h,\tau}||_{H^1}$ $T/\tau$ IM-FEM CN-FEM RE-FEM LCN-FEM 64 0.6590 0.3560 1.9949 1.7527 0.3594 3.9750 1.7272 0.3553 3.9466 1.4095 1.4014 6.6062 128 1.1946 0.3496 2.1666 1.0690 0.3579 2.8822 1.0779 0.3579 2.9085 1.4500 1.2261 8.1410 256 1.1855 0.1893 3.1790 0.8746 0.1714 2.4329 0.8717 0.1707 2.4257 0.5402 0.1431 1.6532 512 0.3416 0.0515 0.9392 0.2362 0.0453 0.6699 0.2373 0.0456 0.6734 0.6026 0.4097 7.1779 1024 0.0847 0.0126 0.2333 0.0575 0.0110 0.1636 0.0599 0.0116 0.1706 0.0400 0.0118 0.2426 2048 0.0192 0.0027 0.0523 0.0110 0.0023 0.0348 0.0150 0.0029 0.0427 0.0111 0.0013 0.0239
•  [1] A. Aftalion, Vortices in Bose-Einstein Condensates, Progress in Nonlinear Differential Equations and their Applications, 67. Birkhäuser Boston, Inc., Boston, MA, 2006. [2] P. L. Christiansen, M. P. Sorensen and A. C. Scott, Nonlinear Science at the Dawn of the 21st Century, Lecture Notes in Physics, 542. Springer-Verlag, Berlin, 2000. doi: 10.1007/3-540-46629-0. [3] G. D. Akrivis, V. A. Dougalis and O. A. Karakashian, On fully discrete Galerkin methods of second-order temporal accuracy for the nonlinear Schrödinger equation, Numer. Math., 59 (1991), 31-53.  doi: 10.1007/BF01385769. [4] T. Aktosun, T. Busse, F. Demontis and C. van der Mee, Exact solutions to the nonlinear Schrödinger equation, in Topics in Operator Theory, Systems and Mathematical Physics, Oper. Theory Adv. Appl., Birkhäuser Verlag, Basel, 2 (2010), 1–12. doi: 10.1007/978-3-0346-0161-0_1. [5] R. Altmann, P. Henning and D. Peterseim, Quantitative Anderson localization of Schrödinger eigenstates under disorder potentials, submitted, 2018, arXiv: 1803.09950. [6] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev., 109 (1958), 1492–1505, URL https://link.aps.org/doi/10.1103/PhysRev.109.1492. doi: 10.1142/9789812567154_0007. [7] X. Antoine, W. Z. Bao and C. Besse, Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations, Comput. Phys. Commun., 184 (2013), 2621-2633.  doi: 10.1016/j.cpc.2013.07.012. [8] D. N. Arnold, G. David, D. Jerison, S. Mayboroda and M. Filoche, Effective confining potential of quantum states in disordered media, Phys. Rev. Lett., 116 (2016), 056602.  doi: 10.1103/PhysRevLett.116.056602. [9] W. Z. Bao and Y. Y. Cai, Mathematical theory and numerical methods for Bose-Einstein condensation, Kinet. Relat. Models, 6 (2013), 1-135.  doi: 10.3934/krm.2013.6.1. [10] W. Z. Bao and Q. Du, Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, SIAM J. Sci. Comput., 25 (2004), 1674-1697.  doi: 10.1137/S1064827503422956. [11] W. Z. Bao, S. Jin and P. A. Markowich, Numerical study of time-splitting spectral discretizations of nonlinear Schrödinger equations in the semiclassical regimes, SIAM J. Sci. Comput., 25 (2003), 27-64.  doi: 10.1137/S1064827501393253. [12] C. Besse, A relaxation scheme for the nonlinear Schrödinger equation, SIAM J. Numer. Anal., 42 (2004), 934-952.  doi: 10.1137/S0036142901396521. [13] C. Besse, S. Descombes, G. Dujardin and I. Lacroix-Violet, Energy Preserving Methods for Nonlinear Schrödinger Equations, 2018, arXiv: 1812.04890. [14] T. Cazenave, Semilinear Schrödinger Equations, Courant Lecture Notes in Mathematics, 10. New York University, Courant Institute of Mathematical Sciences, New York, American Mathematical Society, Providence, RI, 2003. doi: 10.1090/cln/010. [15] Q. S. Chang, E. Jia and W. Sun, Difference schemes for solving the generalized nonlinear Schrödinger equation, J. Comput. Phys., 148 (1999), 397-415.  doi: 10.1006/jcph.1998.6120. [16] D. Clément, A. F. Varón, J. A. Retter, L. Sanchez-Palencia, A Aspect and P. Bouyer, Experimental study of the transport of coherent interacting matter-waves in a 1d random potential induced by laser speckle, New Journal of Physics, 8 (2006), 165, URL http://stacks.iop.org/1367-2630/8/i=8/a=165. [17] L. Erdös, B. Schlein and H.-T. Yau, Derivation of the Gross-Pitaevskii equation for the dynamics of Bose-Einstein condensate, Ann. of Math., 172 (2010), 291-370.  doi: 10.4007/annals.2010.172.291. [18] D. L. Feder, A. A. Svidzinsky, A. L. Fetter and C. W. Clark, Anomalous modes drive vortex dynamics in confined Bose-Einstein condensates, Phys. Rev. Lett., 86 (2001), 564-567.  doi: 10.1103/PhysRevLett.86.564. [19] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature, 415 (2002), 39–44, URL http://dx.doi.org/10.1038/415039a. [20] D. F. Griffiths, A. R. Mitchell and J. L. Morris, A numerical study of the nonlinear Schrödinger equation, Comput. Methods Appl. Mech. Engrg., 45 (1984), 177-215.  doi: 10.1016/0045-7825(84)90156-7. [21] H. Hasimoto and H. Ono, Nonlinear modulation of gravity waves, Journal of the Physical Society of Japan, 33 (1972), 805-811.  doi: 10.1143/JPSJ.33.805. [22] P. Henning and A. Målqvist, The finite element method for the time-dependent Gross-Pitaevskii equation with angular momentum rotation, SIAM J. Numer. Anal., 55 (2017), 923-952.  doi: 10.1137/15M1009172. [23] P. Henning and D. Peterseim, Crank-Nicolson Galerkin approximations to nonlinear Schrödinger equations with rough potentials, Math. Models Methods Appl. Sci., 27 (2017), 2147-2184.  doi: 10.1142/S0218202517500415. [24] E. Jarlebring, S. Kvaal and W. Michiels, An inverse iteration method for eigenvalue problems with eigenvector nonlinearities, SIAM J. Sci. Comput., 36 (2014), A1978-A2001.  doi: 10.1137/130910014. [25] O. Karakashian and C. Makridakis, A space-time finite element method for the nonlinear Schrödinger equation: The continuous Galerkin method, SIAM J. Numer. Anal., 36 (1999), 1779-1807.  doi: 10.1137/S0036142997330111. [26] C. A. Mülle and D. Delande, Disorder and interference: Localization phenomena, arXiv e-prints. [27] J. M. Sanz-Serna, Methods for the numerical solution of the nonlinear Schrödinger equation, Math. Comp., 43 (1984), 21–27, URL http://dx.doi.org/10.2307/2007397. doi: 10.1090/S0025-5718-1984-0744922-X. [28] J. M. Sanz-Serna, Runge-Kutta schemes for Hamiltonian systems, BIT, 28 (1988), 877-883.  doi: 10.1007/BF01954907. [29] J. M. Sanz-Serna and J. G. Verwer, Conservative and nonconservative schemes for the solution of the nonlinear Schrödinger equation, IMA J. Numer. Anal., 6 (1986), 25-42.  doi: 10.1093/imanum/6.1.25. [30] Y. Tourigny, Optimal H1 estimates for two time-discrete Galerkin approximations of a nonlinear Schrödinger equation, IMA J. Numer. Anal., 11 (1991), 509-523.  doi: 10.1093/imanum/11.4.509. [31] J. G. Verwer and J. M. Sanz-Serna, Convergence of method of lines approximations to partial differential equations, Computing, 33 (1984), 297-313.  doi: 10.1007/BF02242274. [32] J. Wang, A new error analysis of Crank-Nicolson Galerkin FEMs for a generalized nonlinear Schrödinger equation, J. Sci. Comput., 60 (2014), 390-407.  doi: 10.1007/s10915-013-9799-4. [33] H. C. Yuen and B. M. Lake, Instabilities of waves on deep water, Annual Review of Fluid Mechanics, 12 (1980), 303-334.  doi: 10.1146/annurev.fl.12.010180.001511. [34] V. E. Zakharov, Stability of periodic waves of finite amplitude on a surface of deep fluid, Journal of Applied Mechanics and Technical Physics, 9 (1968), 190-194.  doi: 10.1007/BF00913182. [35] V. E. Zakharov and A. B. Shabat, Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media, Ž. Éksper. Teoret. Fiz., 61 (1971), 118-134. [36] G. Zhong and J. E. Marsden, Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators, Physics Letters A, 133 (1988), 134-139.  doi: 10.1016/0375-9601(88)90773-6. [37] G. E. Zouraris, On the convergence of a linear two-step finite element method for the nonlinear Schrödinger equation, M2AN Math. Model. Numer. Anal., 35 (2001), 389-405.  doi: 10.1051/m2an:2001121. [38] G. E. Zouraris, Error Estimation of the Besse Relaxation Scheme for a Semilinear Heat Equation, 2018, arXiv: 1812.09273.

Figures(12)

Tables(6)