An Efﬁcient Discrete Model to Approximate the Solutions of a Nonlinear Double-Fractional Two-Component Gross–Pitaevskii-Type System

: In this work, we introduce and theoretically analyze a relatively simple numerical algorithm to solve a double-fractional condensate model. The mathematical system is a generalization of the famous Gross–Pitaevskii equation, which is a model consisting of two nonlinear complex-valued diffusive differential equations. The continuous model studied in this manuscript is a multidimensional system that includes Riesz-type spatial fractional derivatives. We prove here the relevant features of the numerical algorithm, and illustrative simulations will be shown to verify the quadratic order of convergence in both the space and time variables.


Introduction
There have been dramatic developments in the area of fractional calculus in recent decades [1], and many areas in applied and theoretical mathematics have benefited from these developments [2,3]. In particular, there have been substantial developments in the theory and application of numerical methods for fractional partial differential equations. For example, from a theoretical point of view, theoretical analyses of conservative finitedifference schemes to solve the Riesz space-fractional Gross-Pitaevskii system have been proposed in the literature [4], along with convergent three-step numerical methods to solve double-fractional condensates, explicit dissipation-preserving methods for Riesz space-fractional nonlinear wave equations in multiple dimensions [5], energy conservative difference schemes for nonlinear fractional Schrödinger equations [6], conservative difference schemes for the Riesz space-fractional sine-Gordon equation [7], high-order central difference schemes for Caputo fractional derivatives [8], among other examples.
It is important to point out that most of the methods mentioned above refer to discretizations for partial differential equations with fractional derivatives in space. In particular, those models consider fractional partial derivatives of the Riesz type. It is worth noting that the Riesz derivatives are linear combinations of the left and right Riemann-Liouville operators. other approaches strive to provide discretizations of systems of partial differential equations with Caputo-type derivatives in time. As mentioned above, there are reports on high-order central difference schemes for Caputo fractional derivatives [8], modified integral discretization schemes for two-point boundary value problems with Caputo fractional derivatives [9], predictor-corrector schemes for solving nonlinear delay differential equations of fractional order [10], numerical integrators for Caputo fractional differential equations with infinity memory effect at initial conditions [11], second-order schemes for the fast evaluation of the Caputo fractional derivatives [12], and homotopy perturbation methods for solving the Caputo-type fractional-order Volterra-Fredholm integro-differential equations [13].
Various systems mentioned above are capable of preserving some physical quantities, such as mass and energy; the development of numerical schemes that preserve these features is an important area of research. Historically, there have been many reports in this area, including systems of integer-order partial differential equations. For example, there are reports on the numerical solutions of conservative nonlinear Klein-Gordon [14] and sine-Gordon [15,16] equations, symplectic methods for the Schrödinger equation [17], fast and structure-preserving schemes for partial differential equations based on the discrete variational derivative method [18], structure-preserving numerical methods for partial differential equations [19], dissipative or conservative Galerkin methods using discrete partial derivatives for nonlinear evolution equations [20], among other reports [21]. These approaches have been extended to the fractional-case scenario, and have been helpful in designing numerical models that are able to preserve the mass and the energy of nonlinear systems [22][23][24].
Motivated by these developments, the present paper presents an efficient discrete model to approximate the solutions of a nonlinear double-fractional two-component Gross-Pitaevskii system. The system consists of two coupled complex-valued functions, whose dynamics are described by parabolic partial differential equations with nonlinear reactions. We consider here spatial derivatives of fractional order in the Riesz sense. The mathematical model is a complicated system, and the need to approximate its solutions is an interesting task. We propose here an easy-to-implement scheme to approximate the solutions and fully analyze them theoretically. We establish the properties of stability and convergence using a discrete form of Grönwall's inequality. Some simulations will be provided to illustrate the performance of the scheme, and a numerical test of the convergence is provided. Here, it is worth noting that the main advantage of the present discretization is that the computer implementation is easy. Moreover, the scheme is quadratically convergent, and the experiments show that such is the case. Other discretizations [4,25] are more difficult to implement computationally, in view that a fixed point technique should be coded along with the computational algorithm. Moreover, those schemes provide more complicated conditions in order to guarantee the convergence. Finally, one of the limitations of the present methodology is that, to the best of our knowledge, it is not capable of preserving the energy of the system.
Let p ∈ N and T ∈ R + . Define I n = {k ∈ N : k ≤ n}, denote {0} ∪ I n by I n , for each n ∈ N, and let a i , b i ∈ R be such that and Ω T = Ω × (0, T), and let i = √ −1. Agree that ψ 1 : Ω T → C and ψ 2 : Ω T → C, and define x = (x 1 , . . . , x p ) ∈ Ω. Functions defined on Ω will be extended to all of R p by letting them be equal to zero outside Ω.

