Stéphane Balac
email: stephane.balac@univ-rennes1.fr

Arnaud Fernandez
email: afernand@laas.fr

Mathematical analysis of adaptive step-size techniques when solving the nonlinear Schrödinger equation for simulating light-wave propagation in optical fibers

Keywords: Symmetric Split-Step method, adaptive step-size control, nonlinear Schrödinger equation, embedded Runge-Kutta methods, Intercation Picture method, nonlinear optics

come

Mathematical analysis of adaptive step-size techniques when solving the nonlinear Schrödinger equation for simulating light-wave propagation in optical fibers

Introduction

The Nonlinear Schrödinger Equation (NLSE) describes a wide class of physical phenomena among which propagation of light in an optical fiber. Here, we are concerned by the following form of the NLSE referred as the Generalized Nonlinear Schrödinger Equation (GNLSE) (see [START_REF] Fernandez | Nonlinear fiber optics[END_REF][START_REF] Okamoto | Fundamentals of Optical Waveguides[END_REF])

∂ ∂z A(z, t) = - α(z) 2 A(z, t) +        n max n=2 i n+1 β n n! ∂ n ∂t n A(z, t)        + iγ I d + i ω 0 ∂ ∂t • A(z, t) (1 -f R ) A(z, t) 2 + f R R h R (s) |A(z, t -s)| 2 ds , (1)
where the complex valued function A represents the slowly varying pulse envelope of a quasi-monochromatic optical wave at frequency ω 0 in a local frame moving with the pulse at the group velocity v g = c/n g where n g denotes the group index of the fibre, z represents the position along the fiber, t the time in the local frame and I d stands for the identity operator. In the situation considered here, the optical wave is assumed to be an electric field E at frequency ω 0 which is linearly polarized along the vector e x transverse to the propagation's direction e z defined by the fiber axis and expressed as a function of position r = (x, y, z) and time τ in the form E(r, τ) = A(z, t) F(x, y) e -i(ω 0 τ-k z) e x [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF] where F(x, y) is the electric wave transverse representation also called the modal distribution and k is the wavenumber. The relation between the "local" time t in the retarded frame and the absolute time τ is t = τz/v g . The physical effects for wave propagation in fibre taken into account in [START_REF] Fernandez | Nonlinear fiber optics[END_REF] are the following. First, some linear effects are expressed through the linear attenuation/gain coefficient α and the linear dispersion coefficients β n , 2 ≤ n ≤ n max . Some nonlinear effects are involved through the nonlinear parameter γ. Moreover, first order par-tial derivation with respect to time takes into account the dispersion of the nonlinearity through the simplified optical shock parameter τ shock = 1/ω 0 . Instantaneous Kerr effect manifests itself through the term (1f R ) |A| 2 . The delayed Raman contribution in the time domain is taken into account through the convolution product between the instantaneous power |A| 2 and the Raman time response function h R . The constant f R represents the fractional contribution of the delayed Raman response to nonlinear polarization. The most widely used numerical method for solving the NLSE or GNLSE in optics is the Symmetric Split-Step Fourier (S3F) method, see e.g. [START_REF] Deiterding | A reliable Split-Step Fourier method for the propagation equation of ultra-fast pulses in singlemode optical fibers[END_REF][START_REF] Feng | An efficient Split-Step time-domain beam-propagation method for modeling of optical waveguide devices[END_REF][START_REF] Hermansson | Optical coherence calculations with the Split-Step fast Fourier transform method[END_REF][START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF][START_REF] Taha | Analytical and numerical aspects of certain nonlinear evolution equations. II. Numerical, nonlinear Schrödinger equation[END_REF], due to its particular simplicity and efficiency for solving certain type of evolution partial differential equations (PDE). The idea behind the S3F method applied to the GNLSE (1) is to decompose over each subinterval of a given subdivision of the fiber length the PDE problem into a sequence of 3 simpler problems connected to each other. One corresponding to a purely linear PDE problem over the first half of the subinterval, the other over the whole subinterval corresponding to a nonlinear ordinary differential equation (ODE) problem with the time variable as a parameter and last another purely linear PDE problem over the second half of the subinterval, see e.g. [START_REF] Fernandez | Nonlinear fiber optics[END_REF][START_REF] Mclachlan | Splitting methods[END_REF]. Of course solving this sequence of 3 nested problems is not equivalent to solving the GNLSE [START_REF] Fernandez | Nonlinear fiber optics[END_REF] and an error, usually referred to as the "splitting error", is introduced. However, the interest of this numerical approach is that each of the 3 nested problems can be solved much more easily than the GNLSE (1) itself and the approximated solution converges to the solution of the GNLSE (1) when the discretization step-size tends toward 0. The 2 linear PDE problems have an explicit solution obtained by use of the Fourier Transform tool whereas the nonlinear ODE problem is generally solved numerically by standard methods for ODE such as Runge-Kutta (RK) methods [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF]. Among them, the 4th order classical RK method is the most famous due to the good compromise it offers between accuracy and computational cost. Although the S3F method is the most widely used with the 4th order RK method for solving the GNLSE [START_REF] Fernandez | Nonlinear fiber optics[END_REF], alternative ODE solvers such as the Adams multistep methods have been investigated [START_REF] Fornberg | A fast spectral algorithm for nonlinear wave equations with linear dispersion[END_REF][START_REF] Liu | A fast method for nonlinear Schrodinger equation[END_REF]. A second kind of approximation error is therefore introduced by the use of a numerical scheme for solving the ODE problem.

