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

Arnaud Fernandez
email: afernand@laas.fr

SPIP: A computer program implementing the Interaction Picture method for simulation of light-wave propagation in optical fibre

Keywords: generalized non-linear Schrödinger equation, Interaction Picture method, embedded Runge-Kutta method, adaptive step-size control

The computer program spip is aimed at solving the Generalized Non-Linear Schrödinger equation (GNLSE), involved in optics e.g. in the modelling of light-wave propagation in an optical fibre, by the Interaction Picture method, a new efficient alternative method to the Symmetric Split-Step method. In the spip program a dedicated costless adaptive step-size control based on the use of a 4 th order embedded Runge-Kutta method is implemented in order to speed up the resolution.

Introduction

In optics, the non-linear Schrödinger equation occurs for modelling light-wave propagation into a passive optical fibre. This particular form of the Schrödinger equation is obtained from the general set of Maxwell's equations taking advantage of a certain number of assumptions made possible from the very specific characteristics of (quasi-)monochromatic wave propagation in a medium such as a fibre [1,2]. One of the major assumption, referred as the slowly varying envelope approximation, concerns the expression of the electric field in the optical fibre. It assumes that the electric field E is linearly polarized along a direction e x transverse to the direction of propagation e z defined by the fibre and can be represented as a function of time τ and position r = (x, y, z) in a reference frame (O, e x , e y , e z ) as E(r, τ) = A(z, t) F(x, y) e -i(ω 0 τ-kz) e x (1) where the complex valued function A represents the slowly varying optical pulse envelope, F is the electric wave transverse representation, k is the wave number, ω 0 = 2πν 0 is the wave frequency, and t denotes a local time in a moving frame traveling along with the pulse at the group velocity v g = c/n g where n g is the fibre group index and c the speed of light in vacuum. The relationship between the absolute time τ and the local time t is: t = τz/v g . The expression of the electric wave transverse representation F can most of the time be computed explicitly using the classical method of separation of variables for partial differential equations (PDE). For instance, for circular constant transverse section fibres, it is expressed in terms of Bessel functions [1,2].

In the simplest cases of light-wave propagation in an optical fibre, the evolution of the slowly varying pulse envelope A is governed by the following PDE referred in the literature as the Non-Linear Schrödinger Equation (NLSE) [1,2]

∂ ∂z A(z, t) = - α 2 A(z, t) - i 2 β 2 ∂ 2 ∂t 2 A(z, t) + iγA(z, t)|A(z, t)| 2 . (2)
Equation ( 2) describes wave propagation in a single mode fibre taking into account phenomena such as the optical Kerr effect through the non-linear coefficient γ, linear attenuation through the linear attenuation coefficient α and linear dispersion through the chromatic dispersion coefficient β 2 .

In a more accurate model of light-wave propagation in an optical fibre, the evolution of the slowly varying pulse envelope A obeys the so-called Generalized Non-Linear Schrödinger Equation (GNLSE) [1,2]

∂ ∂z A(z, t) = - α 2 A(z, t) +        n max n=2 i n+1 β n n! ∂ n ∂t n A(z, t)        + iγ 1 + i ω 0 ∂ ∂t A(z, t) (1 -f R ) |A(z, t)| 2 (3) + f R +∞ -∞ h R (s) |A(z, t -s)| 2 ds .
In equation (3) linear dispersion is now taken into account through the dispersion coefficients β n , n = 2, . . . , n max whereas non linear dispersion is taken into account 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 non-linear polarization.

Both evolution type PDE (2) and [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] have to be considered together with the initial condition at z = 0 ∀t ∈ R A(0, t) = a 0 (t) [START_REF] Davis | Dynamics in Bose-Einstein condensate[END_REF] where a 0 is a given function and they are solved for all t ∈ R and for all z ∈ [0, L] where L denotes the length of the fibre. The aim of the spip program is to solve the GNLSE (3) and the NLSE (2). For a fixed value of z we denote by A(z) the function t ∈ R → A(z, t) ∈ C. The same notation will be used for functions of 2 variables (z, t) throughout the paper. For convenience, when dealing with the GNLSE (3) we introduce the linear operator

D : A(z) -→ - 1 2 αA(z) + n max n=2 β n i n+1 n! ∂ n t A(z) (5)
and the non-linear operator:

N : A(z) -→ iγ 1 + i ω 0 ∂ ∂t A(z) (1 -f r )|A(z)| 2 + f r (h R * |A(z)| 2 ) ( 6
)
where * stands for the convolution product. With these notations, we are concerned by solving the following PDE problem

         ∂ ∂z A(z) = D A(z) + N(A(z)) ∀z ∈ [0, L] A(0) = a 0 (7)
The NLSE (2) corresponds to n max = 2 in the definition (5) of D and to the non-linear operator

N : A(z) → iγ A(z) |A(z)| 2 (8)
which is formally retrieved from [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF] in the special case when we set f R = 0 and i ω 0 = 0. The spip program solves the PDE problem [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF] with N as given by ( 6) or ( 8) by means of the Interaction Picture (IP) method. The main idea of the IP method is a change of unknown in order to transform the GNLSE for the unknown A into a new equation where only remains an explicit reference to the partial derivation with respect to the space variable z and where the time variable t appears as a parameter. This new equation can be solved numerically using the usual methods for ordinary differential equations (ODE) such as the fourth order Runge-Kutta (RK4) method. Then, by using the inverse transform we obtain approximate values for the unknown A at the grid points of a subdivision of the fibre length interval [0, L]. This numerical approach is referred to as the RK4-IP method. The RK4-IP method has been developed by the Bose-Einstein condensate theory group of R. Ballagh from the Jack Dodd Centre at the University of Otago (New Zealand) in the 90's for solving the Gross-Pitaevskii equation which is ubiquitous in Bose condensation. It was described in the Ph.D. thesis of B. M. Caradoc-Davies [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF] and M. J. Davis [START_REF] Davis | Dynamics in Bose-Einstein condensate[END_REF]. Latter the IP method has been applied for solving the GNLSE in [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibres[END_REF]. A comprehensive mathematical study of the Interaction Picture method for solving the GNLSE has been conducted in [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF]. Recently an efficient embedded RK method based on Dormand and Prince RK4(3)-T formula [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF] and specifically designed for the IP method has been proposed in [START_REF] Balac | Embedded Runge-Kutta scheme for step-size control in the Interaction Picture method[END_REF] to provide a costless adaptive step-size control in the RK4-IP method. It is this method, termed the ERK4(3)-IP method, that is implemented in the spip program.

