1. Introduction

Optical communication systems have evolved to the point where it is now possible to transmit terabits/s for thousands of kilometers through a single optical fiber. This has been accomplished through the use of Wavelength Division Multiplexing (WDM), optical amplifiers, special fiber designs, and dispersion-management techniques. These systems are of such complexity that computer simulations are necessary for their proper design. Because of the large amount of sample points that are often necessary to accurately describe the optical field in these simulations, the computer program has to be highly efficient.

Wave propagation in optical fibers is described by a G-NLSE whose particular form depends on the physical effects considered: higher-order dispersion, Raman scattering, background losses, etc. One of the most efficient algorithms to solve G-NLSEs uses the Split-Step Fourier Method (SSFM) [1]. Four recent publications [2–5] focus on the optimization of this method in modeling optical-fibers communications systems. The approaches of these works are quite different. In Ref. [2] the authors focus on the accuracy and efficiency of different spatial step-size selection criteria. In Ref. [3] is shown that the fourth-order SSFM scheme provides more efficiency and accurate solutions than the first- and second-order schemes. In Ref. [4] the authors study the optimization of the simulations involving non-uniformly time-sampled optical pulses employing a novel numerical algorithm (instead of the well-known SSFM), the Split Step Spline Method. Finally, in Ref. [5] the authors propose a predictor-corrector technique to accelerate SSFM, whose advantage is to retain information on the approximate solution at each of the mesh points z₀,,z₁… z_j before the approximation at z _j+1. These four approaches do not exclude each other and, moreover, they can potentially work in a complementary form.

In this article we focus on the spatial step-size selection criteria, as Ref. [2]. We propose a new user-independent algorithm to the adaptive determination of the step size at each interaction. For reasons later clarified, we call this new method the Uncertainty Principle Method (UPM).

Using the algebra of operators commonly used in Quantum Mechanics (QM), we derive an upper bound to the local error at each numerical step. It is well known that at first order the SSFM error is given by the commutation relation between the dispersive and non-linear operators appearing in the G-NLSE to be solved. Thus, calculating an upper bound to the mean value of this commutation relation, through the mathematical inequality which in QM gives rise to the uncertainty principle, provides an easy way to calculate an upper limit to the local error of the SSFM at an inexpensive computational cost. A step-by-step presentation of this idea and how to implement it is developed in Section 2. In Section 3 we present numerical results for three different optical-fiber communication systems of practical interest, showing the efficiency of the UPM. In Section 4 we discuss some important properties of the UPM in light of a comparative study concerning methods previously presented in the literature. Finally, in Section 5 we draw our conclusions.

2. Presentation of the step-size selection algorithm

Sinkin et al. [2] presented for the first time a robust user- and system-independent step-size selection criterion for solving G-NLSEs, calling it the Local Error Method (LEM). One of the most relevant novelties introduced with the LEM is that the selection algorithm does not rely on the physical intuition or expertise of the user.

In this section, we first present our criterion, which is user- and system independent. Here we consider the standard SSFM and without background loss. Then we show how this criterion can be extended to include background losses and how it can be applied when implementing the Symmetrized-SSFM (S-SSFM) or higher order SSFMs.

2.1 Standard SSFM for G-NLSE without background losses

Wave propagation in lossless optical fibers is described by the G-NLSE [1]

where A=A(z,t) is the slowly varying field amplitude, β₂ and β₃ are, respectively, the second and third order dispersion coefficients, γ is the nonlinear coefficient (optical Kerr effect), z is the position along the fiber, and t is the local time, i.e., in a reference frame that travels with the average group-velocity of the pulse. (Note that we are not considering here higher order dispersion coefficients β₄, β₅,… although high order dispersion and attenuation can be included in the NLSE, see below). Equation (1) can be written as

with D and N defined as

In the SSFM the fiber is divided in steps of length h (h can vary along the fiber) and the field is then propagated at each step assuming first no dispersion (β₂=β₃=0) and then no nonlinearity (γ=0). The G-NLSE with no dispersion is directly solved in the time domain, whereas the G-NLSE without the nonlinear term is solved in the frequency domain using a Fast Fourier Transform (FFT) routine:

