Abstract

The 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. The 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 time-differencing methods.

1. Introduction

The coupled nonlinear Schrödinger equation (CNLSE) can be employed in simulating the propagation of the optical soliton in a birefringent fiber [1–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. The 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–8]. It is also shown in [9, 10] that the anomalous dispersion or diffusion effects are better described by the spatial fractional derivatives. The anomalous effects reflect the Lévy-type particle movement, different from Brownian motion, which depicts the classical random movement of particles. Therefore, 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 conditionsand homogeneous Dirichlet boundary conditions on , where and complex functions and represent the amplitudes of orthogonally polarized waves in a birefringent optical fiber. is the Caputo derivative of in time, is the Riesz derivative of in space, , , 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. Then, 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.

2. Discretization in Space

The spatial Riesz derivative is defined in [10] aswhere . and 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 derivative is calculated using matrix :wherewith , where is a single spatial step.

Similarly, the approximation of the right derivative is calculated using matrix :

Matrices and 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 and is the second-order approximation (8). As been proved by Theorem 11 in [22], compact scheme (10) is fourth-order convergent spatially.

3. The Exponential Time Integrator

We obtain a system of time-fractional equations after discretizing FCNLSE (1) in space:where denotes the Caputo derivative, is the matrix in the Riesz derivative approximation, contains the nonlinear function and the boundary conditions, and with , and the initial condition is .

As been computed using the variation of constant formula in [23], system (11) has an analytical solution:where denotes the inverse function of the Laplace transform to , and can be calculated taking advantage of the Mittag-Leffler (ML) function :

Formula (12) can be written in a discrete form after discretization on with an equal-spaced mesh-grid :

Then, the ETD scheme can be denoted as [19, 23]where is the numerical approximation to , and the convolution weights 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 satisfyingwhere is the temporal step number and is a constant relating to and . This inequality tells us that ETD scheme (15) is first-order convergent temporally.

The two-step ETD scheme is also constructed in [19]:where

Garrappa and Popolizio also proved in [19] that the two-step ETD scheme (18) has the absolute approximation error satisfyingwhere is the temporal step number and is a constant relating to and . This inequality tells us that ETD scheme (18) is -order convergent temporally.

To relief the burden of computing the function , we transform it using the multiplication of eigenvectors and functions of eigenvalues [23]:where is diagonalizable, with ’s to be its eigenvalues, is the composition of ’s eigenvectors, and . 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 to compute the value of the ML function [24, 25]:with coefficients

After simplification, formula (22) becomes

The 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.

4. Approximation Error Analysis

The approximation error of formula (24) is defined as [24]

Then, we havewherein which

Then, we compute the error of approximation as

As stated by Sarumi et al. [24], to make the approximation of reliable for , we need to have . This is why we use to approximate the Mittag-Leffler function.

5. Numerical Experiments

We tested the ETD schemes with Padé approximation on an initial boundary value problem with analytical solutions. The numerical errors in this section are computed as

The rate of convergence in time is computed as

The 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 conditionsand homogeneous Dirichlet boundary conditions on , where the parameters can be chosen as , , , , , and .

The analytical solutions to FCNLSE (32) are given in [12] aswhere and 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 and travel in the same pace and direction. This is due to the fact that and 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.

Figure 1 

The trace of the solution to of FCNLSE (32) with and .

Figure 2 

The trace of the solution to of FCNLSE (32) with and .