Numerical stability of explicit Runge-Kutta finite-difference schemes for the nonlinear Schrödinger equation

R. M. Caplan¹¹1Corresponding author. Present address: Predictive Science Inc. 9990 Mesa Rim Rd, Suite 170, San Diego, CA 92121. email: caplanr@predsci.com, phone: 858-225-2314 and R. Carretero-González
Nonlinear Dynamical System Group²²2URL: http://nlds.sdsu.edu, Computational Science Research Center, and
Department of Mathematics and Statistics, San Diego State University,
San Diego, California 92182-7720, USA

(August 19, 2025)

Abstract

Linearized numerical stability bounds for solving the nonlinear time-dependent Schrödinger equation (NLSE) using explicit finite-differencing are shown. The bounds are computed for the fourth-order Runge-Kutta scheme in time and both second-order and fourth-order central differencing in space. Results are given for Dirichlet, modulus-squared Dirichlet, Laplacian-zero, and periodic boundary conditions for one, two, and three dimensions. Our approach is to use standard Runge-Kutta linear stability theory, treating the nonlinearity of the NLSE as a constant. The required bounds on the eigenvalues of the scheme matrices are found analytically when possible, and otherwise estimated using the Gershgorin circle theorem.

Keywords: Numerical stability, Explicit finite difference schemes, Nonlinear Schrödinger equation.

1 Introduction

The nonlinear Schrödinger equation (NLSE) is used to model a wide variety of physical systems since it describes, to least nonlinear order, modulated wave propagation [10]. The general form of the NLSE can be written as

i\frac{\partial\Psi}{\partial t}+a\nabla^{2}\Psi-V({\bf r})\Psi+s|\Psi|^{2}\Psi=0,

(1)

where $\Psi\in\mathbbm{C}$ is the value of the wavefunction, $\nabla^{2}$ is the Laplacian operator, and where $a>0$ and $s$ are parameters defined by the system being modeled. $V(\bf{r})$ is an external potential term, which when included, makes Eq. (1) known as the Gross-Pitaevskii equation [13].

Often, efficient and easy-to-use numerical methods are employed to simulate the NLSE. One such method is the method of lines where the time-stepping and spatial differencing are treated independently. This transforms the partial differential equation (PDE) into a large number of coupled ordinary differential equations (ODEs). These ODEs can then be solved using a variety of numerical schemes, one of the most common being the fourth order Runge-Kutta (RK4) scheme [2]. Using the RK4 scheme with the NLSE produces a fully explicit scheme where each grid point at time $t$ is only a function of values at time $t-{\bf k}$ where ${\bf k}$ is the time-step. This simplifies computational implementations because no matrices are needed to be formed and stored, and no linear systems are needed to be solved (which in the nonlinear case also require a nonlinear iterative process).

The only drawback to using explicit finite-difference schemes (such as the RK4) for simulating PDEs is that they are conditionally stable. This means that there is an upper bound on the allowed size of the time-step which is dependent on the spatial-step size. If the time-step is larger than this bound, the scheme is unstable and diverges [17]. Although rough estimates of the stability bound can be found through an inefficient educated guess-and-check, for higher dimensional scenarios, as well as long and/or large simulations, a more refined and predictable stability bound is essential for efficient simulations.

In this paper, we formulate linearized stability bounds for simulating the NLSE with the RK4 scheme. The stability bounds depend on the specific spatial differencing scheme being used, as well as on the boundary conditions. We formulate the bounds for both second-order and fourth-order spatial differencing with a variety of boundary conditions (Dirichlet, modulus-squared Dirichlet (MSD), Laplacian-zero (L0), and periodic). Each analysis is done for one, two and three dimensions.

The paper is organized as follows. In Sec. 2 we review the basic RK4 stability properties and apply the results to the NLSE to formulate general stability bounds. Our basic procedure in finding the stability bounds and the linear algebra theorems that we utilize are also discussed. In Sec. 3 we summarize the forms of the boundary conditions we consider. Our main analysis begins in Sec. 4 with the one-dimensional NLSE. Linearized stability bounds are found for each scheme and boundary condition combination. In Sec. 5 and 6, we use the same procedures to formulate the bounds for the two- and three-dimensional NLSE respectively. A few numerical examples showing the accuracy of the bounds are shown in Sec. 7. In Sec. 8, we conclude and summarize all the results from Secs. 4, 5, and 6 into a concise reference.

2 Stability theory

2.1 General Runge-Kutta scheme stability

Given an initial value problem of a set of linear first-order ODEs (in our case, a method-of-lines explicit PDE finite-difference scheme), one can formulate the matrix notation

\frac{\partial\vec{\Psi}}{\partial t}={\cal A}\,\vec{\Psi},

(2)

where ${\cal A}$ contains the coefficients of the right-hand-sides of the ODEs. We now define

\vec{p}={\bf k}\,\vec{\lambda},

(3)

where ${\bf k}$ is the time-step size and $\vec{\lambda}$ contains the eigenvalues of ${\cal A}$ . In our case, the eigenvalues of ${\cal A}$ will have the spatial-step size (denoted $h$ ) included in them, as well as any parameters of the NLSE. As shown in Ref. [16], for the fourth-order Runge-Kutta scheme, if a vector $\vec{R}(\vec{p})$ is defined whose elements are the polynomials

R(p)=1+p+\frac{p^{2}}{2}+\frac{p^{3}}{6}+\frac{p^{4}}{24},

(4)

then the stability of the RK4 scheme is guaranteed if

\lVert\vec{R}(\vec{p})\rVert_{\infty}<1,

(5)

where $\lVert\;\rVert_{\infty}$ denotes the infinity norm defined as $\lVert\vec{x}\rVert_{\infty}=\mbox{max}\{\lvert x_{0}\rvert,\lvert x_{1}\rvert,...,\lvert x_{N-1}\rvert\}$ . Inserting Eq. (3) into Eq. (4) yields

	$\displaystyle\left\|R(\lambda)\right\|^{2}=1$	$\displaystyle+\frac{1}{576}\,{\bf k}^{8}\|\lambda\|^{8}-\frac{1}{72}\,{\bf k}^{6}\|\lambda\|^{6}+\left(\frac{{\bf k}^{6}\|\lambda\|^{6}}{6}-{\bf k}^{4}\|\lambda\|^{4}+24\right)\,\frac{{\bf k}}{12}\,\mbox{Re}(\lambda)$		(6)
		$\displaystyle+\left({\bf k}^{4}\|\lambda\|^{4}+24\right)\,\frac{{\bf k}^{2}}{12}\,(\mbox{Re}(\lambda))^{2}+\left({\bf k}^{2}\|\lambda\|^{2}+4\right)\,\frac{{\bf k}^{3}}{3}\,(\mbox{Re}(\lambda))^{3}+\frac{2{\bf k}^{4}}{3}\,(\mbox{Re}(\lambda))^{4}.$

In Fig. 1 we show the stability region for the RK4 scheme given by Eq. (5) as well as that for lower-order Runge-Kutta schemes (whose $R(p)$ is defined by progressively truncated versions of Eq. (4)) [16].

Refer to caption — Figure 1: (Color online) Left: Stability regions for Runge-Kutta schemes. Schemes from first-order to fourth-order are shown from center outwards. Right: Magnified view of the same plot near the point where $\mbox{Re}(\lambda)=0$ .

