Abstract

In this paper, four compelling numerical approaches, namely, the split-step Fourier transform (SSFT), Fourier pseudospectral method (FPSM), Crank-Nicolson method (CNM), and Hopscotch method (HSM), are exhaustively presented for solving the 1D nonlinear Schrodinger equation (NLSE). The significance of this equation is referred to its notable contribution in modeling wave propagation in a plethora of crucial real-life applications such as the fiber optics field. Although exact solutions can be obtained to solve this equation, these solutions are extremely insufficient because of their limitations to only a unique structure under some limited initial conditions. Therefore, seeking high-performance numerical techniques to manipulate this well-known equation is our fundamental purpose in this study. In this regard, extensive comparisons of the proposed numerical approaches, against the exact solution, are conducted to investigate the benefits of each of them along with their drawbacks, targeting a broad range of temporal and spatial values. Based on the obtained numerical simulations via MATLAB, we extrapolated that the SSFT invariably exhibits the topmost robust potentiality for solving this equation. However, the other suggested schemes are substantiated to be consistently accurate, but they might generate higher errors or even consume more processing time under certain conditions.

1. Introduction

The one-dimensional nonlinear Schrodinger equation (1D NLSE) is a classical field equation. Its most eminent applications are related to the propagation of light waves in optical fibers and planar waveguides along with many others [1]. Specifically, the 1D NLSE is a nonlinear second-order partial differential equation, applicable to both classical and quantum mechanics. It can be written as follows [2]:

More importantly, equation (1) simulates the wave propagation in a lossless fiber optics, which is a nonlinear medium; the unknown function represents a wave. The second-order derivative represents the dispersion, while the nonlinear term represents the nonlinearity of the problem.

Besides, this equation models plenty of major nonlinear effects in an optical fiber such as self-phase modulation, four-wave mixing, second harmonic generation, stimulated Raman scattering, optical solitons, and ultrashort pulses [3, 4], especially by inserting additional terms to its canonical structure. Moreover, the 2D version of the NLSE can be obtained by replacing the second spatial derivative with the Laplacian one; it is frequently used to model the propagation of intense laser beams in a transparent medium [5].

Furthermore, there are two cases to be considered depicted by equation (1). First, when has a negative value, it is called focusing. Hence, it allows for both bright soliton solutions, localized in space and having spatial attenuation towards infinity, and breather solutions as well. Additionally, it can be solved using the inverse scattering transform to obtain its exact solution. Second, when has a positive value, it is called defocusing, which permits dark soliton solutions that have constant amplitude at infinity abreast of a local spatial dip in amplitude [6].

Along with the importance of the NLSE in the optics field, there exists another crucial application for it in water waves. This outstanding equation also describes the evolution of the envelope of the modulated wave groups. In slowly modulated groups, the wave amplitude approximately follows and satisfies the NLSE, while the value of the nonlinearity parameter term depends on the relative water depth and the unknown function represents the amplitude and phase of the water waves [7]. In deep water, when the depth of the water is relatively larger than the wavelength of the water waves, the nonlinear constant is negative and envelope solitons may occur. As a result, the envelope solitons might be expedited under external time-dependent water flow [8]. In shallow water, especially when the wavelengths are longer than 4.6 times the water depth, the nonlinearity parameter is positive and wave groups with envelope solitons do not exist. In other words, the waves of translation exist, but the NLSE does not govern them. Beside this, the NLSE is deemed to be significant in describing the formation of rogue waves, which are unusually large and unpredictable, and suddenly appearing surface waves that can be extremely dangerous to ships, even to large ones. They are distinct from tsunamis, which are caused by the displacement of water due to other phenomena and are often almost unnoticeable in deep waters. A rogue wave appearing at the shore is often referred to as a sneaker wave [9].

In addition to its distinguished applications in both the optical fibers and water waves fields, it also models the Bose-Einstein condensates, which are confined to highly anisotropic cigar-shaped traps in the mean-field regime. This Bose-Einstein condensate (BEC) is a gas of bosons or particles, which are in the same quantum state, so a single wave function can describe them [10].