To be comprehensive we may bring to the attention of the reader that a numerical comparison of the IP method and the Symmetric Split-Step method, the most popular method for solving the GNLSE, has been carried out on benchmark problems in optics in [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibres[END_REF] and a theoretical comparison of the 2 methods has been achieved in [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF]. These investigations demonstrate the supremacy of the IP method over the Symmetric Split-Step method for solving the GNLSE.

The paper is organized as follows. In section 2 we present how the Interaction Picture method is used to solve the GNLSE.

In section 3 we focus on the way the local error is estimated in order to implement the adaptive step-size control. In section 4 we show the algorithm implemented in the spip program to solve the GNLSE and we give information on the way to run the program. Finally, section 5 is devoted to numerical experiments carried out with the spip program.

The Interaction Picture method

Overview of the Interaction Picture method

The interval [0, L] corresponding to the fibre length is divided into K sub-intervals where the spatial grid points are denoted by z k , k = {0, . . . , K} such that ]0

, L] = ∪ K-1 k=0 ]z k , z k+1 ] where 0 = z 0 < z 1 < • • • < z K-1 < z K = L.
The step length between z k and z k+1 is denoted h k and we also set z k+ 1 2 = z k + h k 2 . Solving problem [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF] is equivalent to solving the sequence of connected problems (P k ) k=0,...,K-1 where

(P 0 )          ∂ ∂z A 0 (z) = D A 0 (z) + N(A 0 (z)) ∀z ∈ [z 0 , z 1 ]
A 0 (z 0 ) = a 0 and ∀k ∈ {1, . . . , K -1}

(P k )          ∂ ∂z A k (z) = D A k (z) + N(A k (z)) ∀z ∈ [z k , z k+1 ] A k (z k ) = A k-1 (z k )
Obviously for all k ∈ {0, . . . , K -1} we have

∀z ∈ [z k , z k+1 ] A(z) = A k (z).
Let us consider one of the problems (P k ) for a value of k ∈ {0, . . . , K -1}. Such a problem reads

         ∂ ∂z A k (z) = D A k (z) + N(A k (z)) ∀z ∈ [z k , z k+1 ] A k (z k ) = a k (9)
where a k = A k-1 (z k ) for all k ∈ {1, . . . , K -1}. We introduce as new unknown the mapping

A ip k : (z, t) ∈ [z k , z k+1 ] × R -→ exp(-(z -z k+ 1 2 )D) A k (z, t) (10)
where from a mathematical point of view the exponential term has to be understood in the sense of the continuous group generated by the unbounded linear operator D [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF]. From [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF], one can show [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF] that the new unknown

A ip k is solution to the following problem over each subinterval [z k , z k-1 ] (Q k )            ∂ ∂z A ip k (z) = G k (z, A ip k (z)) ∀z ∈ [z k , z k+1 ] A ip k (z k ) = exp(-(z k -z k+ 1 2 )D) a k (11)
where

G k (z, •) = exp(-(z -z k+ 1 2 )D) • N • exp((z -z k+ 1
2 )D). The major interest in doing such a transformation is that compared to problem [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF], the new problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF] for the unknown A ip k does not involve explicitly partial derivation with respect to the time variable t anymore. Partial derivation with respect to time now occurs through the operator exp(±(zz k+ 1 2 )D). Problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF] can be numerically solved using a standard quadrature scheme for ordinary differential equations such as a Runge-Kutta scheme [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF] as detailed in section 2.2.

Of course, in order to get the initial condition for problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF], one has to compute the mapping exp(-(z kz k+ 1 2 )D) a k . This mapping coincides with the solution for z = z k+ 1 2 to the following linear PDE problem [6]

         ∂ ∂z A k (z) = D A k (z) ∀z ∈ [z k , z k+ 1 2 ] A k (z k ) = a k (12)
Once computed the solution to problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF], one has to use the inverse mapping of (10) to obtain A k (z k+1 ) which is the initial condition for problem (P k+1 ). The function A k (z k+1 ) = exp((z k+1z k+ 1 2 )D) A ip k (z k+1 ) corresponds to the solution for z = z k+1 to the following linear PDE problem [START_REF] Balac | The Interaction Picture method for solving the nonlinear Schrödinger equation in optics[END_REF]

           ∂ ∂z A k (z) = D A k (z) ∀z ∈ [z k+ 1 2 , z k+1 ] A k (z k ) = A ip k (z k+1 ) (13)
The 2 linear PDE problems [START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF] and ( 13) can be numerically solved very efficiently by using the Fourier Integral Transform as shown in section 2.3. In the next section we first study the way of solving the non-linear ODE problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF].

2.2.

Use of an embedded RK scheme for solving problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF] The solution to the non-linear problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF] can be approximated by standard numerical schemes for ordinary differential equations (ODE). Because it represents a good compromise between accuracy and computational cost, we have selected the classical fourth order Runge-Kutta (RK4) scheme. Moreover the values of its internal quadrature nodes (0, 1 2 , 1 2 , 1) imply by symmetry some cancellations in the terms to be computed compared to other 4th order RK schemes when used in conjunction with the IP method, thus reducing the global computational cost of the method. So as to dispose of a costless local error estimator for adaptive step-size purposes, we consider the Embedded Runge-Kutta (ERK) pair schemes of order 3 and 4 defined by the following Butcher extended tableau [START_REF] Balac | Embedded Runge-Kutta scheme for step-size control in the Interaction Picture method[END_REF] where the grey cells correspond to Butcher tableau for the standard RK4 method and the whole array is the Butcher tableau for the RK3 method. This ERK4(3) scheme is referred in literature as the Dormand and Prince Runge-Kutta 4(3) T scheme [START_REF] Dormand | A family of embedded Runge-Kutta formulae[END_REF].

