Long-Time Simulations of Nonlinear Schrödinger-Type Equations using Step Size Exceeding Threshold of Numerical Instability

Abstract

We propose an exponential time-differencing method based on the leapfrog scheme for numerical integration of the generalized nonlinear Schrödinger-type equations. The key advantage of the proposed method over the widely used Fourier split-step method is that in the new method, numerical instability at high wavenumbers is strongly suppressed. This allows one to use time steps that considerably exceed the instability threshold, which leads to a proportional reduction of the computational time. Moreover, we introduce a technique that eliminates numerical instability at low-to-moderate wavenumbers that is common for methods based on the leapfrog scheme. We illustrate the performance of the proposed method with examples from two applications areas: deterministic wave turbulence and solitary waves.

Appendices

Appendix A: Comparison of High-k NI of IF-LF and SSM

We will first illustrate that the high-k NI of these two methods is the same as long as the spectrum of the background solution has negligible content in the vicinity of the resonant wavenumbers $k_{\pi m}$, defined in (18). Then we will present a heuristic argument why the NI of the IF-LF appears to be stronger than that of the SSM when the solution’s spectrum expands to and beyond $k_{\pi m}$.

Figure 17 shows a typical spectrum of a soliton initial condition,

$$\begin{aligned} u(x,0)=\text{ sech }(x) + \xi (x), \end{aligned}$$

(47)

whose evolution is computed by the IF-LF (Fig. 17a) or ETD-LF (Fig. 17b) with a time step exceeding the NI threshold. The spectrum computed by the SSM is indistinguishable from that computed by the IF-LF. The white noise $\xi (x)$ of magnitude of order $10^{-10}$ is added in (47) to facilitate controllable observation of NI; in all simulations $\xi (x)$ was the same. The two spikes in the vicinity of $k_{\pi }$ in Fig. 17a are the unstable modes. In Table 2 we list the locations and magnitudes of these spikes for the SSM as one varies $\Delta t$. The corresponding locations for the IF-LF were found to be exactly the same as for the SSM, while the magnitudes were within $\pm 5$%. This demonstrates the close similarity between the high-k NIs of the two methods.

Fig. 17

Spectra of the numerical solution of (1) with $\gamma =2$, $V(x)\equiv 0$, and the near-soliton initial condition (47). The domain length $L=128\pi$ and the number of grid points $N=2^{12}$ are as in the experiment reported in [23]. Panel a shows a typical spectrum obtained at $t=1000$ with the IF-LF with $\Delta t_\mathrm{thresh}< \Delta t< 2\Delta t_\mathrm{thresh}$. (For some $\Delta t$, there may be more than two unstable modes around $k_{\pi }$; see Table 2.) Panel b shows the spectrum obtained by the ETD-LF at $t=10,000$ (i.e. 10 times longer than in a) with $\Delta t=0.05$. Stabilization, discussed in Sect. 4, did not need to be applied for either the IF- or ETD-LF

Table 2 Locations and magnitudes of the unstable modes for the SSM at $t=1000$; other parameters are listed in the caption to Fig. 17. For reference, $\Delta t_\mathrm{thresh}\approx 0.0031$. Notations $k_1$, etc. are explained in Fig. 17a. For the last two values of $\Delta t$ there are four unstable modes near $k_{\pi }$ [23]; the locations and magnitude for the secondary pair of modes are listed below those for the first pair

Figure 17b shows that the NI remains suppressed in the solution obtained by the ETD-LF over a very long time even when $\Delta t$ exceeds the threshold more than 15-fold. Note that unlike in Fig. 17a, here the spikes do not correspond to numerically unstable modes. These spikes occur in the first few time units of the evolution (probably due to a discretization error) and then remain almost unchanged over a very long time.

We will now present a heuristic explanation of why the unstable modes in the IF-LF method appear to grow much faster than in the SSM (see Fig. 1a). In Fig. 18 we provide evidence that this “extra” growth is not exponential. Rather, it occurs only when the spectrum expands so much that near $\pm k_{\pi m}$ it considerably exceeds the noise floor. Figure 18 shows that the magnitude of unstable modes near $k_{\pi }$ has increased by almost two orders of magnitude precisely at those short time intervals of duration $\sim 0.1$ when the solution’s spectrum has expanded beyond $k_{\pi }$. After the spectrum has shrunk (Fig. 18c), these modes retain their abnormally increased amplitude. They proceed on growing slowly¹² until the next instance of spectrum expansion beyond $k_{\pi }$ occurs.