As we shall show, the eigenvalues of the $\cal A$ matrix are all purely imaginary (or nearly so) in the case of the nonlinear (and linear) Schrödinger equation. Thus, it can be seen from Fig. 1 that the third-order Runge-Kutta is the lowest-order RK scheme that is conditionally stable for the Schrödinger equations (however, as shown in Ref. [8], this is not the case if the real and imaginary parts of the NLSE are computed in a staggered time grid, or, as in Ref. [9], an artificial dissipative term is added to the NLSE), with the RK4 yielding a significantly larger bound on ${\bf k}$ . This is in contrast to similar PDEs such as the heat equation, whose $\cal A$ matrix eigenvalues are typically all real-valued, in which case even forward differencing (RK1) is conditionally stable.

If, as in our case, $\mbox{Re}(\vec{\lambda})=\vec{0}$ , Eq. (6) simplifies greatly and becomes

\left|R(\vec{\lambda})\right|^{2}=1+\frac{1}{576}\,{\bf k}^{8}|\vec{\lambda}|^{8}-\frac{1}{72}\,{\bf k}^{6}|\vec{\lambda}|^{6},

(7)

in which case, Eq. (5) leads to the simple stability bound

{\bf k}<\frac{\sqrt{8}}{\lVert\vec{\lambda}\rVert_{\infty}}.

(8)

2.2 Application to the NLSE

Applying the above stability theory to the NLSE has the obvious problem that the analysis is purely linear, while the NLSE has one (or more) nonlinear terms. A full nonlinear stability analysis is beyond the scope of this paper, so instead we linearize the problem by treating the nonlinearity ( $|\Psi|^{2}$ ) as a constant value $U$ (the external potential term is usually a constant independent of $\Psi$ at each grid point, and so does not need any special treatment). This has been done previously for the one-dimensional coupled NLSE for fourth-order differencing (in the exclusive case where $s<0$ ) in Ref. [12]. Since the value of $|\Psi|^{2}$ changes over time during the simulation, the linearized stability bound will also change over time. This change in many cases is expected to be small (which we have confirmed in numerical simulations, not reported here) and therefore can be ignored, i.e. one may compute the bound using the initial condition of $\Psi$ (and $V(\bf{r})$ ) and just leave a few percent leeway to cover any changes. This is especially true in the repulsive case ( $s<0$ ) where most situations have a constant-density background (or maximum background) and the dynamics do not cause the maximum background value to change significantly (for example, when simulated coherent structures, most of the dynamics are translations of the initial condition with little change in structure). In attractive cases ( $s>0$ ), blow-up can occur which can alter the stability bound greatly, causing the simulation to crash (although in such a case the wavefunction is exploding towards infinity, which most finite-difference schemes cannot handle anyways). Many times, simulations of a steady-state or near-steady-state in the modulus-squared with a constant potential are performed. In such situations, the linearized stability bounds will be (nearly) exact.

It is also useful to formulate stability bounds for the linear Schrödinger equation (LSE) (where $s=0$ and $V({\bf r})=0$ ). In addition to providing bounds for the LSE, as will be discussed below, the results can also be used as practical estimates of the stability bounds for the NLSE (the discrepancy can often be solved by lowering the bound by a few percent).

2.3 Stability analysis procedure

In order to simplify the analysis, we first rewrite Eq. (2) as

\frac{\partial\vec{\Psi}}{\partial t}={\cal A}\vec{\Psi}=\frac{i\,a}{h^{2}}\,A\vec{\Psi},

where $h$ is the step-size of the spatial finite-difference scheme being used. Then, assuming all eigenvalues of $A$ are real-valued, the stability condition of Eq. (8) becomes

{\bf k}<\frac{\sqrt{8}}{\lVert\vec{\lambda}_{A}\rVert_{\infty}}\frac{h^{2}}{a}.

(9)

In order to be able to use the stability bound of Eq. (9), we must first confirm that all eigenvalues of $A$ are purely real (or nearly so) for each scheme/boundary condition combination. In cases where the eigenvalues are not able to be easily computed analytically, we show that the $A$ matrix’s eigenvalues are a set of boundary values with the remaining eigenvalues being those of a symmetric matrix denoted $A^{\prime}$ . Then by Thm. 1, it is known that all the eigenvalues of $A$ are real.

Theorem 1 (Ref. [1]).

The eigenvalues of a real symmetric matrix are real.

Once it has been established that Eq. (9) can be used, in order to get an upper-bound on ${\bf k}$ , we require an upper-bound on the maximum absolute eigenvalue of $A$ . Due to the sparsity and diagonal dominance of $A$ , a good estimate of the upper-bound can be found using the Gershgorin circle theorem (Def. 1 and Thm. 2).

Definition 1 (Ref. [11]).

Let $A$ be a square complex matrix. Around every element $a_{ii}$ on the diagonal of the matrix, a circle with radius equal to the sum of the norms of the other elements in the same row ( $\sum_{j\neq i}|a_{ij}|$ ) is known as a Gershgorin disc.

Theorem 2 (Ref. [11]).

Every eigenvalue of a square complex matrix A lies in one of its Gershgorin discs.

Since every eigenvalue must be contained in a Gershgorin disk, by finding the maximum absolute value of the limits of the disks will yield an upper-bound on the maximum modulus of the eigenvalues of $A$ .

In the one-dimensional LSE case with no external potential and periodic boundary conditions, the $A$ matrix becomes circulant as defined by Def. 2. In this case, the eigenvalues can be computed analytically by Thm. 3. The upper-bound is then taken by finding the limit of the maximum eigenvalue as the size of the matrix goes to infinity.

Definition 2 (Ref. [15]).

A circulant matrix is a square $N\times N$ matrix $C$ that can be fully specified by one vector, $\vec{c}=\{c_{0},c_{1},...,c_{N-1}\}$ , which appears as the first column of $C$ . The remaining columns of $C$ are each cyclic permutations of the vector with the offset equal to the column index.

Theorem 3 (Ref. [15]).

The eigenvalues of a circulant matrix are given by

\lambda_{j}=c_{0}+c_{N-1}\,\omega_{j}+c_{N-2}\,\omega_{j}^{2}+...+c_{1}\,\omega_{j}^{N-1},\qquad j=0,...,{N-1},

where

\omega_{j}=\exp\left(\frac{2\pi\,i\,j}{N}\right).

3 Boundary conditions

Since boundary conditions of the spatial differencing in a PDE like the NLSE have the potential to alter the stability of a scheme, it is necessary to have stability results for each specific boundary condition one would like to use. In this paper we limit ourselves to four boundary conditions which we feel are a good combination of simplicity and usefulness. These boundary conditions are: periodic, Dirichlet, Laplacian-zero, and Modulus-Squared-Dirichlet. As notation, we use the subscript $b$ to represent any boundary point, and $b-1$ to represent the grid position one point inward from the boundary in the normal direction.

For use with the stability analysis, it is desirable to formulate each boundary condition in terms of the temporal derivative in the form

\left.\frac{\partial\Psi}{\partial t}\right|_{b}=\frac{i\,a}{h^{2}}B_{b}\Psi_{b},

(10)

and in terms of the spatial Laplacian in the form

\nabla^{2}\Psi_{b}=\frac{1}{h^{2}}D_{b}\Psi_{b},

(11)

where $B_{b}$ and $D_{b}$ are assumed to be real-valued constants (possibly differing per boundary point) and defined based on the specific boundary condition being used. For periodic boundary conditions (or linear one-sided conditions not discussed here), these forms are not applicable. Writing the boundary conditions in the forms of Eq. (10) and Eq. (11) allows them to be expressed in the $A$ matrix as a single real-valued entry ( $B_{b}$ ), and in the case of the form of the fourth-order differencing chosen here (see Sec. 4.2), the near-boundary interior points will contain $D_{b}$ in their formulation.

3.1 Periodic

