- Research
- Open access
- Published:
Exponential time differencing schemes for the 3-coupled nonlinear fractional Schrödinger equation
Advances in Difference Equations volume 2018, Article number: 476 (2018)
Abstract
Two modified exponential time differencing schemes based on the Fourier spectral method are developed to solve the 3-coupled nonlinear fractional Schrödinger equation. We compare the stability of the schemes by plotting their stability regions. The local truncation errors of the time integrators are proved to be fourth-order. Numerical experiments illustrating the solution to the equations with various parameters and the mass conservation results of the numerical methods are carried out.
1 Introduction
During the last few years, the nonlinear fractional Schrödinger equation (NFSE) has been widely used in modeling physical phenomena such as the propagation of waves in optics and hydrodynamics, see [1,2,3,4] for details. Various numerical methods have been developed and used for solving the fractional Schrödinger equations. For example, a series of difference schemes have been proposed by Wang and Huang [5,6,7] for the one- and two-dimensional space-fractional nonlinear Schrödinger equations. A Fourier spectral method has been developed by Duo and Zhang [8] considering the fractional parabolic equations based on periodic or Neumann boundary conditions. Various numerical methods have been proposed by Dehghan et al. [9,10,11,12,13,14] for solving the systems of one- and multi-dimensional nonlinear Schrödinger equations efficiently.
In this work, we consider the 3-coupled nonlinear fractional Schrödinger equation (3CNFSE) embracing the fractional Laplacian \((-\Delta)^{\alpha/2}\) (\(1<\alpha\leq2\)) as follows:
with the initial conditions
and periodic boundary conditions on \([x_{L}, x_{R}]\), where \(i=\sqrt{-1}\), and γ, σ, ϱ are real-constant parameters.
Although the 2-coupled NFSE was considered in a number of research articles, the 3CNFSE was rarely mentioned. In fact, the 3CNFSE is important in modeling the propagation of periodic solitary waves with perturbation in space, and this kind of model cannot be replaced by the 2-coupled NFSE. This inspired us to find some efficient numerical methods for the 3CNFSE.
The fractional Laplacian in one-dimension has been defined in [15] as
where \(1<\alpha<2\). \(_{-\infty}D_{x}^{\alpha}u(x, t)\) and \(_{x}D_{+\infty}^{\alpha}u(x, t)\) are the Riemann–Liouville derivatives to the left and right, respectively:
and
where \(\varGamma(\cdot)\) denotes the gamma function.
As shown in [16], 3CNFSE (1) conserves mass of the waves, defined as \(Q(t)\) at time t:
and the energy, or Hamiltonian,
where \(\operatorname{Re}(\phi)\) is the real part and \(\phi^{\ast}\) is the complex conjugate of a function ϕ.
The outline of this paper is organized as follows. In Sect. 2, we describe the space discretization for the 3CNFSE, which is the Fourier spectral method. In Sect. 3, we propose the description of the modified exponential time differencing schemes as time integrators for the 3CNFSE. In Sect. 4, the stability issues and rates of truncation errors to the proposed methods are analyzed and proved. In Sect. 5, numerical experiments are demonstrated and the mass conservative property of the proposed schemes is indicated. Finally, we give the conclusion in Sect. 6.
2 The Fourier spectral method
To apply the Fourier spectral method on the 3CNFSE, we follow the procedure given by Weng et al. [17], and we transform \(\psi (x,t)\) into the discrete Fourier space
where \(x_{j}=2\pi j/N\), and N is the number of Fourier nodes. The inverse of formula (5) is
The direct and inverse Fourier transforms in (5)–(6) can be computed efficiently using the fft and ifft algorithms, respectively. As indicated in [16], the fractional Laplacian in (1) can be characterized as
where \(\mathcal{F}\) is the Fourier transform acting on the spatial variable x. For periodic boundary conditions, we can approximate the fractional Laplacian using a truncated series
For Dirichlet boundary conditions, we utilize the discrete sine transform, taking \(x_{l}=a+l\Delta x\) and \(\Delta x=L/(N+1)\), where \(x\in [a, b]\), L represents the range of the spatial domain, and Δx is the space step size. For Neumann boundary conditions, we implement the discrete cosine transform, taking \(x_{l}=a+(l-1)\Delta x+\Delta x/2\) and \(\Delta x=L/N\).
3 The exponential time differential procedure
In this section, we develop two modified exponential time differencing (ETD) integrators for the 3CFNSE. We combine the ETD schemes with the Fourier spectral method discretizing in space. The exponential time differencing scheme was introduced for the stiff semi-linear problems in the form
in which the operator \(\mathcal{L}\) is linear, while the function \(\mathcal{F}\) is nonlinear. After applying the Fourier spectral method (7), we can write equation (9) into the form of a system of ODEs:
According to Duhamel’s principle, we could integrate equation (10) over a one-step interval \(t=t^{n}\) to \(t^{n+1}=t^{n}+\tau\), which gives us
Notice that we use \(u^{n}_{h}\) to represent the discrete approximation to \(u(t^{n})\). Deriving from formulation (11), a fourth-order ETD Runge–Kutta (ETDRK4) scheme was introduced by Cox and Mathew [18]. Also, we can obtain alternative ETD integrators using different methods of approximation to the integral and exponential functions. To improve the performance of the matrix exponential functions calculation, the exponential terms in the ETDRK4 scheme are computed approximately using the fourth-order \((2, 2)\)-Padé approximation or the \((1, 3)\)-Padé approximation described in [19].
The ETDRK4-P22 scheme is modified from the ETDRK4 scheme by approximating the exponential terms with the \((2, 2)\)-Padé approximation:
where
The ETDRK4-P13 scheme is modified from the ETDRK4 scheme by approximating the exponential terms with the \((1, 3)\)-Padé approximation:
where
4 Linear analysis of the ETDRK4-P schemes
4.1 Truncation error
It has been proved by Bueno-Orovio and Kay [20] that the Fourier spectral method applied on the space-fractional derivatives, involved in the 3CNFSE, is of spectral-order. So, we need to compute the overall truncation error in time. First, we investigate the ETDRK4-P22 scheme (12). To accomplish this, we take the following semi-discretized system into consideration:
where A and R denote matrices generated after we linearly discretize the NFSE in space and U represents the solution vector. By applying the ETDRK4-P22 scheme (12) to the linear system (14), we obtain
where
By Taylor expansion, Eq. (15) becomes
The exact solution of (14) is
Hence, the local truncation error of the ETDRK4-P22 scheme (12) is
This result indicates that the ETDRK4-P22 scheme (12) has a local truncation error of fourth order. The local truncation error of the ETDRK4-P13 scheme (13) can be computed in a similar way.
4.2 Stability analysis
4.2.1 Amplification symbol
As defined by Yousuf et al. [21], a rational approximation \(R_{r,s}(z)\) to the exponential \(e^{-z}\) is called A-acceptable when \(|R_{r,s}(-z)|<1\) holds for all −z with negative real part. The approximation is called L-acceptable when it is A-acceptable and it also satisfies \(|R_{r,s}(-z)|\rightarrow0\) as \(\mathfrak{R}(-z)\rightarrow-\infty\).
In Fig. 1, we compare the behavior of functions \(\exp (-z)\), \(R_{2,2}(z)\), and \(R_{1,3}(z)\), as defined in schemes (12) and (13). It can be observed from the traces that \(R_{2,2}(z)\) is A-acceptable and \(R_{1,3}(z)\) is L-acceptable. It can also be seen from Fig. 1 that \(R_{2,2}(z)\) approximates to \(\exp(-z)\) better than \(R_{1,3}(z)\) when z is close to 0.
Figure 2 illustrates the traces of \(\exp(-z)\), \(R_{2,2}(z)\), and \(R_{1,3}(z)\) for \(z=x+iy\in[0,20]\times[-10,10]\). Since the results of the functions are complex, we plot their real parts. It can be seen from the plots that \(R_{2,2}(z)\) approximates to \(\exp(-z)\) well when the absolute value of z is close to 0, while \(R_{1,3}(z)\) converges to 0 for an increasing absolute value of z.
4.3 Stability regions
The stability of the ETD schemes can be observed from the plots of their stability regions (see also [18] and [22]). Let us consider the nonlinear ordinary differential equation
where the function \(F(u)\) is nonlinear. We assume that a fixed point \(u_{0}\) satisfying \(cu_{0}+F(u_{0})=0\) exists, and u is the perturbation of \(u_{0}\). We denote \(\lambda=F'(u_{0})\), and then we linearize equation (17) to obtain
If \(\operatorname{Re}(c+\lambda)<0\), then we can say the fixed point \(u_{0}\) is stable. We denote \(x=\lambda\tau\) and \(y=c \tau\), with τ being a single time step, and then apply the ETDRK4-P22 scheme (12) to Eq. (18). The amplification factor can be calculated in the following way:
where
It can be observed from formula (19) that if \(y = 0\), the amplification factor can be computed as
which corresponds to the fourth-order explicit Runge–Kutta (RK4) method. Notice that we assumed \(r(x, y)<1\) to obtain the stability regions. Suppose that x is complex and y is chosen to be some non-positive values. As can be seen in Fig. 3, the stability regions of the two schemes are plotted. The axes of the plots are real and imaginary parts of x. It was stated by Beylkin et al. [23] that the stability region has to expand as the value of \(|ck|\) grows. It can be observed from Fig. 3 that the stability regions of the ETDRK4-P22 and the ETDRK4-P13 schemes differ in shape, which is due to the different approaches they approximate to the exponential terms.
5 Numerical experiments
In this section, the performance of the ETDRK4-P schemes is tested on a set of initial-boundary value problems. The error in the temporal direction with sufficient Fourier nodes N is computed as follows:
where \(\mathrm{U}(\tau_{s},h)\) represents the “exact” solution. The temporal convergence rates were computed as follows:
All the numerical experiments were undertaken on a Matlab R2013a platform based on an Intel Core i5-6200U 2.40 GHz workstation.
5.1 The split-step scheme
To test the efficiency of the ETDRK4-P schemes, we compare the numerical results with the split-step scheme, which is an explicit time-splitting technique [8]. The second-order split-step scheme utilizes the method of Strang splitting. From \(t=t^{n}\) to \(t^{n+1}=t^{n}+\tau\), the solution of the equation is obtained in the following steps:
where \(0\leq j\leq N-1\), \(n\geq0\), and \(\hat{u_{l}}^{n,1}\) is the lth discrete Fourier transform coefficient of the vector solution \(U^{n,1}\). Combining with the Fourier spectral method, scheme (20) is named the split-step Fourier spectral second-order (SSFS2) scheme.
Furthermore, the split-step Fourier spectral fourth-order (SSFS4) scheme is constructed by combining the SSFS2 schemes using certain weights [24]:
where \(\omega= (2+2^{1/2}+2^{-1/3})/3\), and \(\phi_{2}\) is the SSFS2 operator.
5.2 A one-dimensional 3CNFSE
We consider 3CNFSE (1) together with the initial conditions as suggested in [25]:
and periodic boundary conditions on \([x_{L}, x_{R}]\), where \(a_{0}\), \(b_{0}\), and \(c_{0}\) represent the initial wave amplitudes, the small parameter ε denotes the perturbation strength, and l is the perturbed wave number [26]. In our numerical experiments, the parameters were chosen as follows: \(\gamma_{1}=\gamma_{2}=\gamma_{3}=1\), \(\varrho=1\), \(\sigma=1\), \(a_{0}=0.2\), \(b_{0}=0.3\), \(c_{0}=0.2\), \(\epsilon =0.1\), \(l=0.5\), \(\theta=0\), \(x_{L}=-4\pi\), and \(x_{R}=4\pi\).
In Figs. 4–8, surface plots of the absolute solutions to \(u_{1}\) and \(u_{2}\) of 3CNFSE (1) with initial conditions (22), applying the ETDRK4-P22 method (12) with the Fourier spectral method (\(N=128\), \(\tau=0.01\)), are demonstrated. It can be noticed that, for different α values, the wave periods are different.
The convergence rates in Table 1 were computed using the solutions of \(u_{1}\) by the ETDRK4-P22 scheme (12) with the Fourier spectral method (\(N=256\)) for different α values. Because the analytic solution to this problem has not been discovered, we used the solution of \(\tau_{s} = 1/1600\) as the “exact” solution. It can be observed from Tables 1 and 2 that the ETDRK4-P22 scheme and the ETDRK4-P13 scheme are both fourth-order in time, because the convergence rate p is around 4 for all the α chosen. As we modified the regular ETD scheme with the Padé approximations, we do not have to compute the exponential functions again and again inside the loop, which increases the efficiency of the schemes by reducing the CPU time needed for the computation.
In Tables 3 and 4, the absolute mass error was computed according to the formula \(|(\mathrm{Q}(t^{n})- \mathrm{Q}(t^{0}))/\mathrm{Q}(t^{0})|\). Tables 3 and 4 demonstrate the mass conservative property of the ETDRK4-P22 scheme and the ETDRK4-P13 scheme, respectively. The numerical experiments were undertaken on 3CNFSE (1) with various α values. It can be observed from Tables 3 and 4 that the ETDRK4-P13 scheme (13) conserves mass of the waves better than the ETDRK4-P22 scheme (12), because the absolute mass errors of the former are smaller than the latter.
In Fig. 9, we compare the efficiency of the ETDRK4-P schemes with that of the SSFS4 scheme (21) using the Log-Log plot of CPU time and \(L_{\infty}\) errors. It can be observed from Fig. 9 that to achieve the same level of errors, the ETDRK4-P schemes need less CPU time than the SSFS4 scheme, which indicates the efficiency of the ETDRK4-P schemes.
In Table 5, the ETDRK4-P13 scheme is compared with the explicit fourth-order Runge–Kutta (ERK4) method to illustrate the stability of the scheme. We can observe from Table 5 that the solution of the ERK4 method blows up since \(t=10\), while the ETDRK4-P13 scheme remains stable with \(\tau=0.04\). The reason for this difference is that the ETDRK4-p schemes solve the nonlinear part of the equation implicitly, which improves the stability.
5.3 A two-dimensional 3CNFSE
Consider the 3CNFSE in a two-dimensional space:
where \((x,y,t)\in\varOmega_{T}=\varOmega\times[0,T]\), with initial conditions
and boundary conditions
with \(\gamma=1\), \(\rho=1\), \(d=5\), \(B_{1}=B_{3}=0.6\), \(B_{2}=0.8\), \(\varOmega=[-10, 10]^{2}\).
In Figs. 10 and 11, we plotted the solutions to system (23) with (24) using the spectral ETDRK4-P22 scheme. It can be seen from the surface plots that the three waves propagate towards the origin and they merge into one wave at the interaction. Then, the waves keep moving and separate apart. It can also be noticed from the traces that when \(\alpha=1.8\), the waves travel and disperse faster than \(\alpha=1.5\), which indicates that the value of α affects the speed of wave propagation and dispersion.
Figure 12 depicts the traces of mass errors of solutions to system (23) with (24) using the spectral ETDRK4-P22 scheme. It can be observed from Fig. 12 that the proposed spectral ETDRK4-P22 scheme conserves the mass of the system well, because the mass errors are kept to the level of 10−11 despite the oscillation at the time of the wave interaction.
5.4 A three-dimensional 3CNFSE
Consider the three-dimensional initial-boundary-value problem:
where \((x,y,z,t)\in\varOmega_{T}=\varOmega\times[0,T]\), with initial conditions
with \(\gamma=1\), \(\rho=1\), \(d=5\), \(B_{1}=B_{2}=B_{3}=20\), \(\varOmega =[-20, 20]^{3}\), \(T=3\), and periodic boundary conditions on ∂Ω.
In physics, 3CNFSE (26) can be used to model the dispersion of quantum energy density. In Figs. 13 and 14, the solutions to 3CNFSE (26) using the spectral ETDRK4-P13 scheme are depicted. To better illustrate the merge of the energy density, we set \(z=0\) and plotted the solutions in the x–y plane. It can be concluded from the two figures that the dispersion effect increases as the value of α grows larger.
Figure 15 depicts the trace of energy errors of solutions to system (26) with (27) using the spectral ETDRK4-P13 scheme. It can be observed from Fig. 15 that the proposed spectral ETDRK4-P13 scheme conserves the energy of the system well, because the energy errors are kept to the level of 10−14. We computed the energy error using the formula \((\mathrm{E}(t)-\mathrm {E}(0))/\mathrm{E}(0)\), where \(\mathrm{E}(t)\) is defined in (4).
6 Conclusion
In this work, two modified exponential time differencing Runge–Kutta schemes, using the Padé approximation to the matrix exponential functions, have been developed for the 3-coupled space-fractional nonlinear Schrödinger equation. The stability issues are discussed by computing the amplification factors and plotting the stability regions. Local truncation errors are calculated to indicate the accuracy of the proposed schemes. Numerical experiments are undertaken on the equations with various α values. The results indicate that the proposed schemes conserve mass of the waves well.
References
Longhi, S.: Fractional Schrödinger equation in optics. Opt. Lett. 40, 1117–1120 (2015)
Antoine, X., Tang, Q., Zhang, Y.: On the ground states and dynamics of space fractional nonlinear Schrödinger/Gross–Pitaevskii equations with rotation term and nonlocal nonlinear interactions. J. Comput. Phys. 325, 74–97 (2016)
Fu, C., Lu, C.N., Yang, H.W.: Time-space fractional \((2 + 1)\) dimensional nonlinear Schrödinger equation for envelope gravity waves in baroclinic atmosphere and conservation laws as well as exact solutions. Adv. Differ. Equ. 2018, 56 (2018)
Petroni, N.C., Pusterla, M.: Lévy processes and Schrödinger equation. Physica A 388, 824–836 (2009)
Wang, P., Huang, C.: A conservative linearized difference scheme for the nonlinear fractional Schrödinger equations. Numer. Algorithms 69, 625–641 (2015)
Wang, P., Huang, C.: An energy conservative difference scheme for the nonlinear fractional Schrödinger equations. J. Comput. Phys. 293, 238–251 (2015)
Wang, P., Huang, C.: Split-step alternating direction implicit difference scheme for the fractional Schrödinger equation in two dimensions. Comput. Math. Appl. 71, 1114–1128 (2016)
Duo, S., Zhang, Y.: Mass-conservative Fourier spectral methods for solving the fractional nonlinear Schrödinger equation. Comput. Math. Appl. 71, 2257–2271 (2016)
Dehghan, M., Abbaszadeh, M., Mohebbi, A.: Numerical solution of system of n-coupled nonlinear Schrödinger equations via two variants of the meshless local Petrov–Galerkin (MLPG) method. Comput. Model. Eng. Sci. 100, 399–444 (2014)
Taleei, A., Dehghan, M.: Time-splitting pseudo-spectral domain decomposition method for the soliton solutions of the one- and multi-dimensional nonlinear Schrödinger equations. Comput. Phys. Commun. 185, 1515–1528 (2014)
Mohebbi, A., Abbaszadeh, M., Dehghan, M.: The use of a meshless technique based on collocation and radial basis functions for solving the time fractional nonlinear Schrödinger equation arising in quantum mechanics. Eng. Anal. Bound. Elem. 37, 475–485 (2013)
Dehghan, M., Manafian, J., Saadatmandi, A.: Solving nonlinear fractional partial differential equations using the homotopy analysis method. Numer. Methods Partial Differ. Equ. 26, 448–479 (2010)
Dehghan, M.: Finite difference procedures for solving a problem arising in modeling and design of certain optoelectronic devices. Math. Comput. Simul. 71, 16–30 (2006)
Dehghan, M., Fakhar-Izadi, F.: The spectral collocation method with three different bases for solving a nonlinear partial differential equation arising in modeling of nonlinear waves. Math. Comput. Model. 53, 1865–1877 (2011)
Yang, Q., Liu, F., Turner, I.: Numerical methods for fractional partial differential equations with Riesz space fractional derivatives. Appl. Math. Model. 34, 200–218 (2010)
Wang, D., Xiao, A., Yang, W.: A linearly implicit conservative difference scheme for the space fractional coupled nonlinear Schrödinger equations. J. Comput. Phys. 272, 644–655 (2014)
Weng, Z.F., Zhai, S.Y., Feng, X.L.: A Fourier spectral method for fractional-in-space Cahn–Hilliard equation. Appl. Math. Model. 42, 462–477 (2017)
Cox, S.M., Matthews, P.C.: Exponential time differencing for stiff systems. J. Comput. Phys. 176, 430–455 (2002)
Bhatt, H.P., Khaliq, A.Q.M.: Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equation. Comput. Phys. Commun. 200, 117–138 (2016)
Bueno-Orovio, A., Kay, D.: Fourier spectral methods for fractional-in-space reaction–diffusion equations. BIT Numer. Math. 54, 937–954 (2014)
Yousuf, M., Khaliq, A.Q.M., Kleefeld, B.: The numerical approximation of nonlinear Black–Scholes model for exotic path-dependent American options with transaction cost. Int. J. Comput. Math. 89, 1239–1254 (2012)
Liang, X., Khaliq, A.Q.M., Bhatt, H., Furati, K.M.: The locally extrapolated exponential splitting scheme for multi-dimensional nonlinear space-fractional Schrödinger equations. Numer. Algorithms 76, 939–958 (2017)
Beylkin, G., Keiser, J.M., Vozovoi, L.: A new class of time discretization schemes for the solution of nonlinear PDEs. J. Comput. Phys. 147, 362–387 (1998)
Hederi, M., Islas, A.L., Reger, K., Schober, C.M.: Efficiency of exponential time differencing schemes for nonlinear Schrödinger equations. Math. Comput. Simul. 127, 101–113 (2016)
Aydin, A.: Multisymplectic integration of n-coupled nonlinear Schrödinger equation with destabilized periodic wave solutions. Chaos Solitons Fractals 41, 735–751 (2009)
Qian, X., Song, S., Chen, Y.: A semi-explicit multi-symplectic splitting scheme for a 3-coupled nonlinear Schrödinger equation. Comput. Phys. Commun. 185, 1255–1264 (2014)
Funding
This work was supported by Hubei Provincial Department of Education (B2018158); the Project of Hubei University of Arts and Science (XK2018033).
Author information
Authors and Affiliations
Contributions
The authors declare that the study was realized in collaboration with the same responsibility. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Liang, X., Bhatt, H. Exponential time differencing schemes for the 3-coupled nonlinear fractional Schrödinger equation. Adv Differ Equ 2018, 476 (2018). https://doi.org/10.1186/s13662-018-1936-9
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s13662-018-1936-9