Of course, the practical efficiency of a numerical method such as the S3F method applied to the GNLSE (1) highly depends on the distribution of the discretization grid points along the fibre and the use of an adaptive step-size control strategy is mandatory. The idea behind an adaptive step-size strategy is to introduce the grid points during the progress of the computation taking into account the information available at the current computation stage in order to determine the best suited step size (and therefore the next grid point) so as to maintain a given predefined accuracy of the approximation. Namely, if computations have to be achieved at a given accuracy, the size of the steps have to be chosen small enough so as to attain the desired accuracy but not smaller than required to avoid unnecessary increase of the computation time. Since it is not possible in practice to know the final global error, the step-size is determined so that, at each computational step, a "local error" (or any related quantity) is lower than a prescribed tolerance. A variety of adaptive step-size strategies have been propound to be used in conjunction with the S3F method. We can distinguish 2 types of approaches. The one based on physical concepts (or physical intuition) where at each grid point the step-size is chosen so as a "local error" related quantity estimated from a physical quantity matches the prescribed tolerance value. For instance, in the so-called "nonlinear phase rotation method" [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF] the step-size is chosen so that the phase change due to nonlinearity does not exceed a certain limit. In the "walk-of method", the step-size is chosen to be inversely proportional to the product of the absolute value of the dispersion and the spectral bandwidth of the signal and the method applies to low power, multichannel systems [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF]. In [START_REF] Meirelles | Study on a new Split-Step Fourier algorithm for optical fiber transmission systems simulations[END_REF][START_REF] Rieznik | Uncertainty relation for the optimization of optical-fiber transmission systems simulations[END_REF] a method termed the "uncertainty principle method" is propound where the determination of the stepsize is done from the values of a parameter derived from an inequality which in quantum mechanics gives rise to the uncertainty principle between two non commuting operators. In [START_REF] Blow | Theoretical description of transient stimulated Raman scattering in optical fibers[END_REF] it is shown that when losses in the fibre are neglected the following quantity, referred as the "optical photon number" (OPN), is conserved

P(z) = R n eff A eff ω | A(z, ω)| 2 dω
where A denotes the Fourier transform of the slowly varying pulse envelope of the electric field A, A eff is known as the effective mode area, n eff is the effective index of the fibre, and integration hold over the entire spectrum of the optical wave amplitude assumed to have a bandwidth less than ω 0 3 . In [START_REF] Heidt | Efficient adaptive step size method for the simulation of supercontinuum generation in optical fibers[END_REF] it is made use of the conservation of the OPN to estimate a "local error" related quantity and to define an adaptive step-size control strategy termed the "Conservation Quantity Error (CQE) method". The CQE method applies to low loss fibres.

An other approach for defining adaptive step-size strategies consists in using mathematical concepts from the field of numerical analysis. The "local error" is then defined as the error made by using the numerical scheme when computing an approximation of the solution at the current grid point under the assumption that the value at the previous grid point was exact. In [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF] the authors propound to use the well known step-doubling concept (often referred as Richardson extrapolation) [START_REF] Üeberhuber | Numerical Computation II: Methods, Software, and Analysis[END_REF] to estimate the local error when solving the NLSE by the S3F method. They termed their adaptive step-size strategy the "local error method" (LEM). The main advantage of such a numerical analysis based approach is that since no assumption on the physical quantities involved in the equation is made, the step-size strategy applies to an arbitrary set of parameters in the NLSE. Note that a comprehensive numerical comparison of the LEM and other adaptive step-size methods such as the "nonlinear phase rotation" method and the "walk-off" method is achieved in [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF] on benchmark problems in optics. However a drawback of the step-doubling method is a computational overcost of approximatively 50 % (when compared to the same S3F method with the "optimal" grid points distribution given in advance) and some attempts have been made to find alternative method with the same degree of generality, see e.g. [START_REF] Liu | Adaptive higher-order Split-Step Fourier algorithm for simulating lightwave propagation in optical fiber[END_REF]. How-ever the method proposed in this journal in [START_REF] Liu | Adaptive higher-order Split-Step Fourier algorithm for simulating lightwave propagation in optical fiber[END_REF] suffers from misinterpretations in the mathematical concepts of numerical analysis used to derive the algorithm and this convinces us of the soundness of a clarification on the topic of local error estimate and adaptive step-size control when solving the GNSLE by the S3F method.

The paper is organized as follows. Section 2 is devoted to a succinct presentation of the S3F method aimed at setting the framework for the study of local error estimate approaches made in section 3. We also discuss in section 3 some potential numerical issues and inconsistency between the theory of adaptive step-size control and practical results observed on numerical simulation experiments. In section 4 we detail how an adaptive step-size control strategy for the S3F method can be derived from the knowledge of the local error. Finally, we conclude in section 5 by a reference to an alternative method to the S3F method, termed the "Interaction Picture" (IP) method, that offers the same advantages than the S3F method in terms of ease of implementation and computational cost together with some nicer features for adaptive step-size purposes.

Overview of the Symmetric Split-Step method

Mathematical framework

Let L 2 (R, C) denotes the infinite-dimensional vector space of all complex-valued "square integrable" (also termed "finite energy") functions which are defined on the real line and have

a finite energy norm f = R | f (t)| 2 dt 1 2
. We denote by F the Fourier operator from L 2 (R, C) to L 2 (R, C) defined by continuous extension of the Fourier Transform definition for integrable functions and computed for all f ∈ L 2 (R, C) by an improper integral as

∀ν ∈ R F ( f )(ν) = lim T →+∞ T -T f (t) e 2iπνt dt.
Alternatively we use the notation f (ν) for the Fourier Transform of f and we use the notation F -1 for the inverse Fourier Transform.

For a fixed z ∈ [0, L] we denote by A(z) the first partial function of the slowly varying pulse envelope A in z, i.e. A(z) : t ∈ R → A(z, t). This notation will be used throughout the paper. We also denote by D the linear operator

D : A(z) → n max n=2 β n i n+1 n! ∂ n ∂t n A(z) (3)
and by N the nonlinear operator

N : A(z) → iγ I d + i ω 0 ∂ ∂t • A(z) (1 -f R ) A(z) 2 + f R |A(z)| 2 ⋆ h R - 1 2 αA(z) (4)
where ⋆ stands for the convolution product. With these notations, the GNLSE (1) reads

∂ ∂z A(z) = D A(z) + N(A(z)) ∀z ∈ [0, L] (5)
where L denotes the fibre length, and PDE ( 5) is solved together with the initial condition A(0) = a 0 at the fibre entrance z = 0 where a 0 is a known source term. Note that another splitting is possible: the term -1 2 αA(z) can be added to the linear operator D instead of the nonlinear operator N.

It is worth mentioning that the standard nonlinear Schrödinger equation corresponds to n max = 2 in the definition (3) of D and to the nonlinear operator

N 0 : A(z) → iγ A(z) |A(z)| 2 - 1 2 αA(z) (6)
which is formally retrieved from (4) in the special case when we set f R = 0 and 1/ω 0 = 0.

The splitting approximation scheme

In the S3F method applied to the GNLSE (5), the interval [0, L] is divided into K subintervals where the spatial grid points are denoted z k , k ∈ {0, . . . , K}, where 0