Since the NLSE is only integrable in one dimension [5], Zakharov and Shabat succeeded in solving it in 1972, using the inverse scattering transform method with the aid of a corresponding linear system of equations, which is called the Zakharov-Shabat system. Besides, the NLSE is analogous to the Gross-Pitaevskii equation, which describes the ground state of a quantum system of identical bosons using the Hartree-Fock approximation along with the pseudopotential interaction model. Unlike the NLSE, this equation does not usually have a closed-form analytical solution. However, various numerical solutions exist to solve it. Correspondingly, these numerical methods can also be used to solve the NLSE, for example, the split-step, Crank-Nicolson, and Fourier spectral methods. Software programming is efficiently used to find their solution [11, 12]. Furthermore, copious authors have developed other remarkable approaches to navigate this crucial equation, and for more details, see [1328].

In this paper, the fundamental goal is to seek potential numerical approaches for solving the 1D NLSE in its rudimentary structure, presented in equation (1), which could be employed as a best-fit substitute for its exact solution. More precisely, it is remarkably noticed that the NLSE can be solved using both exact and numerical solutions. Nevertheless, its analytical solution is ineffectual because it is only provided for the canonical formalism abreast of a limited set of initial conditions, for example, the hyperbolic secant and tangent functions. This means that adding extra terms to model additional effects associated with any phenomenon, which shall proportionally change its basic structure, is not supported by an analytical solution [3], and further modifying the initial waveform might not be accompanied by an exact solution as well [29]. For instance, modeling the soliton propagation in a lossy optical fiber requires the inclusion of a fourth term in the NLSE equation to represent the losses effects; the exact solution cannot be reached in this situation because of the equation’s formula modification. Besides, solving this equation in higher dimensions is not guaranteed to have a versatile exact solution as well. As a result, numerous numerical solutions of the NLSE are widely used to opt for appropriate approximation for this equation.

This paper is organized into four sections as follows: The introduction is presented in the first section, while the second section demonstrates the four proposed numerical schemes for solving the 1D NLSE. In the third section, the numerical assessments are extensively conducted, and their simulations and findings are manifested in either graphs or tables. Finally, the overall conclusion of this article is given in the fourth section.

2. Numerical Solutions of the One-Dimensional Schrodinger Equation

Due to the inevitable role that partial differential equations (PDEs) play in modeling and describing numerous significant mathematical physics phenomena, two types of solutions are utilized to manipulate them, either the analytical or numerical methods. Furthermore, analytical solutions are extremely tough or even impossible to be found because of the complexity factor [15]. Therefore, more and more numerical schemes are introduced to offer robust approximate solutions. Moreover, although the analytical solution is usually perplexing to be obtained, it is utterly crucial for conducting extensive research and comparisons with the other numerical solutions [30]. In general, seeking numerical methods for solving the initial value problems for the partial differential equations should often fall within one of two distinct categories, which are either the finite difference methods or the function approximation methods, which contain both the spectral and pseudospectral methods [31]; the detailed organization of the PDEs’ solutions is elaborated in Figure 1.

In this section, we present the mathematical formulation of our four proposed numerical approaches for solving the 1D NLSE, namely, the split-step Fourier transform, Fourier pseudospectral method, Hopscotch method, and Crank-Nicolson method. The first two schemes are pseudospectral methods, while the others are finite difference methods.

2.1. The Pseudospectral Methods

These fascinating methods are a sort of function approximation methods because they expand the solution as a series of orthogonal eigenfunctions of some linear operator, with partial or ordinary derivatives, so that the numerical solution is connected to its spectrum. Usually, the basis function is chosen to be a trigonometric function. The major difference between the spectral methods and the pseudospectral methods is the form of the obtained solution. In the spectral methods, the spectral solutions of the evolution partial differential equations are always established in the frequency space, and therefore the solutions are formulated in terms of spectra. Meanwhile, in the pseudospectral method, the time-dependent partial differential equations are discretely solved at different values of time and space, resembling a finite difference approach. Despite this, the space derivatives are approximated using orthogonal functions such as Fourier integrals and Chebyshev polynomials. Three proposed methods to solve them are either matrix multiplications, fast Fourier transform (FFT), or convolutions. Moreover, plenty of researchers have recently revealed that the pseudospectral methods are much faster than the spectral approach. To recapitulate, the pseudospectral methods are spectral methods, but they occur in a discrete space [32, 33].