For periodic boundary conditions, any element of the scheme that is too small or too large in index (i.e. they are ‘off the grid’) are simply replaced by the grid points on the opposite side of the grid. In the case of the NLSE, periodic boundary conditions can be problematic especially in background-density situations due to the unpredictable phase jump from one side of the grid to the other.

3.2 Dirichlet

Dirichlet boundary conditions are defined as

\Psi_{b}=C,

where $C$ is a constant. In terms of the temporal derivative of the NLSE, this condition is

\left.\frac{\partial\Psi}{\partial t}\right|_{b}=0,

in which case $B_{b}=0$ in Eq. (10). When inserted into the NLSE, this condition in terms of the Laplacian is given by

\nabla^{2}\Psi_{b}=-\frac{1}{a}(s|\Psi_{b}|^{2}-V_{b})\Psi_{b},

and therefore $D_{b}=-h^{2}/a(s|\Psi_{b}|^{2}-V_{b})$ in Eq. (11).

3.3 Modulus-squared Dirichlet

In some situations Dirichlet boundary condition can fail. Such failure typically occurs in simulations with a constant-density background, i.e. a constant value of $|\Psi|^{2}$ at the boundaries. A standard Dirichlet condition will not work in such cases because it does not take into account the phase rotation of $\Psi$ . Instead, one would like to have the modulus-squared of the wavefunction to be constant at the boundaries, i.e.

|\Psi_{b}|^{2}=C,

where $C$ is a constant. We have recently formulated a method for such a boundary condition (which is almost as easy to implement as Dirichlet) called the modulus-squared Dirichlet boundary condition [5]. The MSD boundary condition is given in terms of the temporal derivative of the NLSE as

\Psi_{t,b}\approx i\,\mbox{Im}\left[\frac{\Psi_{t,b-1}}{\Psi_{b-1}}\right]\,\Psi_{b}.

(12)

where $\partial\Psi_{b-1}/\partial t$ is computed by the interior scheme first, and then used to compute the boundary values. Using the MSD boundary condition gives $B_{b}=(h^{2}/a)\mbox{Im}\left[\left.\frac{\partial\Psi}{\partial t}\right|_{b-1}\frac{1}{\Psi_{b-1}}\right]$ , which is nonlinear, and not a constant independent of $\Psi$ . As shown in Ref. [5], due to the underlying assumptions of the MSD boundary condition, Eq. (12) can be viewed as

\left.\frac{\partial\Psi}{\partial t}\right|_{b}\approx i\,\Omega_{b-1}\Psi_{b},

where $\Omega_{b-1}$ is the real-valued frequency of the solution near the boundary. Thus, $B_{b}$ would have the form $B_{b}=(h^{2}/a)\Omega_{b-1}$ . Therefore, we can linearize the MSD boundary condition by treating the $B_{b}$ term as a constant (which can change over the course of the simulation, similar to the nonlinearity of the NLSE).

When inserted into the NLSE, the MSD boundary condition of Eq. (12) yields

\nabla^{2}\Psi_{b}\approx\left[\mbox{Im}\left(i\,\frac{\nabla^{2}\Psi_{b-1}}{\Psi_{b-1}}\right)+\frac{1}{a}\,\left(N_{b-1}-N_{b}\right)\right]\Psi_{b},

(13)

where

N_{b}=s\,|\Psi_{b}|^{2}-V_{b},\qquad N_{b-1}=s\,|\Psi_{b-1}|^{2}-V_{b-1}.

(14)

and therefore $D_{b}=h^{2}\left[\mbox{Im}\left(i\,\frac{\nabla^{2}\Psi_{b-1}}{\Psi_{b-1}}\right)+\frac{1}{a}\,\left(N_{b-1}-N_{b}\right)\right]$ . This too is a nonlinear, non-constant term, and so must be treated as a constant in the same manner as the nonlinearity of the NLSE.

3.4 Laplacian-zero boundary condition

The Laplacian-zero boundary condition is defined by

\nabla^{2}\Psi_{b}=0,

and therefore $D_{b}=0$ . In terms of the time-derivative of the NLSE, the L0 boundary condition is given as

\left.\frac{\partial\Psi}{\partial t}\right|_{b}=i\,(s|\Psi_{b}|^{2}-V_{b})\,\Psi_{b},

making $B_{b}=(h^{2}/a)(s|\Psi_{b}|^{2}-V_{b})$ . This condition is as easy to implement as the Dirichlet, and can be useful in many situations.

To assist the stability analysis, a summary of the values of $B_{b}$ and $D_{b}$ for all the mentioned boundary conditions are given in Table. 1 for future reference.

Table 1: Boundary condition terms for use with stability analysis.

Boundary Condition	$B_{b}$	$D_{b}$
Dirichlet	$0$	$\dfrac{h^{2}}{a}(V_{b}-s\|\Psi_{b}\|^{2})$
Laplacian-zero	$\dfrac{h^{2}}{a}(s\|\Psi_{b}\|^{2}-V_{b})$	$0$
MSD	$\dfrac{h^{2}}{a}\mbox{Im}\left[\dfrac{\Psi_{t,b-1}}{\Psi_{b-1}}\right]$	$h^{2}\left[\mbox{Im}\left(i\,\dfrac{\nabla^{2}\Psi_{b-1}}{\Psi_{b-1}}\right)+\dfrac{1}{a}\,\left(N_{b-1}-N_{b}\right)\right]$

Many other boundary conditions exist for simulating the NLSE, in which case the analysis shown in this paper can be adapted to the other boundary conditions.

4 One-dimensional stability analysis

In the one-dimensional cases we analyze all four boundary conditions mentioned in Sec. 3. As stated, periodic boundary conditions yield a matrix where (in the linear case with $s=0$ and $\vec{V}({\bf r})=0$ ) the eigenvalues can be computed analytically. This allows the results obtained using the upper-bound methods (which we use with other boundary conditions) to be compared with the true eigenvalues giving an idea of how accurate they are.

4.1 Second-order central difference

The second-order central difference is one dimension is given by

\nabla^{2}\Psi_{i}=\left.\frac{\partial^{2}\Psi}{\partial x^{2}}\right|_{i}\approx\frac{\Psi_{i+1}-2\Psi_{i}+\Psi_{i-1}}{h^{2}},

and when implemented into the $A$ matrix, forms a matrix which is tridiagonal (except for the two boundary condition rows).

4.1.1 Periodic boundary conditions

In order to obtain analytic expressions for the eigenvalues of $A$ , we start with the LSE case with no external potential and periodic boundary conditions. This yields the matrix

A=\left[\begin{array}[]{cccccc}-2&1&0&0&1\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&-2&1&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\ddots&\ddots&\ddots&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&1&-2&1\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&0&0&1&-2\end{array}\right],

which, as per Def. 2, is a circulant matrix with $\vec{c}=\{-2,1,0,...,0,1\}$ . Also, since $A$ is a real-valued symmetric matrix, by Thm. 1, all eigenvalues are real and therefore the stability criteria of Eq. (9) can be used. By Thm. 3, the eigenvalues of $A$ are given by

\lambda_{j}=-2+\exp\left[\frac{2\pi ij}{N}\right]+\exp\left[\frac{2\pi ij(N-1)}{N}\right],\qquad j\in\{0,...N-1\}.

The maximum value of $|\lambda_{j}|$ occurs either at $j=N/2$ if $N$ is even, or $j=(N\pm 1)/2$ if $N$ is odd. For $N$ even-valued we have

|\lambda|_{\max}=\left|-2+\exp\left[\pi i\right]+\exp\left[\pi i\right]^{N-1}\right|,

which yields

|\lambda|_{\max}=4.

For $N$ odd-valued we have