= z 0 < z 1 < • • • < z K-1 < z K = L.
We denote by h k = z k+1z k the step-size between grid points z k and z k+1 and we set z k+ 1 2 = z k + h k 2 . The S3F method consists in solving over each subinterval [z k , z k+1 ] for k ∈ {0, . . . , K-1}, the following 3 nested problems with time variable t as a parameter and the operators D and N defined as in ( 3) and (4):

         ∂ ∂z U(z) = D U(z) ∀z ∈ [z k , z k+ 1 2 ] U(z k ) = A k-1 (z k ) (7)
where A k-1 (z k ) represents the approximated solution at grid point z k computed at step k -1;

           ∂ ∂z V(z) = N(V(z)) ∀z ∈ [z k , z k+1 ] V(z k ) = U(z k+ 1 2 ) (8)
where U(z k+ 1 2 ) represents the solution to problem (7) at half grid point z k+ 1 2 ;

           ∂ ∂z W(z) = D W(z) ∀z ∈ [z k+ 1 2 , z k+1 ] W(z k+ 1 2 ) = V(z k+1 ) (9)
where V(z k+1 ) represents the solution to problem (8) at node z k+1 . An approximated solution to the GNLSE (5) at grid node z k+1 is then given by A k (z k+1 ) = W(z k+1 ). The principle of the Symmetric Split-Step scheme is depicted in figure 1. The main interest of such a decomposition is that each one of the 3 problems ( 7), ( 8) and ( 9) can be solved much more easily than the GNLSE (5) considered over the interval [z k , z k+1 ] and for h k small enough, A k (z k+1 ) provides an approximation of the solution to the GNLSE at grid point z k+1 . From a mathematical point of view, it is well known [9] that the linear operator D defined in (3) generates a continuous group of bounded operators on L 2 (R, C), denoted (according to its properties) by exp (zD) with z ∈ R. Moreover, the solution to the nonlinear problem (8) admits the following integral representation form

∀z ∈ [z k , z k+1 ] V(z) = U(z k+ 1 2 ) exp z z k N(V(ζ)) dζ
where

N is such that N(V(ζ)) = N(V(ζ)) V(ζ).
This justifies why the approximated solution given by the S3F method by solving the sequence of nested problems ( 7)-( 8)-( 9) is the most often formally written as

A k (z k+1 ) = exp( h k 2 D) exp z k+1 z k N(V(ζ)) dζ exp( h k 2 D) A k-1 (z k ).
The convergence of Split-Step methods applied to various forms of the Schrödinger equation is widely documented in the litterature, see e.g. [START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF][START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF][START_REF] Zharnitsky | Averaging for Split-Step scheme[END_REF] where the authors prove that the global convergence order of the Symmetric Split-Step method is 2. Namely, in the situation considered here, the following result holds for the local error [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF]. Let A denotes the solution to equation ( 5) and for all k ∈ {0, . . . , K -1} let A k (z k+1 ) denotes its approximation at grid point z k+1 computed by solving the three nested problems ( 7)-( 8)-( 9) with the initial condition A k 1 (z k ) = A(z k ) corresponding to the exact solution (because we are concerned by the local error) assumed to be a sufficiently smooth function. Then, for h k in a neighborhood of 0 we have the following estimate:

A(z k+1 ) = A k (z k+1 ) + O(h 3 k ) ( 10
)
where we recall the significance of the big O notation O(h 3 k ): there exists a positive number C such that for h k in a neighborhood of 0 the local error satisfies

A(z k+1 ) -A k (z k+1 ) ≤ Ch 3
k . From a practical point of view, the quantity of interest is the "global error". The global error at grid point z k+1 takes into account the error accumulation from all the previous k + 1 computational steps and therefore from [START_REF] Deiterding | A reliable Split-Step Fourier method for the propagation equation of ultra-fast pulses in singlemode optical fibers[END_REF] we can deduce that it behaves in O(h 2 ) where h = max i=0,...,k h i . The Symmetric Split-Step method is a 2nd order method.

The 2 linear PDE problems ( 7) and ( 9) admit an explicit solution obtained by use of the Fourier Transform. The solution to problem [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF]

at grid point z k+ 1 2 reads U(z k+ 1 2 ) = F -1 ν → A k-1 (z k , ν) e d(ν) h k 2 where d : ν → i n max n=2 β n n! (2πν
) n and the solution to problem (9) at grid point z k+1 reads

W(z k+1 ) = F -1 ν → V(z k+1 , ν) e d(ν) h k 2 .
When we are concerned by the standard nonlinear Schrödinger equation where the nonlinear operator N has the form given by N 0 in [START_REF] Blow | Theoretical description of transient stimulated Raman scattering in optical fibers[END_REF], there also exists an explicit solution to the nonlinear ODE problem [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF]. In the general case when α 0, the solution at grid point z k+1 reads

V(z k+1 ) = U(z k+ 1 2 ) exp - α 2 h k + iγ α |U(z k+ 1 2 )| 2 (e αh k -1)
and when α = 0

V(z k+1 ) = U(z k+ 1 2 ) exp iγh k |U(z k+ 1 2 )| 2 .
Even if the 3 problems ( 7)-( 8)-( 9) have explicit solutions when the nonlinear Schrödinger equation is considered, it doesn't signify for all that the nonlinear Schrödinger equation itself has an explicit solution. However an approximated solution can be computed by the S3F method and the approximation error solely lies in the use of a Split-Step scheme. Unlike, in the more general case when the GNSLE is considered, a numerical approximation scheme (such as a Runge-Kutta scheme) is mandatory for solving the nonlinear ODE problem [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] resulting in an approximation error made up of the splitting error and the ODE approximation error.

Local error estimate theory

In the field of numerical analysis, an error estimate refers to an approximated value for the error on a numerically calculated quantity. As mentioned in the introduction, the most general way of designing an adaptive step-size control strategy consists in estimating the "local error" at each step of the S3F algorithm, i.e. the error made in computing the approximated solution at a given grid point assuming that the data from the previous grid point was exact. At each grid point, the step-size is then computed so that the local error is lower than a predefined value (termed the tolerance). The tolerance is set depending on the need of accurate results, e.g. to 10 -q with q ranging from 3 to 9. Actually, it is the "global error" (i.e. the total error at a given grid point ensued from the error accumulation from the initial grid point) that we wish to know about but estimating the global error with accuracy is mathematically and numerically out of reach and therefore we must content ourselves with estimating the local error. Note that the global error resulting from the error accumulation from the previous computational steps can be several order of magnitude upper than the tolerance for the local error.