Where F is the FFT (F ^-1 is the IFFT or inverse FFT) and D̃(ω) is the dispersion operator in the frequency domain, given by replacing the differential operator ∂/∂t in Eq. (3) by iω(ω is the frequency). Observe that N is a multiplicative factor in the time domain, whereas D̃ is also a multiplicative factor, but in the Fourier domain.

A good criterion for the choice of the spatial step-size is extremely important for the optimization of the SSFM. The accuracy and execution time of a program implementing the SSFM depends on the choice of the spatial and the temporal step sizes, h and dt, respectively. Here we concentrate on the algorithm for the choice of h, which, as pointed out in reference [2], does not change qualitatively when dt is changed.

It is well known [1] that the error of the SSFM is given by the fact that the dispersion and nonlinear operators do not commute. That is, the error is determined by the fact that it is not the same to propagate first with the nonlinear operator and then with the dispersion, than propagate with these operators but acting in the reverse order. Let us briefly review the mathematical deduction of this fact.

A formal solution of the NLSE, assuming that h is small enough such that variation in N is negligible in an interval from z to z+h, is

But, for any two operators D and N [6]

where [D,N]=DN-ND is the commutator between D and N. Thus, from (4), (5) and (6), to the lowest order in h, the error in the SSFM at each sample time point t is -(1/2)h ²[D,N]A(t).

From a practical point of view, we are not interested in the error at each point in the time domain, but rather in same one-dimensional real number capable of providing enough information to correctly choose the longitudinal step size to be used in the SSFM given a desired precision target. By “providing enough information” we mean that this real number, representing some kind of error averaged over the time domain, must be such that the simulation converges to the correct solution if this number goes to zero. Thus, if the longitudinal step is chosen small enough in order to bound this averaged error by some desired value, one can be sure that the simulation converges. As a good parameter, Reference 2 defines and uses the relative local error, given by

where A_n is the numerical result after one step numerical propagation and A_a is the analytical result, supposed to be somehow known or estimated, and ‖A‖=(∫|A(t)|² dt)^1/2. Using Eqs. (4), (5), and (6) it is easy to show that

Sinkin et al. [4] use an estimative of the relative local error δ value to determine the longitudinal step size at each given desired target precision. In the method we are presenting here, we use another parameter to the adaptive control of the step size, based in the averaged error,

where we used the mean value of the commutator, 〈[D,N]〉, in “state” A as given in Quantum Mechanics (QM). Thus, the field amplitude A plays the analogous role of the wave function in quantum mechanics. We will show trough this article that this error ε, defined with the help of QM, is also a good parameter to determine how much the numerical results differs from the exact result. Since the error at each sample time point t is, as mentioned before, - (1/2)h ²[D,N]A(t), then ε tends to zero when the error at all sample times goes to zero and, moreover, ε gives some kind of “averaged” value of the error. In Section 4.2 we return to discuss theoretically the validity of our quantum mechanically defined error ε as a good parameter to guarantee the convergence of the simulations, comparing it with the relative local error δ (Eqs. (7) and (8)) defined in Ref. [2]. It is not straightforward to show that, as in the case of δ, the simulations converge if the quantum mechanically defined error ε is bounded by some small value. We discuss this point theoretically in Section 4.2, but, before, in Section 3, we perform numerical experiments showing the usefulness of ε as a parameter to guarantee convergence in three cases of practical interest. Now we focus on how the uncertainty principle can be used to establish an upper bound to the value of ε, enabling the optimization of the SSFM through a correct choose of h.

A very useful concept in QM is the uncertainty principle, which gives an upper bound to the average of the commutator of several pairs of physical operators such as position and momentum. This principle arises from the Schwarz inequality [7] and thus applies not only in the context of quantum mechanical operators but to any pair of hermitian operators. In particular, it also applies to the operators D and N (since they are hermitian) and takes the form [8]

