Efficient Exponential Time-Differencing Methods for the Optical SolitonSolutions to the Space-TimeFractionalCoupledNonlinear Schrödinger Equation

+e coupled nonlinear Schrödinger equation is used in simulating the propagation of the optical soliton in a birefringent fiber. Hereditary properties and memory of various materials can be depicted more precisely using the temporal fractional derivatives, and the anomalous dispersion or diffusion effects are better described by the spatial fractional derivatives. In this paper, one-step and two-step exponential time-differencing methods are proposed as time integrators to solve the space-time fractional coupled nonlinear Schrödinger equation numerically to obtain the optical soliton solutions. During this procedure, we take advantage of the global Padé approximation to evaluate the Mittag-Leffler function more efficiently. +e approximation error of the Padé approximation is analyzed. A centered difference method is used for the discretization of the space-fractional derivative. Extensive numerical examples are provided to demonstrate the efficiency and effectiveness of the modified exponential timedifferencing methods.


Introduction
e coupled nonlinear Schrödinger equation (CNLSE) can be employed in simulating the propagation of the optical soliton in a birefringent fiber [1][2][3]. A soliton is a solitary pulse which can travel at a constant speed and keep a stationary shape due to the balancing of the self-phase modulation and the group velocity dispersion effect in fiber optics [4]. According to Agrawal [5], in a fiber communication system, the input pulse may be orthogonally polarized in a birefringent fiber. e polarized components can form solitary waves, which are named as vector solitons. Because of the nonlinear coupling effect, the vector solitons can propagate undistorted even when the components have different widths or peak powers.
During the last few decades, researchers have found that hereditary properties and memory of various materials can be depicted more precisely using the temporal fractional derivatives [6][7][8]. It is also shown in [9,10] that the anomalous dispersion or diffusion effects are better described by the spatial fractional derivatives. e anomalous effects reflect the Lévy-type particle movement, different from Brownian motion, which depicts the classical random movement of particles. erefore, the space-time fractional coupled nonlinear Schrödinger equation (FCNLSE) is useful in modeling solitons in fractional fiber optics.
In this article, we consider the FCNLSE given as follows [11]: with the initial conditions and homogeneous Dirichlet boundary conditions on and complex functions u and v represent the amplitudes of orthogonally polarized waves in a birefringent optical fiber. D α t u � z α u/zt α is the Caputo derivative of u in time, D β x u � z β u/zx β is the Riesz derivative of u in space, 0 < α ≤ 1, 0 < β ≤ 2, and the parameters δ and ρ are some real constants.
Both analytical treatments and numerical methods have been investigated for the fractional Schrödinger equations and some novel types of nonlinear Schrödinger equations. In [12], an extended sinh-Gordon equation expansion method is adopted to solve the space-time fractional Schrödinger equation analytically. In [13], a modified residual power series method is implemented on the fractional Schrödinger equation. In [14], the L1 scheme together with the Fourier-Galerkin spectral method is employed to discretize the time-fractional Schrödinger model. In [15], a Fourier spectral exponential splitting scheme is constructed to solve the space-fractional initial boundary value problems. In [16], a generalized exponential rational function method is applied to a new extension of the nonlinear Schrödinger equation. In [17], a cubic-quartic nonlinear Schrödinger equation is solved analytically for the dark, singular, and bright-singular soliton solutions. In [18], a modified expansion function method and an extended sinh-Gordon method are proposed for the M-fractional paraxial nonlinear Schrödinger equation to obtain soliton solutions.
However, to the best of our current knowledge, numerical methods for the coupled space-time fractional Schrödinger equations are rarely considered. In this paper, we modify the exponential time-differencing (ETD) method for the time-fractional nonlinear PDEs, introduced in [19], by applying the Padé approximation. en, we combine the modified ETD scheme with a fourth-order fractional compact scheme in space. During this procedure, the nonlinear term of the equation is computed explicitly, and the calculation of the fractional exponential time integral is undertaken more efficiently.

Discretization in Space
e spatial Riesz derivative is defined in [10] as x) are the left and right Riemann-Liouville derivatives: in which Γ(·) is the gamma function.
It is stated in [20] that the approximation of the left where Similarly, the approximation of the right derivative Matrices L M are transposes to each other in (5) and (6).
Furthermore, we use the centered difference method for the fractional derivative to approximate the Riesz derivative in the following way [21]: where Noticed that scheme (8) is second-order convergent in space, Ding et al. [22] generated a compact scheme to improve the order of convergence: where x)/h β is the second-order approximation (8). As been proved by eorem 11 in [22], compact scheme (10) is fourth-order convergent spatially.