|\lambda|_{\max}=\left|-2-(-1)^{1/N}+(-1)^{N}\,(-1)^{-1/N}\right|,

which yields

|\lambda|_{\max}=\left|-2-2\,\cos\left(\frac{\pi}{N}\right)\right|.

Taking $N\rightarrow\infty$ , the maximum bound on the maximum absolute eigenvalue becomes

|\lambda|_{\max}<4.

We therefore have an upper bound on the maximum absolute eigenvalue which, for even-valued $N$ , is guaranteed to be one of the eigenvalues. The stability criteria of Eq. (9) is then formulated as

{\bf k}<\frac{\sqrt{8}}{4}\frac{h^{2}}{a}.

(15)

In the general case where $s\neq 0$ and/or $V({\bf r})\neq 0$ , the $A$ matrix is no longer circulant (since the values of the nonlinearity or external potential vary over the diagonal of $A$ ). To get a bound on the maximum absolute eigenvalue, we make use of Thm. 2. The matrix $A$ has $N$ Gershgorin disks, since each diagonal entry of $A$ can be unique, but each disk has the same radius ( $r=2$ ). Also, since the diagonal entries can in theory take on any value, all Gershgorin disk limits must be examined. This yields the stability bound

{\bf k}<\frac{\sqrt{8}}{\max\{\lVert\vec{L}\rVert_{\infty},\lVert\vec{L}-4\rVert_{\infty}\}}\,\frac{h^{2}}{a},

(16)

where we have defined the elements of $\vec{L}$ to be

L_{i}=\frac{h^{2}}{a}(s\,|\Psi_{i}|^{2}-V_{i}),

(17)

where the index $i$ spans over the entire grid. It is important to note that all values of $\vec{L}$ are $O(h^{2})$ . Thus, for $h\ll 1$ , and reasonable values of $|\Psi|^{2}$ and $\vec{V}$ , the linear bound of Eq. (15) should be very close to the true bound of the nonlinear problem.

If we set $\vec{L}=0$ in Eq. (16), we recover the bound in Eq. (15). This shows that (in this case at least), using the Gershgorin circle theorem yields the true bound on the eigenvalues of $A$ .

4.1.2 Dirichlet, MSD, and L0 boundary conditions

As shown in Sec. 3, Dirichlet, Laplacian-zero, and modulus-squared Dirichlet boundary conditions can all be viewed as single entries in the boundary value rows of the $A$ matrix, denoted as $B_{b}$ . As shown there, the values of $B_{b}$ are real-valued and their values for each boundary condition were given in Table 1. Using such a formulation, the $A$ matrix becomes

A=\left[\begin{array}[]{cccccc}B_{0}&0&0&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&L_{1}-2&1&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\ddots&\ddots&\ddots&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&1&L_{N-2}-2&1\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&0&0&B_{N-1}\end{array}\right].

In order to use the simple stability criteria of Eq. (9), we once again need to show that all eigenvalues of $A$ are purely real. The $A$ matrix is no longer symmetric, however it is easy to see that $B_{0}$ and $B_{N-1}$ are eigenvalues of $A$ , and the remaining eigenvalues of $A$ are equivalent to the eigenvalues of the matrix $A^{\prime}$ defined as

A^{\prime}=\left[\begin{array}[]{cccccc}L_{1}-2&1&0&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 1&L_{2}-2&1&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\ddots&\ddots&\ddots&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&1&L_{N-3}-2&1\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&0&1&L_{N-2}-2\end{array}\right].

Since $A^{\prime}$ is real-valued symmetric, we can use the stability bound of Eq. (9).

We now need to find an upper bound on the absolute value of the eigenvalues of $A^{\prime}$ . We use the Gershgorin circle theorem to find all unique Gershgorin disks and take the limits of the disks to find the bounds on the absolute eigenvalues. Many of the Gershgorin disks are similar, differing only in the value of $L_{i}$ of the specific row. Therefore, each disk of different centers and radii has a subset of $L_{i}$ values relevant to it. Although in the current one-dimensional setting it is simple to define the subsets, in higher-dimensional settings, it can become burdensome to separate out each subset of $\vec{L}$ relevant to each Gershgorin disk of the same center and radius. Therefore, for practicality purposes, we define our bounds using all possible values of $L_{i}$ for each Gershgorin disk center and radius. This may make the resulting stability bound slightly higher than necessary in certain cases, but this is outweighed by the ease-of-use of the simplified bounds. The unique forms of the Gershgorin disks of $A^{\prime}$ are shown in Table 2.

Table 2: Unique forms of the Gershgorin disk centers (

a_{ii}

) and radii (

r_{i}

) for the

A^{\prime}

matrix of the one-dimensional second-order central difference scheme.

$a_{ii}$	$r_{i}=\sum_{i\neq j}\|a_{ij}\|$
$L_{i}-2$	$1$
$L_{i}-2$	$2$

The resulting general stability bounds are

{\bf k}<\frac{\sqrt{8}}{\max\{\lVert\vec{B}\rVert_{\infty},\lVert\forall L_{i},L_{i}-\vec{G}\rVert_{\infty}\}}\,\frac{h^{2}}{a},

(18)

where $\vec{B}$ are all boundary condition values (in this case $B_{0}$ and $B_{N-1}$ ), and $\vec{G}$ is defined as

\vec{G}=\left\{4,3,1,0\right\}.

(19)

In general, all possible values of $\vec{G}$ must be taken into consideration since there is no theoretical restriction on what values $\vec{L}$ can take. However, in certain specific circumstances, some of the values of $\vec{G}$ can be ignored (for example, when $s\leq 0$ and $V({\bf r})\geq 0$ , only the largest magnitude value in $\vec{G}$ is needed).

4.2 Fourth order central difference

The standard fourth order central difference scheme is given by

\nabla^{2}\Psi_{i}=\left.\frac{\partial^{2}\Psi}{\partial x^{2}}\right|_{i}\approx\frac{-\Psi_{i+2}+16\Psi_{i+1}-30\Psi_{i}+16\Psi_{i-1}-\Psi_{i-2}}{12\,h^{2}}.

(20)

The stability analysis follows directly from the second-order case. The only major difference is that since the fourth order stencil is five points wide, the grid points near the boundary may need special consideration for the different boundary conditions. For our purposes here, we use the two-step high-order compact (2SHOC) version of the fourth-order scheme as described in Ref. [6], in which case the near-boundary points can be formulated by combining the two steps of the 2SHOC scheme.

4.2.1 Periodic boundary condition

In the periodic case, no special attention is needed near the boundaries, and the $A$ matrix in the LSE case with no external potential is

A=\left[\begin{array}[]{cccccccc}-15/6&4/3&-1/12&0&0&-1/12&4/3\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 4/3&-15/6&4/3&-1/12&0&0&-1/12\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr-1/12&4/3&-15/6&4/3&-1/12&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\ddots&\ddots&\ddots&\ddots&\ddots&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&-1/12&4/3&-15/6&4/3&-1/12\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr-1/12&0&0&-1/12&4/3&-15/6&4/3\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 4/3&-1/12&0&0&-1/12&4/3&-15/6\end{array}\right],

which is a circulant matrix, and its eigenvalues are therefore

	$\displaystyle\lambda_{j}=$	$\displaystyle-\frac{15}{6}+\frac{4}{3}\,\exp\left[\frac{2\pi ij}{N}\right]-\frac{1}{12}\,\exp\left[\frac{4\pi ij}{N}\right]$
		$\displaystyle-\frac{1}{12}\,\exp\left[\frac{2(N-2)\pi ij}{N}\right]+\frac{4}{3}\,\exp\left[\frac{2(N-1)\pi ij}{N}\right].$