where ΔD and ΔN and are standard deviation as defined in QM, calculated as

for any operator B. Combining Eqs. (9) an (10) we obtain

ε=(1/2)h ²|<[D,N]>|≤(1/2)h ² 2ΔDΔN,

and, thus, we can chose the step size h according to

in order to bound the error ε by any desired value chosen to replace it in Eq. (12). This Eq. states our criterion. From the concept in which it is based on, we call it the Uncertainty Principle Method (UPM). Care must be taken to do not misunderstand this nomenclature. We use this name to refer to the mathematical inequality which in QM expresses the Uncertainty Principle, and because we use the algebra of operators and calculate averages in the same way that in QM, and that ultimately inspired our study. The analogy with QM formalism is direct if we think of the field A(z,t) as the wave-function of QM. After this point the analogy breaks down. In physics, the Heisenberg Uncertainty Principle expresses a limitation on accuracy of simultaneous measurement of canonically conjugated observables such as the position and momentum of a particle. There is no such an interpretation for ΔD and ΔN. Dispersion and nonlinearity can be simultaneously measured (or computed) to any degree of accuracy. Our method could be called, perhaps more precisely, the Schwartz inequality method, but for the reasons just exposed, we prefer the UPM name.

A practical implementation of the UPM is straightforward. One needs only to include a routine to calculate the QM average, which is the weighted average over the instantaneous power in the time domain or over the spectral density in the frequency domain (see below). This routine is then used to calculate the standard deviations ΔD and ΔN that are used in Eq. (12) to determine the step size h for a chosen error ε (typically ε=10^-3). When modeling communication systems it is not necessary to re-calculate h on every step. Normally, we call this routine on every 20 steps.

The standard deviation ΔN can be calculated directly in the time domain using Eq. (9), but for the dispersion operator ΔD, which involves a second order time-derivative, this can be computationally slow. Fortunately, we can do better if we realize that we can calculate QM averages (and thus QM variances) in the frequency domain with exactly the same result [9], i.e.

where Ã=Ã(z,ω)=F A(z, t) is the Fourier transform of the field amplitude. For instance, the QM variance of the dispersion operator defined in Eq. (3) is

2.2 Inclusion of background loss

The SSF method can be applied to equations having the general form

where the operators L and N are such that Eq. (15) is easily integrated in the time-domain if N acts alone and in the Fourier space if L acts alone. There are many problems in physics where the SSF is used. Examples from optics are the wave equations describing spatial solitons or linear propagation in gradient index media [10] (where t is substituted by transverse coordinates, and the Fourier transform operates in the bi-dimensional wave-vector space), nonlinear surface waves [11] (where t is the coordinate transverse to the surface), as well as several problems in Beam Propagation Methods. In many of these problems (including optical fibers) the physical system is not conservative on account for losses. In these cases the operators may not be hermitian and the UPM cannot be applied directly. We discuss now how to handle these situations.

Consider for example the case of a single mode with an arbitrary dispersion relation and a frequency dependent attenuation coefficient. The Fourier transform of the linear operator is

Observe that D is hermitian but L is not. If the loss coefficient can be considered frequency independent over the field spectrum, α(z,ω)=α(z), then it commutes with N and the UPM as given by Eq. (12) can be used. If this is not the case, then the quantum mechanically defined error ε is given by

Now, Eq. (10) can be applied to both commutators inside the square root in (17) yielding

Thus, choosing h as

we guarantee that the averaged error ε is smaller than the chosen value to replace it in Eq. (19). In optical fiber communications, normally Δα≪ΔD and Eq. (12) can be used. However, if the field spectrum is broad and close to some resonance (for example, the water peak near 1.4 µm of conventional silica fibers), or in the cases of doped fibers (such as a distributed amplifier) or dispersion flattened fibers with β₂=β₃=0, then Eq. (19) can be safely used.

In the numerical experiments presented in Section 3 we use α=0 and Eq. (12) to determine h. Here, for completeness, we deduced Eq. (19), but we will not use it in the remaining of this article.