Embedded Runge-Kutta schemes [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 RK schemes designed to deliver two approximations of the solution of the ODE problem under consideration, corresponding to 2 different convergence orders. These 2 approximations of the solution can be considered as an accurate approximate solution (here the one computed with the 4 th order RK scheme) and a coarse approximate solution (the one computed with the 3 rd order RK scheme). Combined in a specific way they deliver an estimation of the local error committed while approaching the solution with the lower order method as detailed in section 3.1.

For

k ∈ {0, . . . , K -1}, let u ip,[3] k (resp. u ip,[4] k
) denote the approximate solution to problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF] resulting from the use of the 3 rd order (resp. 4 th order) RK method. The 4 th order accurate solution is given by u ip, [START_REF] Davis | Dynamics in Bose-Einstein condensate[END_REF] k

(z k+1 ) = A k (z k ) + h k Φ 4 (z k , A k (z k ); h k ) ( 15
)
where the mapping Φ 4 is the increment function [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF] of the RK4 method defined according to Butcher Tableau (14) as

Φ 4 (z k , A k (z k ); h k ) = 1 6 (α 1 + 2α 2 + 2α 3 + α 4 ) with α 1 = G k (z k , A k (z k )) = exp( h k 2 D) N(exp(-h k 2 D) A k (z k )) α 2 = G k (z k + h k 2 , A k (z k ) + h k 2 α 1 ) = N(A k (z k ) + h k 2 α 1 ) α 3 = G k (z k + h k 2 , A k (z k ) + h k 2 α 2 ) = N(A k (z k ) + h k 2 α 2 ) α 4 = G k (z k + h k , A k (z k ) + h k α 3 ) = exp(-h k 2 D) N(exp( h k 2 D) [A k (z k ) + h k α 3 ]
). Similarly the 3 rd order accurate solution is

u ip,[3] k (z k+1 ) = A k (z k ) + h k Φ 3 (z k , A k (z k ); h k ) ( 16
)
where the mapping Φ 3 is the increment function of the RK3 method defined according to Butcher Tableau (14) as

Φ 3 (z k , A k (z k ); h k ) = 1 6 (α 1 + 2α 2 + 2α 3 ) + 1 15 α 4 + 1 10 α 5 with α 5 = G k (z k + h, A k (z k ) + h k 6 (α 1 + 2α 2 + 2α 3 + α 4 )) = G k (z k+1 , u ip,[4] (z k+1 )) = exp(-h k 2 D) N(exp( h 2 D) u [4],ip (z k+1 )).
Using the change of unknown (10) back we deduce the following 2 approximations of the solution A k to problem (P k ) at grid point z k+1 u [4] k+1 = exp(

h k 2 D) A k (z k ) + h k 6 (α 1 + 2α 2 + 2α 3 + α 4 ) u [3] k+1 = exp( h k 2 D) A k (z k ) + h k 30 (5α 1 + 10α 2 + 10α 3 + 2α 4 + 3α 5 ) (17)
where u [3] k+1 (resp. u [4] k+1 ) denotes the approximate solution to problem (P k ) at grid point z k+1 resulting from the use of the 3 rd order (resp. 4 th order) RK method.

Actually we are only interested in computing the approximate solutions to problem [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF] given by ( 17) and the use of the new unknown A ip k is a go-between in the computational approach. Therefore we recast the above computational procedure in the following way in order to reduce the cost of the method. For k ∈ {0, . . . , K -1}, step k of the ERK4(3)-IP method reads: [4] k+1 = β + h k 6 α 4 α 5,k+1 = N(u [4] k+1 ) u [3] k+1 = β + h k 30 (2α 4 + 3α 5,k+1 )

u ip k (z k ) = exp( h k 2 D) A k-1 (z k ) α 1 = exp( h k 2 D) α 5,k α 2 = N(u ip k (z k ) + h k 2 α 1 ) α 3 = N(u ip k (z k ) + h k 2 α 2 ) α 4 = N exp( h k 2 D) [u ip k (z k ) + h k α 3 ] β = exp( h k 2 D) u ip k (z k ) + h k 6 (α 1 + 2α 2 + 2α 3 ) u
We refer to section 2.4 for details in the way the non-linear operator N is computed and to section 3.1 for indications in the way the local error at step k is estimated from the knowledge of u [3] k+1 and u [4] k+1 . It is worth mentioning that compared to the standard RK4 scheme the over-cost of using the ERK4(3) scheme as given by ( 18) is low. In both cases, the number of evaluations of the non-linear operator N is 4 and we have 4 evaluations of the exp( h k 2 D) operator. The extra cost per step is limited to 2 additions and 3 multiplications and the need to keep in memory 2 intermediate results.

Solving the linear PDE problems (12) and (13)

Let us now examine how problems ( 12) and ( 13) are solved or actually how the terms involving the operator exp( h k 2 D) in the computational sequence (18) are obtained. The 2 problems (12) and ( 13) are in the form

         ∂ ∂z U(z) = D U(z) ∀z ∈ [a, b] U(a) = ϕ (19)
where U(z) : t ∈ R → U(z, t) ∈ C and the initial condition ϕ is a known data.

The solution to problem (19) can be computed by using the Fourier Integral Transform (FIT) approach.

For a fixed value of the space variable z, let U(z) be the FIT of U(z) with respect to the time variable t:

∀ν ∈ R U(z, ν) = +∞ -∞ U(z, t) e 2iπνt dt.
From (19) we deduce that U satisfies

         ∂ ∂z U(z) = d ν U(z) ∀z ∈ [a, b] U(a) = ϕ (20)
where

d ν = -α 2 + i N n=2 β n n! (2πν) n .
The solution to the linear first order ODE problem (20) is

U(z) = ϕ e d ν (z-a) . ( 21
)
We deduce that

U(b) = F -1 U(b) = F -1 ϕ e d ν (b-a)
where F -1 denotes the inverse Fourier Transform operator. As a consequence, the solution to problem (12) at grid point z

