Solving periodic semilinear stiff PDEs in 1D, 2D and 3D with exponential integrators

Dozens of exponential integration formulas have been proposed for the high-accuracy solution of stiff PDEs such as the Allen-Cahn, Korteweg-de Vries and Ginzburg-Landau equations. We report the results of extensive comparisons in MATLAB and Chebfun of such formulas in 1D, 2D and 3D, focusing on fourth and higher order methods, and periodic semilinear stiff PDEs with constant coefficients. Our conclusion is that it is hard to do much better than one of the simplest of these formulas, the ETDRK4 scheme of Cox and Matthews

High-order numerical method for the nonlinear Helmholtz equation with material discontinuities in one space dimension

The nonlinear Helmholtz equation (NLH) models the propagation of electromagnetic waves in Kerr media, and describes a range of important phenomena in nonlinear optics and in other areas. In our previous work, we developed a fourth order method for its numerical solution that involved an iterative solver based on freezing the nonlinearity. The method enabled a direct simulation of nonlinear self-focusing in the nonparaxial regime, and a quantitative prediction of backscattering. However, our simulations showed that there is a threshold value for the magnitude of the nonlinearity, above which the iterations diverge. In this study, we numerically solve the one-dimensional NLH using a Newton-type nonlinear solver. Because the Kerr nonlinearity contains absolute values of the field, the NLH has to be recast as a system of two real equations in order to apply Newton's method. Our numerical simulations show that Newton's method converges rapidly and, in contradistinction with the iterations based on freezing the nonlinearity, enables computations for very high levels of nonlinearity. In addition, we introduce a novel compact finite-volume fourth order discretization for the NLH with material discontinuities.The one-dimensional results of the current paper create a foundation for the analysis of multi-dimensional problems in the future.Comment: 47 pages, 8 figure

Uniformly accurate time-splitting methods for the semiclassical Schrödinger equationPart 1 : Construction of the schemes and simulations

This article is devoted to the construction of new numerical methods for the semiclassical Schrödinger equation. A phase-amplitude reformulation of the equation is described where the Planck constant epsilon is not a singular parameter. This allows to build splitting schemes whose accuracy is spectral in space, of up to fourth order in time, and independent of epsilon before the caustics. The second-order method additionally preserves the L^2-norm of the solution just as the exact flow does. In this first part of the paper, we introduce the basic splitting scheme in the nonlinear case, reveal our strategy for constructing higher-order methods, and illustrate their properties with simulations. In the second part, we shall prove a uniform convergence result for the first-order splitting scheme applied to the linear Schrödinger equation with a potential

Numerical methods for generalized nonlinear Schrödinger equations

We present and analyze different splitting algorithms for numerical solution of the both classical and generalized nonlinear Schr"odinger equations describing propagation of wave packets with special emphasis on applications to nonlinear fiber-optics. The considered generalizations take into account the higher-order corrections of the linear differential dispersion operator as well as the saturation of nonlinearity and the self-steepening of the field envelope function. For stabilization of the pseudo-spectral splitting schemes for generalized Schr"odinger equations a regularization based on the approximation of the derivatives by the low number of Fourier modes is proposed. To illustrate the theoretically predicted performance of these schemes several numerical experiments have been done

Periodic and solitary waves in systems of coherently coupled nonlinear envelope equations

Exact solutions for two classes of coherently coupled nonlinear envelope equations are derived in terms of products of Jacobi elliptic functions. Physical applications are illustrated in the context of nonlinear optics, namely, polarization of light beams and quadratic (or parametric) solitons. Stabilities of these double-humped solitary pulses are studied by direct numerical simulations. The use of computer is crucial, both in terms of symbolic manipulation in the derivation process and in the implementation of numerical schemes in stability consideration. © 2010 Taylor & Francis.postprin

A time-splitting pseudospectral method for the solution of the Gross-Pitaevskii equations using spherical harmonics with generalised-Laguerre basis functions

We present a method for numerically solving a Gross-Pitaevskii system of equations with a harmonic and a toroidal external potential that governs the dynamics of one- and two-component Bose-Einstein condensates. The method we develop maintains spectral accuracy by employing Fourier or spherical harmonics in the angular coordinates combined with generalised-Laguerre basis functions in the radial direction. Using an error analysis, we show that the method presented leads to more accurate results than one based on a sine transform in the radial direction when combined with a time-splitting method for integrating the equations forward in time. In contrast to a number of previous studies, no assumptions of radial or cylindrical symmetry is assumed allowing the method to be applied to 2D and 3D time-dependent simulations. This is accomplished by developing an efficient algorithm that accurately performs the generalised-Laguerre transforms of rotating Bose-Einstein condensates for different orders of the Laguerre polynomials. Using this spatial discretisation together with a second order Strang time-splitting method, we illustrate the scheme on a number of 2D and 3D computations of the ground state of a non-rotating and rotating condensate. Comparisons between previously derived theoretical results for these ground state solutions and our numerical computations show excellent agreement for these benchmark problems. The method is further applied to simulate a number of time-dependent problems including the Kelvin-Helmholtz instability in a two-component rotating condensate and the motion of quantised vortices in a 3D condensate

Conservation of phase space properties using exponential integrators on the cubic Schrödinger equation

The cubic nonlinear Schrödinger (NLS) equation with periodic boundary conditions is solvable using Inverse Spectral Theory. The nonlinear spectrum of the associated Lax pair reveals topological properties of the NLS phase space that are difficult to assess by other means. In this paper we use the invariance of the nonlinear spectrum to examine the long time behavior of exponential and multisymplectic integrators as compared with the most commonly used split step approach. The initial condition used is a perturbation of the unstable plane wave solution, which is difficult to numerically resolve. Our findings indicate that the exponential integrators from the viewpoint of efficiency and speed have an edge over split step, while a lower order multisymplectic is not as accurate and too slow to compete. © 2006 Elsevier Inc. All rights reserved