2.3 Symmetrized SSFM

The accuracy of the SSFM can be improved if the field propagated with dispersion until the middle of the step, h/2, is used as the input field to perform the non-linear propagation [1]. Mathematically,

The main advantage of this method is that the leading error results from the double commutator in Eq. (6) and is of third order in the step size h, which can be easily shown applying Eq. (6) twice in Eq. (20). Because of the symmetric form of the above equation, this version of the SSFM is called Symetrized SSFM (S-SSFM). It is also known as the second-order SSFM [3] (the first order is the standard SSFM discussed in Section 2.1).

The UPM can be easily extended to the S-SSFM. Since the error of the method is now given by h ³(1/6){(N+2D)[D,N]+[D,N](D+2N)}, the Schwartz inequality can be applied to establish an upper limit to the QM average of this operator. We do not show here explicitly the result because it is rather cumbersome, but its derivation is straightforward using Eq. (10).

For a given target accuracy, the second order or S-SSFM determines larger steps than those given by the standard first-order scheme. Thus, Eq. (12) can be used as a conservative but safe criterion to determine h, even when the S-SSFM is used. This is what we do in our program, i.e., we calculate h using Eq. (12) even though our algorithm uses the symmetrized SSFM. Although this procedure is conservative, we show below that it is more efficient than other methods in a wide range of situations. In principle, we could improve the efficiency of our program using an UPM criterion extended to case of the S-SSFM, but it is not clear that calculating the step in a more complicated expression than Eq. (12) could be faster than just using a few more FFT’s. This is an open question for future works.

3. Numerical results

In this Section we study the efficiency of the UPM, implemented as explained in Section 2, through a set of numerical experiments representing three cases of practical interest: the propagation of a second order soliton, two first-order colliding solitons, and an eight-channel WDM system propagating through 10- and 50- km of a Dispersion-Shifted (DS) fiber.

Following previous work [2, 12], we use the number of FFTs per simulation as a measure of the total computational cost (in section 4.3 we discuss this point). In general, one is interested in the number of FFTs that are necessary to achieve a given accuracy. How to measure the accuracy of a simulation is, then, a fundamental question concerning the comparison of different numerical schemes. In Ref. [3], for instance, particular attention is paid to the conserved quantities as indicators of accuracy. If the fiber is lossless, for instance, the total energy at the end of the simulation should be the same as that at the input. Here, following Ref. [2] and [5], we use the global relative error δ_G as an indicator of the accuracy. This is defined by Eq. (7) but with A_n(t) being the numerical field and A_a(t) the analytical solution at the end of the propagation (and not after a single step, as in the definition of δ). In the case of isolated solitons, for instance, an exact solution actually exist, and this is used for A_a(t). In the general case, however, A_a(t) is computed numerically using the S-SSFM with a very small h. This is done running several simulations with decreasing constant h sizes until the results from two successive simulations are equal at all sample points within machine precision. We have chosen to use δ_G as an accuracy indicator because it has been used in two recent publications on the efficiency of the SSFM in modeling optical communication systems [2,5], and therefore we expect that the community working on this area is familiar with δ_G.

3.1 Second-order soliton

As a first numerical example, we study the propagation of a second-order soliton of the form A(t)=2η(-β₂/γ)^1/2 sech(ηt) with η=0.44 ps^-1, β²=-0.1 ps²/km and γ=2.2 W^-1/km. For these parameters, the peak power is 35 mW and the FWHM pulse duration is 4 ps [2]. We propagated this pulse through ~81.2 km, which correspond to one period of the soliton. We used 2¹⁰=1024 sample points and a simulation time windows of 50 ps. In Fig. 1 we show the results. For completeness, we show the results not only for the UPM and LEM, but also for the Non-Linear Phase-Rotation Method (NPRM), which has been extensively used (and its still used in commercially available software packages). The NPRM is thoroughly discussed in Ref. [2]. As seen in Fig. 1, the UPM works better in the low accuracy regime, while for global errors smaller than ~10^-3 the LEM is more efficient. This happens because the LEM warrants accuracy to third order in h, while the implemented UPM warrants accuracy to second order.

 