k+ 1 2 reads A + k (z k+ 1 2 ) = exp( h 2 D) A k-1 (z k ) = F -1 e d ν h k 2 A k-1 (z k )
and the solution to problem (13) at grid point z k+1 reads

A - k (z k+1 ) = exp( h 2 D) A ip k (z k+1 ) = F -1 e d ν h 2 A ip k (z k+1 ) .
In the computational sequence (18) all the terms involving the exp( h k 2 D) operator can be computed in a similar way.

Computation of the non-linear operator

The non-linear operator N defined in ( 6) can be efficiently evaluated by mean of the FIT since in the frequency domain time derivation of functions amounts for multiplying the FIT of the function by a factor -2iπν. Namely for any given function

g : t ∈ R → g(t) ∈ C, we have ∀ν ∈ R N(g)(ν) = iγ 1 + ν ν 0 × F (1 -f R ) g|g| 2 + f R g (h R * |g| 2 ) (ν)
where F denotes the Fourier Transform operator. Moreover, using the properties of the FIT with respect to convolution we may compute the term h R * |g| 2 as follows

h R * |g| 2 = F -1 [ h R × |g| 2 ] ( 22
)
where h R and |g| 2 denote respectively the FIT of h R and |g| 2 . Therefore, computation of N(g) can be achieved in the following steps:

• compute the FIT h R and |g| 2 of h R and |g| 2 ;

• multiply these 2 mappings and compute the inverse FIT of the result to obtain convolution product h R * |g| 2

• compute the FIT of g (1 -f R )|g| 2 + f R (h R * |g| 2 )
• multiply the result by the mapping ν → iγ 1 + ν/ν 0

• compute the inverse FIT of this last product to get N(g).