The Exponential Time Integrator
We obtain a system of time-fractional equations after discretizing FCNLSE (1) in space: where z α /zt α denotes the Caputo derivative, A is the M x × M x matrix in the Riesz derivative approximation, F: R M x ⟶ R M x contains the nonlinear function and the boundary conditions, and As been computed using the variation of constant formula in [23], system (11) has an analytical solution: (12) where e α,β (t; λ) denotes the inverse function of the Laplace transform s α− β /(s α + λ) to A, and e α,β (t; λ) can be calculated taking advantage of the Mittag-Leffler (ML) function E α,β (z): . (13) Formula (12) can be written in a discrete form after discretization on [0, T] with an equal-spaced mesh-grid t n � nτ, n � 0, 1, . . .: en, the ETD scheme can be denoted as [19,23] U n � e α,1 t n ; where U j is the numerical approximation to U(t j ), and the convolution weights W n,j can be computed as Scheme (15) is called the one-step ETD scheme. Garrappa and Popolizio proved in [19] that the one-step ETD scheme (15) has the absolute approximation error where M is the temporal step number and C is a constant relating to T and α. is inequality tells us that ETD scheme (15) is first-order convergent temporally. e two-step ETD scheme is also constructed in [19]: where W (1) n � e α,α+2 t n− 1 ; A + e α,α+1 t n ; A − e α,α+2 t n ; A , Garrappa and Popolizio also proved in [19] that the twostep ETD scheme (18) has the absolute approximation error Err n � ‖U(t n ) − U n ‖ satisfying where M is the temporal step number and C is a constant relating to T and α. is inequality tells us that ETD scheme (18) is 1 + α { }-order convergent temporally. To relief the burden of computing the function e α,β (t; A), we transform it using the multiplication of eigenvectors and functions of eigenvalues [23]: where A is diagonalizable, with λ k 's to be its eigenvalues, Z is the composition of A's eigenvectors, and D � diag(λ 1 , λ 2 , . . . , λ m ). Using this decomposition, we avoid the computation of the ML function of matrices, which is really time consuming. We only need to calculate the ML function with inputs of numbers and multiply the matrices, which reduces the time of computation significantly.
Moreover, we use the Padé approximation R 3,2 α,β to compute the value of the ML function [24,25]: with coefficients .
After simplification, formula (22) becomes e Padé approximation (24) to the ML function can be applied to the one-step and two-step ETD schemes (15) and (18) to enhance the efficiency.

Approximation Error Analysis
e approximation error of formula (24) is defined as [24] en, we have where en, we compute the error of approximation as As stated by Sarumi et al. [24], to make the approximation of R m,n α,β reliable for β ≠ α, we need to have m ≥ n + 1. is is why we use R 3,2 α,β to approximate the Mittag-Leffler function.

Numerical Experiments
We tested the ETD schemes with Padé approximation on an initial boundary value problem with analytical solutions. e numerical errors in this section are computed as e rate of convergence in time is computed as e experiments were compiled on an Intel Core i5-6200U 2.30 GHz workstation, and MATLAB R2016b was chosen as computation software.
Firstly, we consider the following FCNLSE as suggested by Esen et al. [12]: with initial conditions and homogeneous Dirichlet boundary conditions on [− 20, 20], where the parameters can be chosen as c � 0.25, and ω � − 3. e analytical solutions to FCNLSE (32) are given in [12] as Journal of Mathematics where k � − ����� � μ 2 − ω and μ 2 − ω > 0 for valid solitons. We plot the traces of numerical solutions to FCNLSE (32) with initial conditions (33) using the one-step ETD scheme (15) and the central difference method (8) for different α and β values in Figures 1-4. It can be seen from the plots that |u| 2 and |v| 2 travel in the same pace and direction.
is is due to the fact that u and v model vector solitons in a birefringent fiber. Because of the nonlinear coupling effect, the vector solitons can propagate undistorted even when the components have different widths or peak powers.
In Tables 1 and 2, the temporal convergence rates of the two-step ETD scheme (18) are computed according to formulas (30) and (31). e spatial step size is chosen as h � 0.001 which is relatively small. e experiments are performed for both α � 0.6 and α � 0.8. It can be noticed from the convergence rates that the order of convergence for α � 0.6 is around 1.6, and the order of convergence for α � 0.8 is around 1.8, which means the two-step ETD scheme (18) has a temporal order of 1 + α { }. Secondly, we solve FCNLSE (32) with initial conditions   In Figures 5-12, the evolution traces of solutions to FCNLSE (32) with initial conditions (35) are depicted with different values of α and β, using the two-step ETD scheme (18) in time and the compact scheme (10) in space. It can be observed from the mesh plots that the absolute values of u and v remain the same, which means the magnitudes of the pulses remain identical, while the real parts of u and v remain opposite to each other. It can also be seen from the evolution profiles that different α and β values result in different diffusion effects and time delay effects.
In Table 3, the computation time is recorded solving FCNLSE (32) using the two-step ETD scheme (18) taking advantage of the Padé approximation (24) for different α and β values by counting the CPU time used in MATLAB. As we           took the similar experiments without using the Padé approximation, the CPU time needed for the computation is about 15 times longer. is indicates the efficiency and necessity of the Padé approximation for the ETD schemes.

Conclusion
To solve the space-time fractional coupled nonlinear Schrödinger equation efficiently, we employed exponential time-differencing schemes for the fractional derivative in time. During this process, the Mittag-Leffler function is computed using the Padé approximation. It has been shown in the numerical experiments that the Padé approximation reduces the computational time markedly compared to the original exponential time-differencing scheme. e error of the Padé approximation to the Mittag-Leffler function has been analyzed, and the convergence rates of the schemes have been computed and demonstrated in the Numerical Experiments section. Figures 1-4 express the bright soliton solutions, and Figures 5-12 depict orthogonally polarized optical waves in a birefringent fiber. e main contribution of this paper is the modification of the exponential timedifferencing methods by applying the Padé approximation, as well as obtaining the soliton solutions to the fractional coupled nonlinear Schrödinger equation, which might be applicable in the industry of fiber optics.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.