In this section we carry out an analysis of the local error when a numerical scheme is used for solving the nonlinear ODE problem (8) on the one hand and to an analysis of the Split-Step local error on the other hand.

Local error when solving the nonlinear ODE problem (8)

by a Runge-Kutta method The most effective way for local error estimate when solving a nonlinear ODE problem is by using an Embedded Runge-Kutta (ERK) method. ERK methods [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF] are special Runge-Kutta (RK) methods designed to deliver 2 approximations of the solution of the ODE problem under consideration corresponding to 2 RK schemes of different convergence orders p and p + 1. These 2 approximations of the solution can be considered as an accurate approximated solution (the one computed with the numerical scheme of higher order p + 1) and a coarse approximated solution (the one computed with the one of lower order p). For efficiency, the 2 RK schemes bear several computational stages in common in order to reduce noticeably the computational cost of the local error estimation.

In [START_REF] Liu | Adaptive higher-order Split-Step Fourier algorithm for simulating lightwave propagation in optical fiber[END_REF] the author uses (without mentioning it explicitly) an ERK method of order 1 and 2 (p = 1) for solving problem [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] where the first order approximated solution in z k+1 is given by V [1] k+1 = U(z k+ 1 2 )

+ h k α 1 , α 1 = N(U(z k+ 1 2 )
) and the second order approximated solution is given by V [2] k+1 = U(z k+ 1 2 )

+ h k α 2 , α 2 = N(U(z k+ 1 2 ) + 1 2 h k α 1 )
. Because it represents a good compromise between accuracy and computational cost, the fourth order Runge-Kutta (RK4) scheme [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF] is more commonly used in conjunction with the S3F method. The RK4 method can be formulated as an ERK method of order 3 and 4 (p = 3), referred in the literature as Dormand and Prince RK 4(3) T method [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF], which preserves the ease of implementation of the RK4 formula and provides a local error estimate at no extra computational cost (see [START_REF] Balac | Embedded Runge-Kutta scheme for step-size control in the Interaction Picture method[END_REF] for the implementation details). The 4th order RK method delivers the approximated solution V [4] k+1 = U(z k+ 1 2 ) +

h k 6 (α 1 + 2α 2 + 2α 3 + α 4 )
where

α 3 = N(U(z k+ 1 2 ) + 1 2 h k α 2 ) and α 4 = N(U(z k+ 1 2 ) + h k α 3
) and the 3rd order RK method provides the approximated solution

V [3] k+1 = U(z k+ 1 2 ) + h k 30 (5α 1 + 5α 2 + 5α 3 + 2α 4 + 3α 5 )
where α 5 = N(V [4] k+1 ). When solving problem (8) by the standard RK method of order p, provided all the partial derivatives of z → N(V(z)) exist and are continuous up to the order p + 1, the local error is known to behave as O(h p+1 ), see e.g. [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF]. Let us detail the way the local error can be estimated by using an ERK method with p = 1 or p = 3 as considered above. By use of the Taylor expansion formula, one can show that the local errors for each of the 2 methods are respectively given by [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF]

ℓ RK(p) k+1 = V(z k+1 ) -V [p] k+1 = ϕ p (z k , V [p] k ) h p+1 k + O(h p+2 k ) (11) ℓ RK(p+1) k+1 = V(z k+1 ) -V [p+1] k+1 = ϕ p+1 (z k , V [p+1] k ) h p+2 k + O(h p+3 k ) ( 12
)
where V(z k+1 ) denotes the exact solution to problem (8) at grid point z k+1 and ϕ p (resp. ϕ p+1 ) is a function of the elementary differentials of order p (resp. p + 1) of the function z → N(V(z)). By difference of these 2 relations we obtain

V [p+1] k+1 -V [p] k+1 = ϕ p (z k , V [p] k ) h p+1 k + O(h p+2 k ). (13)
From ( 11) and ( 13) we deduce that the local error for the p-th order RK method at grid point z k+1 can be approximated with an error in O(h p+2 k ) by

∀t ∈ R ℓ RK(p) k+1 (t) ≈ V [p+1] k+1 (t) -V [p] k+1 (t) (14)
and the energy norm of the local error by

L RK(p) k+1 = ℓ RK(p) k+1 ≈ R V [p+1] k+1 (t) -V [p] k+1 (t) 2 dt 1 2 . ( 15
)
It satisfies

L RK(p) k+1 = O(h p+1 k ).

Split-Step local error estimation by step-doubling

As mentioned before, when the nonlinear Schrödinger equation is considered (i.e. N = N 0 as defined in ( 6)) the nonlinear ODE problem (8) admits an analytical solution and the local error in the S3F method only amounts to the splitting error as given by [START_REF] Deiterding | A reliable Split-Step Fourier method for the propagation equation of ultra-fast pulses in singlemode optical fibers[END_REF]. In such a case, although an approach similar to the one performed to estimate the local error when solving a nonlinear ODE by an ERK method could be considered (involving 2 embedded Split-Step schemes of different orders [START_REF] Balac | An Embedded Split-Step method for solving the nonlinear Schrodinger equation in optics[END_REF][START_REF] Koch | Embedded Split-Step formulae for the time integration of nonlinear evolution equations[END_REF]), the usual method for estimating the splitting local error is based on the step-doubling technique. Actually, it is the adaptive step-size control used in conjunction with the S3F method for solving the NLSE presented in [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF]. The idea behind the step doubling method (also known as Richardson extrapolation method and widely used for adaptive step-size control in quadrature methods for integral or ODE [START_REF] Üeberhuber | Numerical Computation II: Methods, Software, and Analysis[END_REF]) is the following. The local error for the S3F method at grid point z k+1 satisfies the following relation [START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF][START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF][START_REF] Zharnitsky | Averaging for Split-Step scheme[END_REF]

ℓ SS k+1 = A(z k+1 ) -A k (z k+1 ) = C h 3 k + O(h 4 k ) (16)
where C denotes a positive number independent of h k . Let u k+1 = A k (z k+1 ) be the solution at grid point z k+1 computed by the Symmetric Split-Step scheme ( 7)-( 8)-( 9) from grid point z k using one step of size h k , and let u k+1 be the solution computed from grid point z k using two half steps of size h k /2, in both cases assuming the solution A k-1 (z k ) at grid point z k to be exact (since we are interested in estimating the local error). Additionally the following assumption is made:

