NUMERICAL METHODS FOR A CLASS OF GENERALIZED NONLINEAR SCHR¨ODINGER EQUATIONS

. We present and analyze diﬀerent splitting algorithms for numerical solution of the both classical and generalized nonlinear Schr¨odinger equations describing propagation of wave packets with special emphasis on applications to nonlinear ﬁber-optics. The considered generalizations take into account the higher-order corrections of the linear diﬀerential dispersion operator as well as the saturation of nonlinearity and the self-steepening of the ﬁeld envelope function. For stabilization of the pseudo-spectral splitting schemes for generalized Schr¨odinger equations a regularization based on the approximation of the derivatives by the low number of Fourier modes is proposed. To illustrate the theoretically predicted performance of these schemes several numerical ex- periments have been done. In particular, we compute real-world examples of extreme pulses propagating in silica ﬁbers.


(Communicated by Zhouping Xin)
Abstract. We present and analyze different splitting algorithms for numerical solution of the both classical and generalized nonlinear Schrödinger equations describing propagation of wave packets with special emphasis on applications to nonlinear fiber-optics. The considered generalizations take into account the higher-order corrections of the linear differential dispersion operator as well as the saturation of nonlinearity and the self-steepening of the field envelope function. For stabilization of the pseudo-spectral splitting schemes for generalized Schrödinger equations a regularization based on the approximation of the derivatives by the low number of Fourier modes is proposed. To illustrate the theoretically predicted performance of these schemes several numerical experiments have been done. In particular, we compute real-world examples of extreme pulses propagating in silica fibers.

Introduction. Nonlinear Schrödinger equation (NLSE)
is widely used to model pulse propagation in nonlinear dispersive media, e.g., along a micro-structured optical fiber [2,22,24]. Optical properties (dispersion) of the fiber are described by a relation between wave vector and frequency, β = β(ω), for a monochromatic linear wave in which the wave field is proportional to e i(βz−ωt) . The propagation coordinate z is directed along the fiber. The dispersion relation is typically encoded in a sequence of the so-called propagation constants β j = d j β/dω j , j = 0, 1, . . ., which are calculated or measured for some reference wave frequency ω 0 . It is convenient to assume that ω 0 equals the pulse carrier frequency. In particular, β 0 is the carrier wave vector and 1/β 1 describes the group velocity at the carrier frequency.
NLSE yields evolution of the complex pulse envelope, ψ(z, τ ), along the fiber. Here the so-called retarded time τ = t − β 1 z corresponds to a coordinate system moving together with the pulse at the group velocity 1/β 1 . The latter generally 216 SHALVA AMIRANASHVILI, RAIMONDASČIEGIS AND MINDAUGAS RADZIUNAS differs from the phase velocity of the carrier wave ω 0 /β 0 . The corresponding electric field E is related to the field envelope ψ by a rapidly oscillating carrier wave e i(β0z−ω0t) : E(z, t) = Re ψ(z, t − β 1 z)e i(β0z−ω0t) .
An important precondition for the NLSE description is that the complex envelope ψ evolves slowly as compared to the carrier wave, such that |∂ τ ψ| |ω 0 ψ|. This slowly varying envelope approximation (SVEA) is used to derive a class of the NLSEs for ψ(z, τ ).
The simplest, classical NLSE is given by the relation [3] i ∂ψ ∂z − β 2 2 where z ∈ [0, L] is the propagation coordinate, L indicates the total length of the fiber, and the retarded time τ plays role of the lateral coordinate of the considered computational domain. The nonlinear term N [ψ] represents an instantaneous cubic (Kerr) nonlinearity, β 2 is referred to as the group-velocity dispersion parameter, and γ characterizes nonlinear properties of the fiber. The initial conditions ψ(z = 0, τ ) = ψ 0 (τ ), (3) are defined by the incoming field envelope ψ 0 , which can be expressed from the corresponding electric field E(0, t) using the SVEA and the so-called sliding average [39]. The lateral boundary conditions are determined either by the natural assumption that field envelope represents a pulse which is localized w.r.t. τ , or by considering a periodic sequence of pulses, such that ψ(z, τ + T ) ≡ ψ(z, τ ), where 1/T is the pulse repetition rate. The cubic NLSE (2) is just one representative of a huge family of the generalized nonlinear Schrödinger equations (GNLSE). The generalizations can be made by using more involved dispersion operator D and/or more general nonlinearity operator N i ∂ψ ∂z + Dψ + N [ψ] = 0.
For example the popular approximation [22], takes into account higher order propagation constants β j , j = 3, . . . , M , which can be important once operating close to the zero dispersion frequency (ZDF) at which β 2 = β (ω 0 ) vanishes, or modeling ultra-short pulses determined on a large frequency range. Note, that in practical situations the dispersion relation β(ω) and the dispersion operator D are known/measured only in an application-relevant interval of frequencies, see e.g. the handbook [14]. In the same way, the Taylor expansion, which is hidden behind the definition (7), has some convergence radius again indicating a finite interval in the frequency domain. In such a case, the operator D should be filtered in the frequency domain. Examples of such problems will be considered in what follows.
Another type of generalizations N = N + Q occurs by assuming a more general form of local nonlinearities N [ψ], as, for example, and nonlocal nonlinearities Q[ψ], as, for example, Here, a real function f can describe a simple saturating nonlinearity e.g., An operator iσ ∂ ∂τ represents a self-steepening of the field envelope function and accounts for non-perfectness of the envelope description, whereas a convolution with the quickly decreasing real function h(τ ) s.t.
∞ 0 h(τ )dτ = 1 represents a delayed medium response Raman effect. Below in this paper we ignore a delayed medium response and consider the GNLSE (6) with the nonlinearity It is worth mentioning that the simplest Eq. (2) possess a very special property: it is completely integrable by means of the inverse-scattering method [41]. Usually this property is destroyed by generalizations, but in some cases the integrability is preserved. For example, the normalized GNLSE preserves integrability for special combinations of coefficients, such as c 1 = 0, c 2 = 1 and either c 3 = 1 or c 3 = 0 (different derivative nonlinear Schrödinger equations), c 1 = 1, c 2 = 6, and c 3 = 0 (Hirota equation), or c 1 = 1, c 2 = 6, and c 3 = 3 (Sasa-Satsuma equation). For a real-valued set of β m , localized in τ pulses, and smooth enough initial function ψ 0 (τ ) the initial value problem (3), (4), (6), (12), can possesses some conservation laws. The classical cubic NLSE (2) has an infinite set of conservation laws [41]. For the general Eq. (6) with the vanishing Q[ψ] in Eq. (12), at least three conservation laws can be found. One can easily derive them by writing (6) as a Hamiltonian equation i ∂ψ ∂z + δH δψ * = 0, where ψ * denotes a complex conjugate of ψ, functional H is defined below in (16), and δH δψ * is a variational derivative. The first two integrals are derived from continuous symmetries of the Hamiltonian H, namely and