It is worth mentioning that the evaluation of the non-linear operator N requires 3 Fourier Transform and 2 inverse Fourier Transform evaluations. The spip program uses the FFTW library for computing the Fourier transforms [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. FFTW for "Fastest Fourier Transform in the West" is a software library for computing discrete Fourier transforms (DFT) developed by Matteo Frigo and Steven G. Johnson at the Massachusetts Institute of Technology. FFTW is known as the fastest free software implementation of the Fast Fourier transform (FFT) algorithm. It can compute DFT of real-valued and complex-valued arrays of arbitrary size n with a complexity in O(n log(n)).

The non-linear operator N defined by [START_REF] Balac | Embedded Runge-Kutta scheme for step-size control in the Interaction Picture method[END_REF] for the NLSE can be computed in a similar but simpler way.

3. Adaptive step-size control in the ERK4(3)-IP method

Local error estimate

Since the IP method is based on the change of unknown defined in [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF], the only approximation in the ERK4(3)-IP method is related to the ERK4(3) scheme used to solve the non-linear ODE problem [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. For k ∈ {0, . . . , K-1}, let u [3] k (resp. u [4] k ) denotes the approximate value of the slowly varying optical pulse envelope A at grid point z k delivered by the 3 rd order (resp. the 4 th order) RK scheme. Assuming that the solution value at grid point z k is regarded as exact (because we are concerned by an estimation of the local error), the local error at grid point z k+1 and time t for the 3 rd order (resp. 4 th order) RK schemes used for solving [START_REF] Frigo | The design and implementation of FFTW3[END_REF] 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] [3] k+1 (t) = A(z k+1 , t)u [3] k+1 (t

) = ψ 3 (t, z k , u [3] k ) h 4 k + O(h 5 k ) [4] k+1 (t) = A(z k+1 , t) -u [4] k+1 (t) = ψ 4 (t, z k , u [4] k ) h 5 k + O(h 6 k ) ( 23
)
where ψ 3 (resp. ψ 4 ) are functions of the elementary differentials of order 3 (resp. 4) of G k . By difference of these 2 relations we obtain u [4] k+1 (t)u [3]

k+1 (t) = ψ 3 (t, z k , u [3] k ) h 4 k + O(h 5 k ).
Thus the local error for the 3 rd order RK method at grid point z k+1 can be approximated, with an error in O(h 5 k ), in the following way [START_REF] Caradoc-Davies | Vortex dynamics in Bose-Einstein condensate[END_REF]

k+1 (t) ≈ ψ 3 (t, z k , u [3] k ) h 4 k + O(h 5 k ) ≈ u [4] k+1 (t) -u [3] k+1 (t). ( 24
)
The L 2 -local error at grid point z k+1 is computed as follows

L [3] k+1 = [3] k+1 L 2 (R,C) = [3] k+1 L 2 (R,C) ≈ R u [4] k+1 (ν) -u [3] k+1 (ν) 2 dν 1 2 ≈         h ν J-1 j=0
u [4] k+1 (ν j )u [3] k+1

(ν j ) 2         1 2 (25)
where (ν j ) j=0,...,J denotes a constant step-size h ν sampling of the observed spectral window and the last approximation results from the use of the left rectangle quadrature rule.

Actually, even if the local error estimate (24) holds only for the 3 rd order method, in practice we use the value given by the 4 th order method, which is more accurate, as the approximation of the solution at grid point z k+1 . This approach, referred in the literature [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: nonstiff problems[END_REF] as the local extrapolation mode for ERK methods, slightly overestimates the actual local error, which is safe but actually not fully optimal.

Step-size control

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] consists in rejecting the current step-size if it gives an estimated local error higher than the specified tolerance; 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. In the spip program the EPS criterion is used. 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. When we assume that the leading term in the asymptotic expansion (24) of the local error dominates, then from ( 24) and (25) there exists C > 0 such that L [3] k+1 = [3] k+1 L 2 (R,C) = C h 4 k . The optimal step-size h opt is the one for which the local error estimate L [3] k+1 is the closest to the prescribed tolerance tol, i.e. C h 4 opt = tol. By eliminating the constant C from these 2 relations we obtain

h opt = h k 4 tol
L [3] k+1 .

In the previous relationship it is common to consider an estimate for the relative local error and a relative tolerance "tol". It is the choice we adopt in the spip program.

For robustness the step-size control has to be designed in order to respond as smoothly as possible to real or apparent abrupt changes in behaviour. This means that the step-size should not vary from one step to the next by an excessive ratio. That is the reason why we impose that the new step-size does not exceed twice the current step-size above and half the current stepsize below. Moreover, in order to avoid situations where the specified tolerance is ever exceeded resulting in rejecting too many steps, a safety factor is introduced. Namely, 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.

Following these requirements, we use the following step-size control formula in the spip program

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

4. The spip program 4.1. Algorithm for the ERK4(3)-IP method with step-size control

Since the evaluation of the non-linear operator N and the exp( h k 2 D) operator both involve the use of Fourier Transforms, it is convenient to slightly modify the computational sequence (18) in order to save some Fourier Transforms by optimizing the occurrence of successive FT and inverse FT. This leads to the following algorithm for the ERK4(3)-IP method.

ERK4(3)-IP algorithm

Input: -Array [t j ] j=1,...,J contains the time sampling points -Array u contains the incident slowly varying pulse envelope sampled in time -Array [ν j ] j=1,...,J contains the frequency sampling points -Real number α : linear attenuation coefficient -Array [β n ] n=2,...,n max contains the dispersion coefficients -Real number tol: the tolerance value for the local error -Real number h initial : the initial step size value Output: Array u contains the outgoing slowly varying pulse envelope at fibre end sampled in time 1: {Initialisation} 2: for j = 1, . . . , J do 3:

d[ j] ← -1 2 α + i n max n=2 β n
n! (2πν j ) n 4: end for 5: u ← FFT(u, f orward) 6: Nu ← COMPUTE FTN(u, u) 7: z k = 0, h = h initial 8: {Loop over the propagation sub-interval} 9: while z k < L do 10:

for j = 1, . . . , J do 11: e[ j] ← exp( h 2 d[ j]) 12: u ip [ j] ← e[ j] × u[ j] 13:
end for 14:

for j = 1, . . . , J do 15:

α 1 [ j] ← e[ j] × Nu[ j] 16: u 2 [ j] ← u ip [ j] + h 2 α 1 [ j] 17:
end for 18:

u 2 ← FFT( u 2 , backward) 19: α 2 ← COMPUTE FTN(u 2 , u 2 )
20:

for j = 1, . . . , J do 21:

u 3 [ j] ← u ip [ j] + h 2 α 2 [ j] 22:
end for 23:

u 3 ← FFT( u 3 , backward) 24: α 3 ← COMPUTE FTN(u 3 , u 3 )
25:

for j = 1, . . . , J do 26: for j = 1, . . . , J do 31:

u 4 [ j] ← e[ j] × ( u ip [ j] + h α 3 [ j])
r[ j] ← e[ j] × ( u ip [ j] + h 6 α 1 [ j] + h 3 α 2 [ j] + h 3 α 3 [ j]) 32: u 1 [ j] ← r[ j] + h 6 α 4 [ j] {RK4 solution} 33:
end for 34:

u 1 ← FFT( u 1 , backward) 35: α 5 ← COMPUTE FTN(u 1 , u 1 )
36:

for j = 1, . . . , J do 37:

v 1 [ j] ← r[ j] + h 30 (2 α 4 [ j] + 3 α 5 [ j]) {RK3 solution} 38:
end for 39:

{Step-size control} 40:

err ← 0 {L 2 error} 41:

for j = 1, . . . , J -1 do 42: err ← err + | u 1 [ j] -v 1 [ j]| 2 43:
end for 44:

err ← √ h t err {see Formula (25)} 45:

h opt = max 0.5 , min 2.0 , 0.9 4 tol err h {Optimal stepsize for the given prescribed tolerance, see Formula (26)} 46:

if err ≤ tol then end if 57: end while FFT(u, f orward) stands for a call to the Fast Fourier Transform (FFT) algorithm to compute the Discrete Fourier Transform (DFT) of array u, whereas FFT(u, backward) stands for a call to FFT algorithm to compute the inverse DFT of array u. COMPUTE FTN refers to the following function.

FUNCTION g = COMPUTE FTN( f , f ) Compute the Fourier Transform of g : t → N( f )(z, t)
for a given z where N refers to the non-linear operator defined by ( 6) for the GNLSE and by [START_REF] Balac | Embedded Runge-Kutta scheme for step-size control in the Interaction Picture method[END_REF] for j = 1, . . . , J do 4:

op 1 [ j] ← | f [ j]| 2 5:
end for 6:

op 1 ← FFT(op 1 , f orward) 7:
for j = 1, . . . , J do 8:

op 2 [ j] ← op 1 [ j] × h R [ j] 9:
end for 10:

op 2 ← FFT( op 2 , backward) {Array op 2 contains the convolution product h R * | f | 2 } 11:
for j = 1, . . . , N do 12:

op 3 [ j] ← f [ j] × (1 -f R ) op 1 [ j] + f R op 2 [ j] 13:
end for 14:

op 3 ← FFT(op 3 , f orward) 15:
for j = 1, . . . , J do 16: for j = 1, . . . , J do 21:

g[ j] ← iγ(1 + 2πν j ω 0 ) op 3 [ j]
op 1 [ j] ← f [ j] × | f [ j]| 2 22: end for 23: g ← FFT(op 1 , f orward) 24: end if
When the solution is required at every grid points it suffices to modify line 50 of the ERK4(3)-IP algorithm to record or draw the solution u.

As mentioned earlier, the main cost of the algorithm lies in the number of evaluations of the non-linear operator N which is 4 per step and to a lower degree in the 4 evaluations of the exp( h k 2 D) operator. In total, it amounts for 16 FFT per step.

Installation instructions

The spip program requires 2 external librairies: It is also possible to record the time evolution of the solution along the fibre by setting the variable REC to the value 1 in the spip.h header file. Note that this option is not set by default since it increases the computational time and produces a very large file on the disc (e.g. around 78 Mo for the present simulation).

• FFTW a C subroutine
A(L) 2 = R |A(L, t)| 2 dt 1 2 , A(L) 1 = R |A(L, t)| dt, A(L) ∞ = sup t∈R |A(L, t)|. (27)
At the end of the simulation, the directory specified to store the results files contains the following files (xxx stands for the key label of the current execution):

• data_xxx.txt: this ASCII file contains the values of the physical and numerical parameters used for the current simulation.

• pulse_in_xxx.txt: this ASCII file contains the incident slowly varying optical pulse envelope sampled over the time window.

• pulse_out_xxx.txt: this ASCII file contains the slowly varying optical pulse envelope at the fibre end sampled over the time window.

• zstep_xxx.txt: this ASCII file contains the step-size values along the propagation direction computed by the adaptive step-size method.

When the REC variable in the spip.h header file has been set to 1 to record the solution at every computational step along the fibre, the directory contains and additional binary file named pulse_all_xxx.spip.

Matlab/Octave tools

Together with the spip program are provided several Matlab/Octave scripts to handle spip results files. File plotsol.m contains a script to draw under Matlab/Octave the slowly varying optical pulse envelope A, from the results files pulse_in_ xxx.txt and pulse_out_xxx.txt generated by the spip program. It also plots the Fourier Transform of the slowly varying optical pulse envelope A. The script allows to plot again the solution when the Gnuplot graphical windows of the spip program have been closed. It also allows image manipulations thanks to the Matlab/Octave graphics interface. The file plotspip.m contains a script to draw the solution at every computational step along the fibre when the REC variable in the spip.h header file has been set to 1 and a binary file named pulse_all_xxx.spip has been created. Namely, the script draws the slowly varying pulse envelope power |A| 2 as a function of time and space (two different views are propound) and the power spectral density in dBm as a function of wavelength and position along the fibre.

The file plotstep.m can be used to draw the variation of the step-size along the fibre length resulting from the adaptive step-size strategy used in the spip program.

Numerical experiments

We present in this section some numerical experiments carried out with the spip program. All the simulations presented here were achieved under Linux Ubuntu 14.04 LTS on a desktop computer equipped with an Intel Core i5-4200M processor and 8 GO RAM. The figures were obtained under Matlab with the various scripts provided with the spip program.

Soliton solution to the NLSE

As a first test example, let us consider the case of the NLSE (2) where α = 0, f R = 0, n max = 2. When β 2 < 0, there exists an exact solution to the NLSE known as the optical Soliton [1]. Namely, if the source term is given by

a 0 : t → N s √ γL D cosh(t/T 0 ) (28)
where N s is the soliton order, T 0 is the pulse half-width and L D = -T 2 0 /β 2 is the dispersion length, then the solution to the NLSE at the soliton period z p = π 2 L D is given by

∀t ∈ R A(z p , t) = N s √ γL D e iz p /(2L D ) cosh(t/T 0 ) . (29)
We consider a 3 rd order soliton (N s = 3) and the physical parameters of the numerical experiment are λ = 1550 nm, T 0 = 5.673 ps, L = 8 km, γ = 4.3 W -1 km -1 , β 2 = -19.83 ps 2 km -1 (SMF-28 fibre). For the simulation, the number of sampling points for the FFT computations was set to 2 12 and the time window to 200 ps. For adaptive step-size control purposes, the tolerance was set to 10 -6 and the initial step-size was 10 m. The number of discretisation steps along the fibre was found to be 1340 and the computation time was 2.04 s (this time doesn't include the time devoted to the recording of the solution at every step-size along the fibre). We have depicted in Fig. 3 the time and spectral variations of the slowly varying pulse envelope A along the fibre length; namely we have drawn the power |A| 2 as a function of the time and space variables and the power spectral density expressed in dBm (defined as 10 log 10 ( |A|)) as a function of the wavelength and position along the fibre. The figures were obtained using the Matlab/Octave tools provided with the spip program. Since the exact solution is known, the computational error has been measured. After one period (z = z p = 2.5493 km) the relative global error evaluated with the L 2 -norm, defined in (27), is 7.77 10 -5 whereas the maximum relative error, defined in (27), is 1.19 10 -4 . After 3 periods (z = 3z p = 7.6479 km), the quadratic relative global error is 8.01 10 -4 whereas the maximum relative error is 1.42 10 -3 .

The variation of the step-size resulting from the adaptive step-size strategy used in the spip program is given at Fig. 4. Since the step-size is adapted so that the local error matches the user defined tolerance at each computational step, small stepsizes indicate an area in the fibre where the slowly varying pulse envelope A varies a lot whereas larger step-sizes indicate areas where the solution varies in a very smooth way.

We have also tested the propagation of a 10 th order Soliton in the same fibre. For this experiment, the number of sampling points for the FFT computations was set to 2 18 and the time window to 200 ps. The tolerance was set to 10 -9 and the initial step-size was 0.1 m. After one period (z = zp = 2.5493 km), the quadratic relative error is 4.64 10 -4 whereas the maximum relative error is 6.33 10 -6 . The number of steps was 33968 and the simulation time 41 mn. After 3 periods (z = 3z p = 7.6479 km), the quadratic relative global error is 1.38 10 -3 whereas the maximum relative error is 1.81 10 -3 . The number of steps was 101192. We have depicted in Fig. 5 the time variation along the fibre of the slowly varying pulse envelope power |A| 2 (top figure) and the power spectral density expressed in dBm as a function of the wavelength along the fibre (bottom figure) for the 10 th order Soliton over one period. We have also depicted in Fig. 6 the variation of the step-size resulting from the adaptive step-size strategy for the 10 th order Soliton over one period.

Solving the GNLSE (I)

We consider the case of the GNLSE (3) with the following set of physical parameters : λ = 1064 nm, γ = 4.3 W -1 km -1 , n max = 3, β 2 = 19.83 ps 2 km -1 , β 3 = 0.031 ps 3 km -1 , α = 0.046 km -1 , L = 96, 77 m, f R = 0.245. The Hyperbolic Secant pulse at the fibre entrance (z = 0) is expressed as

∀t ∈ R a 0 (t) = P 0 exp(-1 2 iC t 2 /T 2 0 ) cosh(t/T 0 ) ( 30
)
where T 0 = 2.8365 ps is the pulse half-width, C = 1 is the chirp parameter and P 0 = 100 W is the pulse peak power. The number of sampling points for the FFT computations was set to 2 14 and the time window width to 100 ps. The tolerance for the adaptive step-size control was set to 10 -6 and the initial step-size was 1 m. The number of discretisation steps along the fibre was found to be 301 and the computation time was 3.72 s.

We have depicted in Fig. 7 the time variation along the fibre of the slowly varying pulse envelope power |A| 2 (top figure) and the power spectral density expressed in dBm as a function of the wavelength along the fibre (bottom figure).

Solving the GNLSE (II)

We use the spip program to solve the GNLSE 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 infra-red spectrum (λ = 1550 nm). The following set of fibres parameters were used for the simulation: α = 0.046 km -1 , γ = 4.3 W -1 km -1 , f R = 0.245, n max = 3, β 2 = -19.83 ps 2 km -1 , β 3 = 0.031 ps 3 km -1 . 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 where T 0 = 6.8 ps is the pulse half-width at 1/e intensity point and P 0 = 25 mW is the pulse peak power.

The number of sampling points for the FFT computations was set to 2 14 and the time window size to 500 ps. The tolerance for the adaptive step-size control was set to 10 -9 and the initial step-size was 0.1 m. The number of discretisation steps along the fibre was found to be 99 and the computation time was 1.87 s. We have depicted in Fig. 8 the time variation along the fibre of the slowly varying pulse envelope power |A| 2 (top figure) and the power spectral density expressed in dBm as a function of the wavelength along the fibre (bottom figure).

Soliton collisions

We now present numerical simulation results for the collision of 2 first order Solitons [1]. It is known that when two neighbouring Solitons are launched with the same phase, they are initially attracted towards each other and then the two pulses periodically coalesce to form one pulse and separate [START_REF] Desem | Reducing soliton interaction in single-mode optical fibres[END_REF]. The source term was

a 0 : t ∈ R → 1 √ γL D 1 cosh((t -T 1 )/T 0 ) + Re iφ cosh(R(t + T 1 )/T 0 )
where T 0 is the pulse half-width, L D = -T 2 0 /β 2 is the dispersion length, R accounts for the relative amplitude, φ for the relative phase shift and T 1 for the initial separation time. To solve the NLSE, the following physical parameters were taken for the numerical experiment: L = 5000 km, λ = 1550 nm, γ = 2.2 W -1 km -1 , β 2 = -0.1 ps 2 km -1 , T 0 = 4 ps, T 1 = 100 ps, R = 1 and φ = 0. The time windows was 400 ps wide and the number of FFT nodes was 2 14 . The initial step-size was set to 1 km and the tolerance to 10 -6 . The number of discretisation steps along the fibre was found to be 480 and the computation time was 3.80 s.

We have depicted in Fig. 9 the evolution of the lowly varying pulse envelope power |A| 2 as a function of time and position along the fibre and the evolution of the power spectral density in dBm as a function of wavelength and position along the fibre. We can easily identify in the two figures the position of the Solitons collision. With the values considered for this simulation, the collision of the two Solitons is predicted to happen at a distance of 4161 km [1]. This is confirmed by the plot of the variation of the step-size resulting from the adaptive step-size strategy, see Fig. 10. When the two Solitons catch up, the numerical solution varies a lot over a small distance and therefore the step-size decreases so that the local error meets the prescribed tolerance. The minimal step-size value in the area of collision is found to occur at a distance of 4165 km.

Super-continuum generation

Finally, we consider a typical problem for super-continuum generation in Photonic Crystal Fiber (PCF) employing femtosecond lasers and pumping in the anomalous dispersion regime [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibres[END_REF]. The initial pulse was a 6 th order Soliton as defined in (28) with T 0 = 28.4 fs, λ = 850 nm and the fibre parameters were: L = 10 cm, α = 0.046 km -1 , f R = 0.18, γ = 45 W -1 km -1 , β 2 = -12.76 ps 2 km -1 , β 3 = 8.119 10 -2 ps 3 km -1 , β 4 = -1.321 10 -4 ps 4 km -1 , β 5 = 3.032 10 -7 ps 5 km -1 , β 6 = -4.196 10 -10 ps 6 km -1 and β 7 = 2.57 10 -13 ps 7 km -1 . The simulation parameters were the following: time windows width 6 ps, number of FFT nodes 2 14 , initial step-size 10 -4 m and tolerance 10 -6 . The computation time was 14.9 s and the number of steps in the adaptive stepsize strategy was found to be 1475. We have depicted in Fig. 11 the evolution of the lowly varying pulse envelope power |A| 2 as a function of time and position along the fibre and the evolution of the power spectral density in dBm as a function of wavelength and position along the fibre. We also provide in Fig. 12 the variation of the step-size along the fibre. On can observe how the step-size is adapted in the area where the spectral broadening caused by the interaction between self-phase modulation and group velocity dispersion takes place in the first centimetre of the fibre [START_REF] Hult | A fourth-order Runge-Kutta in the Interaction Picture method for simulating supercontinuum generation in optical fibres[END_REF].

Conclusion

We have presented the open-source software spip that solves the Generalized Non-Linear Schrödinger Equation (GNLSE) as well as the Non-Linear Schrödinger Equation (NLSE) by the Interaction Picture method combined with an adaptive stepsize control strategy based on the use of a dedicated Embedded Runge-Kutta method. We have shown the results provided by the spip program on typical problems involving the GNLSE/NLSE in optics.

  47:{the current local error matches the tolerance} 48:z k = z k + h {New grid point is confirmed} 49: h = min(h opt , Lz k ) {New step-size value} 50: u ← u 1 {Array u contains the time sampled values [A k (z k , t j )]j=1,...,J of the pulse envelope at grid point z k } {The current local error does not match the tolerance} 55: h = h opt {New computation from z k with smaller stepsize h opt is necessary} 56:

Dispersion length = 4 .

 4 057354e-01 km Non-linear length = 2.325581e-03 km Number of spatial steps = 300 L^2-norm of the solution at fibre end = 2

Figure 1 :

 1 Figure 1: Screen-shot of the Gnuplot graphics windows where are depicted the time variation of the modulus, real and imaginary parts of the slowly varying pulse envelope A(0, t) at the fibre entrance (top figure) and the Fourier Transform of these quantities (bottom figure, zoom on the area of interest).

Figure 2 :

 2 Figure 2: Screen-shot of the Gnuplot graphics windows where are depicted the time variation of the modulus, real and imaginary parts of the slowly varying pulse envelope A(L, t) at the fibre end (top figure) and the Fourier Transform of these quantity (bottom figure, zoom on the area of interest).

Figure 3 :

 3 Figure 3: Slowly varying pulse envelope power |A| 2 in W (top figure) as a function of time and position along the fibre and power spectral density in dBm (bottom figure) as a function of wavelength and position along the fibre for the propagation of a 3rd order Soliton along a 8 km long fibre.

Figure 4 :

 4 Figure 4: Variation of the step-size resulting from the adaptive step-size strategy used in the spip program for the 3 rd order Soliton.

Figure 5 :

 5 Figure 5: Slowly varying pulse envelope power |A| 2 in W (top figure) as a function of time and position along the fibre and power spectral density in dBm (bottom figure) as a function of wavelength and position along the fibre for the propagation of a 10 th order Soliton along one period.

Figure 6 :

 6 Figure 6: Variation of the step-size resulting from the adaptive step-size strategy used in the spip program for the 10 th order Soliton over one period.

Figure 7 :

 7 Figure 7: Slowly varying pulse envelope power |A| 2 in W (top figure) as a function of time and position along the fibre and power spectral density in dBm (bottom figure) as a function of wavelength and position along the fibre for the propagation of a Hyperbolic Secant pulse along a 96.77 m long fibre.

Figure 8 :

 8 Figure 8: Slowly varying pulse envelope power |A| 2 in W (top figure) as a function of time and position along the fibre and power spectral density in dBm (bottom figure) as a function of wavelength and position along the fibre for the propagation of a Gaussian pulse along a 20 km long fibre.

Figure 9 :

 9 Figure 9: Slowly varying pulse envelope power |A| 2 in W (top figure) as a function of time and position along the fibre and power spectral density in dBm (bottom figure) as a function of wavelength and position along the fibre for the propagation of two neighbouring Solitons in a 5000 km long fibre.

Figure 10 :

 10 Figure 10: Variation of the step-size resulting from the adaptive step-size strategy (top figure) and zoom in the area of collision (bottom figure).

Figure 11 :

 11 Figure 11: Slowly varying pulse envelope power |A| 2 in W (top figure) as a function of time and position along the fibre and power spectral density in dBm (bottom figure) as a function of wavelength and position along the fibre for the propagation of a 6th order Soliton pulse into a 10 cm PCF.

Figure 12 :

 12 Figure 12: Variation of the step-size resulting from the adaptive step-size strategy for the super-continuum generation.

  Display requests to a Gnuplot session from spip is achieved thanks to the gnuplot_i module developed by Nicolas Devillard, see http://ndevilla.free.fr/gnuplot/.

	4.2.1. Under Linux
	Unzip the spip archive spip-1.1.zip in a suitable location
	of the home directory by typing in a shell window:
	$ unzip spip-1.1.zip
	Go to the Linux directory:
	$ cd spip-1.1/linux
	Compile the C program thanks to the Makefile available in the
	directory:

library for computing the discrete Fourier transform (DFT) developed by Matteo Frigo and Steven G. Johnson, see http://www.fftw.org • Gnuplot a portable command-line driven graphing utility, see http://www.gnuplot.info.

Acknowledgment

The present work has been undertaken under the framework of the Green Laser Project supported by Région Bretagne, France. The authors would like to thank F. Mahé, F. Méhat and R. Texier-Picard from the Research Institute of Mathematics in Rennes (IRMAR CNRS UMR 6625), France, for their involvement in the mathematical study of the IP method as well as T. Chartier from Foton laboratory in Lannion (CNRS UMR 6082), France, for his support.

$ make

The executable is named spipxx. Optionally, you can delete object files by typing: $ make clean

Under Windows

You can use the Unix-like environment and command-line interface Cygwin for Microsoft Windows (http://cygwin. com) and proceed pretty much as for the Linux installation. Alternatively, you can use an C language IDE (e.g. Code::Blocks, http://www.codeblocks.org) to compile the sources. A Code::Blocks project spip_windows.cbp is available in the archive under the directory windows for convenience.

Note that display requests to a Gnuplot session from spip is not available under Windows. It is however possible to display the slowly varying optical pulse envelope using the Matlab/Octave scripts provided in the archive, see section 4.4 for details.

Running the program

Under Linux, the spip program is launched by typing in a shell window:

Under Windows, double-click on the spipxx program icon located in the Windows directory of the spip archive to launch the program.

In the command line user interface, the user is first asked to provide the name and path of the folder where the result files will be stored.

[?] Directory where result files will be stored: res Typing a carriage return key will set this directory to the current directory by default. If the folder doesn't exist, the program will end to let the user create it.

Then the user is asked to provide a keyword for the present simulation. The result files generated by the spip program will automatically bear this keyword in their name with an additional unique key corresponding to the elapsed time in second from Epoch.

[?] Keyword for the present simulation: green Then the type of equation to solve must be specified: Non-Linear Schrödinger Equation or Generalized Non-Linear Schrödinger Equation.

[?] Equation type [0 for NLS or 1 for GNLS] = 1 Four pre-defined incident shapes for the slowly varying electric field pulse envelope are available: Soliton, Gaussian and Super-Gaussian, Hyperbolic Secant, see [1]. The mathematical expression for these pulses are given in the spip documentation available in the spip archive. Additionally, it is possible to provide any other incident pulse shape in the program. To do so, one must append the expression of the incident pulse in file inpulse shape.c in the space provided and one must recompile the program.

[?] Shape of the incident pulse envelope [1: Soliton, 2: Gaussian,3: Hyperbolic Secant, 4: Super-Gaussian, 5: user-defined] = 2

The next stage consists in providing the incident pulse envelope features. The items may vary depending on the selected shape. For the Gaussian pulse they are : the wavelength of the pulse, its Peak Power, the value of the chirp constant and the half-width at 1/e-intensity point of the pulse.