# Efficient Exponential Time-Differencing Methods for the Optical Soliton Solutions to the Space-Time Fractional Coupled Nonlinear Schrödinger Equation

Authors: Xiao Liang; Bo Tang
Journal: Journal of Mathematics
 (2021)

---

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


---

## Body

## 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]:(1)iDtαu+Dxβu+δu2+γv2u=0,x∈ℝ,0<t≤T,iDtαv+Dxβv+δv2+γu2v=0,x∈ℝ,0<t≤T,with the initial conditions(2)u0,x=u0x,v0,x=v0x,and homogeneous Dirichlet boundary conditions on xL,xR, where i=−1 and complex functions u and v represent the amplitudes of orthogonally polarized waves in a birefringent optical fiber. Dtαu=∂αu/∂tα is the Caputo derivative of u in time, Dxβu=∂βu/∂xβ 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. 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] as(3)Dxβut,x=∂β∂xβut,x−12cosπβ/2Dxβ−∞ut,x+D+∞βxut,x,where 1<β<2. Dxβ−∞ut,x and D+∞βxut,x are the left and right Riemann–Liouville derivatives:(4)Dxβ−∞ut,x=1Γ2−β∂2∂x2∫−∞xut,ξx−ξβ−1dξ,D+∞βxut,x=1Γ2−β∂2∂x2∫x+∞ut,ξξ−xβ−1dξ,in which Γ⋅ is the gamma function.It is stated in [20] that the approximation of the left derivative Dxβ−∞vt,x is calculated using matrix BMβ:(5)vMβvM−1β…v1βv0β⊤=BMβvMvM−1…v1v0⊤,where(6)BMβ=1hβω0βω1β⋱⋱ωM−1βωMβ0ω0βω1β⋱⋱ωM−1β00ω0βω1β⋱⋱⋯⋯⋯⋱⋱⋱0⋯00ω0βω1β00⋯00ω0β,ωjβ=−1jβj,with x=jhj=0,1,…,M, where h is a single spatial step.Similarly, the approximation of the right derivativeD+∞βxvt,x is calculated using matrix LMβ:(7)vMβvM−1β⋯v1βv0β⊤=LMβvMvM−1⋯v1v0⊤.MatricesLMβ and BMβ 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]:(8)vMβvM−1β⋯v1βv0β⊤=HMβvMvM−1⋯v1v0⊤,where(9)HMβ=1hβω0βω1βω2βω3β⋯ωMβω1βω0βω1βω2β⋯ωM−1βω2βω1βω0βω1β⋯ωM−2β⋱⋱⋱⋱⋯⋯ωM−1β⋱ω2βω1βω0βω1βωMβωM−1β⋱ω2βω1βω0β,ωjβ=−1jΓβ+1cosβπ/2Γβ/2−j+1Γβ/2+j+1,j=0,1,…,M.Noticed that scheme (8) is second-order convergent in space, Ding et al. [22] generated a compact scheme to improve the order of convergence:(10)∂βvt,x∂xβ=1hββ24Δhβvt,x−h−1+β12Δhβvt,x+β24Δhβvt,x+h+Oh4=−1hβ1−β24δx2Δhβvt,x+Oh4=−1hβ1+β24δx2−1Δhβvt,x+Oh4,where δx2vt,x=vt,x−h−2vt,x+vt,x+h and −Δhβvt,x/hβ 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:(11)∂α∂tαUt+AUt=FUt,where ∂α/∂tα denotes the Caputo derivative, A is the Mx×Mx matrix in the Riesz derivative approximation, F:ℝMx⟶ℝMx contains the nonlinear function and the boundary conditions, and Ut=U1t,U2t,…,UMxtT with Ujt=uxj,t,j=1,…,Mx, and the initial condition is U0=U0.As been computed using the variation of constant formula in [23], system (11) has an analytical solution:(12)Ut=eα,1t;AU0+∫0teα,αt−s;AFUsds,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)eα,βt;λ=tβ−1Eα,β−tαλ,Eα,βz=∑k=0∞zkΓαk+β.Formula (12) can be written in a discrete form after discretization on 0,T with an equal-spaced mesh-grid tn=nτ,n=0,1,…:(14)Utn=eα,1tn;AU0+∑j=0n−1∫tjtj+1eα,αtn−s;AFUsds.Then, the ETD scheme can be denoted as [19, 23](15)Un=eα,1tn;AU0+∑j=0n−1Wn,jFUj,where Uj is the numerical approximation to Utj, and the convolution weights Wn,j can be computed as(16)Wn,j=eα,α+1tn−tj;A−eα,α+1tn−tj+1;A.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 Errn=Utn−Un satisfying(17)Utn−Un≤Cτ,n=1,…,M,where M is the temporal step number and C is a constant relating to T and α. This inequality tells us that ETD scheme (15) is first-order convergent temporally.The two-step ETD scheme is also constructed in [19]:(18)Un=eα,1tn;AU0+Wn1FU0+∑j=0n−1Wn,j2FUj−Wn,n2FUn−2+2Wn,n2FUn−1,where(19)Wn1=eα,α+2tn−1;A+eα,α+1tn;A−eα,α+2tn;A,Wn,j2=eα,α+2t1;A,n=j,eα,α+2tn−tj+1;A−2eα,α+2tn−tj;A+eα,α+2tn−tj−1;A,n>j.Garrappa and Popolizio also proved in [19] that the two-step ETD scheme (18) has the absolute approximation error Errn=Utn−Un satisfying(20)Utn−Un≤Cτ1+α,n=1,…,M,where M is the temporal step number and C is a constant relating to T and α. This inequality tells us that ETD scheme (18) is 1+α-order convergent temporally.To relief the burden of computing the functioneα,βt;A, we transform it using the multiplication of eigenvectors and functions of eigenvalues [23]:(21)fA=ZfDZ−1=Zfλ1fλ1⋱fλ1Z−1,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é approximationRα,β3,2 to compute the value of the ML function [24, 25]:(22)Eα,β−x≈Rα,β3,2x=1Γβ−αx⋅p1+xq0+q1x+x2,with coefficients(23)p1=cα,βΓβΓβ+α−Γβ+αΓ2β−αΓβ−2α,q0=cα,βΓ2βΓβ+αΓβ−α−ΓβΓβ+αΓβ−αΓβ−2α,q1=cα,βΓβΓβ+α−Γβ−αΓ2βΓβ−2α,cα,β=1Γβ+αΓβ−α−Γ2β.After simplification, formula (22) becomes(24)Eα,β−x≈Rα,β3,2x=αΓ1+α+2Γ1−α2/Γ1−2αx+Γ1−αx2.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](25)eα,β3,2x≔Eα,β−x−Rα,β3,2x,x>0.Then, we have(26)Eα,β−x=1sα,βxEα,βx,where(27)sα,βx=Γβ−αx,β>α,−Γ−αx2,β=α,Eα,βx=Eα,β0x+Oxm,asx⟶0,m≥2,β>α,3,β=α,Eα,β∞x−1+Ox−n,asx⟶∞,n≥1,β>α,2,β=α,in which(28)Eα,β0x=Γβ−αx∑k=0m−2−xkΓβ+αk,β>α,−Γ−αx2∑k=0m−3−xkΓα+αk,β=α,Eα,β∞x−1=−Γβ−αx∑k=1n−xkΓβ−αk,β>α,Γ−αx2∑k=1n−x−k+1Γ−αk,β=α.Then, we compute the error of approximation as(29)eα,β3,2x=Eα,β−x−Rα,β3,2x=1sα,βxEα,βx−pxqx=1sα,βxOx3+Ox=O1,β>α,Ox−1,β=α,asx⟶0,eα,β3,2x=Eα,β−x−Rα,β3,2x=1sα,βxEα,βx−pxqx=1sα,βxOx−2=Ox−3,β>α,Ox−4,β=α,asx⟶∞.As stated by Sarumi et al. [24], to make the approximation of Rα,βm,n reliable for β≠α, we need to have m≥n+1. This is why we use Rα,β3,2 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(30)errτ=Utn−UnL2.The rate of convergence in time is computed as(31)p=logerrτk/errτk+1logτk/τk+1.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]:(32)iDtαu+Dxβu+δu2+γv2u=0,iDtαv+Dxβv+δv2+γu2v=0,with initial conditions(33)ux,0=μ−2δ1+γsechμxβ/2β/2ei−kxβ/2β/2+p,vx,0=−μ−2δ1+γsechμxβ/2β/2ei−kxβ/2β/2+p,and homogeneous Dirichlet boundary conditions on −20,20, where the parameters can be chosen as γ=0.25, μ=0.45, δ=−0.35, p=1.5, k=−μ2−ω, and ω=−3.The analytical solutions to FCNLSE (32) are given in [12] as(34)ux,t=μ−2δ1+γsechμxβ/2β/2+2ktααei−kxβ/2/β/2+ωtα/α+p,vx,t=−μ−2δ1+γsechμxβ/2β/2+2ktααei−kxβ/2/β/2+ωtα/α+p,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 u2 and v2 travel in the same pace and direction. This 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.Figure 1
The trace of the solution tou2 of FCNLSE (32) with α=0.6 and β=1.4.Figure 2
The trace of the solution tov2 of FCNLSE (32) with α=0.6 and β=1.4.Figure 3
The trace of the solution tou2 of FCNLSE (32) with α=0.8 and β=1.6.Figure 4
The trace of the solution tov2 of FCNLSE (32) with α=0.8 and β=1.6.In Tables1 and 2, the temporal convergence rates of the two-step ETD scheme (18) are computed according to formulas (30) and (31). The spatial step size is chosen as h=0.001 which is relatively small. The 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+α.Table 1
Temporal convergence rates of the two-step ETD scheme with Padé approximation for FCNLSE (32) with α=0.6.
τerrτ for β=1.4Orderperrτ for α=1.8Orderp0.12.3309e − 3—2.2265e − 3—0.057.6231e − 41.61247.2183e − 41.62500.0252.4983e − 41.60942.3503e − 41.61880.01258.1755e − 51.61167.6287e − 51.6233Table 2
Temporal convergence rates of the two-step ETD scheme with Padé approximation for FCNLSE (32) with α=0.8.
τerrτ for β=1.4Orderperrτ for α=1.8Orderp0.11.8527e − 3—1.7532e − 3—0.055.3631e − 41.78854.9810e − 41.81550.0251.5159e − 41.82291.3967e − 41.83440.01254.2812e − 51.82413.9562e − 51.8198Secondly, we solve FCNLSE (32) with initial conditions(35)ux,0=sechx,vx,0=−sechx,and homogeneous Dirichlet boundary conditions on −10,10, where the parameters are chosen as δ=−1 and γ=1.In Figures5–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.Figure 5
The trace of the solution tou of FCNLSE (32) with α=0.6 and β=1.8.Figure 6
The trace of the solution to the real part ofu with α=0.6 and β=1.8.Figure 7
The trace of the solution tov of FCNLSE (32) with α=0.6 and β=1.8.Figure 8
The trace of the solution to the real part ofv with α=0.6 and β=1.8.Figure 9
The trace of the solution tou of FCNLSE (32) with α=0.8 and β=1.2.Figure 10
The trace of the solution to the real part ofu with α=0.8 and β=1.2.Figure 11
The trace of the solution tov of FCNLSE (32) with α=0.8 and β=1.2.Figure 12
The trace of the solution to the real part ofv with α=0.8 and β=1.2.In Table3, 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. This indicates the efficiency and necessity of the Padé approximation for the ETD schemes.Table 3
Computation time (CPU time in s) needed for solving FCNLSE (32) with initial conditions (35) via the two-step ETD scheme (18) utilizing the Padé approximation (24), taking h=0.05 and τ=0.01.
α0.30.30.30.60.60.6β1.21.51.81.21.51.8t = 0.252.4522.4452.3602.3132.3882.426t = 0.54.7364.8234.6334.5774.5794.752t = 0.757.3057.5377.3457.2927.3347.436t = 19.2389.4269.3329.4439.5639.573


## 6. 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. The 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. Figures1–4 express the bright soliton solutions, and Figures 5–12 depict orthogonally polarized optical waves in a birefringent fiber. The main contribution of this paper is the modification of the exponential time-differencing 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.



---

*Source: 5575128-2021-04-24.xml*