Numerical Algorithm
For the remainder, we will let N and M i be natural numbers, with i ∈ I p . Define ). Finally, let ∂J represent the collection of all j ∈ J with the property that x j ∈ ∂Ω. Definition 4. If w = u, v and α ∈ (0, 1) ∪ (1, 2], we define the averages and the differences Moreover, we agree that x p w n j ). (16) With this nomenclature, we introduce the following numerical algorithm to approximate the solution of (4) on Ω × [0, T], for each (j, n) ∈ J × I N−1 : It is clear that, if u n , u n−1 , v n and v n−1 are known for some n ∈ I N−1 , then the resulting iterative formulas wield expressions where the only unknowns are u n+1 and v n+1 . On the other hand, u 0 and v 0 are prescribed by the initial data.
In order to determine the approximations when the time is equal to t 1 , we employ the initial conditions µ (1) As a consequence, we obtain and For the sake of convenience, we let h = (h 1 , . . . , h p ), and denote with V h the complexvalued functions whose domains are {x j : j ∈ J}. Moreover, set w j = w(x j ). Theorem 1. For all initial conditions, the system (17) has a unique solution.
Proof. We provide the proof for the one-dimensional case only, the higher dimensional scenario being similar. Beforehand, notice that u 0 and v 0 are defined by the initial conditions, and so are u 1 and v 1 , as shown by the identities (20) and (21), respectively. Suppose that u n , u n−1 , v n and v n−1 are obtained. The first numerical model is expressed alternatively in the following form: for each (j, n) ∈ J × I N−1 . Here, As a consequence, the first equation of the system is rewritten as where For the existence of u n+1 , it suffices show that A is invertible [30]. Observe that Also, the following hold: By the previous lemma, A is non-singular, which means that there exists a unique solution u n+1 of (25). The existence of v n+1 is proved similarly, noting that the second recursive equation of (17) is rewritten by where c(u n j , v n j ) = λu n j + (β 12 |u n j | 2 + β 22 |v n j | 2 )v n j , Moreover, B is the square complex matrix defined by which is also strictly diagonally dominant. The conclusion follows now by induction.
The discretization proposed in this work is similar to a linear implicit discretization of (4) in various senses. Indeed, notice that our approach hinges on approximating the nonlinear term at the time t n , while the linear terms are approximated at the time t n+1 . The difference is that the linear term of the numerical model (17) is approximated by the average of the numerical solutions at the levels n + 1 and n − 1 through µ (1) t u n j and µ (1) t v n j . In that sense, the present discretization would seem computationally more complex than the linear implicit scheme. In this point, we would like to clarify that the linear implicit scheme has the advantage of being a two-step method, but the computational implementation would require solving systems of linear equations similar to those associated to the discrete model (17). On the other hand, as we will see in the following section, our current discretization has convergence of the second order in space, while the corresponding order of the linear implicit scheme is known to be linear.
Proof. Suppose that the functions u, v, and V are sufficiently smooth. Then there exist real numbers C 1,k , C 2,k,i , and C 3,k , such that The conclusion follows from the triangle inequality and Taylor's theorem.
In the following, it is worth recalling that the square-root operator of − (α) h is the discrete fractional operator It is easy to check that, if α ∈ (1, 2] and w = (w n ) N n=1 ⊆ V h , then Im iδ (1) (1) These identities will be employed in the proofs of Theorems 3 and 4. Similarly, the following result will be crucial in those proofs.
Proof. It is easy to check that Im n , ζ n+1 + n−1 , ζ n + ζ n , n+1 + ζ n−1 , n which is what we wanted to prove.

Lemma 2.
For each n ∈ I N , let n = u n −ũ n and ζ n = v n −ṽ n . Agree that for each (j, n) ∈ J × I N−1 . Then there exists a constant C ≥ 0, which dependents on τ and h, such that max{|P n j |, |Q n j |} ≤ C(| n j | + |ζ n j |), for each (j, n) ∈ J × I N−1 .
Algebraic calculations establish that t ζ n , ∀n ∈ I N−1 .

Proof.
The proof is similar to that of the stability property of (17).
Proof. The proof of this result is similar to that of Theorem 3 and, for that reason, we provide here a shortened proof only. Beforehand, let n = u n − U n and ζ n = v n − V n , for each n ∈ I N . It is obvious then that the pair ( , ζ) satisfies the discrete problem where for each (j, n) ∈ J × I N−1 . Proceeding now as in the proof of Theorem 3, there exists a constant C ≥ 0, such that (50) is satisfied. Following the same steps used to obtain (51) and (52), we may readily obtain the following identities, valid for each n ∈ I N−1 : Again, as in the proof of Theorem 3, we add these last two equations and take the sum on both sides over all n ∈ I m , for some m ∈ I N−1 . Using telescoping sums, rearranging terms, taking then absolute values, using Young's inequality and (50), and applying the inequality of the hypotheses, it follows that there is a constant C ≥ 0, such that for sufficiently small values of τ. In this expression, ω n is the constant used in the proof of Theorem 3, for each n ∈ I N−1 . Moreover, in this case, Letting now C 1 = 8e C T , using the discrete Grönwall inequality [31] for sufficiently small values of τ and substituting the initial conditions of (55), we readily reach that As a consequence, there exists a constant C ≥ 0 with the property that m 2 , ζ m 2 ≤ C(τ 2 + h 2 2 ), for each m ∈ I N . We conclude that the scheme (17) converges to the exact solution of (4), with the quadratic order of convergence.

Conclusions
Before concluding this paper, we must note that our discretization of the double fractional Gross-Pitaevskii system follows a standard approach [4,25]. In particular, this means that the numerical discretization uses local approximations, estimating each of the terms of the mathematical model. A natural question would be whether a non-standard approach could be used to discretize our mathematical model. Here, we understand the adjective 'non-standard' in the sense of R.E. Mickens [32]. Such topic of investigation is outside the scope of this work, but it is an interesting avenue of research for a future study. On the other hand, the constant C ≥ 0 in Lemma 2 are not necessarily independent of h and τ. This is due to the fact that the nonlinear terms in (4) are not globally Lipschitz continuous. This is one of the limitations of our study, and it remains an open problem of research. To solve this shortcoming, it would be desirable to establish the uniform boundedness of the numerical approximations obtained through (4). Following recent reports available in the literature [24], the determination of conserved positive quantities for the finite-difference method (17) would be helpful in bounding the numerical solutions [7]. The authors have attempted to employ the discrete energy method to derive such quantities [33,34]. Unfortunately, their efforts have not yielded conserved positive quantities to this day.