Here, we present two pseudospectral approaches for solving 1D NLSE, which are the split-step Fourier transform and Fourier pseudospectral method.

2.1.1. The Split-Step Fourier Transform

This spectacular numerical method is essential for understanding the fiber optics nonlinearity because both the dispersion and the nonlinear effects are represented alone in this process. Precisely, it can efficiently be used to model the light pulses’ propagation in an optical fiber over a short distance of . In addition to this, it is beneficial to address its advantage as a faster speed approach, especially when compared to the finite difference approach. More importantly, although this transform is deemed to be an extremely compelling numerical technique, it is an easily implemented one [34]. Besides, it is an unconditionally stable approach.

Let us start with the 1D NLSE presented in equation (1), while setting the parameter , to indicate a focusing equation associated with the initial condition , as shown below:

Then, we define the time-independent linear operator as and the nonlinear operator as , and split it into two parts as shown below:First, introduce the nonlinear step: , where . Therefore, the analytical solution will be given by the following equation:Second, introduce the linear step: . Applying Fourier transform to both sides converts the PDE into an ODE in the frequency domain, which can easily be solved as follows:The partial derivates are converted into a multiplication process in the frequency domain, and thus, the analytical solution can be computed using the relation in the following equation:Finally, the final equation can be written as follows [35]:

2.1.2. Fourier Pseudospectral Method

This method fosters the implementation of the Fourier transform for the space component, where the unknown function and its derivatives are transformed into the Fourier space with respect to . Consequently, the transformed variables become algebraic and easy to be solved, along with the finite difference discretization for the time component.

It is incredibly noticed that the second-order linear derivative can be converted into , which is a multiplication process in the Fourier space instead of derivative in the spatial space [36]. Additionally, this approach is only applicable for the periodic functions over an interval .

Solving the NLSE presented in equation (2), using the mentioned initial condition, requires the following steps.

Firstly, replace the temporal first derivative with the following difference relation:

Secondly, substitute the above relation into equation (2) to get the following equation:

Or

The previous equation provides a solution, which is only stable for values of .

However, adjusting equation (9), by using the Fornberg-Whitham principles [31], leads to an unconditionally stable solution, as shown in the equation below:

2.2. The Finite Difference Methods

Roughly speaking, these outstanding methods embrace a three-step procedure to find the solution of the nonlinear PDEs. The first step is to generate the mesh, while the second step is to apply the initial and boundary conditions. The last one is to replace the authentic function and its partial derivatives with the corresponding averages and difference relations, respectively. Furthermore, the finite difference methods could be derived either explicitly or implicitly [3741].

Here, two finite difference techniques are presented in detail for solving the 1D NLSE: the Hopscotch method, which is an explicit finite difference method, and the Crank-Nicolson method, which is an implicit scheme.

2.2.1. The Hopscotch Method

This approach is a fast explicit type of the finite difference method that is unconditionally stable; the distinguished feature about it is to replace the nonlinear term by an average formula computed at row as illustrated below [31, 42]:

Substitute equation (11) and the other appropriate difference relations in the following focusing 1D NLSE, which was obtained after reordering equation (2), subject to the initial and boundary conditions:

Following the implementation of the first two FDM fundamental steps, the previous substitutions yield the following:

or

Thus, equation (15) provides the final explicit formula, where :

2.2.2. The Crank-Nicolson Method