The maximum absolute value once again occurs occurs at either $j=N/2$ if $N$ is even, or $j=(N\pm 1)/2$ if $N$ is odd. For $N$ even-valued we have

\lambda_{N/2}=-\frac{15}{6}-\frac{4}{3}-\frac{1}{12}-\frac{1}{12}(-1)^{N-2}+\frac{4}{3}(-1)^{N-1}=-\frac{16}{3}.

For $N$ odd, we have

\lambda_{(N+1)/2}=-\frac{15}{6}-\frac{4}{3}\,\left((-1)^{1/N}+(-1)^{-1/N}\right)-\frac{1}{12}\left((-1)^{2/N}+(-1)^{-2/N}\right),

which yields

\lambda_{(N+1)/2}=-\frac{15}{6}-\frac{4}{3}\,\left(2\,\cos\left(\frac{\pi}{N}\right)\right)-\frac{1}{12}\,\left(2\,\cos\left(\frac{2\pi}{N}\right)\right).

As $N\rightarrow\infty$ , $|\lambda|\rightarrow\frac{16}{3}$ , which is the same bound as the $N$ -even case. Thus, the stability bound is given by

{\bf k}<\left(\frac{3}{4}\right)\frac{\sqrt{8}}{4}\frac{h^{2}}{a},

(21)

which we note is only a $25\%$ reduction of the second-order bound of Eq. (15). In the general case where $\vec{L}\neq 0$ , the bound becomes

{\bf k}<\frac{6\sqrt{2}}{\max\{\lVert 3\vec{L}-16\rVert_{\infty},\lVert 3\vec{L}+1\rVert_{\infty}\}}\frac{h^{2}}{a}.

(22)

If $s\leq 0$ and $V({\bf r})\geq 0$ , the first term in the denominator is the maximum of the two terms, and the resulting stability bound is equivalent to that found in Ref. [12] for using the RK4 scheme and fourth-order spatial differencing with the coupled NLSE.

4.2.2 Dirichlet, MSD, and L0 boundary conditions

As per Sec. 3, we formulate all three boundary conditions in terms of a $B_{b}$ entry in the $A$ matrix. As discussed, an important issue is that we need to handle the grid points near the boundary due to the width of the scheme. A common way of dealing with the closest-interior points is to compute the Laplacian at those points using second-order differencing, however this can lead to the overall scheme becoming second-order. However, since we are using the 2SHOC version of the fourth-order differencing, we can derive the closest-interior points which, if the assumptions of the chosen boundary conditions hold, should maintain fourth-order accuracy. In one dimension, the 2SHOC scheme is defined as [6]

	$\displaystyle 1)\qquad$	$\displaystyle D_{i}$	$\displaystyle=\frac{1}{h^{2}}\left(\Psi_{i+1}-2\Psi_{i}+\Psi_{i-1}\right),$		(23)
	$\displaystyle 2)\qquad$	$\displaystyle\nabla^{2}\Psi_{i}$	$\displaystyle\approx\frac{7}{6}D_{i}-\frac{1}{12}\left(D_{i+1}+D_{i-1}\right).$		(24)

In the first step, the second-order Laplacian is computed with the chosen boundary condition applied to it. Next, the result is used to compute the fourth-order Laplacian. As mentioned in Ref. [6], this two-step scheme is equivalent to the standard wide-stencil of Eq. (20) for the interior points. We use the form of Eq. (11) for the boundary conditions on the Laplacian, and after combining the steps of Eq. (23) and Eq. (24), we get the matrix

A=\left[\begin{array}[]{cccccccc}B_{0}&0&0&0&0&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\frac{14-D_{0}}{12}&L_{1}-\frac{29}{12}&\frac{4}{3}&-\frac{1}{12}&0&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr-\frac{1}{12}&\frac{4}{3}&L_{2}-\frac{15}{6}&\frac{4}{3}&-\frac{1}{12}&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\ddots&\ddots&\ddots&\ddots&\ddots&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&-\frac{1}{12}&\frac{4}{3}&L_{N-3}-\frac{15}{6}&\frac{4}{3}&-\frac{1}{12}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&0&-\frac{1}{12}&\frac{4}{3}&L_{N-2}-\frac{29}{12}&\frac{14-D_{N-1}}{12}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&0&0&0&0&B_{N-1}\end{array}\right],

where the $B_{b}$ and $D_{b}$ ( $b\in\{0,N-1\}$ ) terms for each different boundary condition are once again given in Table 1. As in Sec. 4.1.2, the $A$ matrix is not symmetric and has eigenvalues equal to $\vec{B}$ . The remaining eigenvalues are those of the matrix $A^{\prime}$ defined as

A^{\prime}=\left[\begin{array}[]{ccccccc}L_{1}-\frac{29}{12}&\frac{4}{3}&-\frac{1}{12}&0&0&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr\frac{4}{3}&L_{2}-\frac{15}{6}&\frac{4}{3}&-\frac{1}{12}&0&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr-\frac{1}{12}&\frac{4}{3}&L_{3}-\frac{15}{6}&\frac{4}{3}&-\frac{1}{12}&0&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\ddots&\ddots&\ddots&\ddots&\ddots&0\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&-\frac{1}{12}&\frac{4}{3}&L_{N-4}-\frac{15}{6}&\frac{4}{3}&-\frac{1}{12}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&0&-\frac{1}{12}&\frac{4}{3}&L_{N-3}-\frac{15}{6}&\frac{4}{3}\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&0&0&-\frac{1}{12}&\frac{4}{3}&L_{N-2}-\frac{29}{12}\end{array}\right],

which is once again real-symmetric so the bounds of Eq. (9) can be used. It is interesting to note that the values of $D_{b}$ do not appear in any of the eigenvalues of $A^{\prime}$ .

Table 3: Unique forms of the Gershgorin disk centers (

a_{ii}

) and radii (

r_{i}

) for the

A^{\prime}

matrix of the one-dimensional fourth-order 2SHOC scheme.

$a_{ii}$	$r_{i}=\sum_{i\neq j}\|a_{ij}\|$
$L_{i}-5/2$	$11/4$
$L_{i}-5/2$	$17/6$
$L_{i}-29/12$	$17/12$

The unique forms (see the discussion in Sec. 4.1.2) of the Gershgorin disk centers and radii of $A^{\prime}$ are shown in Table 3. The full stability bound is the same form as Eq. (18), but with $\vec{G}$ defined as

\vec{G}=\frac{1}{12}\times\left\{64,63,46,12,-3,-4\right\}.

(25)

Once again, in the most general case, all values of Eq. (25) must be considered in finding the maximum allowed time-step value.

5 Two-dimensional stability analysis

In higher dimensions, the $A$ matrix is formed by unwrapping the solution into a one-dimensional vector and then formulating the scheme matrix accordingly.

In Secs. 4.1.1 and 4.2.1 we noted that the stability bounds given using the Gershgorin circle theorem were equivalent to those obtained analytically for the linear case with periodic boundary conditions. We therefore justify relying exclusively on the Gershgorin theorem for higher dimensions, and focus on the stability bounds for Dirichlet, MSD, and Laplacian-zero boundary conditions (since the periodic boundary condition bounds will be a subset of the bounds computed for the other boundary conditions).

5.1 Second-order central difference

The second-order central difference scheme in two dimensions is given by

\qquad\nabla^{2}\Psi_{i,j}=\left.\dfrac{\partial^{2}\Psi}{\partial x^{2}}\right|_{i,j}+\left.\dfrac{\partial^{2}\Psi}{\partial y^{2}}\right|_{i,j}\approx\dfrac{1}{h^{2}}

	$1$
