Comparison of Split-Step and Hamiltonian Integration Methods for Simulation of the Nonlinear Schrödinger Equation
Anastassiya Semenova, Sergey A. Dyachenko, Alexander O. Korotkevich, Pavel M. Lushnikov
CComparison of Split-Step and Hamiltonian IntegrationMethods for Simulation of the Nonlinear Schr¨odingerEquation
Anastassiya Semenova a, ∗ , Sergey A. Dyachenko c , Alexander O. Korotkevich a,b ,Pavel M. Lushnikov a,b a Department of Mathematics & Statistics, The University of New Mexico, MSC01 1115, 1University of New Mexico, Albuquerque, New Mexico, 87131-0001, USA b L. D. Landau Institute for Theoretical Physics, 2 Kosygin Str., Moscow, 119334, RussianFederation c Department of Applied Mathematics, University of Washington, Lewis Hall 201, Box353925, Seattle, Washington 98195-3925, USA
Abstract
We provide a systematic comparison of two numerical methods to solve thewidely used nonlinear Schr¨odinger equation. The first one is the standard secondorder split-step (SS2) method based on operator splitting approach. The secondone is the Hamiltonian integration method (HIM), originally proposed in Ref.[1]. It allows the exact conservation of the Hamiltonian at the cost of requiringthe implicit time stepping. We found that numerical error for HIM methodis systematically smaller than the SS2 solution for the same time step. At thesame time, one can take orders of magnitude larger time steps in HIM comparedwith SS2 still ensuring numerical stability. In contrast, SS2 time step is limitedby the numerical stability threshold.
Keywords: nonlinear Schr¨odinger equation, numerical methods,pseudospectral methods, computational physics
1. Introduction
A nonlinear Schr¨odinger equation (NLSE) is one of the most generic non-linear partial differential equation in numerous branches of mathematical andtheoretical physics is [2]. NLSE naturally appears if one considers envelope dy-namics of a quasi–monochromatic nonlinear wave [3]. In quantum mechanicsa version of NLSE is called a Gross–Pitaevskii equation [4] which describes aBose-Einstein condensate for short-range interactions of particles. ∗ Corresponding author
Email addresses: [email protected] (Anastassiya Semenova), [email protected] (Sergey A. Dyachenko), [email protected] (Alexander O. Korotkevich), [email protected] (Pavel M. Lushnikov)
Preprint submitted to Journal of Computational Physics August 11, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] A ug typical NLSE application is the dynamics of optical pulses in an opticalfiber. The time evolution of the envelope of an optical pulse in a fiber is wellapproximated by NLSE, including the description of very long, transoceanic op-tical communication lines, see e.g. [5, 6]. The Langmuir waves in plasmas are de-scribed by NLSE as well, see e.g. Ref. [7, 8]. Dynamics of quasi-monochromaticoceanic waves (which is typical e.g. for ocean swell) is reduced to NLSE orits modifications [9]. For example, the analysis of NLSE offers a possible ex-planation to the mystery of appearance of the rogue waves [10]. All these andnumerous other applications of NLSE and its modification require efficient nu-merical simulation.Many techniques can be used in simulation of NLSE: the Crank-Nicholsonscheme, the hopscotch method, the Ablowitz–Ladik scheme, the pseudo–spectralsplit-step method, the Hamiltonian preserving method, and many others (seeRefs [11], [12], [13]). One of the most popular method of integration of NLSE wasproposed by F. Tappert [14], and its performance was studied in the article [11].The split-step method can be considered as a version of the Strang’s operatorsplitting approach [15] combined with pseudo-spectral method. The split-stepmethod can be constructed to any order of accuracy, in this work we consider thesecond order symmetrized split-step (SS2) method. A recent study of stabilityof the split-step method can be found in the work [16] and references therein.In 1992 a novel method for simulation of NLSE has been proposed in thepaper [1]. It has been successful to study turbulence in two–dimensional NLSE,however it passed largely unnoticed by a wide audience. Perhaps, that is thereason why it was not mentioned in the recent papers such as [17], [18], thatdescribe somewhat similar numerical methods. This numerical method, whichwe call the Hamiltonian integration method (HIM), conserves the numericalHamiltonian and the optical power exactly (in exact arithmetic), and it is basedon discrete Hamilton’s equations. In finite digit arithmetic, the error in con-servation of Hamiltonian is due to round-off errors inherent to floating pointrepresentation.By using the discrete Hamilton’s equations in other systems one may deriveHamiltonian–preserving numerical schemes. As an example, we refer the readerto the recent work [19] on numerical simulations of nonlinear water waves. Onecan trace similarities with the symplectic methods [20], whereas HIM is a com-pletely self–contained ad hoc method which can be derived for other Hamiltoniansystems having canonical symplectic structure.We compare the two numerical methods by performing a set of simulationswith various initial conditions. In these experiments we observe that in somescenarios HIM method can outperforms SS2 when very high accuracy is notessential. The SS2 method requires a stringent condition on time step for sta-bility, whereas HIM is an implicit method and as such allows the time step tobe a hundred times larger. Our observations illustrate that HIM method can bethe method of choice for efficient simulations of interaction of solitons, where atight balance between nonlinearity and dispersion occurs.The paper is organized as follows: we describe the mathematical problem insection 2; the description of numerical methods is given in section 3; the section 42iscusses the relation between the dimensionless NLSE and the physical unitsrelevant to optical fibers communications; the section 5 describes the set ofsimulations and discusses obtained results; and in section 6 we summarize ourobservations and discuss the applicability of both methods. The derivationof HIM method is placed in Appendix A and the convergence conditions arediscussed in Appendix B.
2. Problem Formulation
Let us consider NLSE in its simplest form (rescaling of coordinate, time, andamplitude can bring NLSE into this form without loss of generality): i Φ t + Φ xx + γ | Φ | Φ = 0 , (1)where Φ( x, t ) is a complex function, γ = ± x and t . The latter equation is solved on an interval x ∈ [ − L, L ] subject to periodicboundary conditions, and t ∈ [0 , T ]. For the sake of simplicity, we considerNLSE in one spatial dimension, although both methods are applicable to anydimensions (for example, HIM was originally formulated for 2D problem [1]). The Hamiltonian, H , and the number of particles, N , given by: H = (cid:90) (cid:16) | Φ x | − γ | Φ | (cid:17) dx and N = (cid:90) | Φ | dx, (2)are conserved quantities for (1). Here and further we integrate over one spatialperiod [ − L, L ] and drop the integration limits for brevity. The NLSE is anintegrable system [21], and it has infinitely many nontrivial integrals of motion,that may be used to track accuracy of numerical simulation. We consider firsttwo nontrivial integrals of motion, that are given by [21], [22]: C = (cid:90) (cid:20) Φ ¯Φ xxx + 3 γ x | Φ | (cid:21) dx, (3) C = (cid:90) (cid:20) | Φ xx | + γ | Φ | − γ (cid:0) | Φ | x (cid:1) − γ | Φ | | Φ x | (cid:21) dx. (4)We denote them C and C because the first three are so called trivial integralsof motion: the number of particles N and the Hamiltonian H (2), and themomentum. The equation (1), has soliton solutions [21], and when NLSE is consideredon infinite spatial interval, it may be solved by means of the inverse scattering3ransform (IST). Some solutions of NLSE that decay at x → ±∞ , such as N -soliton solutions may be used on a periodic interval when the magnitude of | Φ | is close enough to zero at the endpoints x = ± L .The one-soliton solution (here and further we use γ = 1) is given by theformula: Φ = √ λe i ( vx + ( λ − v ) t +Φ )cosh (cid:104) √ λ ( x − vt − x ) (cid:105) , (5)where x , and v are the constants that determine initial position and the prop-agation speed of the soliton, and the constants λ and Φ determine the solitonamplitude and the initial phase respectively.Another exact solution of (1) on infinite line is the two-soliton solution whichcan be obtained by dressing method [23], given by the formula:Φ = (cid:20) e η +¯ η ( p − p ) p + ¯ p ) ( p + ¯ p ) (cid:21) e η + (cid:20) e η +¯ η ( p − p ) p + p ) ( p + ¯ p ) (cid:21) e η D , (6)where D is the following expression: D = 1 + e η +¯ η p + ¯ p ) + e η +¯ η p + ¯ p ) + e η +¯ η p + ¯ p ) + e ¯ η + η p + p ) ++ e η +¯ η + η +¯ η | p − p | p + ¯ p ) ( p + ¯ p ) | p + ¯ p | (7)and η , η are determined by the expression: η , = p , x + ip , t + a , , (8)here p , and a , are complex constants. The width and the propagation speedof solitons are defined by the real and the imaginary parts of p , respectively.The initial positions of each soliton are defined by a , . It is natural to use Fourier series to approximate Φ( x, t ) on the periodicinterval x ∈ [ − L, L ] using a pseudo spectral approach by the means of the dis-crete Fourier transform (DFT) that is computed using the fast Fourier transformlibrary FFTW [24]. In physical space we use a uniform grid, x j = 2 LN j − L where j = 0 , . . . N − − L, L ]. We introduce a grid function, Φ nj = Φ( x j , n ∆ t ),where ∆ t is an elementary time step. 4 . Description of Numerical Methods In the SS2 method, the linear and nonlinear terms of (1) are treated sepa-rately in a style of Strang splitting [15].Let ˆ L = i∂ /∂x represent the operator for the linear term and ˆ N = iγ | Φ | represent the operator for the nonlinear term of (1), then Φ t ( x, t ) = ( ˆ L +ˆ N )Φ( x, t ). This equation has the formal solution Φ( x, t +∆ t ) = e (ˆ L + ˆ N )∆ t Φ( x, t )on a time step ∆ t . In the SS2 method [11] we approximate the exponential termby the product of separate exponents: e (ˆ L + ˆ N )∆ t = e ˆ L ∆ t e ˆ N ∆ t e ˆ L ∆ t + ∆ t { [ ˆ L, [ ˆ N , ˆ L ]] + 12 [ ˆ N , [ ˆ
N , ˆ L ]] } + . . . , (10)that is accurate up to third order in time. This is a special case of applicationof Campbell-Baker-Hausdorff formula [25]. By doing this, the evolution of thelinear part and nonlinear part on the step ∆ t can be carried out separately.In the context of NLSE this is particularly attractive because both evolutionscan be carried out analytically. Note that the linear PDE i Φ t = − Φ xx , can besolved exactly in the Fourier domain:Φ k ( t + ∆ t ) = e − ik ∆ t Φ k ( t ) , (11)where Φ k ( t ) denotes the Fourier coefficient, corresponding to wavenumber k , ofΦ( x, t ). The nonlinear part of (1) given by i Φ t = − γ | Φ | Φ is an ODE, and canbe solved exactly: Φ( x, t + ∆ t ) = e iγ | Φ | ∆ t Φ( x, t ) . (12)Equations (11) and (12) give us explicit expressions for e ˆ L and e ˆ N correspond-ingly. The only complexity is that these two exact solutions are given in Fourierand coordinate spaces which requires switching between them in order to rep-resent e ˆ L ∆ t e ˆ N ∆ t e ˆ L ∆ t in (10) consecutively.In a similar manner one may construct higher order split step methods, byalternating linear and nonlinear steps. The SS2 method is stable if the condition,∆ t ≤ ∆ x π (13)described in [26] is satisfied. However, one can violate this condition when thehighest Fourier coefficients are small enough.One can note that both steps (linear and nonlinear one) in SS2 methods areperforming only rotation of phase, so conservation of number of particle N isan intrinsic property of the method. 5 .2. Hamiltonian Integration Method The main feature of the HIM method (introduced in [1]) is its exact conser-vation of the Hamiltonian, H , and number of particles, N . This is achieved byrequiring that the difference in H (and N ) on subsequent time steps vanishes,the details of derivation of HIM are given in the Appendix A. HIM is an implicitscheme: i Φ n +1 j − Φ nj ∆ t = − (cid:2) Φ n +1 j + Φ nj (cid:3) xx − (Φ n +1 j + Φ nj )( | Φ n +1 j | + | Φ nj | )4 . (14)that is solved by means of fixed point iterations on every time step Equation (14)implicitly defines the solution at the subsequent time steps. In the Fourier spacethe formula (14) becomes the following:ˆΦ n +1 k − ˆΦ nk = − ik ∆ t n +1 k + ˆΦ nk ) + i ∆ t F (cid:2) (Φ n +1 + Φ n )( | Φ n +1 | + | Φ n | ) (cid:3) , (15)where ˆΦ nk = ˆ F [Φ n ] is the k -th Fourier coefficient of the grid function Φ nj . Fol-lowing the work [19], the linear part of the equation (15) can be resolved forˆΦ n +1 k which yields:ˆΦ n +1 k = 1 − i k ∆ t i k ∆ t ˆΦ nk + i ∆ t i k ∆ t ) ˆ F (cid:2) (Φ n +1 + Φ n )( | Φ n +1 | + | Φ n | ) (cid:3) . (16)The equation (16) can be solved by a fixed point iterations:ˆΦ n +1 ,s +1 k = 1 − i k ∆ t i k ∆ t ˆΦ nk + i ∆ t i k ∆ t ) ˆ F (cid:2) (Φ n +1 ,s + Φ n )( | Φ n +1 ,s | + | Φ n | ) (cid:3) , (17)where s denotes the iteration number and ˆΦ n +1 , k = ˆΦ nk . We iterate (17) untilthe residual condition is satisfied: (cid:13)(cid:13)(cid:13) ˆΦ n +1 ,s +1 k − ˆΦ n +1 ,sk (cid:13)(cid:13)(cid:13) = (cid:115)(cid:88) k (cid:12)(cid:12)(cid:12) ˆΦ n +1 ,s +1 k − ˆΦ n +1 ,sk (cid:12)(cid:12)(cid:12) ≤ ε , (18)where (cid:107)·(cid:107) denotes the l norm on [ − L, L ], and ε is the tolerance for fixed pointiterations. The initial values Φ n, are computed by using one step of ForwardEuler. Following [1], the fixed point iterations of HIM converge for∆ t <
23 max j ( | Φ nj | ) . (19)Derivation of this condition is given in Appendix B.For the time step that satisfies the above condition, the fixed point iterationstypically converge in 4 to 6 steps with the tolerance ε ≤ − .6 . Physical Units Relevant to Optical Fiber Before the investigation of the performance of the two methods on a longtime scale, we would like to estimate the characteristic time of simulation thatcorresponds to the dynamics of a pulse in a physically realistic fiber. In orderto do so we consider a trans–Atlantic fiber described in the reference paper [6]subject to: iA z − β A ττ + σ | A | A = 0 . (20)We use the values for β = −
20 ps km − , the group velocity dispersion (GVD),and σ = 1 . × − km − mW − , the strength of nonlinearity for a fiber, pro-vided therein.The dimensionless NLSE given by (1) must be rewritten in the originaldimensional units. We transform the dimensionless NLSE to dimensional unitsas follows: z = lt, τ = xω , and A = A Φ (21)The derivatives with respect to t and x are given by: ∂ t = l∂ z and ∂ x = 1 ω ∂ τ . (22)The resulting equation transforms into: iA z + 1 ω l A ττ + | A | AA l = 0 . (23)Comparison of two equations (20) and (23) reveals that: β (cid:20) ps km (cid:21) = − ω l , (24) σ (cid:20) km mW (cid:21) = 1 A l , (25)where A = 1 mW / . By using the parameters β and σ from the referencepaper [6], we find that l ≈
769 km, ω = 1 . × − ps − from the equations (24)–(25). We find that it is necessary to simulate the fiber until the dimensionlesstime t max ≈
13 in order to mimic a 10 km fiber. The nonlinear time is thengiven by t NL = π | Φ | = π | λ | which in physical units corresponds to z NL = t NL l .
5. Numerical Methods Performance
In this simulation we check the convergence rate of HIM and SS2 by runninga sequence of simulations with various time steps. As the initial condition weconsider a one-soliton solution (5) with the following parameters: λ = 2 , and Φ = x = v = 0 . (26)7 -8 -7 -6 -5 -4 -3 -2 -1 -4 -3 -2 M a x | E rr | ∆ tHIMSS2 -14 -12 -10 -8 -6 -4 -4 -3 -2 ∆ t ∆ N HIM ∆ N SS2 ∆ H HIM ∆ H SS2 ∆ C5 HIM ∆ C5 SS2
Figure 1: (Stationary one-soliton solution on a fully resolved grid) (Left) Convergence rate ofnumerical methods, HIM (green) and SS2 (red). Both methods have second order convergence,but L ∞ error in solution is about one order smaller for HIM compared to SS2 for the sametime steps. (Right) Error in conserved quantities: number of particles N (solid), Hamiltonian H (dotted), and C (dash-dotted) for various time steps. When time step is larger than thestability condition of SS2, errors in H and C start to grow. For HIM, the error is dominatedby accumulation of round-off errors and is smaller by several orders of magnitude comparedwith SS2. We run the simulation on a fully resolved (highest harmonics are of round-offlevel) uniform grid of N = 2048 grid points, and L = 25 π . The tolerancefor HIM iterations is set to ε = 10 − and simulation time is T = 5. Theconvergence of both methods is demonstrated in Figure 1. We omit the C inthe Figure 1 because this quantity is identically zero for a stationary one-solitonsolution. The error in the integrals of motion for SS2 method is dominatedby accumulation of round-off errors for small ∆ t , and by the order of methodfor large ∆ t as shown in the Figure 1. The critical value of ∆ t for which thetransition occurs is close to the stability condition of SS2 method. In these simulations we investigate how the traveling speed v of the one-soliton solution (5) affects the accuracy of both numerical methods. It is knownthat dispersion of waves by SS2 method is identical to the dispersion of NLSE,while from (16) it follows that the dispersion of HIM is only accurate up to thirdorder in k ∆ t . We expect that for sufficiently large time step the travel speed ofsoliton will deviate from its true value. We show the results of the simulationswith various travel speeds in Figure 2. The initial data for these simulations isgiven by (5) with parameters: λ = 2 , and Φ = x = 0 , and v ∈ [0 , . (27)The computational box size is L = 25 π and the number of grid points is N =2048. The tolerance for HIM iterations is ε = 10 − and the simulation time is T = 100. Time step for both methods is set to be ∆ t = . x π .It should be noted, that soliton velocity is given in dimensionless units.In the left panel of Figure 2, we observe that the error in the solution has no8 -4 -3 -2
0 1 2 3 4 5 M a x | E rr | v, speedHIMSS2 -12 -10 -8 -6
0 1 2 3 4 5v, speed ∆ N HIM ∆ N SS2 ∆ H HIM ∆ H SS2 -12 -10 -8 -6
0 1 2 3 4 5v, speed ∆ C HIM ∆ C SS2 ∆ C HIM ∆ C SS2
Figure 2: (Moving one-soliton solution on a fully resolved grid) (Left) The maximum absoluteerror of the solution at time T = 100 as a function of propagation speed of the soliton. TheSS2 method (red) has no dependence of the error on travel speed of the soliton because itnaturally captures the dispersion relation of NLSE, while HIM (green) has dispersion relationaccurate up to ∆ t . (Center) The error in integral quantities, N (solid), and H (dotted) isabout seven orders of magnitude smaller than the error in the solution. (Right) The errorin integral quantities, C (solid), and C (dotted) is about seven orders of magnitude smallerthan the error in the solution. For travel speed v ≤ C and C , but HIM behaves worse as soon as v is larger than 3. dependence on travel speed of the soliton for SS2 method. For HIM, the error inthe solution depends on travel speed which is due to inexact dispersion relationof HIM method: ω HIM ( k ) = i ∆ t ln 1 − i k ∆ t i k ∆ t = k (cid:18) − k ∆ t
12 + . . . (cid:19) , (28)where ω HIM ( k ) is the angular frequency of the k –th Fourier harmonic.In the center panel, we look at the absolute error in integral quantities, N and H . It is about seven orders of magnitude smaller than the error in the solution.On the right panel, we consider the absolute error in integral quantities, C and C . Similarly to N and H , it is about seven orders of magnitude smaller thanthe error in the solution. We notice that for travel speed v ≤ C and C , but the error in HIM becomes larger assoon as v gets larger than 3. Note that the error in the solution does not alwayscorrelate with the error in integral quantities. In this simulation we demonstrate the difference between SS2 and HIM whenthe initial data is a two-soliton solution with the following set of parameters: p = 2 . p = 1 . a = 60 + i = − ¯ a . (29)The simulation time is T = 5, the solution is underresolved on a grid with N = 1024 points. The computation box is x ∈ [ − L, L ] where L = 25 π . The time9 igure 3: (Stationary Two-Soliton Solution) Solution of NLSE with the initial data (29) withHIM method (left), and SS2 method (right). The SS2 method radiates waves continuouslyover the course of the simulation, while the HIM emits localized small amplitude perturbationsthat travel in the computational box and are reflected and transmitted through the stationarysolitons. At the time T = 5, the background radiation around the stationary solitons emittedin SS2 is several orders of magnitude larger than for HIM. step is ∆ t = . x π . It is typical to have solution not resolved to round–off errorin long and/or multichannel simulations of light pulses propagating in opticalfibers. A smaller number of Fourier harmonics implies faster computations. Forthis experiment, the smallest amplitudes were of the order 10 − . We presentthe results of the simulation in Figures 3 - 4.In the course of simulation we observe that the error in H and C is one totwo orders of magnitude smaller in HIM than in SS2. The number of particlesis better conserved by SS2 and the error is two orders of magnitude smaller. In this section, we study the dynamics of the two-soliton solution (6). Wepresent parameters of simulations in sections 5.4.1 - 5.4.3, and discuss resultsof simulations in the section 5.4.4. We use periodic box with L = 25 π and N = 4096 grid points for fully resolved simulations and N = 1024 for unre-solved simulations. The time step is ∆ t = . x π < ∆ x π to satisfy the sta-bility condition (13) in all three simulations. The HIM iterations tolerance is (cid:15) = 10 − . The initial condition is given by the two-soliton solution formula (6) whereone of the solitons is moving towards the other soliton which is at rest. Thesimulation time is T = 50, and over the course of simulation two solitons interactonce. We present the results of the simulation in the Figures 5 - 7.The parameters for this two–soliton solution are given by: p = 1 . p = 1 . ia = 2 . i and a = 65 + i. (30)10 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6
0 1 2 3 4 5t, time ∆ N HIM ∆ N SS2 ∆ H HIM ∆ H SS2 -12 -11 -10 -9 -8 -7 -6 -5
0 1 2 3 4 5t, time ∆ C HIM ∆ C SS2 ∆ C HIM ∆ C SS2
Figure 4: (Stationary two-soliton solution on underresolved grid) Conserved integrals in asimulation with initial data (29). (Left) The number of particles (solid) and the Hamiltonian(dotted) computed via SS2 (red) and HIM (green). (Right) The integrals C (solid) and C (dotted) via SS2 (red) and HIM (green).Figure 5: (Collision with stationary soliton) (Top) Numerical solution for HIM (left) and SS2(right) methods on a fully resolved grid N = 4096. (Bottom) Numerical solution for HIM(left) and SS2 (right) methods on an underresolved grid with N = 1024. -4 -3 -2 -1 M a x | E rr | ∆ t, timestepHIMSS2 -5 -5
0 10 20 30 40 50 M a x | E rr o r | t, timeHIMSS2 Figure 6: (Collision with stationary soliton on a fully resolved grid) (Left) Error in the solu-tion in L ∞ -norm as a function of time step in double-logarithmic scale shows second orderconvergence in ∆ t . (Right) Absolute error as a function of time, the solitons interact at ap-proximately t = 25. The error vs time is close to a straight line before and after the collision.Its slope, m , changes from m = 6 . × − to m = 7 . × − for HIM method, and from m = 8 . × − to m = 1 . × − for SS2. -14 -12 -10 -8 -6
0 10 20 30 40 50t, time ∆ N HIM ∆ N SS2 ∆ H HIM ∆ H SS2 -14 -12 -10 -8 -6
0 10 20 30 40 50t, time ∆ C HIM ∆ C SS2 ∆ C HIM ∆ C SS2
Figure 7: (Collision with stationary soliton on a fully resolved grid) The conserved quantitiesplotted as a function of time over the course of the simulation, note that SS2 demonstrates astrong peak in error in H at the time of solitons interaction. After the moment of interactionthe C , and the C exhibit jump and increase in error with in SS2 and HIM. igure 8: (Headon collision of solitons) (Top) Numerical solution for HIM (left) and SS2(right) methods on a fully resolved grid N = 4096. (Bottom) Numerical solution for HIM(left) and SS2 (right) methods on an underresolved grid with N = 1024. The initial condition is given by the two-soliton solution formula (6) withsolitons moving toward each other. The final time of simulation is T = 45, andtwo solitons interact once. We present the results of the simulation in Figures 8- 10. The parameters for this simulation are the following: p = 1 . − . i and p = 1 . ia = −
20 + i and a = 60 + i. (31) The initial condition is given by the two-soliton solution formula (6) withone soliton pursuing another soliton. The final time of simulation is T = 54.The pursuing soliton overtakes and interacts with the slower soliton once. Theresults of this simulation are presented in the Figures 11 - 13, and parametersof the initial condition are as follows: p = 1 . . i and p = 1 . ia = 50 + i and a = 110 + i. (32) In the latter sequence of three simulations involving two-soliton collision, wefound that the radiation level in SS2 simulation has been consistently higher13 -5 -5
0 10 20 30 40 M a x | E rr o r | t, timeHIMSS2 Figure 9: (Headon collision of solitons on a fully resolved grid) Error in L ∞ -norm of thesolution vs time computed for SS2 (red) and HIM (green) methods. The collision occurs atthe time approximately t = 21 where we observe a spike in the error. The error vs time isclose to a straight line before and after collision. Its slope, m changes from m = 8 . × − to m = 1 . × − for SS2 method, and from m = 6 . × − to m = 8 . × − for HIMmethod. -14 -12 -10 -8 -6 -4
0 5 10 15 20 25 30 35 40 45t, time ∆ N HIM ∆ N SS2 ∆ H HIM ∆ H SS2 -14 -12 -10 -8 -6 -4
0 5 10 15 20 25 30 35 40 45t, time ∆ C HIM ∆ C SS2 ∆ C HIM ∆ C SS2
Figure 10: (Headon collision of solitons on a fully resolved grid) The conserved quantities (left)∆ N (solid), ∆ H (dotted), and (right) ∆ C (solid), and ∆ C (dotted) as a function of time overthe course of the simulation with HIM (green) and SS2 (red). Note that SS2 demonstrates astrong peak in error in H at the time of soliton interaction. After the interaction time the C ,and the C exhibit large error with both SS2 and HIM. igure 11: (Collision with pursuing soliton) (Top) Numerical solution for HIM (left) and SS2(right) methods on a fully resolved grid with N = 4096 points. (Bottom) Numerical solutionfor HIM (left) and SS2 (right) methods on an underresolved grid with N = 1024 points. -4 -4
0 10 20 30 40 50 M a x | E rr o r | t, timeHIMSS2 Figure 12: (Collision with pursuing soliton on a fully resolved grid) Error in the solution vstime for SS2(red) and HIM(green) methods in the simulation with one soliton pursuing theother. The time of collision is approximately t = 28. We observe that the slope, m of thestraight line of error vs time changes at the collision for both methods. In SS2 it changes from m = 1 . × − to m = 1 . × − , and in HIM the slope changes from m = 6 . × − to m = 7 . × − . -14 -12 -10 -8 -6
0 10 20 30 40 50t, time ∆ N HIM ∆ N SS2 ∆ H HIM ∆ H SS2 -14 -12 -10 -8 -6 -4
0 10 20 30 40 50t, time ∆ C HIM ∆ C SS2 ∆ C HIM ∆ C SS2
Figure 13: (Collision with pursuing soliton on a fully resolved grid) The error in conservedquantities (left) ∆ N (solid), ∆ H (dotted), and (right) ∆ C (solid), and ∆ C (dotted) as afunction of time over the course of the simulation with HIM (green) and SS2 (red). Note thatSS2 demonstrates a strong peak in error in H at the time of soliton interaction. After theinteraction time the C , and the C exhibit large error with both SS2 and HIM. than in simulations with HIM method. In both methods we observe that con-servation of integrals of motion H , N , C and C does not imply highly accuratesolution in L ∞ -norm. In all the cases we found that HIM method gives smaller L ∞ error in the solution by a factor of at least 1 . L ∞ error we use the exact solution given by the for-mula (6). The simulation time is chosen so that there is a single collision inthe periodic box [ − L, L ]. The formula (6) gives a solution on an infinite line,whereas the simulation is performed on a periodic box and thus the simulationtime must not exceed the time it takes the solitons to reach the boundary ofthe box. Moreover, the soliton must still be exponentially small near the end ofthe box, so that the comparison with the exact formula is applicable.Despite the L ∞ error of the solution not being smaller than 10 − , we observethat the integrals of motion H , N are conserved up to 5 × − . Nevertheless,at the time of collision we found that ∆ H experiences a jump up to 5 orders ofmagnitude in SS2 method, while in HIM it is conserved by construction of themethod. Both methods exactly conserve N aside from accumulation of round-off errors over the course of simulations. The two nontrivial integrals of motion, C and C are not conserved exactly, nevertheless we observe that until the timeof collision these quantities vary only in 9-th decimal place. After the collisionthese values demonstrate a large jump (up to four orders of magnitude) in bothmethods. Unlike the Hamiltonian, H , in SS2 method, these integrals do notrevert to their original values after the collision.16 .5. Three Solitons Interactions Simulation It is known that solitons of the NLSE interact as particles, and interchangemomenta during collision [22]. The details of the process can be complicated,but once the solitons move sufficiently far from each other, they behave likeseparate pulses propagating without change of shape.In dimensionless units the one–soliton solution is given by (5). For thissimulation, the initial condition is the sum of three distinct one–soliton solutions:Φ( x, t = 0) = Φ + Φ + Φ , (33)where Φ , , are given by (5) with the following set of parameters: λ = 2 . , λ = 2 . , λ = 3 . , (34) v = 0 , v = 0 , v = 23 , (35) x , = 40 , x , = − , x , = − , (36)and zero initial phases. This set of parameters gives us two stationary solitonsand one moving. To make sure that we use approximation of a three-solitonsolution on a periodic boundary, we make the overlap between solitons is about10 − and at the boundary | Φ( x, t = 0) | ≈ − .After using the formulas (21), we translate this initial data to dimensionalunits. In the dimensional units the characteristic widths, τ c , and amplitudes, A , are given by: τ c = 1 ω √ λ ≈
50 ps A = √ λA ≈ . / and the value of λ varies from approximately is 2 . .
2. Whereas in theoriginal paper [6] the parameters of Gaussian pulses at the end of the fiber varyin amplitude from approximately 1 . − . / and have characteristic widths10 −
20 ps.The nonlinear time is given by t NL = π | Φ | = π | λ | ≈ . z NL = t NL l ≈
377 km. If transatlantic fiber is considered,this amounts to approximately 26 t NL . We will illustrate the performance ofHIM and SS2, on time scale of 400 t NL ≈
200 which is still physically relevant.The solution is computed on a grid of N = 4096 points (which correspondsto fully resolved spectrum of solution) with L = 25 π . The fixed point iterationstolerance is ε = 10 − for HIM method. The time step for the split step methodis chosen to be ∆ t SS = . x π . During simulation time 200 the numericalsolitons interact two times.In this simulation the results are presented in the Figure 14, we take ∆ t HIM =64∆ t SS , and due to larger time step HIM computation time is approximately5 .
76 times smaller. It takes 27 .
15 seconds for HIM, and 156 .
45 seconds for SS217 -8 -6 -4 -2 -80 -40 0 40 80x HIMSS2 10 -8 -6 -4 -2 -80 -40 0 40 80x HIMSS2 Figure 14: (Left) Soliton solution for SS2 (red) with ∆ t SS = . x π and HIM (green)∆ t HIM = 64∆ t SS . (Right) Soliton solution for SS2 (red) with ∆ t SS = . x π and HIM(green) ∆ t HIM = 128∆ t SS . to complete the computation on Intel Core I7-6700HQ CPU with frequency 2 . − for SS2,and 10 − for HIM while the time step for HIM is 64 times larger than for SS2.This time step allows HIM to accurately depict the positions of the interactingsolitons: at the final time the discrepancy in the location of stationary solitonswas less than ∆ x . Moreover, if the time step for HIM is increased to 128∆ t SS then the discrepancy in the location is still below 2∆ x and CPU time is 21 . .
35 times faster than SS2). We note that theamplitude of radiation in the tails of solitons scales as ∆ t for both methods. Inexact arithmetic and infinitely small ∆ t the magnitude of the solution in theseregions is exponentially small.In the Figure 15, we illustrate the conservation of integrals of motion byshowing the difference between the Hamiltonian, the number of particles, andthe integrals C and C at time t and its value at initial time. We note that thenumber of particles varies no more than 10 − for HIM, and less than 10 − forSS2. The value of Hamiltonian varies no larger than 10 − for HIM, however forSS2 it varies significantly at the time of soliton interaction. We note however,that the accuracy of actual solution is not representative of these number, andthe pointwise error of the numerical solution can be much larger. The integral C is equal to zero in this example, and is not presented in the figure, but theintegral C is not zero. It experiences jumps at the time of soliton interactions,and is conserved up to 10 − in HIM method due to the much larger time step,∆ t HIM = 64∆ t SS . 18 -14 -12 -10 -8 -6 -4 -2
0 50 100 150 200t, time ∆ N HIM ∆ N SS2 ∆ H HIM ∆ H SS2 -14 -12 -10 -8 -6 -4 -2
0 50 100 150 200t, time ∆ C HIM ∆ C SS2
Figure 15: (Left) Error in number of particles, ∆ N , and Hamiltonian, ∆ H for SS2 (red) with∆ t SS = . x π and HIM (green) with ∆ t HIM = 64∆ t SS . (Right) Error ∆ C for SS2 (red)with ∆ t SS = . x π and HIM (green) with ∆ t HIM = 128∆ t SS .
6. Conclusion
We performed detailed comparison of two algorithms for simulation of NLSE:Hamiltonian integration, proposed in [1] and the widely used split-step method.In all cases Hamiltonian integration demonstrates better conservation of Hamil-tonian at the time of soliton collision even for very large time steps. The otherconstants of motion N , C and C are conserved better by HIM when the timestep is the same or slightly larger than for split-step method. However, if thetime step is increased several orders of magnitude, the accuracy of conservationof integarls of motion in HIM may become lower. On the other hand, the point-wise error between the numerical solution and analytic formula is significantlylarger than the variation of conserved quantities, and thus reflects the qualityof the solution rather poorly. In experiments we observed this error to be about10 − –10 − in the maximum norm. For this reason a criterion on convergenceof fixed iterations by the number of particles that was used in the original pa-per [1] is suboptimal, and it is more accurate to control convergence by theresidual (18) as it was proposed in this paper.However, if the primary goal is to accurately portray the interaction of soli-tons over the physically relevant time such as propagation distance in opticalfiber, it is significantly more advantageous to use the HIM method with largetime step rather than SS2 method which requires smaller time step to satisfy thestability criterion. In our simulations for 400 nonlinear times, the time step forHIM is about 64–128 times larger than the instability criterion for SS2. How-ever, in a simulation for significantly longer time it may lead to accumulation19f errors in positioning of the solitons (jitter). For example if one simulates for4000 nonlinear times, the inaccuracy in the soliton position is about 10∆ x , andin order to keep the accuracy at ∆ x one would need to decrease the time stepfor HIM which results in smaller gains in computation time.The accurate portrayal of soliton interactions is crucial for the simulationof interactions in soliton gas [27, 28, 29, 30], or the fast–developing field ofintegrable turbulence [3], both methods are well–suited for this. At the sametime one should mention that split-step is simpler to implement and is moreefficient memory–wise. The split-step method is explicit, whereas HIM is animplicit method.As a summary, we would recommend to use Hamiltonian integration methodfor simulations requiring accurate description of soliton-soliton interactions orother subtle phenomena in Hamiltonian systems especially when computationtime is of the essence. Relevance of fast computational algorithms for opticalproblems can be illustrated by paper [6], [31] where massively parallel algorithmfor modification of NLSE was proposed and implemented. For multidimensionalturbulence (see for instance [32]), or for high accuracy short term dynamics, thesplit-step scheme of the order two, and higher order split step method [20], [33]can be the methods of choice. Acknowledgments
The authors wish to acknowledge gratefully the following contributions: SAAduring her graduate study was partially supported by NSF grants DMS-1554456and DMS-1500704; SAD was supported by NSF grant DMS-1716822; KAOwas supported by the Simons Collaboration on Wave Turbulence; PML thankssupport of Russian Ministry of Science and Higher education. The work ofPML was supported by the National Science Foundation, grant DMS-1814619.Simulations were performed at the Texas Advanced Computing Center using theExtreme Science and Engineering Discovery Environment (XSEDE), supportedby NSF Grant ACI-1053575. Also authors would like to thank developers ofFFTW [24] for this useful and free software.
Appendix A. Derivation of HIM method
Let H n = (cid:82) (cid:0) | Φ nx | − γ | Φ n | (cid:1) be the discretized in time Hamiltonian at the n -th time step. We consider change of Hamiltonian after one time step ∆ t :∆ H = H n +1 − H n = I + I , (A.1)where I := (cid:82) (cid:0) | Φ n +1 x | − | Φ nx | (cid:1) dx and I := γ (cid:82) (cid:0) | Φ n | − | Φ n +1 | (cid:1) dx . Weconsider I and I separately.By addition and subtraction to I of the following terms, Φ nx ¯Φ n +1 x and Φ n +1 x ¯Φ nx , under the integral sign, combining terms and using integration byparts, one gets: I = − (cid:90) (cid:0) ¯Φ n +1 xx ∆Φ + Φ nxx ∆ ¯Φ + Φ n +1 xx ∆ ¯Φ + ¯Φ nxx ∆Φ (cid:1) dx, n +1 − Φ n .By addition and subtraction to I of the four following terms, γ | Φ n +1 | Φ n ¯Φ n +1 , γ | Φ n +1 | Φ n +1 ¯Φ n , γ | Φ n | Φ n +1 ¯Φ n and γ | Φ n | Φ n ¯Φ n +1 , under the integral signand combining terms, we arrive at I = − γ (cid:90) (cid:0) ∆Φ( ¯Φ n +1 + ¯Φ n )( | Φ n +1 | + | Φ n | )+∆ ¯Φ(Φ n +1 +Φ n )( | Φ n +1 | + | Φ n | ) (cid:1) dx. After combining the like terms, we arrive at the formula∆ H = 12 (cid:90) [∆Φ (cid:16) − ¯Φ n +1 xx − ¯Φ nxx − γ n +1 + ¯Φ n )( | Φ n +1 | + | Φ n | ) (cid:17) +∆ ¯Φ (cid:16) − Φ n +1 xx − Φ nxx − γ n +1 + Φ n )( | Φ n +1 | + | Φ n | ) (cid:17) ] dx. (A.2)If we divide (A.2) by ∆ t and require that the first and second expressions insquare brackets are equal to i ∆¯Φ∆ t and − i ∆Φ∆ t correspondingly, then ∆ H vanishes.We note that: i Φ t = δ H δ ¯Φ i ¯Φ t = − δ H δ Φ . We get the following numerical scheme in time: i Φ n +1 − Φ n ∆ t = − (cid:2) Φ n +1 + Φ n (cid:3) xx − γ (Φ n +1 + Φ n )( | Φ n +1 | + | Φ n | )4 . (A.3) Appendix B. Derivation of the stability condition
We consider Φ n +1 ,s +1 = Φ n +10 + δ Φ s +1 and Φ n +1 ,s = Φ n +10 + δ Φ s whereΦ n +10 is the exact solution at the ( n + 1)-st time step. The iteration scheme isthe following:Φ n +1 ,s +1 k = 1 − ik ∆ t ik ∆ t Φ nk + i ∆ tγ ik ∆ t ˆ F (cid:2) ( | Φ n +1 ,s | + | Φ n | )(Φ n +1 ,s + Φ n ) (cid:3) . (B.1)Let’s keep only terms linear in δ Φ s +1 and neglect terms with small scale per-turbations δ Φ s : δ Φ s +1 k = i ∆ tγ ik ∆ t (cid:2) | Φ n +10 | + | Φ n | + Φ n ¯Φ n +10 (cid:3) δ Φ sk ++ i ∆ tγ ik ∆ t (cid:2) (Φ n +10 ) + Φ n Φ n +10 (cid:3) ¯ δ Φ sk (B.2)21herefore, we can compose the following system of linear equations: (cid:20) δ Φ s +1 k δ ¯Φ s +1 k (cid:21) = (cid:18) c (cid:2) | Φ n +10 | + | Φ n | + Φ n ¯Φ n +10 (cid:3) c (cid:2) (Φ n +10 ) + Φ n Φ n +10 (cid:3) ¯ c (cid:2) ( ¯Φ n +10 ) + ¯Φ n ¯Φ n +10 (cid:3) ¯ c (cid:2) | Φ n +10 | + | Φ n | + ¯Φ n Φ n +10 (cid:3)(cid:19) (cid:20) δ Φ sk δ ¯Φ sk (cid:21) (B.3)where c = iγ ∆ t ik t We need the matrix on the right hand side of (B.3) (letsname it A ) to be a contracting map. As a result, we require its determinant tobe smaller than 1. From | det( A ) | <
1, we can get the condition for iterationsconvergence of HIM: ∆ t < | γ |√ | Φ n | ) (B.4) References [1] S. Dyachenko, A. C. Newell, A. Pushkarev, V. E. Zakharov, Optical turbu-lence: weak turbulence, condensates and collapsing filaments in the nonlin-ear Schr¨odinger equation, Physica D: Nonlinear Phenomena 57 (1) (1992)96–160.[2] C. Sulem, P. Sulem, The Nonlinear Schr¨odinger Equation: Self-Focusingand Wave Collapse, Applied Mathematical Sciences, Springer New York,1999.URL https://books.google.com/books?id=uTxYaEztjzgChttps://books.google.com/books?id=uTxYaEztjzgC