Fig. 1. Number of FFTs versus global relative error for a second order soliton propagating through ~81.2 km of fiber. Results form the UPM, LEM, and NPRM are shown. The lines joining the points are just aids to the eyes in identifying the general behavior of the methods.

Download Full Size | PDF

3.2 Soliton collision

We simulate the collision of two first-order soliton of the form A(t)=η(-β₂/γ)^1/2sech(ηt), where the values of η, β₂, and γ are the same as in the previous section, giving a pulse duration of 4 ps and a peak power of 8.8 mW [2]. The central-frequency difference between the solitons is 800 GHz, and at the input, they are separated in time by 100 ps. The fiber length is 400 km. We used 3072 sample points and a simulation time windows of 400 ps. The results are shown in Fig. 2. Clearly, the UPM behaves better for global errors larger than ~10^-4, while the LEM is more efficient if extreme accuracy is desired.

 

Fig. 2. Number of FFTs versus global relative error for a two first-order soliton collision propagating through 400 km of fiber. Results from the UPM, LEM, and NPRM are shown.

Download Full Size | PDF

3.3 WDM system

It has been shown [2] that the Walk-Off Method (WOM), in which the step size is chosen to be inversely proportional to the fiber dispersion, performs even better than the LEM for multichannel systems. Thus, in this section we compare the UPM against both the LEM and the WOM. We simulate a 10 Gb/s eight-channels RZ (return-to-zero) system propagating through 10- and 50- km of a Dispersion Shifted (DS) fiber. Other simulation parameters are shown in Table 1. In each channel, the bits were modeled by a fourth-order super-gaussian function and the bit sequence is random and independent of the other channels. The number of FFTs performed to obtain a given accuracy for the simulations implementing the UPM, the LEM, and the WOM and for the propagation distances of 10 and 50 km are shown in Fig. 3(a) and Fig. 3(b), respectively. We stress that in these simulations the WOM determines a constant step size since we are using only one type of fiber.

Table 1. Main parameters used in the simulation of an eight-channel WDM system

View Table

 

Fig. 3. Number of FFTs versus global relative error for a WDM eight-channel system simulated using the UPM, LEM and WOM for propagating distances of (a) 10 km and (b) 50 km.

Download Full Size | PDF

It is clear from Fig. 3(a) that the UPM and the WOM behaviors are very similar and they are more efficient that the LEM over the entire range of values simulated here. Fitting the curves in Fig. 3(a), we find that the LEM curve would cross the other curves, becoming the more efficient method, at global accuracies of ~10^-5, far from the region of practical interest in system design (which is between 10^-1 and 10^-3 [2]). With a propagation distance of 10 km, the number of FFTs performed to obtain a given accuracy using the LEM is in all cases ~5 times greater than when using the UPM or WOM methods. In Fig. 3(b), however, when the propagation distance is, more realistically, 50 km, the UPM and LEM are, in average,~two times more efficient than the LEM for global accuracies between 10^-1 and 10^-3. Moreover, the LEM curve crosses the others, becoming more efficient, at global accuracies of ~10^-4, closer to the region of practical interest. This results can be generalized so that, for any given input field, the longer the distance to be simulated, the greater the region of accuracies in which the LEM becomes the most efficient method.

4. Discussion

4.1 Variation of method parameters

We now address two important questions concerning the method parameters (ε or δ) used in the UPM and the LEM. These questions are, literally from reference [2, Section III.E]: “First, for a given global error, how much does the method parameter depend on the particular system? Second, by what factor should the method parameter be decreased to halve the global error?” To address these questions we show in Fig. 4 the dependence of the relative global error on the method parameter (ε in the UPM case and δ in the LEM case) for different systems.