It is an implicit finite difference method that is widely used to solve the heat equation and some other partial differential equations. More specifically, it resembles the Runge-Kutta method. In addition, it is unconditionally numerically stable with a truncation error of order . However, unstable behavior may appear in the actual simulation due to the dominance of the nonlinear term in the dynamics [29]. Subject to the initial and boundary conditions, the implicit difference scheme stands on replacing the second partial derivatives by an average of two central difference quotients, where one is evaluated at node and the other at node :

Substituting equation (11), equation (16), and the temporal first-order difference relation in equation (12),

Hence, equation (17) can be written as

Similarly, the above equation can be rewritten in the following form as well:

Knowing that , the implicit nature of the difference method that is expressed in equation (19) implies a system of equations in the unknowns, which can smoothly be solved by any approximate method or software program such as MATLAB [31].

3. Numerical Tests

In this section, numerical experiments are conducted, as followed elsewhere [4347], to elaborate on the accuracy, speed, and stability of the proposed numerical solutions, compared to the exact analytical solution. Therefore, two tests are presented to reach a cogent milestone in measuring the efficacy among various numerical techniques, estimating the error, and deciding the most robust and fastest approach for solving the 1D NLSE. To conclude, our main purpose is to benchmark plenty of numerical techniques against the suggested analytical one for solving the focusing 1D NLSE.

The 1D NLSE is an evolution equation that depends on time. It is also widely known to be a general envelope equation. The standard form of this equation is fundamentally presented in equation (1); meanwhile, it can be reordered and rewritten as shown below:

As mentioned earlier, is the complex wave amplitude, is the spatial second-order derivative that governs the nonlinear effect, and is the cubic nonlinear term, which retains the dispersion phenomena. It has been found that this equation provides a solution if and only if the nonlinear effect is balanced with the dispersive phenomena. When this balance occurs, the one soliton solution, the multisoliton solution, and the boundary soliton solution are all acceptable and valid. For performing this test, we shall select the bright optical single soliton solution, shown below, among all the other approved solutions [19, 48, 49]:

For the numerical assessments, we discretize the space and time domains, setting parameters and to one unit and , , and to zero, respectively. This leads to the initial condition:

It also leads to the boundary conditions, at , and , shown below:

Moreover, for better characterization of the outcome of these experiments, Figure 2(a) exhibits a 3D graph of our proposed exact analytical solution, which is the hyperbolic secant soliton solution, of the focusing 1D NLSE under study. On the other hand, a 3D graph of its approximate solution at a spatial step and temporal step , using the split-step Fourier transform approach, is presented in Figure 2(b).

3.1. Numerical Simulations over a Short Range of Temporal and Spatial Values

The key pillars of this experiment stand on utilizing adaptive spatial and temporal step sizes over a predefined short range of space and time values. More precisely, to verify the efficacy of the proposed numerical algorithms, we have plotted the exact solution against the other numerical techniques over space domain and time domain , using various space steps at time step , computed at time . First, we have compared between the pseudospectral methods, which are the split-step Fourier transform and the Fourier pseudospectral method, as shown in Figure 3. Second, we have repeated the same comparison in Figure 4 between the finite difference methods, which are the Hopscotch scheme and the Crank-Nicolson method. Thereafter, the same operation is performed but by changing the values of the used time step instead, while keeping the space step fixed at 0.2, computed at the same time . First, run the test on the proposed pseudospectral methods. Then, redo the test to compare between the proposed finite difference approaches. These comparisons are presented in Figures 5 and 6, respectively. Additionally, Table 1 is generated, which contains the sum of squares error (SSE) of each of the four mentioned numerical techniques over a range of values from , at a given space step and time step , while selecting different times. As asserted by the numerical simulations, all proposed approaches are corroborated to be accurate and compatible, when compared to the exact solution. Despite this, some of them demonstrate a more reliable attitude than others under certain conditions. Specifically, the split-step Fourier transform has achieved the best performance, by obtaining the least possible sum of squares errors, followed by the Hopscotch method and then the Crank-Nicolson method, while the Fourier pseudospectral method has come in the final rank, with the most produced errors. These results are clearly elaborated in Figures 36 and Table 1.

