Multiple scales analysis of a nonlinear timestepping instability in simulations of solitons

TL;DR

The paper identifies a nonlinear timestepping instability in pseudospectral simulations of solitons that cannot be explained by linear stability analysis or CFL arguments. A novel multi-scale asymptotic framework is developed to capture the slow, nonlinear accumulation of discretization error across multiple time scales, yielding predictions for blow-up or decay of solitons depending on the IMEX scheme. The framework accurately predicts slow-time evolution of the soliton parameter and corresponding blow-up or decay times, aligning with numerical experiments and highlighting the limitations of linear analyses for long-time nonlinear wave simulations. The results guide the selection of timestepping schemes for robust long-time soliton simulations and suggest extensions to fully implicit methods and dissipative systems.

Abstract

The susceptibility of timestepping algorithms to numerical instabilities is an important consideration when simulating partial differential equations (PDEs). Here we identify and analyze a pernicious numerical instability arising in pseudospectral simulations of nonlinear wave propagation resulting in finite-time blow-up. The blow-up time scale is independent of the spatial resolution and spectral basis but sensitive to the timestepping scheme and the timestep size. The instability appears in multi-step and multi-stage implicit-explicit (IMEX) timestepping schemes of different orders of accuracy and has been found to manifest in simulations of soliton solutions of the Korteweg-de Vries (KdV) equation and traveling wave solutions of a nonlinear generalized Klein-Gordon equation. Focusing on the case of KdV solitons, we show that modal predictions from linear stability theory are unable to explain the instability because the spurious growth from linear dispersion is small and nonlinear sources of error growth converge too slowly in the limit of small timestep size. We then develop a novel multi-scale asymptotic framework that captures the slow, nonlinear accumulation of timestepping errors. The framework allows the solution to vary with respect to multiple time scales related to the timestep size and thus recovers the instability as a function of a slow time scale dictated by the order of accuracy of the timestepping scheme. We show that this approach correctly describes our simulations of solitons by making accurate predictions of the blow-up time scale and transient features of the instability. Our work demonstrates that studies of long-time simulations of nonlinear waves should exercise caution when validating their timestepping schemes.

Figures (13)

• Figure 1: Evolution of the nonlinear timestepping instability. The numerical solution experiences unstable growth in the $L_2$ norm while maintaining the general shape of a soliton. The instability consists of the soliton increasing in amplitude, shrinking in width, and gaining speed. The left panel depicts three snapshots of the soliton when $\alpha = 0.00697$ in a simulation with $\Delta t = 0.00493$ using the SBDF2 timestepping scheme and a Fourier basis with $N = 512$. Snapshots are recorded at $T_0 = 0$, $T_{300} = 5657$, and $T_{600} = 10160$, which denote the times at which the soliton has traversed the domain 0, 300, and 600 times, respectively. The right panel tracks the evolution of the $L_2$ norm of the numerical solution for three different spatial discretizations---including Fourier and Chebyshev bases with different numbers of modes---and shows that they are indistinguishable.
• Figure 2: Error in simulations (Sim) of solitons compared to predictions from von Neumann (VN) analysis. We find linear instabilities in SBDF3 and SBDF4 in good agreement with the VN predictions. By contrast, SBDF1, SBDF2, and RK222 exhibit fast non-modal growth at early times as the numerical solution quickly becomes out of phase with the exact solution. RK443 exhibits similar behavior but, unlike the preceding schemes, is stable and does not result in blow-up. Shown are the $L_2$ norms of the error in pseudospectral simulations evolved using a timestep size of $\Delta t = 0.00324$ and dispersion parameter $\alpha = 0.00697$, alongside the norms of the fastest growing modes predicted by the von Neumann analysis. Each panel corresponds to a different timestepping scheme. Errors are plotted from $t = 0$ to $t = T_{\rm end}$ where $T_{\rm end}$ is the time at which the simulation terminated. The cause of termination is either due to blow-up (SBDF1: $T_{\rm end} = 17.35$, SBDF2: $T_{\rm end} = 36770$, SBDF3: $T_{\rm end} = 4.510$, SBDF4: $T_{\rm end} = 0.4367$, and RK222: $T_{\rm end} = 56680$) or reaching a user-specified time-limit (RK443: $T_{\rm end} = 3(30030/77069) \Delta t^{-3} \alpha^2 c_0^{-6} = 107400$). The errors in SBDF2, RK222, and RK443 momentarily "dip" whenever the numerical solution is in phase with the exact solution; since the numerical solution gains speed over time, it repeatedly laps the exact solution.
• Figure 3: Linear instability error growth rates measured in simulations (Sim) of solitons with SBDF3 and SBDF4 compared to rates predicted by von Neumann (VN) analysis. Measured growth rates for both schemes exhibit good agreement with the rates predicted by taking the largest VN growth rates, with better agreement as $\Delta t \rightarrow 0$. Shown are (crosses) the measured exponential growth rates of the $L_2$ norm of the error in simulations with SBDF3 and SBDF4 and (circles) the real parts of the growth rates $\lambda$ inferred from the largest eigenvalues in the von Neumann eigenspectra. The former values were obtained from least squares fits of the $L_2$ data. Each point represents a measurement/computation for a particular ($\alpha, \Delta t$). The growth rate $\lambda$ is defined such that $L_2(u_{\rm sim} - u_{\rm ex}) \propto e^{{\rm Re}(\lambda) t}$. Both schemes exhibit ${\rm Re}(\lambda) \sim {\cal O}(\Delta t^{-1})$, indicating that reduction of the timestep size leads to faster unstable growth.
• Figure 4: Measured and predicted soliton blow-up times $T_{\rm blowup}$ in SBDF1, SBDF2, and RK222 simulations. The accuracy of the blow-up time prediction improves as $\epsilon \rightarrow 0$. Blow-up times $T_{\rm blowup}$ in SBDF1, SBDF2, and RK222 are plotted against $\epsilon \propto \alpha^{-1/2} \Delta t$ as defined in \ref{['eq:ms-epsilon']}. We plot $T_{\rm blowup}$ values measured in simulations (Sim) and those predicted by the multiple scales analysis (MS) assuming the solitons are evolved over a finite periodic domain.
• Figure 5: Relative errors between measured blow-up times $T_{\rm blowup}$ and predicted blow-up times $T_{\rm blowup, MS}$ in SBDF1, SBDF2, and RK222 simulations. The relative error in the blow-up time prediction converges in the $\epsilon \rightarrow 0$ limit. As $\epsilon \rightarrow 0$, the errors decay as $\epsilon^2$, except for the SBDF1 timestepping scheme. In that case, we find errors instead decrease like $\epsilon$ for $\epsilon \lesssim 10^{-2}$.