(H) the local error after 2 half steps is twice the local error after one half step (which consists in an approximation since the initial data for the second half-step is actually the approximated solution computed after the first half-step and not the exact one). Then, from [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibers[END_REF] we deduce that

A(z k+1 ) -u k+1 = C h 3 k + O(h 4 k ) (17) A(z k+1 ) -u k+1 = 2 C h k 2 3 + O(h 4 k ) (18)
and therefore by difference between these 2 relations we get

u k+1 -u k+1 = 3 4 C h 3 k + O(h 4 k ). ( 19
)
Thus the local error related to the Split-Step scheme at grid point z k+1 can be approximated, with an error in O(h 4 k ), by

∀t ∈ R ℓ SS k+1 (t) ≈ C h 3 k ≈ 4 3 ( u k+1 (t) -u k+1 (t)). ( 20
)
The energy norm of the local error at grid point z k+1 is then expressed as

L SS k+1 = ℓ SS k+1 ≈ 4 3 R | u k+1 (t) -u k+1 (t)| 2 dt 1 2 (21)
and from [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibers[END_REF] we have L SS k+1 = O(h 3 k ). We have to point out that relation [START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF] gives an approximation of the local error corresponding to the solution u k+1 computed by the Symmetric Split-Step scheme over the coarse grid with step-size h k . However the fine mesh grid solution u k+1 is a better approximation and in practice it is kept as the approximated solution and propagated along the fibre. This process is referred in the literature as the "local extrapolation mode". Moreover, by difference of equation [START_REF] Koch | Embedded Split-Step formulae for the time integration of nonlinear evolution equations[END_REF] and equation [START_REF] Liu | Adaptive higher-order Split-Step Fourier algorithm for simulating lightwave propagation in optical fiber[END_REF] pre-multiplied by 4 we obtained under assumption (H) an even more accurate approximation (referred in the literature as Richardson extrapolation): A(z k+1 ) can be approximated with an error in O(h 4 k ) (instead of O(h 3 k ) as given by ( 18)) as

A(z k+1 ) ≈ 1 3 4 u k+1 (t) -u k+1 (t) .
The over-cost of estimating the local error in this way is the cost of the computation of the coarse mesh grid solution and this cost is approximately half the cost of the computation of the fine mesh grid solution since the step-size is twice larger. Thus, estimating the local error using the step doubling approach is liable of an extra computational cost of 50% more than the cost of the computation of the approximated solution itself. In fact, the additional cost is slightly less than 50% since some computations are shared by each of the 2 methods, and needs to be carried out only once. This approach is the one of the LEM [START_REF] Rieznik | Uncertainty relation for the optimization of optical-fiber transmission systems simulations[END_REF]. Finally, we would like to mention that the step-doubling approach could be used to estimate the local error when solving the nonlinear ODE problem [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] in the S3F method by the RK4 method. (In the above reasoning, equations ( 17) and ( 18) would be replaced by [START_REF] Feng | An efficient Split-Step time-domain beam-propagation method for modeling of optical waveguide devices[END_REF] and [START_REF] Fornberg | A fast spectral algorithm for nonlinear wave equations with linear dispersion[END_REF].) However the computational cost of this way of estimating the local error of the RK4 method is prohibitive compared to the use of the ERK4(3) method.

Local error analysis when solving the NLSE and GNLSE

by the S3F method

Local error analysis for the GNLSE

In the case when the GNLSE (1) is solved by the S3F method, a dedicated numerical method such as a Runge-Kutta method is to be used to solve the nonlinear ODE problem [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF]. As a consequence, the local error consists in the error produced by the Split-Step scheme and the error produced by the (p + 1)-th order Runge-Kutta scheme (p = 1 or p = 3 for the 2 above mentioned ERK methods). Since the 2 linear PDE problems [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF] and ( 9) are solved exactly by use of the Fourier Transform, from ( 12) and ( 17) we deduce that the local error behaves for h k in a neighborhood of 0 as

A(z k+1 ) -A k (z k+1 ) ≤ C h 3 k + K h p+2 k ( 22
)
where C and K denote 2 positive numbers. The constant number K depends on the elementary differentials of order 3 of the mapping A(z) → N(A(z)) whereas C depends on the elementary differentials of order p + 1 of the mapping F : A(z) → DA(z) + N(A(z)). Therefore the 2 constant numbers C and K depend in a different way on the physical parameters involved in the GNLSE (1) and it can happen that they have values different by several order of magnitude. Note that the local error can not be estimated only for the nonlinear ODE problem [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] as it is done in [START_REF] Liu | Adaptive higher-order Split-Step Fourier algorithm for simulating lightwave propagation in optical fiber[END_REF] and in no case when p = 1 the whole numerical approximation scheme can be 4th order accurate as the author states. When p = 1, relation [START_REF] Meirelles | Study on a new Split-Step Fourier algorithm for optical fiber transmission systems simulations[END_REF] shows that the global convergence order of the S3F method is 2. From a theoretical point of view, we deduce from ( 22) that when p = 3 the S3F method with problem (8) solved by the RK4 scheme (referred as the S3F-RK4 method in the sequel) has global convergence order 2 since for h k in a neighborhood of 0 we have K h 5 k = O(h 3 k ). Thus, the local error ( 22) could be estimated by the step doubling approach just as detailed in section 3.2. However in practice it can happen that under special circumstances (i.e. for a special set of values for the physical parameters involved in the GNLSE (1)) the value of the constant C is much lower than the value of K. In such a case, when the convergence order of the S3F-RK4 method is explored on numerical examples, the experimental convergence curve can behave as a function of the step-size h proportional to h 4 instead of h 2 . If so, the local error estimate presented in section 3.2 is likely to provide wrong results since it was assumed that the convergence order is 2. Actually, it is not difficult to convince oneself that whenever the nonlinear effects predominate in the GNLSE (1) the convergence order observed experimentally will be close to 4 and C << K whereas when it is the linear effects that predominate the convergence order observed experimentally will be close to 2 and K << C. (Note that when dealing with the convergence behavior, the 2 constants K and C we are referring to may differ from the constants appearing in the local error estimate [START_REF] Meirelles | Study on a new Split-Step Fourier algorithm for optical fiber transmission systems simulations[END_REF] but we have maintained the same notation for simplicity.)

In order to numerically illustrate the mathematical analysis presented in this section we have solved the GNLSE by the S3F-RK4 method on a test example chosen to match with a typical case of high speed data propagation through a L = 20 km single mode fibre in optical telecommunication with a data's carrier frequency located in the C band of the infrared spectrum ( f 0 = 193 Thz). The following set of fibre's parameters were used for the simulation: α = 0.046 km -1 , γ = 4.3 W -1 km -1 , f R = 0.245, β 2 = -19.83 ps 2 km -1 , β 3 = 0.031 ps 3 km -1 and β n = 0 for n ≥ 4. An expression for the Raman time response function for silica core fibre is given in [START_REF] Fernandez | Nonlinear fiber optics[END_REF]. The source term a 0 = A(z = 0) was represented as a first order Gaussian pulse:

a 0 : t → P 0 e -1 2 (t/T 0 ) 2 (23)
where T 0 is the pulse half-width at 1/e intensity point and P 0 is the pulse peak power. Simulations were carried out for a pulse-width T 0 = 6.8 ps and for three different peak power values: P 0 = 5 mW, 25 mW and 50 mW. As explained in [START_REF] Fernandez | Nonlinear fiber optics[END_REF] the dispersive length L D and nonlinear length L NL defined as

L D = T 2 0 |β 2 | and L NL = 1 γP 0 ( 24
)
provide the length scales over which dispersive or nonlinear effects become important for pulse evolution. In the test example considered here the dispersion length L D = 2.26 km was kept constant and the nonlinear length L NL was decreasing from 46.5 km to 4.65 km while P 0 was increasing from 5 mW to 50 mW. Thus when the pulse peak power was increasing the nonlinear Kerr effect was becoming more important through propagation. The convergence curves representing the global error at the fiber end for the energy norm when using a constant step-size S3F-RK4 method for solving the GNLSE versus the step-size are depicted in Figure 2. The global error at fiber end takes into account the error accumulation from the fiber entrance and therefore its asymptotic behavior is proportional to the local error by a factor 1/h. In Figure 2, one can observe that for a 50 mW peak power pulse (where L 4 L NL ) the experimental convergence curve (solid line with square symbols) drawn in a logarithmic scale has a slope coefficient close to 4 (dotted line labeled 4). In this specific case the nonlinear effects predominate in the GNLSE and the error behavior is mainly governed for the step-size range considered in the simulation by the error committed when solving the nonlinear ODE problem by the 4th order RK scheme. On the contrary, for a 5 mW peak power pulse (where L 1/2 L NL ) the experimental convergence curve (solid line with star symbols) has a slope coefficient close to 2 (dotted line labeled 2). The linear dispersion effects predominate in the GNLSE and the error behavior is mainly governed for the step-size range considered in the simulation by the error inherent to the use of the Symmetric Split-Step scheme (the error committed when solving the nonlinear ODE problem by the 4th order RK scheme is much lower). The convergence curve corresponding to a 25 mW peak power pulse (where L L NL and both linear dispersion and nonlinear effects have a comparable strength) shows a very interesting phenomenon (solid line with circle symbols). The first part of the convergence curve has a slope coefficient close to 4 whereas it is close to 2 on the second part. This can be understood as follows. The value of the constant C involved in the expression of the error related to the use of the Symmetric Split-Step scheme is a little smaller than the value of the constant K involved in the expression of the error of the RK scheme. As a consequence, for the largest values of the step-size h in the range considered in the simulation we have K h 4 > C h 2 whereas K h 4 < C h 2 when h becomes small enough, i.e. for the lowest values in the range considered in the simulation.

To be comprehensive in the analysis of the results depicted on Figure 2, one can notice the somehow erratic behavior of the curves for the largest values of the step-size (upper right corner). It can be explained by the large values of the global error for such step-sizes which are not really compatible with an asymptotic behavior analysis. As well one can notice a distinguishing behavior of the curves for the smallest values of the step-size where the slope of the curves decreases (lower left corner). This can be explained by the round-off error due to the finite precision (floating-point) arithmetic used on the computer. Round-off error increases in proportion to the total number of arithmetic operations (related to the number of steps here). For the smallest values of the step-size, it grows up to the point of being comparable in magnitude to the error of the method and becomes apparent on the curves.

To conclude this section we would like to emphasis the fact that the simulation results we obtained do not reveal any mathematical contradiction on the convergence order of the S3F-RK4 method: it is 2 and not 4. The apparent contradiction can be understood as follows. The theoretical behavior is proved for a step-size h arbitrarily close from 0 (in other words the property is true asymptoticly when h tends to 0) whereas the experimental curves are obtained certainly with small values of h but inevitably larger than a limit imposed by the performances of the computer. When C is several order of magnitude lower than K, it is not possible experimentally to reach values of h small enough so that K h 4 becomes much smaller than C h 2 . On the contrary, on the range of experimental possible values for h we may have C h 2 much smaller than K h 4 and what it is observed on Figure 2 for P 0 = 50 mW is this latter behavior.

Local error analysis for the NLSE

The above mentioned difficulty in estimating the asymptotic behavior of the local error when solving the GNLSE by the S3F method disappears with the NLSE (where the nonlinear operator is defined in ( 6)) since the nonlinear problem ( 8) is solved analytically and only remains the splitting error. The Symmetric Split-Step method is known to be second order accurate see e.g. [START_REF] Besse | Order estimates in time of splitting methods for the nonlinear Schrödinger equation[END_REF][START_REF] Lubich | On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations[END_REF][START_REF] Zharnitsky | Averaging for Split-Step scheme[END_REF], i.e. the global error behavior at grid point 2 where h = max i=1,...,k h i and C denotes a positive constant. Normalized global error (a.u.)

z k is E k = A k (z k ) -A(z k ) ≤ Ch
Step-size (km) 2 Figure 3: Experimental convergence curve in logarithmic scale (solid line with square symbols) for the S3F method applied for solving the NLSE on the test example. The theoretical convergence order 2 is depicted in dotted line.

In order to illustrate this behavior of the local error we have solved the NLSE by the S3F method, neglecting the Raman induced frequency shift and self-steepening effects (i.e. considering the simplified expression (6) of the nonlinear operator (4)), for the same set of parameters as in section 4.1 when the peak power value P 0 = 50 mW. In Figure [START_REF] Balac | Embedded Runge-Kutta scheme for step-size control in the Interaction Picture method[END_REF] we have depicted the experimental convergence curve in logarithmic scale, i.e. the global error at fiber end versus the step length. We can observe that the experimental convergence curve fits well the theoretical convergence order curve (the dotted line labeled 2 has a slope coefficient equal to 2). Since the global error asymptotic behavior (which takes into account the error accumulation from the fiber entrance) is proportional to the local error by a factor 1/h, this result is in full accordance with relation [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibers[END_REF]. On the contrary to the GNLSE case presented in section 4.1 where the global error behavior is a mixture of the error contributions of the splitting approximation and RK approximation (with different convergence orders), for the NLSE the error is only due to the splitting approximation.

Step-size control

Step-size control for the NLSE solved by the S3F method

For step-size control, a tolerance "tol" is given as bound on the local error estimate. A step-size control strategy [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF] consists in rejecting the current step-size if it gives an estimated local error higher than the specified tolerance and in accepting the solution computed with this step-size otherwise. There are 2 criteria usually employed for step-size control purposes. The criterion of error per step (EPS) selects the step size h k at each step so that the local error is lower than the tolerance tol whereas the criterion of error per unit step (EPUS) selects the step size h k at each step so that the local error is lower than tol × h k . It is clear that for sufficiently small tolerance value EPUS criterion selects a smaller step-size than EPS criterion. When the current step-size is rejected, a new smaller step-size has to be chosen to recompute the solution over the current step. As well, when the current step-size meets the tolerance requirement for the local error it has to be scaled up for the next step computations. In both case, the new step-size has to be estimated using the available information on the previous step computations. On the basis of the results set in the previous section, according to [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibers[END_REF] we assume that there exists a positive real number C such that the local error has the form

L S S k+1 = ℓ S S k+1 = C h 3 k .
The optimal step-size h opt is the one for which the local error estimate L k+1 is the closest to the prescribed tolerance tol, i.e. C h 3 opt = tol. By eliminating the constant C from these 2 relations we obtain that the optimal step-size is

h opt = h k 3 tol L k+1
.

For robustness the step-size control has to be designed in order to respond as smoothly as possible with real or apparent abrupt changes in behavior. This means that the step-size should not vary from one step to the other by an excessive ratio. That is the reason why it is imposed that the new step-size does not exceed twice the current step-size above and half the current step-size below. As well, in order to avoid situations where the specified tolerance is ever exceeded resulting in rejecting too many steps, a safety factor is sometimes introduced: if h opt is the value of the step-size estimated to give a predicted truncation error equal to the tolerance, then the smaller value 0.9h opt for instance is used instead. These considerations lead the following step-size control formula

h new = max        0.5 , min        2.0 , 0.9 3 tol err               h k ( 25
)
where "tol" denotes the tolerance value specified by the user as a bound on the local error and "err" denotes the estimation of the local error for the current step. The 3 constants with values 0.9, 0.5 and 2.0 are somewhat arbitrary and have to be regarded as design parameters.

In order to numerically illustrate the above described stepsize control strategy, we have solved the NLSE for the set of parameters given in section 4.2 by the S3F method with local error estimations based on the Local Error method (LEM) [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF] (i.e. the Step-Doubling method presented in section 3.2) and on the CQE method [START_REF] Heidt | Efficient adaptive step size method for the simulation of supercontinuum generation in optical fibers[END_REF]. For an initial step-size of 10 m, a tolerance of 10 -6 and a number of FFT grid points of 2 14 , the LEM method computes an approximated solution with a relative quadratic error of 1.24 10 -5 using 124 steps. The computation (CPU) time was 7.17s on a AMD A8 Desktop Computer with 8 GO RAM. With the same tolerance value, the CQE method provides a solution with a relative quadratic error of 2.46 10 -6 using 763 steps and a CPU time of 14.7s. Actually, a relative quadratic error similar to the one of the LEM method is obtained with a tolerance of 5 10 -6 . The number of steps is then 343 and the CPU time of 6.77 s.

On the difficulty in designing a step-size control for the GNLSE solved by the S3F method

The difficulty in designing an efficient step-size control strategy for the GNLSE solved by the S3F method is due to the fact that the numerical error is compound of 2 kinds of approximation errors with different asymptotic behavior: the error arising from the Split-Step scheme itself and the error coming from the Runge-Kutta scheme used to solve the nonlinear ODE problem resulting from the Split-Step approach. A step-size control formula such as formula [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF] requires the knowledge of the convergence order of the method but we have shown in section 4.1 that the convergence order of the S3F method applied to the GNLSE is hard to predict. It can vary from 2 to 4 depending on the balance of linear and nonlinear effects. When we known in advance the convergence order p of the method and when we are able to estimate the local error, a step-size control formula similar to formula [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF] with the cube root replaced by the (p + 1) th root (p = 2, 3 or 4) can be used.

For instance, when solving the GNLSE with the set of parameters given in section 4.1 for a peak-power value of 0.5 W the nonlinear effects predominate and it can be shown experimentally that the convergence order of the method is 4. We are then in a position to compare the efficiency of various methods for step-size control on this specific case: Step-Doubling or Local Error (SD/LEM) method, ERK4(3) method and Conservation Quantity Error (CQE) method. For comparison purposes, we have given in Table 1 the CPU time (in s.), the relative quadratic error A ref (L) -A K-1 (L) / A ref (L) at fiber end (where A ref denotes a reference solution computed with a very small constant step-size h = 1 m and A K-1 (L) the approximated solution at fiber end) and the number of steps for the 3 methods.

If we wish a relative quadratic error at fiber end for the ERK4(3) method similar to the one of the SD/LEM method (≈ 7 10 -4 ), a tolerance of 5 10 -10 has to be chosen. CPU time is then 118 s and the number of steps 3773.

A comparison of the evolution of the step-size along the fiber from an initial step-size of 10 m and a tolerance of 10 -6 for the 3 methods is depicted on Figure 4. When compared to the SD/LEM method on this particular test example, we can observe that the ERK4(3) method underestimates the local error and therefore select higher step-sizes than desirable. The reason is that the ERK4(3) method neglects the Splitting error. However the global shape of the step-size curves are similar for the 2 methods (in particular the position of the humps are the same).

We can also observe that on this test example the CQE method provides a step-size curve which is almost constant. It underestimates the local error as well and therefore select higher stepsize than desirable. The reason is that in the CQE method the step-size is adjusted so as to guarantee the conservation of the optical photon number along the propagation and not directly on an estimation of the local error.

3) method and the SD when solving the GNLSE.

Discussion

When implementing a step-size control method, one must be aware of the following issues. First, all the estimates for the local error in the previous sections have been derived for a step-size h k in a neighborhood of 0 and they are used for stepsize control where the goal is to determine step-size as large as possible to reduce the computational cost. It is one of the reason why a good step-size control strategy must include safeguard so that the step-size does not increase in an inconsiderate way. Additionally, local error estimate is based on relation [START_REF] Heidt | Efficient adaptive step size method for the simulation of supercontinuum generation in optical fibers[END_REF] (resp. ( 20)) corresponding to an approximation of the principal part of the Taylor expansion (11) (resp. ( 17)) of the local error for the Runge-Kutta scheme (resp. the Split-Step scheme). For h k small enough, the principal term in the local error expansion is usually large in comparison with the other terms involved in the expansion, which justifies its use to set the step length. However, one must be aware that under special circumstances this can not be anymore the case and the other terms in the expression of the local error expansion can overwhelm the principal term. Here again this is likely to happen when the step-size h k is not small enough.

As we have seen it, implementing a relevant step-size control strategy based on an estimation of the local error for the S3F method when solving the GNLSE (1) is not so easy due to the accumulation of 2 kinds of approximation errors with different asymptotic behavior: the error arising from the Split-Step scheme itself and the error coming from the Runge-Kutta scheme used to solve the nonlinear ODE problem resulting from the Split-Step approach. Recently, an alternative method to the S3F method for solving the GNLSE (and other similar evolution type PDE such as the Gross-Pitaevskii equation) has been propound. The method is termed the "Interaction Picture" (IP) method. In a very similar way to the S3F method, the IP method consists in decoupling the linear and nonlinear terms in the GNLSE in order to make resolution simpler and the approximated solution is identically obtained by solving a sequence of linear PDE problems and nonlinear ODE problem in a given sequential order. However, while in the S3F method the decoupling of the linear and nonlinear terms is obtained at the cost of an approximation (the splitting approximation), the IP method can be considered as exact since it amounts from a mathematical point of view to a change of unknown. The IP method has been developed in the 90's by the Bose-Einstein Condensate Theory Group of R. Ballagh from the Jack Dodd Centre at the University of Otago for solving the Gross-Pitaevskii equation (GPE) in the context of Bose condensation. It was first described in the Ph.D. thesis of B.M. Caradoc-Davies [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] for solving the GPE. Latter it has been applied for solving the GNLSE in optics by J. Hult in [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibers[END_REF] where an experimental comparison of the IP method to Split-Step methods based on their numerical efficiency on benchmark problems in optics is presented. A comprehensive mathematical analysis of the properties of the IP method is presented in [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF]. Since the IP method is based on a change of unknown rather than on a splitting approximation formula, implementing a step-size control strategy based on an estimation of the local error for the GNLSE is straightforward. Indeed, in the IP method the only approximation error results from the use of a Runge-Kutta scheme to solve the nonlinear ODE problem (which is slightly different to the nonlinear ODE problem [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] in the S3F method). In [START_REF] Balac | Embedded Runge-Kutta scheme for step-size control in the Interaction Picture method[END_REF] a costless step-size control strategy for the IP method based on the use of a 4th order embedded Runge-Kutta method, an termed the ERK4(3)-IP method, is propound. Indubitably the IP method offers an elegant and efficient method for solving the GNLSE and its vicinity with the S3F method makes its computer implementation easy from a S3F program. When solving the test example of section 5.2 by the ERK4(3)-IP method with an initial step-size of 10 m, a tolerance of 10 -6 and a number of FFT grid points of 2 14 , we obtain an approximated solution with a relative quadratic error of 3.7 10 -4 using 605 steps and a CPU time of 17 s. Thus the ERK4(3)-IP method provides an accuracy very similar to the one provided by the SD/LEM method, for one-tenth of the CPU time.

Conclusion

We have presented an overview of the mathematical concepts underlying adaptive step-size control techniques for solving the NLSE or GNLSE. We have shown that when solving the NLSE (where the nonlinear problem (8) admits an explicit solution), the Step-Doubling method (also referred as the Local Error method [START_REF] Sinkin | Optimization of the Split-Step Fourier method in modeling optical-fiber communications systems[END_REF]) is the method of choice for adaptive step-size control. One alternative approach would be to use an Embedded Split-Step (ESS) scheme, see [START_REF] Balac | An Embedded Split-Step method for solving the nonlinear Schrodinger equation in optics[END_REF] for a presentation of the ESS method applied to the NLSE and a comparison of the 2 methods. We have also highlighted the difficulty in designing an adaptive step-size control strategy for solving the GNLSE by the S3F method. The reason is that when solving the GNLSE by the S3F method, 2 kinds of approximation errors with different asymptotic behavior are involved: the error arising from the Split-Step scheme itself and the error coming from the Runge-Kutta scheme used to solve the nonlinear ODE problem resulting from the Split-Step approach. This difficulty can be circumvented by using the Interaction Picture (IP) method instead of the Symmetric Split-Step method. The algorithms of the 2 methods are very similar but the IP method amounts from a mathematical point of view to a change of unknown and therefore the local error has only a contribution from the Runge-Kutta approximation scheme.

Figure 1 :

 1 Figure 1: Principle of the Symmetric Split-Step scheme (computational step k is displayed in the dashed box).

Figure 2 :

 2 Figure 2: Experimental convergence curve (logarithmic scale) for the S3F-RK4 method applied for solving the GNLSE on the test example with different peakpower values (solid lines). The convergence order is ranging from 2 to 4 depending on the importance of Kerr nonlinearity (theoretical curves in dotted lines).

Figure 4 :

 4 Figure4: Evolution of the step-size along the fiber for the CQE method, ERK4(3) method and the SD when solving the GNLSE.

Table 1 :

 1 Comparison ofStep-Doubling or Local Error (SD/LEM) method, ERK4(3) method and Conservation Quantity Error (CQE) method for solving the GNLSE by the S3F method with tolerance for the local error set to 10 -6 .

	Method	CPU time (s)	rel. quad. error	nb of steps
	SD/LEM	186 s	6.6 10 -4	2211
	ERK4(3)	17 s	2.3 10 -2	569
	CQE	22 s	5.3 10 -2	850