On the Numerical Simulation of Time-Space Fractional Coupled Nonlinear Schrödinger Equations Utilizing Wendland's Compactly Supported Function Collocation Method

This research describes an efficient numerical method based on Wendland’s compactly supported functions to simulate the time-space fractional coupled nonlinear Schrödinger (TSFCNLS) equations. Here, the time and space fractional derivatives are considered in terms of Caputo and Conformable derivatives, respectively. The present numerical discussion is based on the following ways: we first approximate the Caputo fractional derivative of the proposed equation by a scheme order O(∆t2−α), 0 < α < 1 and then the Crank-Nicolson scheme is employed in the mentioned equation to discretize the equations. Second, applying a linear difference scheme to avoid solving nonlinear systems. In this way, we have a linear, suitable calculation scheme. Then, the conformable fractional derivatives of the Wendland’s compactly supported functions are established for the scheme. The stability analysis of the suggested scheme is also examined in a similar way to the classic Von-Neumann technique for the governing equations. The efficiency and accuracy of the present method are verified by solving two examples.


Introduction
Our main intention of this current investigation is to simulate the TSFCNLS equations where 0 < α ≤ 1, 1 < β ≤ 2, i = √ −1, γ, ρ and η are known constants. u(x, t) and v(x, t) are complex valued functions. α and β are the fractional order of the time and spatial derivatives, respectively. In case α = 1 and β = 2, this system represents classical coupled nonlinear Schrödinger (CNLS) equations. When η = 0, the system (1.1)-(1.2) reduces to single time and space fractional nonlinear Schrödinger (TSFNLS) equation. If we choose ρ = 0, it is decoupled and turns into the time and space fractional Schrödinger (TSFS) equation of free particles. In the current paper, we choose the parameters η = γ = ρ = 1 in the equation system.
To investigate several physical situations especially in optics, hydrodynamics, CNLS equations have been used by many researchers [5,6,14,15]. For instance, Manakov [15] used integrable or nearly integrable systems of two CNLS equations to solve the nonlinear phenomenon about the optical solitons that store and transfer information in polyfibrous optical media. Nonlinear Schrödinger equations can be considered to verify the modulation instability of gravity waves in fluids at finite or great depths in hydrodynamics. Furthermore, numerous computational and effective analytical and numerical methods have been proposed to investigates the nonlinear Schrödinger equation and CNLS equation [1,10,11,12,13,33].
In recent years, fractional partial differential equations play an important role in physics. The growing number of investigation for the fractional NLS and CNLS equations has attracted the researchers. That's why these equations have been studied and obtained important achievements. Laskin generated that the fractional NLS via using a new path integral over Lévy-like quantum mechanical paths. The fractional Schrödinger equation include space fractional derivative instead of the second-order derivative in the classical Schrödinger equation. Space fractional NLS or CNLS can modify the shape of the wave, while the nonlinearity and dispersion effects remain unchanged. Thus, there are numerous works on this topic recently. Later, some physical applications of the fractional Schrödinger equations with time-space fractional derivatives are studied by the authors [8,30]. Besides, the finite difference method is constructed for time-space fractional Schrödinger equations by Liu, Zheng and Li [19]. Another study is done by [22]. They used a linearized Crank-Nicolson method to obtain the numerical solution of the time-space fractional NLS equations. Also, Bhrawy et al. [3] used a Jacobi spectral collocation method for solving fractional Schrödinger equations and coupled fractional Schrödinger systems. Meanwhile, the papers [21,29] used a conservative difference method to obtain the numerical solutions of space fractional CNLS equations. The Crank-Nicolson difference scheme for the CNLS with the Riesz space fractional derivative is proposed by Wang et. al. [28]. Another study is given by Mei et. al. [31]. They used a spectral Galerkin method for the space-fractional CNLS equations.
The main contribution of this study is to display the proposed method as an alternative approach to solve numerically these type fractional partial differential equations. To our knowledge, TSFCNLS equations are not solved by using the collocation method based on Wendland's compactly supported functions. Therefore, we use the proposed method estimates the solution of the mentioned equation system. Meanwhile, we recognize that Riesz space fractional derivative has been generally used for the space fractional CNLS equations in the literature. In this study, the conformable fractional derivative is preferred for space fractional derivative in TSFCNLS equations. Because a relationship between the conformable derivative and classical derivative provides ease of implementation.
The layout of this paper is as follows: in Section 2, the most popular fractional derivatives are given. The next Section 3 describes the Wendland's compactly supported functions and its properties. In Section 4, we firstly construct a Crank-Nicolson difference scheme for the discretization of the equation system. Secondly, the equation system is linearized. In Section 5, we construct the implementation of the suggested method for the equation system. We investigate the stability of the current scheme in Section 6. Finally, numerical results are illustrated in Section 7 and end the study with a short conclusion given.