Fig. 18

Spectra (on the $\log _{10}$-scale) of the solution of the pure NLS obtained by the IF-LF with $\Delta t=3.2\cdot 10^{-4} > \Delta t_\mathrm{thresh}$. Parameters are similar to those in Sect. 5.1, but $L=80\pi$, $N=2^{13}$. Shown are successive stages of the spectrum expansion beyond $k_{\pi }$ (a, b) and subsequent shrinking (c)

In contrast, the SSM does not produce such quickly growing modes; numerically unstable modes in its solution grow at an approximately constant rate and independently of the behavior of the spectrum. This indicates that a possible reason for the quickly growing modes in the IF-LF is the presence of two solutions of its three-level difference equation (8).¹³ Indeed, the SSM, being a two-level method, has only one solution of its difference scheme. An indirect confirmation of our conjecture comes from the solution by the ETD-LF, which, like the IF-LF, is a three-level method but, unlike the IF-LF, has no or a greatly reduced high-k NI. Figure 1 shows that the ETD-LF solution also develops spectral peaks around $|k|=k_{\pi m}$, and these peaks grow only when the spectrum expands beyond their locations. However, we have been unable to quantitatively explain the formation of these peaks around $|k|=k_{\pi m}$ during spectrum expansions.

Appendix B: Reduction of Spectral Distortions Near $|k|=k_{\pi }$ and $k_{2\pi }$

We have observed that the distortions in questions appeared only near $k_{\pi m}$ but not, for example, near $k_{\pi /2}$ (where $k_{\pi /2}^2\Delta t=\pi /2$). Therefore, one expects that distortions for $|k|\approx k_{\pi }$ will be reduced if one somehow could use a twice as small time step, $\Delta t_\mathrm{new}=\Delta t/2$, to resolve the dynamics of those specific Fourier modes (because then $k_{\pi }^2\Delta t_\mathrm{new}=\pi /2$). This can be done by using, instead of (8), a formula

$$\begin{aligned} e^{ik^2\Delta t/2 }\,\widehat{u}(t_{n+1}) - e^{-ik^2\Delta t/2}\,\widehat{u}(t_{n}) = i\Delta t\,\mathrm{sinc}\,(k^2\Delta t/2)\,\widehat{\mathcal N}(t_{n+1/2})\,, \end{aligned}$$

(48a)

where the term on the r.h.s. is found by the standard linear extrapolation:

$$\begin{aligned} \widehat{\mathcal N}(t_{n+1/2}) \equiv \frac{3}{2} \widehat{\mathcal N}(t_{n}) - \frac{1}{2} \widehat{\mathcal N}(t_{n-1}). \end{aligned}$$

(48b)

For the result shown in Fig. 1(b) with the green solid line, this technique was applied to harmonics with $|k|\in [k_{\pi }-12,\,k_{\pi }+10]$ (for reference, $k_{\pi }=41.8$ there).

We have empirically found that this idea does not work if extended to $|k|\approx k_{2\pi }$, where one would use $\Delta t_\mathrm{new}=\Delta t/4$. The problem appeared to be in a linear NI that sets in when the term $\widehat{\mathcal N}(t_{n+3/4})$, required for an extension of (48a), was computed by an extrapolation analogous to (48b).

Therefore we used a different idea for $|k|\approx k_{2\pi }$. Those harmonics do not need to be resolved accurately, because they contribute little to the solution. Rather, one simply needs to arrest their abnormal growth (see the spectrum with “no corrections” in Fig. 1b). A well-known way to gain stability at the expense of accuracy is to use a linearly implicit method:

$$\begin{aligned} \widehat{u}(t_{n+1}) - \widehat{u}(t_{n}) = (-ik^2\Delta t/2)\big ( \widehat{u}(t_{n+1}) + \widehat{u}(t_{n}) \big ) + i\Delta t\,\widehat{\mathcal N}(t_{n+1/2})\,, \end{aligned}$$

(49)

where $\widehat{\mathcal N}(t_{n+1/2})$ is computed by (48b). We have empirically found that the best results are obtained when (49) is applied every third step. The solid green line in Fig. 1b shows the result of doing so for harmonics with $|k|\in (k_{\pi }+10,\,k_{2\pi }+9]$.