$1$	$-4$	$1$
	$1$

\Psi_{i,j}.

(26)

The corresponding $A$ matrix has a tri-banded structure, with diagonal sub-sections corresponding to the boundary values. The form of the $A$ matrix is shown in Fig. 2 (we do not show the values of the entries of the matrix due to space considerations, but they can be obtained through symbolic math codes). As in the one-dimensional case, all diagonal entries (which are the boundary value entries $B_{b}$ ) of $A$ are eigenvalues, and the remaining eigenvalues are real and equivalent to those of a matrix $A^{\prime}$ which is real-symmetric, thus allowing the use of the bounds in Eq. (9). The form of $A^{\prime}$ is also shown in Fig. 2.

The unique forms of the Gershgorin disk centers and radii for $A^{\prime}$ are shown in Table 4.

Table 4: Gershgorin disk centers and radii for the

A^{\prime}

matrix of the two-dimensional central-difference scheme.

$a_{ii}$	$r_{i}=\sum_{i\neq j}\|a_{ij}\|$
$L_{i}-4$	$2$
$L_{i}-4$	$3$
$L_{i}-4$	$4$

The stability bounds are once again the same as in Eq. (18) with $\vec{B}$ being the set of all boundary values $B_{b}$ and with $\vec{G}$ now defined as

\vec{G}=\{0,1,2,6,7,8\}.

(27)

In the linear case ( $s=0$ ) with no external potential and periodic, Dirichlet or Laplacian-zero boundary conditions, we get the linear stability bound

{\bf k}<\frac{\sqrt{8}}{8}\,\frac{h^{2}}{a}.

(28)

As before, since within the $A$ matrix, all boundary, potential, and nonlinear terms are $O(h^{2})$ the simple bound of Eq. (28) with slight adjustment can be used in many applications.

5.2 Fourth-order central difference

The fourth-order central difference scheme in two dimensions is given by

\nabla^{2}\Psi_{i,j}\approx-\dfrac{1}{12\,h^{2}}

		$1$
		$-16$
$1$	$-16$	$60$	$-16$	$1$
		$-16$
		$1$

\Psi_{i,j}.

(29)

The low-storage version of the 2SHOC equivalent scheme is defined as [6]

			(34)

			(42)

The corresponding $A$ matrix has a five-banded structure. The structure of the $A$ matrix and its corresponding $A^{\prime}$ matrix are shown in Fig. 3.

The resulting unique forms of the Gershgorin disk centers and radii for $A^{\prime}$ are shown in Table 5.

Table 5: Gershgorin disk centers (

a_{ii}

) and radii (

r_{i}

) for the

A^{\prime}

matrix of the two-dimensional 2SHOC scheme.

$a_{ii}$	$r_{i}=\sum_{i\neq j}\|a_{ij}\|$
$L_{i}-5$	$11/2$
$L_{i}-5$	$67/12$
$L_{i}-5$	$17/3$
$L_{i}-29/6$	$17/6$
$L_{i}-59/12$	$25/6$
$L_{i}-59/12$	$17/4$

The linearized stability bound is once again Eq. (18) but with $\vec{G}$ defined as

\vec{G}=\dfrac{1}{12}\times\left\{128,127,126,110,109,92,24,9,8,-6,-7,-8\right\}.

(43)

The linear bound (with $s=0$ and $V({\bf r})=0$ ) is then given by

{\bf k}<\frac{3\sqrt{8}}{32}\,\frac{h^{2}}{a},

(44)

which, as in the one-dimensional case, is $25\%$ lower than the second-order linear bound given in Eq. (28).

6 Three-dimensional stability analysis

For the stability analysis in three dimensions, the same procedure utilized in the two-dimensional case of Sec. 5 is used.

6.1 Second-order central difference

The second-order central difference scheme in three dimensions is given by

$\nabla^{2}\Psi_{i,j,k}$	$=\left.\dfrac{\partial^{2}\Psi}{\partial x^{2}}\right\|_{i,j,k}+\left.\dfrac{\partial^{2}\Psi}{\partial y^{2}}\right\|_{i,j,k}+\left.\dfrac{\partial^{2}\Psi}{\partial z^{2}}\right\|_{i,j,k}$
	$\approx\dfrac{1}{h^{2}}$ $\left(\;\begin{tabular}[]{\|c\|c\|c\|}\hline\cr&&\\ \hline\cr&$1$&\\ \hline\cr&&\\ \hline\cr\end{tabular}\,\Psi_{i,j+1,k}+\begin{tabular}[]{\|c\|c\|c\|}\hline\cr&$1$&\\ \hline\cr$1$&$-6$&$1$\\ \hline\cr&$1$&\\ \hline\cr\end{tabular}\,\Psi_{i,j,k}+\begin{tabular}[]{\|c\|c\|c\|}\hline\cr&&\\ \hline\cr&$1$&\\ \hline\cr&&\\ \hline\cr\end{tabular}\,\Psi_{i,j-1,k}\;\right),$

(45)

and the structure of the corresponding $A$ and $A^{\prime}$ matrix are given in Fig. 4.

The unique forms of the Gershgorin disk centers ( $a_{ii}$ ) and radii ( $r_{i}$ ) for $A^{\prime}$ are shown in Table 6.

Table 6: Gershgorin disk centers (

a_{ii}

) and radii (

r_{i}

) for the

A^{\prime}

matrix of the three-dimensional central difference scheme.

$a_{ii}$	$r_{i}=\sum_{i\neq j}\|a_{ij}\|$
$L_{i}-6$	$3$
$L_{i}-6$	$4$
$L_{i}-6$	$5$
$L_{i}-6$	$6$

The stability bounds of Eq. (18) in this case has $\vec{G}$ defined as

\vec{G}=\{12,11,10,9,3,2,1,0\}.

(46)

In the linear case ( $s=0$ ) with no external potential and periodic, Dirichlet or Laplacian-zero boundary conditions, the linear stability bound becomes

{\bf k}<\frac{\sqrt{8}}{12}\,\frac{h^{2}}{a}.

(47)

6.2 Fourth-order central difference

The fourth-order central difference scheme in three dimensions is given by

	$\displaystyle\nabla^{2}\Psi$	$\displaystyle\approx\frac{1}{12\,h^{2}}\left[\Psi_{i+2,j,k}+\Psi_{i-2,j,k}+\Psi_{i,j+2,k}+\Psi_{i,j-2,k}+\Psi_{i,j,k+2}+\Psi_{i,j,k-2}\right.$		(48)
		$\displaystyle\left.-16\,(\Psi_{i+1,j,k}+\Psi_{i-1,j,k}+\Psi_{i,j+1,k}+\Psi_{i,j-1,k}+\Psi_{i,j,k+1}+\Psi_{i,j,k-1})+90\,\Psi_{i,j,k}\right].$

The single-storage version of the 2SHOC equivalent scheme is defined as [6]

		$\displaystyle 1)\;D_{i,j,k}=$		(49)
				(60)

		$\displaystyle 2)\;\nabla^{2}\Psi_{i,j,k}\approx$		(61)
				(72)
				(83)

and the structure of the corresponding $A$ and $A^{\prime}$ matrix are given in Fig. 5, while the unique forms of the Gershgorin disk centers and radii for $A^{\prime}$ are shown in Table 7.

Table 7: Gershgorin disk centers and radii for the

A^{\prime}

matrix of the three-dimensional 2SHOC scheme.