Preliminaries of the fractional derivatives
A fractional derivative of a function u(x, t) can be described in several ways. Several definitions of fractional derivatives due to Riemann-Liouville, Caputo, Grünwald-Letnikov and Conformable can be found in the literature [16,17,18,20]. Since the applicability of the standard initial and boundary conditions to be involved in the formulation of the problems, we consider Caputo and Conformable derivative for the time and spatial derivatives, respectively.
where the function Γ (.) is called a Gamma function.
This fractional derivative was introduced by the Italian mathematician Caputo in 1967 [4]. Sun and Wu [25] proposed a useful approximation to compute this fractional derivative to t. This approximation is given in their paper [25]. Let us state in the following lemma.
We will use the following expression for the discretization of the Caputo type derivative.
. Now, let us introduce the Conformable fractional derivative. The Conformable derivative of order α, 0 < α < 1, which is based on the main limit definition of derivative, is stated in the following definition.
6. If f is a differentiable function, there is a relationship between conformable derivative and classical derivative such as T α (f )(t) = t 1−α df dt (t). As the equation system has a fractional derivative of order α, 1 < α < 2 in our study, we need to use a fractional derivative definition whose order α, α ∈ (1, 2]. Khail et. al. [20] proposed the following definition and remark in their paper.

Definition 4.
[20] Let α ∈ (n, n + 1] and f be an n-times differentiable function for t > 0. The conformable derivative of f order α is defined by where α is the smallest integer greater than or equal to α.
Remark 2. The truncation error of the approximation for the space derivative by using the conformable derivative is given as

Wendland's compactly supported functions
In the literature, there are a lot of basis functions. In this investigation, we use the Wendland functions, which were originally described in Wendland [32]. The Wendland functions are piecewise polynomial compactly supported. They have minimum polynomial degree for any level of smoothness and are positive definite because of being strictly positive Fourier transform. They are given in the following as: φ l,k (r) = (1 − r) n + p l,k (r), with the following conditions where p is described polynomial for k ≥ 1 and l is the dimension number. This condition displays that all the Wendland functions are equal to zero outside [0, 1]. Also, the choice of l ensures that the obtained functions are positive definite. In our algorithms, we use the Wendland functions (see the details in [32]) which are defined as follows: , W 6,4 (r) = (1 − r) 10 + (5 + 50r + 210r 2 + 450r 3 + 429r 4 ), W 7,5 (r) = (1 − r) 12 + (9 + 108r + 566r 2 + 1644r 3 + 2697r 4 + 2048r 5 ).
For ease notation in the rest of the paper, φ l,k will be used as W l,k . Since the Wendland functions have compact support, the interpolation matrices are sparse and for the evaluation of the interpolants, only a few terms have to be considered. This provides efficient algorithms for the computation [23].
That's why, the Wendland functions are first preferred to obtain approximate solutions of the TSFCNLS equations. As we have already noted, Conformable fractional is used for the space-fractional derivative of the equation system in this paper. To our knowledge, for computing the Conformable derivative of a function f (t), the function must be defined for all t > 0. Meanwhile, the basis function W l,k (r) can be scaled to have compact support on [0, δ] by replacing r with r δ for δ > 0. Here δ is called a scaling factor. Therefore, the calculations are made with the Wendland's compactly supported functions in our work.
Note that finding the optimal size of the scaling factor δ is still an unsolved problem. In this investigation, to determine the optimal sizes of the scaling factor δ condition number estimation of the kernel matrix is used [24]. This procedure is based on two parts: first, we can guess scale, and then a loop is created for a good scaling number by using a condition number estimation of the kernel matrix.