To answer the first question we observe that, from Fig. 4, for a given global relative error of 10^-3 the method parameter varies in both cases (UPM and LEM) over one order of magnitude for the systems considered. Thus both methods are equally robust in the sense that no matter the system to be simulated, the method parameter value to achieve this desired target accuracy are not significantly different. As a comparison, the WOM exhibits variations of 2 to 3 order of magnitude larger than the LEM [2]. On the other hand, it is clearly seen in Fig. 4 that for achieve a global relative errors of 10^-4 the UPM becomes more sensitive to the specific system. Also, it is clear that the curves corresponding to the 10- and 50- km WDM propagation are more separated in the UPM than in the LEM graphs. We attribute this observed larger variance on the method parameter to the average error ε used in the UPM.

Since ε performs the integration and then takes the absolute value (Eq. (9)), it is quite broad in its measuring range and because of that we expect that, in general, the UPM will show a larger sensitivity to the system condition to achieve a given target accuracy than the LEM. We return to this particular point in Section 4.2.

To answer the second question we observe that the fitted curve slopes, shown in the labels of Fig. 4, are in the three cases smaller than 1 for the UPM, while in the LEM case the curve for the second-order soliton has a slope of ~1.2, greater than 1. Ideally, these slopes should be equal to one, thus ensuring that a reduction of the method parameter by a certain factor will cause a reduction in the global error by the same factor. As discussed in reference [2, Section III.E], the value of ~1.2 for the LEM simulating a second order soliton may be related to the fact that the true local relative error is actually estimated in the LEM.

 

Fig. 4. Global relative error as a function of the method parameter (ε or δ), for the UPM and the LEM, and linear fit of each curve. We show the results from the systems simulated in Section 3. The slopes given by the linear fits of the curves are shown in the labels.

Download Full Size | PDF

4.2 Quantum mechanically defined error versus relative local error (ε vs. δ)

In comparing the LEM and UPM, we first note that h in the LEM is determined by a recursive algorithm that is time consuming, whereas in the UPM h is determined using an analytical expression, Eq. (12), which is very fast to be calculated. This explains why the UPM is more efficient for moderate accuracy.

On the negative side, there might be cases where ΔN or ΔD are very small and Eq. (12) gives an excessively large step, leading to errors at some particular time points. This might happen because Eq. (12) warrants only that the time integrated error is bounded. In the LEM, δ is defined by an integral of the squared absolute values at each sample point, thus warranting that if δ goes to zero, then the error is bounded at all time points. This is not assured in the UPM. Since in Eq. (9) we first integrate the error and then take the absolute value of the integral, if the QM averaged error goes to zero this does not ensure that the error is bounded at all time points. Because of that the UPM will show a large sensitivity to the system condition to achieve a given target accuracy, as said in Section 4.1 and shown in Fig. 4.

Although ε do not ensures point-to-point convergence, two comments are relevant and enable the use of ε in a safe manner: First, the numerical experiments which results were shown in Section 3 show that for these representative situations, ε works even better than δ as an indicator of accuracy and as a parameter to decide the step size to be taken in a wide range of situations of practical interest. Second, the cases in which the real and imaginary parts of the argument in the integral in Eq. (9) are perfectly canceled during the integration seem to be very rare. During one step propagation the error at half of the points must have exactly the same and opposite value than in the other half. But generally the errors at different sample points do not cancel each other and, moreover, in many cases they have the same sign. In a lossy fiber, for instance, if the step is too large, the effects of the attenuation will be over-estimated, reducing the field magnitude at all sample points.

A final relevant comment on this issue: The UPM was thought to solve actual communication systems, where statistically sound results are more desirable than detailed accuracy at all time points, and where computational time should be minimized. In communication systems design, one is more interested in eye diagrams than in the actual shape of individual bits. Furthermore, in actual communications systems, ΔN and ΔD do not vanish unless there is no nonlinearity or no dispersion at all. In both cases the UPM gives quickly the correct answer: you can take h as being the total fiber length (in both cases there is an exact solution either in the time domain or in the frequency domain). The UPM was not intended to solve problems such as supercontinuum generation, self-steppening, or the propagation of a cw laser; although we feel that it could be easily extended to deal accurately with these problems.