$a_{ii}$	$r_{i}=\sum_{i\neq j}\|a_{ij}\|$
$L_{i}-15/2$	$33/4$
$L_{i}-15/2$	$25/3$
$L_{i}-15/2$	$101/12$
$L_{i}-15/2$	$17/2$
$L_{i}-22/3$	$67/12$
$L_{i}-22/3$	$17/3$
$L_{i}-29/4$	$17/4$
$L_{i}-89/12$	$83/12$
$L_{i}-89/12$	$7$
$L_{i}-89/12$	$85/12$

The stability bounds are the same as in Eq. (18), but with $\vec{G}$ now defined as

\vec{G}=\frac{1}{12}\times\left\{192,191,190,189,174,173,172,156,155,138,36,21,20,6,5,4,-9,-10,-11,-12\right\}.

(84)

In the linear case ( $s=0$ ) with no external potential and periodic, Dirichlet or Laplacian-zero boundary conditions, we get the linear stability bound

{\bf k}<\frac{\sqrt{8}}{16}\,\frac{h^{2}}{a},

(85)

which, as in the one- and two-dimensional cases, is simply $3/4$ of the second-order bound given in Eq. (47).

7 Numerical examples

Here we show some numerical examples to demonstrate the accuracy of the predicted stability bounds. Since the accuracy of the bounds are highly dependent on the problem including the values of $s$ , $a$ , and $V({\bf r})$ , the tests given here are not exhaustive, but serve as a good indication of the general accuracy of the bounds.

We choose to use three initial conditions, one for each dimensionality case, and integrate them using both the CD and 2SHOC schemes. In one dimension, we use the exact steady-state bright soliton solution to the NLSE with $V(x)=0$ and $s>0$ given as [14]

\Psi(x,t)=\sqrt{\frac{2\,\Omega}{s}}\,\mbox{sech}\left[\sqrt{\frac{\Omega}{a}}\,x\right]\,\mbox{exp}\left(i\,\Omega\,t\right),

(86)

where $\Omega$ is the frequency, and we set $V(x)=0$ , $\Omega=1$ , $s=1$ , and $a=1$ and use Dirichlet boundary conditions ( $\Psi=0$ ). In two dimensions, we use an approximation to a co-rotating dark vortex pair solution to the NLSE. Each vortex is given by [7]

\Psi(r,\theta,t)=f(r)\,\mbox{exp}[i\,(m\,\theta+\Omega\,t)],

(87)

where $m$ is the topological charge of the vortex (in our case, we use $m=1$ ), $\Omega=-1$ , and we use $s=1$ and $a=1$ . The term $f(r)$ is the real-valued radial profile centered at the vortex core which can be found numerically [4]. The pair of vortices are then combined to yield the initial condition

\Psi(x,y,t=0)=f(r_{1})f(r_{2})\,\mbox{exp}[i\,m\,(\theta_{1}+\theta_{2})],

where $r_{1}=\sqrt{(x-x_{0})^{2}+y^{2}}$ , $r_{2}=\sqrt{(x+x_{0})^{2}+y^{2}}$ , $\theta_{1}=\tan^{-1}(y/(x-x_{0}))$ , and $\theta_{2}=\tan^{-1}(y/(x+x_{0}))$ , which approximates the true initial condition of a co-rotating steady-state vortex pair solution located at $(-x_{0},0)$ and $(x_{0},0)$ . Here we choose $x_{0}=4$ . Since $|\Psi|^{2}$ does not decay at the boundaries, we use the modulus-squared Dirichlet boundary condition $|\Psi|^{2}=1$ . In three dimensions, we use a steady-state bright Gaussian solution of the LSE in a potential trap with an added initial ‘kick’ in the $x$ -direction which causes the structure to oscillate in the $x$ -direction. The initial condition is given by

\Psi(x,y,z,t=0)=\mbox{exp}\left(-\frac{x^{2}+y^{2}+z^{2}}{2\,a}\right)\mbox{exp}\left(-i\,\frac{x}{2}\right),

(88)

with external potential

V(x,y,z)=\frac{x^{2}+y^{2}+z^{2}}{a},

(89)

and we use the Dirichlet boundary condition $\Psi=0$ . For all simulations, we set the spatial step-size of the grid to $h=1/5$ . The three initial conditions are shown in the left column of Fig. 6.

To test the stability bounds, each solution is integrated to an ending time of $t=100$ and it is observed if the solution remains stable. We increase the time-step $k$ until the solution becomes unstable within the $t=100$ simulation time, at which point the largest time-step that was stable is denoted $k_{\mbox{\scriptsize num}}$ . This is then compared to the computed linear [Eqs. (90) and (91) denoted $k_{\mbox{\scriptsize lin}}$ ] and fully linearized [Eq. (92) denoted $k_{\mbox{\scriptsize linz}}$ ] stability bounds formulated in Secs. 4–6. The time-step is incremented to yield the numerical stability limit to within four significant figures. All of the simulations are performed using the NLSEmagic software package [3].

Example:	$k_{\mbox{\scriptsize lin}}$	$k_{\mbox{\scriptsize linz}}$	$k_{\mbox{\scriptsize num}}$	%-diff $k_{\mbox{\scriptsize lin}}$	%-diff $k_{\mbox{\scriptsize linz}}$
1D CD	0.02828	0.02828	0.02832	0.14	0.14
1D 2SHOC	0.02121	0.02121	0.02124	0.14	0.14
2D CD	0.01414	0.01407	0.01402	-0.85	-0.36
2D 2SHOC	0.01061	0.01057	0.01054	-0.66	-0.28
3D CD	0.009428	0.008650	0.009213	-2.28	6.51
3D 2SHOC	0.007071	0.006624	0.006992	-1.12	5.56

Table 8: Numerical test results of finding the numerical stability bound (

k_{\mbox{\scriptsize num}}

) for the example problems described in Sec. 7 compared to the predicted linear (

k_{\mbox{\scriptsize lin}}

) and linearized (

k_{\mbox{\scriptsize linz}}

) bounds.

Before displaying the results, we point out that there are some sources of error to consider. First, the predicted stability bounds are linearized and therefore will not be the same as the corresponding true nonlinear stability bounds. Second, in our analysis, we chose to use every possible combination of $\vec{L}$ which may lead to predictions of the bounds which are stricter than the true bound. Finally, it is sometimes difficult to determine the true stability bound numerically, as some unstable time-steps may only exhibit their instability after a very long simulation time. For our test, we choose a moderately long simulation time, but the exact bound may be slightly higher than the given result.

The results are shown in Table 8, while Fig. 6 shows the solutions before and after the recorded numerical stability bounds for three chosen examples. We see that overall, the numerical results match the predicted stability values quite well (especially in one and two dimensions) with a maximum percent difference of $6.5\%$ when $V({\bf r})\not=0$ in the three-dimensional example, but with a typical percent difference less that $1\%$ when $V({\bf r})=0$ . It is noted that in some cases the predicted bounds are stricter than the numerical result, while in other cases, they are too lenient, noting that the examples with $s>0$ were all too strict, while those with $s<0$ were all too lenient. However, due to the small number of tests, no conclusions about the effect of the sign and presence of the parameters and external potential of the LSE and NLSE on the stability bound predictions can be drawn from these observations. In terms of choosing a stable time-step for LSE and NLSE simulations, the results given are well within a tolerable range, and in practice one would use a time-step some percentage (say $10$ – $20\%$ ) lower than the predicted bound to ensure stability over long integration times.

8 Conclusion and summary of results

In this paper we have formulated linearized stability bounds for using second- and fourth-order spatial finite-differencing with fourth-order Runge-Kutta time-stepping for the multi-dimensional nonlinear Schrödinger equation (NLSE) with Dirichlet, modulus-squared Dirichlet, Laplacian-zero, and periodic boundary conditions.