Solution of the TSFCNLS equations
This section first presents the discretization of the equation system. Since the Equations (1.1)-(1.2) include complex valued functions u and v, they are composed in its real and imaginary parts by writing (4.1) By using (4.1), we obtain following equalities Substituting these above equations into the Equations (1.1)-(1.2), we have 2) The equation system (4.2) can be considered as four partial differential equations: The initial and boundary conditions of these equations are given The matrix-vector form of the equation system (4.3) can be written as follows: ∂t α denotes the Caputo fractional derivative of U (x, t) with respect to t. Using the Lemma 2, we compute the approximation of this fractional derivative as the following: . Let's apply first Caputo fractional derivative which are given Equation (2.1), and then Crank-Nicolson method to the equation (4.4)-(4.6), we get where ∂ β Ui ∂|x| β n denotes the n-level of discretization and (u (v Substituting Equations (4.10)-(4.11) into the Equation (4.9), we obtain In order to linearize the nonlinear terms in the mentioned equation system, we first apply the Caputo fractional derivative to ∂ α ur ∂t α , ∂ α us ∂t α , ∂ α vr ∂t α ve ∂ α vs ∂t α . Then we have As in [7], we suppose that and substituting Equation (4.12) into Equation (4.8), we get a following linear equation as

Governing of the proposed method to the TSFCNLS equations
This part illustrates the application of Wendland's compactly supported function to TSFCNLS equations which are given in Equations (1.1)-(1.2). Now, we approximate U (x i , t n ) and fractional derivative with respect to space variables at each time step can be expressed as: where T β is a conformable fractional operator with respect to the variable x. Substituting the Equations (5.1)into the Equation (4.13) and rearrange, we get where A k and B k are both N × N matrices and The first and last row of these matrices consist entirely of φ functions, the N −2 rows consist (N − 2) × (N − 2) matrices C k and D k . For k = 1, 2, the matrices C k and D k are given as: and G k = [0 H k 0] T , where for k = 1, 2, 3, 4, H k are given as follows: 4. The initial value of λ n is obtained by using the initial condition (4.5).
Then substituting this value in the Equation (5.1), the value of U 0 can be computed.
7. The approximate value of U n+1 can be computed by using the Equation (5.1).

By calculating
, the envelope solutions of u and v at the time step (n + 1). 9. Iteration is continued until n∆t = T .

Von-Neumann stability technique
The stability analysis of the present scheme is examined by using a technique similar to the traditional Von-Neumann analysis method. Since this stability technique applies only linear differential equations, nonlinear term s in Equation (4.4) is changed by a constant q. The Equation (4.4) can be written as: Here ∂ α ∂t α represents Caputo type derivative according to t, (4.7) is rearranged and we get . . , n − 1. By substituting Equation (6.2) into Equation (6.1) and then Crank-Nicolson method is applied to Equation (6.1), we get We have after applying collocation method based on Wendland's compactly supported functions Equation (6.3) As in the von Neumann stability analysis method, suppose that where i 2 = −1, θ ∈ R, Υ ∈ R 4×1 and ξ ∈ R 4×4 is the amplification matrix. Substituting (6.5) into (6.4), we have the following expression: After some manipulations we get where I is the unit matrix and The necessary condition is |λ j | ≤ 1, j = 1, 2, 3, 4, where λ j are the eigenvalues of the amplification matrix ξ. We consider the time independent limit value ξ = I as in [26], then So, Equation (6.6) can be written as The eigenvalues of the amplification matrix ξ are |λ 1,3 | = 1 and |λ 2,4 | = 1, so the scheme is unconditionally stable.

Numerical results
In this section, two test problems are used for TSFCNLS equations to indicate the efficiency of the proposed method. In order to examine the temporal convergence order, we use the following formula Order = log (||error(∆t 1 )||/||error(∆t 2 )||) log (∆t 1 /∆t 2 ) .