4.3 Number of FFTs as a measure of computational cost

A natural concern when studying the efficiency of the UPM is if the calculation of the QM average for the calculus of ΔD and ΔN are not mayor time-consuming computer tasks and, then, if this operation take longer than performing an FFT. In such case, the number of FFTs would not be a good parameter to quantify the computational cost of the UPM.

This concern could also be raised in relation to the LEM, because it performs the calculus of h at each step, which, as shown by Eq. (7), involves the calculus of two integrals. Roughly speaking, 5N arithmetical operations (N is the number of sample points) are necessary to calculate δ (N subtractions in the numerator, 2N multiplication to calculate the square values of the field at each sample point at the numerator and at the denominator and 2N sums to integrate them). Similarly, 15N arithmetical operations are necessary to compute ΔD and ΔN in the UPM. On the other hand, it is well-known that each FFT performs between 4Nlog₂N and 5Nlog₂N operations, depending on the implemented algorithm.

Now, to simulate a real optical-fiber communications system, typically not less than 2¹⁵ sample points are necessary (here, for instance, we used 2¹⁶ to simulate the simple WDM system of Section 3.3). Using these numbers in the estimates above, the number of arithmetical operations performed to the calculus of h in the UPM is less than 10% of that used by the FFT. Of course, care must be taken when simulating systems with very small number of sample points because in this case the time consumed to calculate h can be of the order of the time consumed by the FFT algorithm.

5. Conclusions

We presented a user- and system-independent spatial step-size selection algorithm for solving G-NLSEs using the SSFM. We called this method the UPM. We compared our method against other existent methods and showed that the UPM is the most efficient method in a wide range of situations of practical interest, including WDM systems and second-order solitons simulations. We also showed that the longer the propagation length to be simulated, the larger the region of accuracies in which the LEM [2] becomes the more efficient method. Since the UPM, as the LEM, does not rely on the physical intuition of the user, it can be, in principle, applied to solve any G-NLSE concerning, for instance, fluid dynamics, plasma physics, and wave propagation in atmosphere [13].

Future improvements of the method here presented can follow two paths, because the UPM here implemented is conservative in two senses. First because we used Eq. (12) to determine h but we implemented the S-SSFM, which enables the use of an even larger h to obtain the same accuracy. As said in Section 2.3, how much improvement it would be obtained if we used the UPM extended to the S-SSFM case is left as an open question to be addressed in future works. Second, because we use inequality (10) in a conservative way, using always the worst case, i.e., the equal case, with gives the smallest h. If the right hand side of inequality (10) is much larger than the left hand side, the UPM is being conservative. It is well known that the equal sign holds only for gaussian wave-packets [7]. So, unless the field is gaussian and retains a gaussian shape during the propagation, the UPM is being conservative. The challenge to improve the UPM is to find situations in which one can be sure that (10) is far from the equality establishing how much the step can be increased in relation to the equal sign case.

Several methods where developed in recent years [2–5] for the optimization of the SSFM. Although different, the approaches to improve the SSFM presented in these articles can work in a complementary form. Some “hybrid” schemes employing at the same time two or more of these proposed methods could be easily implemented. Future trends in this area include the study of algorithms capable of identifying the most efficient “hybrid” implementation of the SSFM in modeling any given system. We end this article, then, giving three examples easily imaginable on how these approaches could work together. First, in Ref. [3] it has been shown that, in a wide range of situations, the fourth-order SSFM performs better than the first and second order schemes (this last called here S-SSFM). Thus, a future improvement of the UPM, the LEM or the method presented in Ref. [5] can be expected if implemented employing the fourth-order SSFM. Second, the criterion to determine the step size used in the UPM can be employed to determine the initial step size of the LEM in situations in which this last performs better. Third and last, both the UPM and the LEM can be implemented to determine the spatial step-size when implementing the predictor-corrector SSFM proposed in Ref. [5].