SHALVA AMIRANASHVILI, RAIMONDASČIEGIS AND MINDAUGAS RADZIUNAS
Finally, the Hamiltonian as such yields another conservation law, where F (ξ) = ξ 0 f (s) ds. The same conservation laws (14), (15), and (16) with integration limits 0 and T are valid in the case of T -periodic boundary conditions (5). Generalizations of these integrals can be derived also for non-vanishing Q[ψ], see Ref. [8,28,31] and Section 3.2 of this paper.
2. Numerical approximation of GNLSE. There are many numerical algorithms for the solution of nonlinear Schrödinger type problem (for various types of local nonlinearity operators) . They are based on finite-difference schemes (see, e.g., [10,17,18,19,33,37]), finite-element and Galerkin approaches (see, [6,10]), or spectral and pseudo-spectral methods (see the well-known pioneer paper [23], then these methods where developed in many papers [12,15,31,38], see also references given therein). Various splitting techniques are used to construct efficient integrators of the obtained semi-discrete problems [27,29,38]. Recently, very good reviews on numerical methods for simulation of various types of solitons in the nonlinear Schrödinger equation are presented in [11,13]. In our paper we are following the classification and notation of numerical methods proposed in these papers. Our main goal is to select appropriate numerical solvers and to compare their efficiency when they used to solve the generalized NLSEs. with the grid points denoted by τ j , the grid boundary ∂Ω h = {τ 0 = a, τ J = b}, and the inner part of the grid Ω h =Ω h \∂Ω h . The constant step h is assumed in order to apply the same grid for both finite difference and spectral methods. In the case of finite difference approximations non-uniform grids can be used and approximations can be adapted to the behavior of solutions of GNLSE (see, e.g. [32], where highorder finite difference schemes are constructed on non-uniform grids).
Let ω k be a "time" grid where k n is the discretization step. For simplicity of notation, in most definitions the step-size k n is taken constant but in computational experiments we use adaptive strategies to define k n . We consider numerical approximations U n j to the exact solution values ψ n j = ψ(z n , τ j ) at the grid points (z n , τ j ) ∈ ω k ×Ω h . The following notation for difference and averaging operators is used All finite difference numerical schemes for the GNLSE (6) without higher order dispersion constants β j , j > 2, and with initial condition (3) presented below can be written as where the selection of the grid Ω h ⊂Ω h depends on the assumed conditions at the boundaries ω k × ∂Ω h of the computational domain. These conditions depend on the considered problem and/or preferred numerical scheme. For example, in the case of localized pulses satisfying (4) we can set We use homogeneous Dirichlet condition µ = 0 in (19) if the evolving envelope function does not reach lateral bounds of the computational domain, or some more sophisticated µ minimizing envelope back-reflection to the computational domain once the envelope function reaches the lateral bounds [10,18]. For (b − a)-periodic sequence of envelope pulses satisfying (5) assume that J is even and set We note, that due to periodicity the consideration of the schemes on the above defined grid Ω h is equivalent to the study of the problem on the shifted grid With these periodic boundary conditions we avoid artificial envelope reflections from the lateral boundaries at large z, but feel an influence of the envelope pulses from the neighboring periods, instead. In many cases the usage of these conditions is presupposed by the fast Fourier transform driven split-step numerical schemes resolving the dispersion operator in the wave-vector domain. In these cases we approximate the wave function ψ with the trigonometrical interpolant of J/2-degree where U are the Fourier coefficients defined as and the grid function U j is written as a discrete Fourier sum: Here, [F(U )] and [F −1 ( U )] j are corresponding components of the vectors U =F(U ) and U =F −1 ( U ) determined by the forward discrete Fourier transform F and backward Fourier transform F −1 , respectively. The spectral approximations are very accurate in the case when solutions are smooth functions. The estimate for the interpolation error of the trigonometric interpolant I J u is given by [15] .