3.2. Numerical Simulations over a Broad Range of Temporal and Spatial Values

In this section, we have plotted the exact solution against the other numerical techniques over a wider space domain and time domain , using numerous space steps at time step 0.001, computed at time . First, we have compared between the pseudospectral methods, which are the split-step Fourier transform and the Fourier pseudospectral method, as shown in Figure 7. Second, we have repeated the same comparison in Figure 8 between the finite difference methods, which are the Hopscotch scheme and the Crank-Nicolson method, but at time instead of . The reason for this special selection is that the Crank-Nicolson method does not properly function well for higher time values, and, further, it produces large errors.

Subsequently, the same operation is performed but by changing the values of the used time step instead, while keeping the space step fixed at 0.2, computed at the same time First, run the test on the proposed pseudospectral methods. The graphs are presented in Figure 9. Then, redo the test to compare between the proposed finite difference approaches at time instead of , while using small values of the time step to guarantee and achieve better accuracy for the Crank-Nicolson scheme in Figure 10.

Two tables are generated in this section. Table 2 comprises an error comparison between the four numerical methods over a wide spatial domain from −40 to 40 at a given space step and time step , while selecting different times from 1 to 10. Likewise, the generated Table 3 contains the sum of squares errors of all the mentioned numerical techniques except the Crank-Nicolson approach, at the same space step and time step and wide space domain from −40 to 40, while opting for higher time values from 20 to 80.

As bolstered by the simulation findings in Figures 710 and Tables 2 and 3, we get the following:(i)The split-step Fourier transform has retained its supreme performance over the other methods, achieving the least produced SSEs, along with the highest processing speed.(ii)The Hopscotch method has also preserved its reasonable accuracy abreast of its moderate processing time.(iii)The Fourier pseudospectral method has enhanced its performance, over the broad range of spatial values, producing smaller sum of squares errors than the errors that have been achieved at the short spatial values range, presented in Table 1 under the first experiment.(iv)The Crank-Nicolson method’s performance has diminished because when time increased, this technique’s error dramatically increased, establishing a direct proportionality correlation. However, employing smaller values of the temporal steps sizes might mitigate the situation and thus improve the accuracy, but, unfortunately, this attempt inversely constitutes a trade-off with the computational speed. In addition, for computing higher temporal values, the processing time has been observed to be extremely high, especially when compared to the other suggested approaches.

4. Conclusion

The prominent 1D NLSE, which has a major application in the fast-growing field of soliton propagation in optical fibers, was intensively investigated in this paper. It is remarkably noted that this crucial equation can be solved using different closed-form analytical solutions. However, these exact solutions can only be obtained under certain restrictions. Therefore, a plenty of studies have been dedicated to seeking plausible equivalent numerical solutions for this essential problem. In this regard, we have presented four different reliable numerical approaches, which are the SSFT, FPSM, HSM, and CNM, to solve this equation and further conducted a systematic comparison among them with the exact solution, through carrying out numerous numerical simulations, with the aid of MATLAB, to disclose the advantages and disadvantages of each of them. Consequently, the obtained numerical results endorsed that the SSFT has demonstrated a superlative capability in handling this equation by attaining the least sum of squares error alongside fastest computational speed; meanwhile, the HSM has accomplished the second rank of reasonable accuracy over the short and wide ranges of temporal values from 0 to 80. Moreover, The FPSM has improved its accuracy when manipulating larger ranges of spatial values from −40 to 40. In contrast, over this broad range of spatial values, the CNM has been an impractical method when dealing with temporal values higher than 15. Additionally, it has consumed long processing time, along with generating remarkably high errors upon reaching this certain limit. To sum up, the four proposed approaches have been corroborated to have compatible approximate solutions, when compared to the exact solution. However, some limitations might appear when adjusting the spatial and temporal ranges.

Data Availability

All the relevant data are available upon request from the authors.

Conflicts of Interest

The authors declare that they have no potential conflicts of interest.