Test Problem I
In this section, we present some numerical results. Let's consider an initial and boundary value problem in the form When α = 1 and β = 2, the exact solution of this problem which is derived in [27] is given by On the contrary, the numerical outcomes and exact solutions are compatible for classical order, in case of α ∈ (0, 1) and β ∈ (1, 2), there are no exact solutions in the literature. For this reason, an error estimates method in [9] is used and the computed numerical solutions for a very small time and space step size ∆t, ∆x are accepted as reference solutions. These reference solutions are obtained by taking the time and space step size ∆t = 0.00125, ∆x = 0.09375, and then the error estimates are obtained different choices of ∆t, ∆x and the fractional derivative orders α and β. Since numerical results of v(x, t) is similar to numerical results of u(x, t), the results of u(x, t) are only given in this paper. First, the discrete conservation laws of the suggested scheme are discussed. Table 1 illustrates the invariants of u(x, t) for different choices Wendland functions with α = 1, β = 2 and T = 1. The analytical values of invariants are calculated as 0.8944271909 and 0.3279382901. It is easy to find that physical invariants are preserved at very high accuracy. Also, the values of mass and energy for different selected of α and β are depicted at T = 1 in Table 2. We can see that the scheme preserves the discrete masses. On the other hand, Table 3 indicates the changes in amplitudes and peak positions for different values of α and β. It is seen that the changes of amplitudes are 0.1787%, 0.8048% and 3.35% for the values of order α = 1, 0.95, 0.75 and β = 2, 1.95, 1.75 respectively. It has to be noted that when α and β become smaller, the motion of the wave falls back and the amplitudes increase quickly but the peak position of the wave shifts to left. This situation is expected due to the use of the fractional Schrödinger equations in physics to adapt the shape of a wave without the switch of the nonlinearity and dispersion effects. Table 4 illustrates the error norms in the sense of L ∞ for different selections of α and β. It is observed that the maximal value of the error norm is about 10 −2 and minimal value is about 10 −4 . This means that the accuracy of the method is very reliable. Errors, convergence orders for the scheme is listed in Table 5 for different α with β = 1.5. The numerical results indicate that the order of temporal convergence is about O(∆t 2−α ), 0 < α < 1. In Figure 1, the kernel matrix is presented.

Test Problem II
In this test problem, we take the initial conditions for interaction of two waves in the form                   When α = 1 and β = 2, the system is completely integrable and so the collision of the soliton waves is elastic. In the rest of the values of α and β, we can conclude that the shape of the solitons will change more quickly as α, β increases and the waveforms become similar to the case of α = 1 and β = 2. Also, it is observed that the waves go out after collision in Figures 16-25.
It should be noted that the collision will occur earlier with the growing of α and β. We would also like to say that the time fractional order α greatly affects the collision. In addition, the authors in [28] have emphasized that the collision of the solitons are not elastic any more for the values α = η = 1 and β = 2. It is reverified by Figures 12,14,16,18. Finally, the invariants of numerical solution of u(x, t) are listed in Table 6.

Conclusions
We have proposed a method based on Wendland's compactly supported functions to solve numerically the TSFCNLS equations. Caputo and Comformable fractional derivative are used for the time and space fractional derivative, respectively. Numerical outcomes are presented Tables 1-6 and Figures 1-25 for different choices of α and β. We would also like to say that the approximation solutions are in the ideal agreement exact solution in case of α = 1 and β = 2. Another point that is worthy of being emphasized is that when αandβ tends to 1 and 2, respectively, the numerical solutions of the mentioned equation are convergent to the solutions of the coupled NLS. Besides, the stability analysis of the scheme is examined by a similar procedure to the classic von Neumann analysis. It is seen that the proposed method is unconditionally stable. Some of the major advantages of the numerical technique are geometrical flexibility, high-order computational accuracy. As a final remark, we would like to mention that the suggested method is an efficient alternative way to solve these type fractional differential equations because of the very simple and easy implementation of the scheme.