SHALVA AMIRANASHVILI, RAIMONDASČIEGIS AND MINDAUGAS RADZIUNAS
Thus for functions u ∈ C ∞ we have the exponential convergence rate of I J u to u. The differential operators ∂ s ∂τ s , s = 1, 2, . . . can be approximated by the exact derivative of the J-degree trigonometric interpolant: Consequently, a general differential dispersion operator Dψ defined in (6) is approximated as We should take into account that on the uniform space grid we can only represent Fourier modes with ∈ J. The energy of higher waves is moved by the discrete Fourier transform to the low number waves.
It is noteworthy, that both transformations can be computed very efficiently by the fast Fourier transform (FFT) algorithm, which reduces the number of arithmetic operations to 5J log J multiplications.
In this paper we are also interested to check how accurately the discrete schemes constructed below preserve the conservation laws (14) and (16). In the case of standard NLSE (2) the discrete analogs of these laws are given by where the discrete inner products and norms in the expressions above are defined as follows: Our aim is to investigate the optimality of the following alternatives: • Finite difference method versus spectral method; • Full approximation algorithms versus splitting algorithms. Crank-Nicolson scheme. One of the most popular full-approximation finite difference scheme is the standard Crank-Nicolson scheme (CNFD). Here we restrict to the approximation of NLSE (2) allowing the local nonlinearities of the form (8) [11,17] For sufficiently smooth solutions it is easy to check that the approximation error of the discrete scheme (26), (18) with the corresponding boundary conditions (19) or (20) is O(k 2 + h 2 ). The convergence of a discrete solution of the Crank-Nicolson scheme (26) has been investigated in many papers, and the second-order accuracy is proved in various norms, see [6,33] and references therein.
The implementation of the nonlinear discrete scheme (26), (18), (19) is done by using iterative algorithms, then a sequence of simple linear systems with tridiagonal matrices are solved. For periodic boundary conditions (20) the tridiagonal structure of the corresponding matrices is slightly violated, and approximately doubled number of the arithmetic operations is required for the effective implementation of the scheme.
The template of splitting algorithms. Here we formulate the symmetrical "time" splitting (TS) template of numerical algorithms for solving GNLSE. We split the diffraction and nonlinear interaction processes by using the symmetrical Strang splitting approach (see [5,11,25,27,34,36,38]): Now depending on the method which is used to approximate "space" operators D and Q we get the "time"-splitting finite difference (TSFD) or exponent pseudospectral (TSEP) methods.
Different forms of the GNLSE problem (6), (12) are solved numerically in papers [28,31], see also references therein. Pathria and Morris [31] have investigated spectral split step methods, they have restricted to the first order splitting techniques. The splitting algorithms treat separately the linear diffraction D, the local nonlinear interaction N and Q. All algorithms differ only in the way how the subproblem Q is treated. Four different algorithms are proposed to solve the subproblem Q. The high-order spectral splitting algorithms are investigated in [28], where only the linear diffraction D and nonlinear reaction N + Q processes are split in the schemes.
"Time"-splitting finite difference method. Assuming that n is even, we formulate the corresponding TSFD scheme as follows (here the standard second order dispersion operator is considered): The important property of the splitting algorithm is that the nonlinear part of the problem is solved exactly by the explicit formula. For sufficiently smooth solutions it is easy to check that the approximation error of this discrete scheme is O(k 2 + h 2 ). The TSFD scheme satisfies the first discrete conservation law (24). This scheme is also "time" transverse invariant.
"Time"-splitting pseudospectral method. We use the same template (27), but the diffraction and non-local nonlinearity terms are approximated by using the pseudospectral method [12,15,31,38]. In general, this algorithm assumes periodicity of functions ψ (5), i.e., the related grid functions U satisfy the periodic boundary conditions (20). The "time"-splitting exponent pseudospectral (TSEP) method can be written as where β(ω ) defines the dispersion relation in the frequency domain. Mapping D(U ) is approximated by the exact mapping of J-degree trigonometric interpolant. For example, the second derivative of a discrete solution is approximated by D 2 J U , defined in (22).
The subproblem Q is solved in a different way than in [31]. We define the following righthand side of the system for the vector U of spectral coefficients The first order derivatives in the definition of Q[U ] are approximated by the exact derivatives of the J 1 -degree trigonometric interpolant D 1 J1 defined in (22). In order to control the possible ill-posedness of the numerical differentiation, we introduce the regularization by restricting the degree of the trigonometric interpolant J 1 < J.
The obtained system of ordinary differential equations is solved by using one of the Runge-Kutta type adaptive algorithms. More details on such solvers will be given in the next subsection.
If results at intermediate levels z n are not required, we can combine two neigbouring steps to solve dispersion D or local nonlinear interaction N subproblems into one explicit step, this trick twice reduces the number of required operations.
Full approximation pseudospectral algorithm. We formulate a more general pseudospectral algorithm for solving GNLSE (6). The discrete exponent Fourier transformation is applied for the GNLSE in order to get a system of nonlinear ODEs for the Fourier spectral coefficients: This system can be rewritten in the operator form as where A = diag β(ω −J/2 ), . . . , β(ω J/2−1 ) , coefficients ω l are defined in (22), and is a column vector with T denoting a transpose. It is solved by using one of the Runge-Kutta type adaptive algorithms and requires performance of forward and backward Fourier transforms each time when the local nonlinearities should be estimated at a certain "time" moment. The described algorithm is called the full-approximation exponent pseudospectral (FAEP) method.
2.2. "Time" integration. For time integration of the generalized nonlinear process Q in both FAEP and TSEP algorithms the adaptive explicit Runge-Kutta type Dormand-Prince (DP45) method is applied [25]. The fifth order method is used as a basic integration method, and the fourth order solution is used to estimate the error of the discrete solution. Such a selection is based on the fact, that the A-stability region of the fifth-order DP45 method includes an interval in the imaginary axis [−0.9972i, 0.9972i]. Also it is known, that the coefficients of DP45 method are chosen to minimize the error of the fifth-order solution. This RK algorithm has seven stages, but it uses only six function evaluations per step because it has the FSAL (First Same As Last) property.
Here we should put two important remarks dealing with the computational efficiency of the pseudo-spectral algorithm (31). First, in order to step forward from U n to new approximation U n+1 by using the DP45 algorithm we should compute the right-hand side of (31) six times, thus the forward and backward discrete Fourier transformations needed to evaluate the nonlinear function f (|U | 2 )U should be computed six times.
Second, the specific stability properties of the explicit Runge-Kutta methods should be taken into account. Let us consider the standard test problem Here the parameter λ defines the the spectral coefficient and it satisfies the asymptotic relation where p is the corresponding order of the derivative in the dispersion term, j defines the jth spectral mode. By introducing z = λz, we put equation (32) into nondimensional form We solve this problem by using the classical explicit fourth-order Runge-Kutta (RK4) method or by the DP45 method [25]. Let us denote the error of the obtained discrete solution U (z) by E(z) = u(z) − U (z). It satisfies the estimate where k = λk and k is the discretization step. The linear dependence of the error on L follows from the fact that the exact solution is a periodical function and the error depends linearly on the number of full cycles of the solution. Returning back to the physical parameters we get the error estimate Thus in order to guarantee some specified accuracy ε, the "time" discretization step k must be related to λ as Thus the restriction on k is weaker for the DP45 method. The A-stability region of the classical fourth-order Runge-Kutta (RK4) method includes an interval in the imaginary axis −i2 √ 2, i2 √ 2 [25]. Due to the stability requirements, we get the restrictions on the integration step they guarantee the boundedness of the solution. Comparison of estimates (34) and (35) shows that the A-stability requirement is slightly weaker than the approximation requirement (34). But here we should take into account that the approximation estimate must be guaranteed only for spectral modes having a sufficient amount of energy, while the stability estimate must be fulfilled for all modes. Thus the stability estimate will be the most restrictive requirement if the number of modes is large. For illustration of the obtained theoretical results we present results of computational experiments. The test problem (32) is solved in domain 0 < z ≤ 5 with the dispersion parameters We have listed in Table 1 the numbers of the discretization points N 1 RK4 , N 1 DP 45 , which are sufficient to get a bounded solution |U (z)| ≤ 1 to solve the test problem with the classical RK4 and the adaptive DP45 methods respectively. We also have listed the numbers of the discretization points N 2 RK4 , N 2 DP 45 , which are sufficient to solve the test problem with the accuracy 1 · 10 −4 .
The results of computational experiments have confirmed the theoretical conclusion that the classical RK4 solver guarantees the boundedness of the discrete solution for weaker requirements on the discrete "time" step k, but the adaptive DP45 solver gives a more accurate approximation for the same fixed values of k. Fully in agreement with the theoretical conclusions, the asymptotic of N 1 RK4 , N 1 DP 45 follow the estimate (35), while N 2 RK4 and N 2 DP 45 satisfy the estimate (34). We note that the Runge-Kutta type methods are not "time" symmetric. The latter symmetry is an important feature of discrete algorithms (see e.g. [11]) and is significant for some applications, such as delicate investigations of chaos onset and periodic trajectories in nonlinear systems [26]. On the other hand, it is well known that the time reversibility is neither necessary nor sufficient for convergence of the Runge-Kutta methods. Moreover, this kind of symmetry is immediately destroyed by dissipative effects, such as Raman term (10) or the imaginary part of the parameter β 2 which leads to diffusion in Eq. (2). On the other hand, an important advantage of the higher-order Runge-Kutta methods in the demanding applications is a natural possibility to implement an adaptive selection of the integration step. That is why we prefer to use the DP45 solver.
3. Numerical experiments. In this section we solve a few GNLSE in order to test the efficiency of time splitting and full approximation methods.
3.1. NLSE with cubic nonlinearity. Let us consider the normalized nonlinear Schrödinger equation The single soliton solution of (36) is given by [25] ψ(z, τ ) = e −i( 1 For the numerical tests we consider the superposition of two solitons. The initial profile is given with constants c 1 = 1, c 2 = 1/10 and δ = 25. We solve the problem in the domain −20 ≤ τ ≤ 80 for 0 ≤ z ≤ 40. For these values of z the solution is negligibly small outside the τ -domain, thus the choice of boundary conditions is not important for this test problem. The most difficult part of the test is to simulate accurately the collision of the solitons, when the faster soliton catches the slower and passes through it. The shape and velocity of solitons remains unchanged after the collision. This equation was used as a test problem in many papers. Very good reviews on extensive comparison of different computational methods are presented in recent papers [11,13] (see also references given therein). We note, that a standard approach is to investigate the influence of time and space discretization errors separately: • To test the spatial discretization errors, the time step is fixed such that the time discretization errors are negligible.
• To test the time discretization errors, the space step is fixed such that the space discretization errors are negligible. Our aim is to investigate the situation when both sources of the error are of the same order. In all computations the errors have been estimated by comparing the results with a numerical reference solution U ref calculated by the pseudo-spectral splitting algorithm (29) on a fine grid with J = 8192, N = 200000.
First, for illustration we have investigated the error introduced mainly by the finite difference approximation of the diffraction operator. We have listed in Table 2  The presented results show that the error changes as O(h 2 ) in accordance with the theoretical results.
Next we have investigated the interaction of both discretization errors. The Crank-Nicolson scheme CNFD (26) is used with a fixed "space" step J = 8192 and various "time" steps k. Results are presented in the first row of Table 3. Table 3. Convergence analysis of CNFD scheme (26).  (26) in the maximum norm. The behavior of the error is quite intricate. These results can be explained assuming that the errors due to the "time" and "space" approximations can cancel out each other. Our hypothesis is that the global error can be represented in the following form: where Z CN F D (L; N ) describes the "space" discretization error. By taking from Table 2 that C 2 h 2 = 0.0011 for J = 8192 and solving the equation (37), we compute the modified errors Z CN (L; N ) ∞ . They are presented in the second row of Table 3. It follows from presented results that Z CN (L; N ) ∞ is of order O(k 2 ). The indicated interaction of two types of errors can be used to minimize the global error of the discrete solution. We keep the linear relation for discrete steps k = Ch, in order to maximize the cancellation effect of the two different sources of approximation error. In the third row of Table 3 we present errors Z CN (L; N, J) ∞ of CNFD scheme when J = 8192N/3200. The results of numerical experiments show that the obtained errors are smaller even though the computational grid is more sparse.
Next we consider very briefly the results of simulations done by using the TSEP scheme (29) and FAEP scheme (31) (for more results see [11,13]). The accuracy of spectral approximation of the solitons is very high even for a quite small number of modes. The exponential convergence rate can be illustrated by the following results: where N = 200000 is used in numerical experiments. Thus if we take moderate values of J then the influence of spectral approximation error is negligible and we can investigate the accuracy of time splitting approximation. The results are presented in Table 4 and they confirm the second order accuracy. Table 4. Convergence analysis of the discrete solutions of the TSPE scheme (29) with respect to the integration step k. Errors are computed in the maximum norm for fixed h = 100/512 and various k = L/N at L = 40. In the last experiment of this section we have investigated the stability and approximation accuracy properties of the explicit FAEP solver (31), which is based on the pseudo-spectral approximation and the adaptive DP45 algorithm. In Table 5 we list the number of integration steps N for different numbers of spectral modes J, the integration accuracy is fixed to ε = 1 · 10 −8 . Table 5. Stability and accuracy analysis of the discrete solution of FAEP scheme (31). N is the number of integration steps of the explicit adaptive DP45 solver for different numbers of spectral modes J. The presented results have confirmed the theoretically established results, that for smaller J the more restrictive is the approximation requirement (34), which can be written as k ≤ CJ −2.4 , while for larger numbers of spectral modes the stability requirement (35) starts to be a dominant one.
Our main conclusion of this section is that for the given benchmark problem the most efficient is the TSEP splitting scheme (29). The high-order approximation accuracy of the full-approximation FAEP scheme (31) is not sufficient to compensate two drawbacks -the need to compute seven forward and backward FFT per one step and to keep integration step sufficiently small due to the stability requirements of the explicit DP45 solver.
3.2. Generalized NLSE. In this section we compare the numerical performance and accuracy of different numerical algorithms for a generalized nonlinear Schrödinger equation (6) with a nonlinear term defined in (12). We note, that for simple cubic nonlinear function f = γ|ψ| 2 ψ, vanishing higher order dispersion terms, β j , j ≥ 4, and β 2 γ < 1 this problem is equivalent to Eq. (13). It is known that for certain values of the coefficients and certain initial conditions, solutions of this problem experience finite-time blow up [30]. Therefore, the development of robust numerical schemes for such problems is important and challenging task.
We consider a generalized nonlinear Schrödinger problem (6), (12) with a nonlinear function f (|ψ| 2 ) = γ 1 |ψ| 2 + γ 2 |ψ| 4 , an additional higher-order (e.g., saturating) term in the nonlinear part, and a standard second-order dispersion operator Assuming that the solution ψ and all it derivatives converge to zero as τ → ±∞, one can show that the solution of this problem conserves the following three quantities [30], which can be used as convenient aposteriori measures to test the accuracy of the discrete solution [28,31]: It is noted in [31] that the quantities M, S and H provide aposteriori check on the numerical results and that the energy integral H is the most sensitive of the three in indicating numerical difficulties.
To test the performance of the numerical schemes for GNLSE with a standard second order dispersion operator we have solved numerically the test problem from [28,31] i ∂ψ ∂z This equation has a traveling solitary wave solution: The purpose of the numerical experiments is to verify numerically (i) that the proposed split-step and full approximation schemes exhibit the second-order and fifth-order convergence in "time", respectively and (ii) to compare the stability regions of both schemes. In all computations the regularization parameter J 1 = 120 is used to compute the approximation of the derivative D 1 J1 U . In Table 6 we presents errors of the solution of TSEP scheme (29) with respect to the space integration step. The errors are computed at L = 2, the computations are done in the transversal domain [0,35].
The presented results show that the error changes as O(k 2 ) in accordance with the theoretical results. The regularized algorithm is stable with respect to the Table 6. Convergence analysis of the discrete solutions of the TSEP scheme (29)  Next we present results of computational experiments with the full approximation FAEP scheme (31). They have confirmed the main theoretical conclusions. Due to conditional stability of the explicit scheme the restriction on integration step k leads to the need to integrate the system of ODEs with small "time" steps k, much smaller than required by the approximation accuracy. The efficiency of FAEP-DP45 solver can be improved if the regularization algorithm is applied to approximate the first and second derivatives of the solution. Results of computational experiments are presented in Table 7. Here N denotes the total number of integration steps, including steps rejected by the adaptive algorithm. In all computations the integration accuracy is fixed to 1 · 10 −8 . Table 7. The numbers of integration steps N for FAEP-DP45 scheme with respect to different numbers of modes J. J 1 and J 2 are the regularization parameters used to compute derivatives D 1 J1 U and D 2 J2 U . The integration accuracy is fixed to 1 · 10 −8 . It follows from the presented results, that full approximation FAEP-DP45 scheme (31) computes in stable way for any integration steps if derivatives of the solution are approximated by using the regularization technique. The approximation accuracy of order O(10 −6 ) is achieved for all presented results.
As the next example, we consider a generalized NLSE with an additional thirdorder dispersion (TOD) term (see [4,40]) where the TOD operator is included as a perturbation of the classical nonlinear Schrödinger equation (2). It is well-known that for the unperturbed (ε = 0) equation (38) the fundamental soliton solution reads where k = A 2 /2 is the soliton wave number. To explain the TOD effect, let us use the ansatz ψ(z, τ ) = A ϕ(s) exp(iA 2 z/2), s = Aτ motivated by Eq. (39) with the yet unknown function ϕ(s). The latter should be derived from the equation Now, for = 0 the operator T is positive-defined and Eq. (40) yields the solitary solution (39). The soliton can be derived even for a more general nonlinearity (12) by properly rescaled successive iterations [1]. For = 0 this is not the case, because T becomes singular. In terms of pulse propagation that means that the initial pulse (39) emits radiation. Two different cases should be distinguished here. If Aε < 0.04 the emitted radiation is extremely weak [40]. To first order of ε we have a quasi-solitary solution with the shape where τ = τ −εA 2 z. If Aε > 0.04, the radiation is strong and the carrier wavelength of the pulse is quickly shifted into the red. The underlying physical mechanism is similar to that of Cherenkov radiation created by superluminal objects in dispersive media [4].
Our main aim is to analyze the influence of the TOD to the stability and accuracy of the split step TSEP and the full approximation FAEP-DP45 methods. It follows from the computational experiments that for the accurate approximation of the soliton (41) it is sufficient to use J = 128 spectral modes. Thus if we use more spectral modes, then it is important to check how different algorithms are resolving these excessive modes.
As it was proved in the previous section the DP54 method used is only conditionally A-stable, thus in the case of TOD the integration step k automatically will be restricted to k ≤ Ch 3 . This conclusion is confirmed by results of computational experiments. The problem (38) was solved in the domain −20 ≤ τ ≤ 20 for 0 ≤ z ≤ 100, and A = 1, ε = 0.02. In Table 8 we present the total CPU times T in seconds and the number of integration steps N for different numbers of spectral modes J, the accuracy of DP45 algorithm was fixed to 10 −8 . Computations were performed on VGTU cluster of computers "Vilkas", consisting of nodes with Intel R Core TM processor i7-860 @ 2.80 GHz and 4 GB RAM. Table 8. Results of numerical experiments: TOD problem with ε = 0.02 is solved by the FAEP-DP45 scheme. T is the total CPU time in seconds and N is the total number of integration steps of the explicit adaptive DP45 solver for different numbers of spectral modes J. The accuracy of integration is 1 · 10 −8 . In the case of TSEP method the unconditional stability of the discrete scheme guarantees a more flexible control of the integration step k. It is sufficient to resolve accurately the active spectral modes and the remaining modes will be kept bounded.
The adaptive TSEP integration algorithm is constructed by using the Runge rule, the aposteriori error estimator is defined as , U (2k, z n ) are two solutions of the discrete TSEP scheme, computed starting from the same initial condition U n−2 and taking two different "time" steps k and 2k. Results of computational experiments are presented in Table 9. Table 9. Results of numerical experiments: TOD problem with ε = 0.02 is solved by the TSEP scheme. T is the total CPU time in seconds and N is the number of integration steps for different numbers of spectral modes J, the monitoring accuracy is 1 · 10 −8 . It follows from the presented results, that the number of integration steps does not depend on J, while the total CPU time depends linearly on the number of modes J.
As an interesting remark we note, that our computations have confirmed the theoretical accuracy of the corrected form of the soliton (41). For this reason we have integrated numerically the TOD problem for A = 1 and different values of the small parameter ε. The error Z (L; N, J) = U − ψ ε (L) ∞ is estimated at L = 500, here U is the numerical solution computed by TSEP method. These results are presented in Table 10. Table 10. Results of numerical experiments: problem (38) is solved by TSEP method until L = 500 using J = 512, and different ε. ε = 0.02 ε = 0.01 ε = 0.005 Z ε (500; J) 0.0619 0.0164 0.0048 The obtained computational results agree well with the theoretical prediction that these errors should be of order O(ε 2 ).
3.3. Dispersion operator defined in the frequency domain. In this section we follow [8,9] and consider the following model equation for the complex electric field E(z, t) where c is the speed of light, n 2 is a material constant [16], and the real-valued electric field E = Re [E]. The propagation constant β(ω) is related to the index of refraction n(ω) via β(ω) = ωn(ω)/c. Note, that iβ(ω) is invariant under ω → −ω transformation plus complex conjugation, this behavior is typical for all response functions. The factor "−" in the definition of β is chosen such that the linearized Eq. (42) still has a monochromatic-wave solution of the standard e iβ(ω)z−iωt form. It is convenient to change to a moving frame by replacing E(z, t) with E(z, τ ), where τ = t − z/V for some velocity V , such that Equation (43) can be related to Eq. (2) if the field in question has a pronounced carrier frequency ω 0 . One takes V −1 = β 1 = β (ω 0 ), defines ψ(z, τ ) such that E(z, τ ) = ψ(z, τ )e i(β0z−ω0t) , and applies the SVEA assuming that ψ is slow on the time scale 1/ω 0 . Finally β(ω) is approximated by a Taylor series calculated at ω = ω 0 . Then a generalized NLSE of the form (6) with the nonlinear term appears. Equation (43) has an advantage over the NLSE because it does beyond the SVEA. Recently Eq. (43) was applied to the non-SVEA regimes of pulse switching and pulse compression (see, e.g. [20,21]). Here we examine a test problem of shock formation for a 5µm silica strand fiber [35]. The effective index of refraction was approximated using approach of [7] n(ω) = 1+ 2.09315−23.98991ω 2 +2.13808ω 4 −0.0495842ω 6 +0.0000626444ω 8 1 − 52.7265ω 2 + 4.81505ω 4 − 0.117971ω 6 + 0.000323795ω 8 .
where the transparency region is 0.67 < ω < 4.7 but the rational approximation provides reasonable values on a wider interval 0.2 < ω < 16. The angular frequency ω is measured in PHz, ω 0 = 1.3. In Table 11 we present the total CPU times T in seconds and the number of integration steps N for different numbers of spectral modes J, the accuracy of DP45 algorithm was fixed to 10 −8 . Table 11. Results of numerical experiments: problem (43) is solved by the FAEP-DP45 and TSEP-DP45 schemes. The total CPU time in seconds is presented for different values of z. The accuracy of integration is 1 · 10 −8 . The presented results again show that the adaptive TSEP-DP45 method is more efficient than the FAEP-DP45. We note that both methods are robust and flexible to resolve the given generalized nonlinearity and dispersion processes.

Conclusions.
We have considered one full approximation and several splitting algorithms for numerical integration of GNLSE. The splitting algorithms are based on separate treatment of local and/or nonlocal nonlinearities and linear dispersion processes. We have shown, that the numerical schemes for GNLSE obtained by the straightforward generalization of the related pseudo-spectral schemes of NLSE are only conditionally stable. The stability of these schemes can be improved by neglecting higher-order Fourier modes in the approximation of the derivative terms of GNLSE. The efficiency of proposed methods are tested for three representative GNLSE problems, including one real-world application.