A summary of the stability results for easy reference is given presently. For the nonlinear Schrödinger equation defined as

i\,\frac{\partial\Psi}{\partial t}+a\,\nabla^{2}\Psi-V({\bf r})\Psi+s\,|\Psi|^{2}\Psi=0,

where $a>0$ and $s$ are parameters of the system and $V({\bf r})$ is an external potential, the numerical stability bounds on the time-step when using the fourth-order Runge-Kutta time-stepping scheme is as follows:

In the linear case where $s=0$ and with no external potential ( $V({\bf r})=0$ ), utilizing periodic, Dirichlet, or Laplacian-zero boundary conditions, the stability bound on the time-step ${\bf k}$ when using the second-order central difference (CD) scheme in a $d$ -dimensional setting is

{\bf k_{\mbox{\scriptsize CD}}}<\frac{h^{2}}{d\,\sqrt{2}\,a},

(90)

while that of using a fourth-order central difference scheme (with interior points computed in the two-step high-order compact (2SHOC) methodology of Ref. [6]) is

{\bf k_{\mbox{\scriptsize 2SHOC}}}<\left(\frac{3}{4}\right)\frac{h^{2}}{d\,\sqrt{2}\,a}.

(91)

The linearized stability bounds for the general NLSE are

{\bf k}<\frac{\sqrt{8}}{\max\{\lVert\vec{B}\rVert_{\infty},\|\forall L_{i},L_{i}-\vec{G}\rVert_{\infty}\}}\,\frac{h^{2}}{a},

(92)

where $\vec{B}$ are the boundary points as defined by Table 9 (or in the periodic case is ignored), the elements of $\vec{L}$ is defined as

L_{i}=\frac{h^{2}}{a}\left(s|\Psi_{i}|^{2}-V_{i}\right),

where the index $i$ spans the entire grid, and $\vec{G}$ is a set of values defined in Table 10, determined by the dimension and method being used.

Table 9: Values of

\vec{B}

in Eq. (92).

	Dirichlet ( $\Psi_{b}=\mbox{const}$ )	Laplacian-zero ( $\nabla^{2}\Psi_{b}=0$ )	MSD ( $\|\Psi_{b}\|^{2}=\mbox{const}$ )
$B_{b}$	$0$	$\dfrac{h^{2}}{a}\left(s\|\Psi_{b}\|^{2}-V_{b}\right)$	$\dfrac{h^{2}}{a}\mbox{Im}\left[\dfrac{\Psi_{t,b-1}}{\Psi_{b-1}}\right]$

Table 10: Values of

\vec{G}

in Eq. (92).

Scheme $\rightarrow$	CD $O(h^{2})$	2SHOC $O(h^{4})$
1D	$\{4,3,1,0\}$	$\dfrac{1}{12}\times\left\{64,63,46,\right.$
		$\left.\qquad 12,-3,-4\right\}$
2D	$\{8,7,6,2,1,0\}$	$\dfrac{1}{12}\times\left\{128,127,126,110,109,\right.$
		$\left.\qquad 92,24,9,8,-6,-7,-8\right\}$
3D	$\{12,11,10,9,3,2,1,0\}$	$\dfrac{1}{12}\times\left\{192,191,190,189,174,\right.$
		$\qquad 173,172,156,155,138,36,21,$
		$\qquad\left.20,6,5,4,-9,-10,-11,-12\right\}$

We have found through numerical testing (those of Sec. 7, as well as others not reported here) that to ensure stability in all dimensions for typical problems, the bounds must be lowered by about $10\%$ – $20\%$ (most likely due to nonlinear effects). Also, we note that the reduced linear results are often similar to the full linearized bounds and can therefore be used as a good quick estimate of the stability bound.

Acknowledgments

This research was supported by NSF-DMS-0806762 and the Computational Science Research Center (CSRC) at SDSU. We gratefully acknowledge insightful discussions with Peter Blomgren.

References

[1] R. Bronson, Schaum’s Outline of Matrix Operations, 1st ed., McGraw-Hill, 1988.
[2] J. Butcher, A history of Runge-Kutta methods, Appl. Num. Math. 20 (3) (1996) 247–260.
[3] R. M. Caplan, NLSEmagic: Nonlinear Schrödinger equation multi-dimensional matlab-based GPU-accelerated integrators using compact high-order schemes, Comput. Phys. Commun. (2012) 10.1016/j.cpc.2012.12.010.
[4] R. M. Caplan, Study of vortex ring dynamics in the nonlinear Schrödinger equation utilizing GPU-accelerated high-order compact numerical integrators, Ph.D. thesis, Claremont Graduate University and San Diego State University (2012).
[5] R. M. Caplan, R. Carretero-González, A modulus-squared Dirichlet boundary condition for time-dependant complex partial differential equations and its application to the nonlinear Schrödinger equation, arXiv:1110.0569.
[6] R. M. Caplan, R. Carretero-González, A two-step high order compact scheme for the Laplacian operator and its implementation in a fully explicit method for integrating the nonlinear Schrödinger equation, arXiv:1109.1027.
[7] R. Carretero-González, D. Frantzeskakis, P. Kevrekidis, Nonlinear waves in Bose-Einstein condensates: physical relevance and mathematical techniques, Nonlinearity 21 (2008) R139–R202.
[8] M. M. Cerimele, M. L. Chiofalo, F. Pistella, S. Succi, M. P. Tosi, Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates, Phys. Rev. E. 62 (1) (2000) 1382–1389.
[9] W. Dai, An unconditionally stable three-level explicit difference scheme for the Schrödinger equation with a variable coefficient, SIAM J. Num. Anal. 29 (1) (1992) 174–181.
[10] L. Debnath, Nonlinear Partial Differential Equations for Scientists and Engineers, 2nd ed., Birkhauser Boston, 2005.
[11] G. H. Golub, Matrix Computations, 3rd ed., Johns Hopkins University Press, 1996.
[12] M. S. Ismail, A fourth-order explicit schemes for the coupled nonlinear Schrödinger equation, Appli. Math. Comp. 196 (2008) 273–284.
[13] P. Kevrekidis, D. Frantzeskakis, R. Carretero-González, Emergent Nonlinear Phenomena in Bose-Einstein Condensates: Theory and Experiment, vol. 45, Springer Series on Atomic, Optical, and Plasma Physics, 2008.
[14] Y. S. Kivshar, B. Luther-Davies, Dark optical solitons: physics and applications, Phys. Rep. 298 (2-3) (1998) 81–197.
[15] I. Kra, S. R. Simanca, On circulant matrices, Not. AMS 59 (3) (2012) 368–377.
[16] J. D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, 1991.
[17] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd ed., SIAM, 2004.

	$\displaystyle\left\|R(\lambda)\right\|^{2}=1$	$\displaystyle+\frac{1}{576}\,{\bf k}^{8}\|\lambda\|^{8}-\frac{1}{72}\,{\bf k}^{6}\|\lambda\|^{6}+\left(\frac{{\bf k}^{6}\|\lambda\|^{6}}{6}-{\bf k}^{4}\|\lambda\|^{4}+24\right)\,\frac{{\bf k}}{12}\,\mbox{Re}(\lambda)$		(6)
		$\displaystyle+\left({\bf k}^{4}\|\lambda\|^{4}+24\right)\,\frac{{\bf k}^{2}}{12}\,(\mbox{Re}(\lambda))^{2}+\left({\bf k}^{2}\|\lambda\|^{2}+4\right)\,\frac{{\bf k}^{3}}{3}\,(\mbox{Re}(\lambda))^{3}+\frac{2{\bf k}^{4}}{3}\,(\mbox{Re}(\lambda))^{4}.$