# American Institute of Mathematical Sciences

December  2019, 12(6): 1247-1271. doi: 10.3934/krm.2019048

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

 Department of Mathematics, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden

Received  October 2018 Revised  March 2019 Published  September 2019

Fund Project: 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).

Citation: Patrick Henning, Johan Wärnegård. Numerical comparison of mass-conservative schemes for the Gross-Pitaevskii equation. Kinetic & Related Models, 2019, 12 (6) : 1247-1271. doi: 10.3934/krm.2019048
##### References:

show all references

##### References:
$L^1$-errors of density and $H^1$-errors for all five methods for the single soliton test case
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
Figures (a)-(c) show the solution to the two stationary soliton test case
Long time behaviour for Model Problem 3.2
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$
Solution plots of density $|u_{h,\tau}|^2$ in Model Problem 4.1 for $\tau = 2^{-8}$ and $h = 0.06$
Solution plots of density $|u_{h,\tau}|^2$ in Model Problem 4.1 for $\tau = 2^{-9}$ and $h = 0.06$
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$
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
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
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}$
Density plots, $|u_{h,\tau}|^2$, for $\tau = 0.515\cdot2^{-11}$ and $h = 0.01$
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
 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
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
 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
$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
 $||\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
$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
 $||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
$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
 $||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
$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
 $||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] Roy H. Goodman, Jeremy L. Marzuola, Michael I. Weinstein. Self-trapping and Josephson tunneling solutions to the nonlinear Schrödinger / Gross-Pitaevskii equation. Discrete & Continuous Dynamical Systems - A, 2015, 35 (1) : 225-246. doi: 10.3934/dcds.2015.35.225 [2] Jeremy L. Marzuola, Michael I. Weinstein. Long time dynamics near the symmetry breaking bifurcation for nonlinear Schrödinger/Gross-Pitaevskii equations. Discrete & Continuous Dynamical Systems - A, 2010, 28 (4) : 1505-1554. doi: 10.3934/dcds.2010.28.1505 [3] Xiaoyu Zeng, Yimin Zhang. Asymptotic behaviors of ground states for a modified Gross-Pitaevskii equation. Discrete & Continuous Dynamical Systems - A, 2019, 39 (9) : 5263-5273. doi: 10.3934/dcds.2019214 [4] Georgy L. Alfimov, Pavel P. Kizin, Dmitry A. Zezyulin. Gap solitons for the repulsive Gross-Pitaevskii equation with periodic potential: Coding and method for computation. Discrete & Continuous Dynamical Systems - B, 2017, 22 (4) : 1207-1229. doi: 10.3934/dcdsb.2017059 [5] Gianluca Frasca-Caccia, Peter E. Hydon. Locally conservative finite difference schemes for the modified KdV equation. Journal of Computational Dynamics, 2019, 6 (2) : 307-323. doi: 10.3934/jcd.2019015 [6] M. D. Todorov, C. I. Christov. Conservative numerical scheme in complex arithmetic for coupled nonlinear Schrödinger equations. Conference Publications, 2007, 2007 (Special) : 982-992. doi: 10.3934/proc.2007.2007.982 [7] Walid K. Abou Salem, Xiao Liu, Catherine Sulem. Numerical simulation of resonant tunneling of fast solitons for the nonlinear Schrödinger equation. Discrete & Continuous Dynamical Systems - A, 2011, 29 (4) : 1637-1649. doi: 10.3934/dcds.2011.29.1637 [8] Bernard Ducomet, Alexander Zlotnik, Ilya Zlotnik. On a family of finite-difference schemes with approximate transparent boundary conditions for a generalized 1D Schrödinger equation. Kinetic & Related Models, 2009, 2 (1) : 151-179. doi: 10.3934/krm.2009.2.151 [9] Thierry Cazenave, Yvan Martel, Lifeng Zhao. Finite-time blowup for a Schrödinger equation with nonlinear source term. Discrete & Continuous Dynamical Systems - A, 2019, 39 (2) : 1171-1183. doi: 10.3934/dcds.2019050 [10] Norman E. Dancer. On the converse problem for the Gross-Pitaevskii equations with a large parameter. Discrete & Continuous Dynamical Systems - A, 2014, 34 (6) : 2481-2493. doi: 10.3934/dcds.2014.34.2481 [11] E. Norman Dancer. On a degree associated with the Gross-Pitaevskii system with a large parameter. Discrete & Continuous Dynamical Systems - S, 2019, 12 (7) : 1835-1839. doi: 10.3934/dcdss.2019120 [12] Philipp Bader, Sergio Blanes, Fernando Casas, Mechthild Thalhammer. Efficient time integration methods for Gross-Pitaevskii equations with rotation term. Journal of Computational Dynamics, 2019, 6 (2) : 147-169. doi: 10.3934/jcd.2019008 [13] Ko-Shin Chen, Peter Sternberg. Dynamics of Ginzburg-Landau and Gross-Pitaevskii vortices on manifolds. Discrete & Continuous Dynamical Systems - A, 2014, 34 (5) : 1905-1931. doi: 10.3934/dcds.2014.34.1905 [14] Thomas Chen, Nataša Pavlović. On the Cauchy problem for focusing and defocusing Gross-Pitaevskii hierarchies. Discrete & Continuous Dynamical Systems - A, 2010, 27 (2) : 715-739. doi: 10.3934/dcds.2010.27.715 [15] D.G. deFigueiredo, Yanheng Ding. Solutions of a nonlinear Schrödinger equation. Discrete & Continuous Dynamical Systems - A, 2002, 8 (3) : 563-584. doi: 10.3934/dcds.2002.8.563 [16] Pavel I. Naumkin, Isahi Sánchez-Suárez. On the critical nongauge invariant nonlinear Schrödinger equation. Discrete & Continuous Dynamical Systems - A, 2011, 30 (3) : 807-834. doi: 10.3934/dcds.2011.30.807 [17] Tarek Saanouni. Remarks on the damped nonlinear Schrödinger equation. Evolution Equations & Control Theory, 2020, 9 (3) : 721-732. doi: 10.3934/eect.2020030 [18] Younghun Hong. Scattering for a nonlinear Schrödinger equation with a potential. Communications on Pure & Applied Analysis, 2016, 15 (5) : 1571-1601. doi: 10.3934/cpaa.2016003 [19] Alexander Komech, Elena Kopylova, David Stuart. On asymptotic stability of solitons in a nonlinear Schrödinger equation. Communications on Pure & Applied Analysis, 2012, 11 (3) : 1063-1079. doi: 10.3934/cpaa.2012.11.1063 [20] Dario Bambusi, A. Carati, A. Ponno. The nonlinear Schrödinger equation as a resonant normal form. Discrete & Continuous Dynamical Systems - B, 2002, 2 (1) : 109-128. doi: 10.3934/dcdsb.2002.2.109

2019 Impact Factor: 1.311