Next Article in Journal
Impulsive Fractional Differential Inclusions and Almost Periodic Waves
Next Article in Special Issue
Collocation Methods for High-Order Well-Balanced Methods for Systems of Balance Laws
Previous Article in Journal
Finite-Time Projective Synchronization of Caputo Type Fractional Complex-Valued Delayed Neural Networks
Previous Article in Special Issue
Integrating Semilinear Wave Problems with Time-Dependent Boundary Values Using Arbitrarily High-Order Splitting Methods
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Convergent Three-Step Numerical Method to Solve a Double-Fractional Two-Component Bose–Einstein Condensate

by
Adán J. Serna-Reyes
1,
Jorge E. Macías-Díaz
2,3,* and
Nuria Reguera
4
1
Centro de Ciencias Básicas, Universidad Autónoma de Aguascalientes, Aguascalientes 20131, Mexico
2
Department of Mathematics and Didactics of Mathematics, School of Digital Technologies, Tallinn University, 10120 Tallinn, Estonia
3
Departamento de Matemáticas y Física, Universidad Autónoma de Aguascalientes, Aguascalientes 20131, Mexico
4
Departamento de Matemáticas y Computación, Universidad de Burgos, IMUVA, 09001 Burgos, Spain
*
Author to whom correspondence should be addressed.
Submission received: 19 May 2021 / Revised: 1 June 2021 / Accepted: 2 June 2021 / Published: 18 June 2021
(This article belongs to the Special Issue Numerical Methods for Evolutionary Problems)

Abstract

:
This manuscript introduces a discrete technique to estimate the solution of a double-fractional two-component Bose–Einstein condensate. The system consists of two coupled nonlinear parabolic partial differential equations whose solutions are two complex functions, and the spatial fractional derivatives are interpreted in the Riesz sense. Initial and homogeneous Dirichlet boundary data are imposed on a multidimensional spatial domain. To approximate the solutions, we employ a finite difference methodology. We rigorously establish the existence of numerical solutions along with the main numerical properties. Concretely, we show that the scheme is consistent in both space and time as well as stable and convergent. Numerical simulations in the one-dimensional scenario are presented in order to show the performance of the scheme. For the sake of convenience, A MATLAB code of the numerical model is provided in the appendix at the end of this work.

1. Introduction

In this work, p N represents the number of spatial dimensions, and T R + defines a period of time. We agree that I n = { 1 , 2 , . . . , n } and I ¯ n = I n { 0 } , for each n N , and assume that a i and b i are real numbers with the property that a i < b i , for all i I p . For convenience, we use the notations Ω = Π i = 1 p ( a i , b i ) R p and Ω T = Ω × ( 0 , T ) and agree that Ω ¯ and Ω ¯ T denote the closures of Ω and Ω T , respectively, in the Euclidean metric topology of R p . We assume that ψ 1 : Ω ¯ T C and ψ 2 : Ω ¯ T C are sufficiently smooth functions and employ the conventions x = ( x 1 , x 2 , , x p ) Ω and i = 1 . Throughout this manuscript, Γ will represent the usual Gamma function that extends factorials.
Definition 1
(Podlubny [1]). Suppose that f is a real function defined in all of R and assume that n N { 0 } and α R are such that n 1 < α < n . We define the fractional derivative in the Riesz sense of order α of f at the point x R as
d α f ( x ) d | x | α = 1 2 cos ( π α 2 ) Γ ( n α ) d n d x n f ( ξ ) d ξ | x ξ | α + 1 n .
Definition 2.
Let ψ : Ω ¯ T C be a function and let i I p . Suppose that α and n are a real and a natural number, respectively, satisfying the properties in the previous definition. Whenever it exists, the space-fractional partial derivative in the sense of Riesz of order α of ψ with respect to the variable x i at x , t Ω T is defined as
α ψ ( x , t ) | x i | α = 1 2 cos ( π α 2 ) Γ ( n α ) n x i n ψ ( x 1 , , x i 1 , ξ , x i + 1 , , x p , t ) | x i ξ | α 1 d ξ .
Here, we extend ψ by allowing ψ 0 outside of Ω ¯ . If these fractional derivatives exist for all i I p , then we define
α ψ ( x , t ) = i = 1 p α ψ | x i | α ( x , t ) .
In this manuscript, we will use D, β 11 , β 12 , β 22 , and λ to represent constant real numbers, and we employ P : Ω ¯ R to denote a real function. Throughout, we will assume that α 1 , α 2 R satisfy the properties 1 < α 1 2 and 1 < α 2 2 . We let ϕ 1 : Ω ¯ C and ϕ 2 : Ω ¯ C be sufficiently regular functions. In this work, we tackle the numerical solution of the following double-fractional initial-boundary value problem ruled by two nonlinear coupled Gross–Pitaevskii-type equations:
i ψ 1 t = P ( x ) + D + β 11 ψ 1 2 + β 12 ψ 2 2 1 2 α 1 ψ 1 + λ ψ 2 , ( x , t ) Ω T , i ψ 2 t = P ( x ) + β 12 ψ 1 2 + β 22 ψ 2 2 1 2 α 2 ψ 2 + λ ψ 1 , ( x , t ) Ω T , satisfying ψ 1 ( x , 0 ) = ϕ 1 ( x ) , x Ω ¯ , ψ 2 ( x , 0 ) = ϕ 2 ( x ) , x Ω ¯ , ψ 1 ( x , t ) = ψ 2 ( x , t ) = 0 , ( x , t ) Ω × ( 0 , T ) .
This system is a two-component Bose–Einstein condensate that considers the presence of an internal atomic Josephson junction. In this physical context, ψ 1 and ψ 2 denote the two-component condensate wave functions.
There are systems that describe single-component Bose–Einstein condensates. As examples, there are works that investigate round, symmetric, and central vortex states in rotating Bose–Einstein condensates [2], considering only a parameter that characterizes the interaction between particles and a constant that corresponds to the angular speed of the laser beam. As in the present work, a real-valued external trapping potential is employed in those models. From the numerical point of view, there are already works that tackle the numerical solution of single-component fractional Nose–Einstein condensates.
As an example, there are papers that computationally studied the ground states of space-fractinoal nonlinear Schrödinger–Gross–Pitaevskii equations with a rotation term and non-local nonlinear interactions [3]. However, we must clarify that the single-component non-fractional case has been more studied that the double-component fully fractional scenario [4,5], possibly due to the mathematical difficulties in the double-component and fractional scenario. Clearly, this case is physically more relevant as it describes the dynamics of mixtures Bose–Einstein condensates that have been experimentally observed in the laboratory [6,7,8,9,10].
The following are crucial definitions and properties to discretize the fractional derivatives in the initial-boundary value problem in Equation (4).
Definition 3
(Ortigueira [11]). Let f be any real function defined on R . Suppose that h R + and that α is a real number satisfying α > 1 . We define the centered difference of fractional order α of f at the number x as
Δ h ( α ) f ( x ) = k = f ( x k h ) g k ( α ) , x R ,
with
g k ( α ) = ( 1 ) k Γ ( α + 1 ) Γ ( α 2 + k + 1 ) Γ ( α 2 k + 1 ) , k Z .
Lemma 1
(Wang et al. [12]). Let 0 < α 2 and α 1 .
(a) 
If k 0 , then g k ( α ) = g k ( α ) < 0 .
(b) 
g 0 ( α ) R + .
(c) 
k = g k ( α ) = 0 .
Lemma 2
(Wang et al. [12]). Let f : R R be a function in the set C 5 ( R ) , which has integrable derivatives to the fifth order on R . If α ( 0 , 1 ) ( 1 , 2 ] , then
Δ h α f ( x ) h α = α f ( x ) | x | α + O ( h 2 ) ,
for almost all x R .
The field of fractional calculus has witnessed a vertiginous development in recent years, partly due to the vast amount of potential applications in the physical sciences [13]. Today, many different fractional derivatives and integrals have been proposed. For example, some of the first fractional derivatives introduced historically in mathematics were the Riemann–Liouville fractional derivatives [14], which generalized the classical integer-order derivatives with respect to specific analytical properties [1].
The fractional derivatives in the senses of Caputo, Riesz, and Grünwald–Letnikov are also extensions of the traditional derivatives of integer order. These fractional operators are nonequivalent in general, and various applications of all of them have been proposed in science and engineering [15,16]. For instance, some reports have provided theoretical foundations for the application of fractional calculus to the theory of viscoelasticity [17], while others have proposed possible applications of fractional calculus to dynamic problems of solid mechanics [18], continuous-time financial economics [19,20], Earth system dynamics [21], mathematical modeling of biological phenomena [22], and the modeling of two-phase gas/liquid flow systems [23], to mention some potential applications.
However, it is important to recall that Riesz-type derivatives may be the only fractional derivatives that have real physical applications. This is due to a well-known result by Tarasov, which establishes that Riesz fractional derivatives result from systems with long-range interactions in a continuum-limit case [24]. In turn, systems consisting of particles with long-range interactions are useful in statistical mechanics, thermostatics [25,26], and the theory of biological oscillator networks [27].
From the mathematical point of view, many interesting avenues of investigation have been opened by the progress in fractional calculus. Indeed, the different fractional derivatives have found discrete analogues, which have been used extensively in the literature. As examples, Riesz fractional derivatives have been discretized consistently in various fashions using fractional-order centered differences [11,28] and weighted-shifted Grünwald differences [29,30].
Clearly, those discrete approaches have been studied to determine their analytical properties, and they have been used extensively to provide discrete models to solve Riesz space-fractional conservative/dissipative space-fractional wave equations [31], a Hamiltonian fractional nonlinear elastic string equation [32], an energy-preserving double fractional Klein–Gordon–Zakharov system [33], and even a Riesz space-fractional generalization with generalized time-dependent diffusion coefficient and potential of the Higgs boson equation in the de Sitter space-time [32], among other complex systems.
On the other hand, Caputo fractional derivatives have been discretized consistently using various criteria. For instance, some high-order L 2 -compact difference approaches have used to that end [34], as well as L 1 formulas [35,36] and L 1 –2 methodologies [37]. Using those approaches, various numerical schemes have been proposed to efficiently solve Caputo time-fractional diffusive and wave differential equations [38,39,40]. Various potential applications of these systems have been reported in the sciences [25,41].
From the analytical point of view, the literature offers a wide range of reports that focus on the extension of integer-order methods and results for the fractional case. For example, there are various articles that tackle the existence, uniqueness, regularity, and asymptotic behavior of the solution for the fractional porous medium equation [42], nonlinear fractional diffusion equations [43], nonlinear fractional heat equations [44], the Fisher–Kolmogorov–Petrovskii–Piscounov equation with nonlinear fractional diffusion [45], fractional thin-film equations [46], and the fractional Schrödinger equation with general nonnegative potentials [47].
From a more particular point of view, the fractional generalization of the classical vector calculus operators (that is, the gradient, divergence, curl, and Laplacian operators) has also been an active topic of research developed from different approaches. Some of the first attempts to extend those operators to the fractional scenario were proposed in [48,49] using the Nishimoto fractional derivative.
These operators were used later on in [50] to provide a physical interpretation for the fractional advection-dispersion equation for flow in heterogeneous porous media (see [51] and the references therein for a historic account of the efforts to formulate a fractional form of vector calculus). More recently, a new generalization of the Helmholtz decomposition theorem for both fractional time and space was proposed in [52,53] using the discrete Grünwald–Letnikov fractional derivative. As in the present manuscript, the authors of [52,53] considered different derivative orders assuming non-homogeneous models and non-isotropic spaces.
This manuscript is organized as follows. Section 2 is devoted to present the discrete nomenclature and the numerical model to solve the continuous problem under investigation in this work. We introduce suitable discrete norms and recall a useful theorem by Desplanques. With that result at hand, we rigorously establish the properties of existence and uniqueness of numerical solutions. Section 3 is devoted to establishing the most important computational features of our numerical scheme. More concretely, we show that the discrete model proposed in this work has second order consistency under suitable conditions on the smoothness of the solutions.
We also rigorously prove the stability property of the scheme along with its quadratic order of convergence in both time and space. In turn, Section 4 is devoted to show some illustrative simulations that exhibit the stability and the convergence of the method proposed in this work. We will also provide numerical evidence that the finite-difference scheme proposed here has a quadratic order of convergence in both space and time. Finally, we close this manuscript with some concluding remarks. For the sake of convenience to the reader, an appendix is provided at the end of this paper, in which we present the numerical scheme used to produce the simulations.

2. Numerical Method

Throughout this manuscript, we let M i , N N , for each i I p , and fix the spatial and temporal step-sizes h i = b i a i / M i and τ = T / N , for each i I p . We will consider regular partitions of the intervals [ a i , b i ] and [ 0 , T ] represented as
x i , j i = a i + j i h i , i I p , j i I ¯ M i ,
t n = n τ , n I ¯ N ,
respectively. Let J = i = 1 p I M i 1 and J ¯ = i = 1 p I ¯ M i , and convey that x j = ( x 1 , j 1 , , x p , j p ) , for each j = ( j 1 , , j p ) J ¯ . Throughout this work, we will employ the symbol ( u j n , v j n ) to denote a numerical approximation to ( U j n , V j n ) , where U j n = ψ 1 ( x j , t n ) and V j n = ψ 2 ( x j , t n ) ) , for each ( j , n ) J ¯ × I ¯ N . Let us agree that P j = P ( x j ) , for each j J ¯ . We will employ the symbol J to represent the boundary of J ¯ , which is defined as the set of all multi-indexes j J that satisfy x j Ω .
Definition 4.
We introduce the following average operators, for each α ( 0 , 1 ) ( 1 , 2 ] , ( j , n ) J × I N 1 and w = u , v :
μ t w j n = w j n + 1 + w j n 2 ,
μ t ( 1 ) w j n = 3 w j n w j n 1 2 .
Under the same conditions, we introduce the difference operators
δ t w j n = w j n + 1 w j n τ ,
δ x i ( α ) w j n = 1 h i α k = 0 M i g j i k ( α ) w ( x 1 , j 1 , , x i 1 , j i 1 , x i , k , x i + 1 , j i + 1 , , x p , j p , t n ) .
In the present manuscript, we convey that
h ( α ) w j n = δ x 1 ( α ) w j n + δ x 2 ( α ) w j n + + δ x p ( α ) w j n .
Moreover, for the sake of convenience, we also introduce the discrete fractional gradient vector
h ( α ) = ( δ x 1 ( α ) w j n , δ x 2 ( α ) w j n , , δ x p ( α ) w j n ) .
With this nomenclature at hand, the following numerical method will be employed in this work to approximate the solution of the initial-boundary value problem Equation (4) on the set Ω ¯ × [ 0 , T ] , for each ( j , n ) J × I ¯ N 1 :   
i δ t u j n = 1 2 h ( α 1 ) + P j + D μ t u j n + β 11 μ t ( 1 ) u j n 2 + β 12 μ t ( 1 ) v j n 2 μ t ( 1 ) u j n + λ μ t ( 1 ) v j n , i δ t v j n = 1 2 h ( α 2 ) + P j μ t v j n + β 12 μ t ( 1 ) u j n 2 + β 22 μ t ( 1 ) v j n 2 μ t ( 1 ) v j n + λ μ t ( 1 ) u j n , such that u j 0 = μ t u j 0 = μ t ( 1 ) u j 0 = ϕ 1 ( x j ) , j J , v j 0 = μ t v j 0 = μ t ( 1 ) v j 0 = ϕ 2 ( x j ) , j J , u j n = v j n = 0 , ( j , n ) J × I ¯ N .
One may readily check that the general iterative formula of this numerical model makes it a three-step implicit model. On the other hand, to determine the approximations at the time t 1 , we need to employ the initial conditions u j 0 = μ t u j 0 = μ t ( 1 ) u j 0 = ϕ 1 ( x j ) and v j 0 = μ t v j 0 = μ t ( 1 ) v j 0 = ϕ 2 ( x j ) , for all j J . Substituting these formulas into the general iterative equations at time t 0 , we obtain
u j 1 = ϕ 1 ( x j ) + i τ 1 2 h ( α 1 ) P j D ϕ 1 ( x j ) i τ λ ϕ 2 ( x j ) i τ β 11 ϕ 1 ( x j ) 2 + β 12 ϕ 2 ( x j ) 2 ϕ 1 ( x j ) , j J ,
and
v j 1 = ϕ 2 ( x j ) + i τ 1 2 h ( α 2 ) P j ϕ 2 ( x j ) i τ λ ϕ 1 ( x j ) i τ β 12 ϕ 1 ( x j ) 2 + β 22 ϕ 2 ( x j ) 2 ϕ 2 ( x j ) , j J .
For the remainder of his manuscript, we define h = ( h 1 , , h p ) and h * = i = 1 p h i . We will employ V h to represent the collection of all complex functions on the spatial mesh { x j : j J ¯ } R p . Finally, for each w V h and all j J ¯ , let w j = w ( x j ) .
Definition 5.
Let us introduce the functions · , · : V h × V h C and · 1 : V h R , defined respectively by
u , v = h * j I u j v j ¯ ,
u 1 = h * j I | u j | ,
for any u , v V h . We define the function · 2 : V h R as the Euclidean norm obtained from · , · . If u , v V h , then we agree that u v V h represents the pointwise multiplication of u and v. Moreover, | u | V h will denote the absolute value of the complex function u.
Lemma 3
(Desplanques [54]). Let A be a square complex matrix. If A is strictly diagonally dominant, then A is invertible.
Theorem 1
(Unique solvability). If ( ϕ 1 , ϕ 2 ) is any set of initial conditions, then the finite-difference model Equation (16) has a unique solution.
Proof. 
The proof will be provided for the one-dimensional scenario only in view that the case for higher dimensions is similar. Note that, beforehand, the approximations u 0 , v 0 , u 1 and v 1 are provided explicitly by the initial conditions and the formulas Equations (17) and (18). Proceeding by induction, suppose that u n , v n , u n 1 , and v n 1 were calculated, for some n I N 1 . It is important to observe now that the first equation of the numerical method can be equivalently rewritten as
i τ 1 4 h α 1 g 0 α 1 1 2 P j + D u j n + 1 1 4 h α 1 l = 0 l j M 1 g l j ( α 1 ) u l n + 1 = b ( u j n ) + a ( u j n 1 , u j n , v j n 1 , v j n ) , ( j , n ) J × I N 1 ,
where
a ( u j n 1 , u j n , v j n 1 , v j n ) = ( β 11 | μ t ( 1 ) u j n | 2 + β 12 | μ t ( 1 ) v j n | 2 ) μ t ( 1 ) u j n + λ μ t ( 1 ) v j n ,
b ( u j n ) = i τ + 1 2 1 2 h ( α 1 ) + P j + D u j n ,
for each ( j , n ) J × I N 1 . It follows that the first identity of the numerical scheme Equation (16) can be alternatively expressed in vector form as
A u n + 1 = b ( u n ) a ( u n 1 , u n , v n 1 , v n ) , n I N 1 .
In this expression, the square complex matrix A = ( a i j ) is defined by -4.6cm0cm
A = g 0 ( α 1 ) 4 h 1 α 1 + P 0 + D 2 i τ g 1 ( α 1 ) 4 h 1 α 1 g 2 M 1 ( α 1 ) 4 h 1 α 1 g 1 ( α 1 ) 4 h 1 α 1 g 0 ( α 1 ) 4 h 1 α 1 + P 1 + D 2 i τ g 3 M 1 ( α 1 ) 4 h 1 α 1 g M 1 2 ( α 1 ) 4 h 1 α 1 g M 1 3 ( α 1 ) 4 h 1 α 1 g 0 ( α 1 ) 4 h 1 α 1 + P M 1 + D 2 i τ .
Notice now that the diagonal entries of A satisfy
a i i = g 0 ( α 1 ) 4 h 1 α 1 + 1 2 P i 1 + D 2 + 1 τ 2 .
On the other hand, the numbers in the real sequence ( g k ( α 1 ) ) k = guarantee that
j i M 1 a i j = j i M 1 g l j ( α 1 ) 4 h 1 α 1 = j i M 1 g l j α 1 4 h α 1 < l = l j g l j α 1 4 h α 1 = g 0 α 1 4 h α 1 a i i .
It follows that the matrix A is non-singular by the previous lemma, which implies then that there exists a unique solution u n + 1 of the vector equation Equation (24). In similar fashion, we can readily prove the existence of the vector v n + 1 after noticing that the second difference equation of the discrete system Equation (16) may be rewritten in vector form as
B v n + 1 = d ( v n ) c ( u n 1 , u n , v n 1 , v n ) ,
where, for each ( j , n ) J × I N 1 ,
c ( u n 1 , u n , v n 1 , v n ) = ( β 12 | μ t ( 1 ) u j n | 2 + β 22 | μ t ( 1 ) v j n | 2 ) μ t ( 1 ) v j n + λ μ t ( 1 ) u j n ,
d ( v j n ) = i τ + 1 2 1 2 h ( α 2 ) + P j v j n .
In this expression, B is the square complex matrix given by
B = g 0 ( α 1 ) 4 h 1 α 2 + P 0 2 i τ g 1 ( α 2 ) 4 h 1 α 2 g 2 M 1 ( α 2 ) 4 h 1 α 2 g 1 ( α 2 ) 4 h 1 α 2 g 0 ( α 2 ) 4 h 1 α 2 + P 1 2 i τ g 3 M 1 ( α 2 ) 4 h 1 α 2 g M 1 2 ( α 2 ) 4 h 1 α 2 g M 1 3 ( α 2 ) 4 h 1 α 2 g 0 ( α 2 ) 4 h 1 α 2 + P M 1 2 i τ .
One may readily check that this matrix is strictly diagonally dominant; therefore, Equation (28) has a unique solution v n + 1 . The conclusion of this theorem readily follows now by mathematical induction.    □

3. Numerical Properties

The aim of the present section is to theoretically demonstrate the most relevant features of the discrete model Equation (16). Concretely, we, herein, prove properties, including the consistency, stability, and convergence of our numerical model. In a first stage, we will require to define the following continuous operators for each ( x , t ) Ω T :
L 1 ( ψ 1 ( x , t ) , ψ 2 ( x , t ) ) = i ψ 1 ( x , t ) t + 1 2 α 1 P ( x ) D β 11 ψ 1 ( x , t ) 2 β 12 ψ 2 ( x , t ) 2 ψ 1 ( x , t ) λ ψ 2 ( x , t ) ,
and
L 2 ( ψ 1 ( x , t ) , ψ 2 ( x , t ) ) = i ψ 2 ( x , t ) t + 1 2 α 2 P ( x ) β 12 ψ 1 ( x , t ) 2 β 22 ψ 2 ( x , t ) 2 ψ 2 ( x , t ) λ ψ 1 ( x , t ) .
For each ( j , n ) J × I N 1 , we also introduce the discrete operators
L 1 ( U j n , V j n ) = i δ t U j n + 1 2 h ( α 1 ) P j D μ t U j n β 11 μ ( 1 ) U j n 2 + β 12 μ ( 1 ) V j n 2 μ ( 1 ) U j n λ μ ( 1 ) V j n ,
and
L 2 ( U j n , V j n ) = i δ t V j n + 1 2 h ( α 2 ) P j μ t V j n β 12 μ ( 1 ) U j n 2 + β 22 μ ( 1 ) V j n 2 μ ( 1 ) V j n λ μ ( 1 ) U j n .
Here, recall that we previously agreed that U j n = ψ 1 ( x j , t n ) and V j n = ψ 2 ( x j , t n ) , for each ( j , n ) J ¯ × I ¯ N . For the sake of convenience, when ( x , t ) Ω T and ( j , n ) J × I N 1 , then we agree that
L ( ψ 1 ( x , t ) , ψ 2 ( x , t ) ) = L 1 ( ψ 1 ( x , t ) , ψ 2 ( x , t ) ) , L 2 ( ψ 1 ( x , t ) , ψ 2 ( x , t ) ) ,
L ( ψ 1 ( x j , t n ) , ψ 2 ( x j , t n ) ) = ( L 1 ( ψ 1 ( x j , t n ) , ψ 2 ( x j , t n ) ) , L 2 ( ψ 1 ( x j , t n ) , ψ 2 ( x j , t n ) ) ) .
Theorem 2
(Consistency). Suppose that ψ 1 , ψ 2 C x , t 5 , 4 ( Ω T ¯ ) and P C ( Ω ¯ ) . There are positive numbers C 1 and C 2 independent of h and τ, which satisfy the inequalities
L ( ψ 1 ( x j , t n ) , ψ 2 ( x j , t n ) ) L ( ψ 1 ( x j , t n + 1 2 ) , ψ 2 ( x j , t n + 1 2 ) ) C 1 ( τ 2 + h 2 2 ) ,
L ( ψ 1 ( x j , t n ) , ψ 2 ( x j , t n ) ) L ( ψ 1 ( x j , t n ) , ψ 2 ( x j , t n ) ) C 2 ( τ + h 2 2 ) ,
for each ( j , n ) J × I N 1 .
Proof. 
To prove these results, we require Taylor’s theorem, Lemma 2, and the smoothness assumptions on the functions ψ 1 and ψ 2 . There exist nonnegative constants C 1 , k , C 2 , k , C 3 , k , i , C 4 , k , and C 5 , k . l , which are independent of τ and h, for each i I p and k , l I 2 , with the property that
μ t ( 1 ) ψ k ( x j , t n ) ψ k ( x j , t n + 1 2 ) C 1 , k τ 2 ,
δ t ψ k ( x j , t n ) ψ k t ( x j , t n + 1 2 ) C 2 , k τ 2 ,
μ t δ x i ( α 1 ) ψ k ( x j , t n ) α 1 ψ k | x | α 1 ( x j , t n + 1 2 ) C 3 , k , i ( τ 2 + h i 2 ) ,
P ( x j ) μ t ψ k ( x j , t n ) P ( x j ) ψ k ( x j , t n + 1 2 ) C 4 , k τ 2 ,
μ t ( 1 ) ψ k ( x j , t n ) 2 μ t ( 1 ) ψ l ( x j , t n ) ψ k ( x j , t n + 1 2 ) 2 ψ l ( x j , t n + 1 2 ) C 5 , k , l τ 2 ,
for each ( j , n ) J × I N 1 . The inequality in the conclusion of this proposition follows now applying the triangle inequality for an appropriate number C 1 , which depends only on C 1 , k , C 2 , k , C 3 , k , i , C 4 , k , and C 5 , k . l for each i I p and k , l I 2 , and the model parameters. This implies, in particular, that C 2 is independent of h and τ , as desired. The constant C 2 is obtained applying the same arguments around the node ( x j , t n ) instead of ( x j , t n + 1 2 ) .    □
We tackle now the stability and the convergence of the numerical scheme Equation (16). We require the next technical results.
Lemma 4
(Pen-Yu [55]). Suppose that ( ω n ) n = 0 N and ( ρ n ) n = 0 N are sequences of non-negative numbers. Let C 0 satisfy the inequality
ω n ρ n + C τ k = 0 n 1 ω k .
Then, ω n ρ n e C n τ for each 0 n N .
Lemma 5
(Macías-Díaz [56]). Suppose that α is a real number with α ( 0 , 1 ) ( 1 , 2 ] . If  u , v V h , then h ( α ) u , v = h ( α / 2 ) u , h ( α / 2 ) v .
Lemma 6.
If α ( 1 , 2 ] and w = ( w n ) n = 1 N V h , then the following holds, for each n I N 1 :
(a) 
2 δ t w n , μ t w n = δ t w n 2 2 + i τ 2 Im w n + 1 , w n ,
(b) 
2 h ( α ) μ t w n , μ t w n = 2 h ( α / 2 ) μ t w n 2 2 ,
(c) 
2 V μ t w n , μ t w n = 2 V , | μ t w n | 2 .
Lemma 7.
Let ϵ = ( ϵ n ) n = 0 N and ζ = ( ζ n ) n = 0 N be sequences in V h . If m I N 1 , then
2 n = 1 m Im μ t ( 1 ) ϵ n , μ t ζ n + Im μ t ( 1 ) ζ n , μ t ϵ n 1 4 ϵ 0 2 2 + ζ 0 2 2 + 3 2 n = 1 m ϵ n 2 2 + ζ n 2 2 + ϵ m + 1 2 2 + ζ m + 1 2 2 .
Proof. 
The proof follows after applications of Young’s inequality and some algebraic regrouping.    □
The stability of the discrete method Equation (16) will require us to consider two sets of initial conditions, namely ( ϕ 1 , ϕ 2 ) and ( ϕ ˜ 1 , ϕ ˜ 2 ) , where all the initial data are complex functions, which are defined on Ω ¯ . The respective solutions will be denoted by ( u , v ) and ( u ˜ , v ˜ ) . This means that ( u , v ) satisfies Equation (16), while ( u ˜ , v ˜ ) satisfies the problem
i δ t u ˜ j n = 1 2 h ( α 1 ) + P j + D μ t u ˜ j n + β 11 μ t ( 1 ) u ˜ j n 2 + β 12 μ t ( 1 ) v ˜ j n 2 μ t ( 1 ) u ˜ j n + λ μ t ( 1 ) v ˜ j n , i δ t v ˜ j n = 1 2 h ( α 2 ) + P j μ t v ˜ j n + β 12 μ t ( 1 ) u ˜ j n 2 + β 22 μ t ( 1 ) v ˜ j n 2 μ t ( 1 ) v ˜ j n + λ μ t ( 1 ) u ˜ j n , such that u ˜ j 0 = μ t ( 1 ) u ˜ j 0 = ϕ ˜ 1 ( x j ) , j J , v ˜ j 0 = μ t ( 1 ) v ˜ j 0 = ϕ ˜ 2 ( x j ) , j J , u ˜ j n = v ˜ j n = 0 , ( j , n ) J × I ¯ N .
for each ( j , n ) J × I N 1 .
Lemma 8.
Let ( u , v ) and ( u ˜ , v ˜ ) be solutions of Equation (16), respectively, and let ϵ n = u n u ˜ n and ζ n = v n v ˜ n , for each n I ¯ N . For each ( j , n ) J × I N 1 , let us define
P j n = β 11 μ t ( 1 ) u j n 2 + β 12 μ t ( 1 ) v j n 2 μ t ( 1 ) u j n β 11 μ t ( 1 ) u ˜ j n 2 + β 12 μ t ( 1 ) v ˜ j n 2 μ t ( 1 ) u ˜ j n ,
Q j n = β 22 μ t ( 1 ) v j n 2 + β 12 μ t ( 1 ) u j n 2 μ t ( 1 ) v j n β 22 μ t ( 1 ) v ˜ j n 2 + β 12 μ t ( 1 ) u ˜ j n 2 μ t ( 1 ) v ˜ j n .
Then, there exists a constant C 0 , which is independent of h and τ, with the property that
max { | P j n | , | Q j n | } C ( | ϵ j n 1 | + | ζ j n 1 | + | ϵ j n | + | ζ j n | ) , ( j , n ) J × I N 1 .
Proof. 
The result is straightforward.    □
Theorem 3
(Nonlinear stability). Consider the solutions ( u , v ) and ( u ˜ , v ˜ ) of Equation (16) corresponding to initial conditions ( ϕ 1 , ϕ 2 ) and ( ϕ ˜ 1 , ϕ ˜ 2 ) , respectively. For each n I ¯ N , define ϵ n = u n u ˜ n and ζ n = v n v ˜ n . If ( C + λ ) τ < 1 2 , then there is a nonnegative constant C that is not dependent on τ and h, such that
ϵ n 2 2 + ζ n 2 2 ( 1 + ( C + λ 4 ) τ ) ( ϵ 0 2 2 + ζ 0 2 2 ) e 2 C T , n I ¯ N .
Proof. 
Beforehand, take C as in the last result. Note that the existence of ( u , v ) and ( u ˜ , v ˜ ) is guaranteed by Theorem 1. It is clear that ( ϵ , η ) satisfies the following discrete problem, for each ( j , n ) J × I N 1 :
i δ t ϵ j n = 1 2 h ( α 1 ) + P j + D μ t ϵ j n + P j n + λ μ t ( 1 ) ζ j n , i δ t ζ j n = 1 2 h ( α 2 ) + P j μ t ζ j n + Q j n + λ μ t ( 1 ) ϵ j n , such that ϵ j 0 = μ t ( 1 ) ϵ j 0 = ϕ 1 ( x j ) ϕ ˜ 1 ( x j ) , j J ¯ , ζ j 0 = μ t ( 1 ) ζ j 0 = ϕ 2 ( x j ) ϕ ˜ 2 ( x j ) , j J ¯ , ϵ j n = ζ j n = 0 , ( j , n ) J × I ¯ N .
Here, P j n and Q j n are as in Lemma 8. Then, there is a nonnegative real number C such that, for each n I N 1 ,
max { P n 2 2 , Q n 2 2 } C ϵ n 1 2 2 + ζ n 1 2 2 + ϵ n 2 2 + ζ n 2 2 .
The number C does not depend on h and τ . Take the first identity in Equation (52) and calculate its inner product with 2 μ t ( 1 ) ϵ n . On both sides of the result, next calculate the imaginary part. Similarly, consider the second identity of Equation (52) and take the inner product with 2 μ t ( 1 ) ζ n . Calculate the imaginary part on both sides. After applying the identities after Theorem 6, we readily obtain that
δ t ϵ n 2 2 = 2 Im P n , μ t ϵ n + 2 λ Im μ t ( 1 ) ζ n , μ t ϵ n , n I N 1 ,
δ t ζ n 2 2 = 2 Im Q n , μ t ζ n + 2 λ Im μ t ( 1 ) ϵ n , μ t ζ n , n I N 1 .
After adding those two results side by side and bounding from above, we notice that there exists a nonnegative number C that is not dependent on τ and h, with
ϵ m + 1 2 2 + ζ m + 1 2 2 = ϵ 0 2 2 + ζ 0 2 2 + 2 τ n = 1 m Im ( P n , μ t ϵ n + Q n , μ t ζ n ) + λ τ 4 ϵ 0 2 2 + ζ 0 2 2 + 4 ϵ m + 1 2 2 + ζ m + 1 2 2 + 6 n = 1 m ϵ n 2 2 + ζ n 2 2 ( 1 + λ τ 4 ) ( ϵ 0 2 2 + ζ 0 2 2 ) + τ 2 n = 1 m 2 P n 2 2 + 2 Q n 2 2 + ϵ n 2 2 + ζ n 2 2 + ϵ n + 1 2 2 + ζ n + 1 2 2 + λ τ 2 2 ϵ m + 1 2 2 + ζ m + 1 2 2 + 3 n = 1 m ϵ n 2 2 + ζ n 2 2 ( 1 + ( C + λ 4 ) τ ) ( ϵ 0 2 2 + ζ 0 2 2 ) + 3 τ C + λ 2 n = 1 m ϵ n 2 2 + ζ n 2 2 + τ ( C + λ ) ( ϵ m + 1 2 2 + ζ m + 1 2 2 ) , m I M 1 .
Subtract the term τ ( C + λ ) ( ϵ m + 1 2 2 + ζ m + 1 2 2 ) from both ends of this chain of inequalities, group the common terms, and we notice that ( C + λ ) τ < 1 2 . As a consequence,
ϵ m + 1 2 2 + ζ m + 1 2 2 ( 1 + ( C + λ 4 ) τ ) ( ϵ 0 2 2 + ζ 0 2 2 ) + 3 τ C + λ 2 n = 1 m ϵ n 2 2 + ζ n 2 2 ,
for each m I M 1 . The conclusion is reached then using Lemma 4 with the constants ω n = ϵ n 2 2 + ζ n 2 2 and ρ n = ( 1 + ( C + λ 4 ) τ ) ϵ 0 2 2 + ζ 0 2 2 , for each n I ¯ N 1 .    □
Finally, we investigate now the convergence of the finite-difference method Equation (16). As mentioned previously, ( u , v ) represents a solution of the numerical model, while ( U , V ) = ( ψ 1 , ψ 2 ) denotes a solution of the continuous problem. It is clear that ( U , V ) satisfies the discrete initial-boundary value problem
i δ t U j n = 1 2 h ( α 1 ) + P j + D μ t U j n + β 11 μ t ( 1 ) U j n 2 + β 12 μ t ( 1 ) V j n 2 μ t ( 1 ) U j n + λ μ t ( 1 ) V j n + ρ j n , i δ t V j n = 1 2 h ( α 2 ) + P j μ t V j n + β 12 μ t ( 1 ) U j n 2 + β 22 μ t ( 1 ) V j n 2 μ t ( 1 ) V j n + λ μ t ( 1 ) U j n + σ j n , such that U j 0 = μ t ( 1 ) U j 0 = ϕ 1 ( x j ) , j J , V j 0 = μ t ( 1 ) V j 0 = ϕ 2 ( x j ) , j J , U j n = V j n = 0 , ( j , n ) J × I ¯ N .
Here, ρ j n and σ j n denote local truncation errors. Under the regularity assumptions on the solutions U and V stated in Theorem 2, there is a nonnegative number C 0 , such that ρ n 2 , σ n 2 C 0 ( τ 2 + h 2 2 ) . This number is independent of τ and h. This fact will be employed in the proof of the following theorem.
Theorem 4
(Convergence). Let ( U , V ) be a solution of Equation (4), such that U , V C x , t ( Ω T ¯ ) and P C ( Ω ¯ ) . If τ is sufficiently small, then the corresponding solution of the finite-difference scheme Equation (16) converges to ( U , V ) with order O ( τ 2 + h 2 2 ) in the L 2 -norm.
Proof. 
The proof is analogous to the previous result. Beforehand, let ϵ n = u n U n and ζ n = v n V n , for each n I ¯ N . It is clear then that the pair ( ϵ , ζ ) satisfies
i δ t ϵ j n = 1 2 h ( α 1 ) + P j + D μ t ϵ j n + P j n + λ μ t ( 1 ) ζ j n + ρ j n , i δ t ζ j n = 1 2 h ( α 2 ) + P j μ t ζ j n + Q j n + λ μ t ( 1 ) ϵ j n + σ j n , such that ϵ j 0 = μ t ( 1 ) ϵ j 0 = 0 , j J ¯ , ζ j 0 = μ t ( 1 ) ζ j 0 = 0 , j J ¯ , ϵ j n = ζ j n = 0 , ( j , n ) J × I ¯ N ,
where, for each ( j , n ) J × I N 1 ,
P j n = β 11 μ t ( 1 ) u j n 2 + β 12 μ t ( 1 ) v j n 2 μ t ( 1 ) u j n β 11 μ t ( 1 ) U j n 2 + β 12 μ t ( 1 ) V j n 2 μ t ( 1 ) U j n
and
Q j n = β 22 μ t ( 1 ) v j n 2 + β 12 μ t ( 1 ) u j n 2 μ t ( 1 ) v j n β 22 μ t ( 1 ) V j n 2 + β 12 μ t ( 1 ) U j n 2 μ t ( 1 ) V j n .
Proceeding now as in the proof of Theorem 3, we obtain a nonnegative number C , which does not depend on τ and h, with the property that Equation (53) is satisfied. Applying the same steps used to obtain Equations (54) and (), we may readily obtain the following identities, valid for each n I N 1 :
δ t ϵ n 2 2 = 2 Im P n , μ t ϵ n + 2 λ Im μ t ( 1 ) ζ n , μ t ϵ n + 2 Im ρ n , μ ϵ n ,
δ t ζ n 2 2 = 2 Im Q n , μ t ζ n + 2 λ Im μ t ( 1 ) ϵ n , μ t ζ n + 2 Im σ n , μ t ζ n .
From here, it is possible to check that there are numbers C 3 , C 4 0 independent of τ and h, such that
ϵ m + 1 2 2 + ζ m + 1 2 2 2 τ n = 0 m Im P n , μ t ϵ n + Im Q n , μ t ζ n + Im ρ n , μ t ϵ n + Im σ n , μ t ζ n ( 1 + ( C + λ 4 ) τ ) ( ϵ 0 2 2 + ζ 0 2 2 ) + 3 τ C + λ 2 n = 1 m ϵ n 2 2 + ζ n 2 2 + τ ( C + λ ) ( ϵ m + 1 2 2 + ζ m + 1 2 2 ) + τ 2 n = 1 m 2 P n 2 2 + 2 Q n 2 2 + ϵ n 2 2 + ζ n 2 2 + ϵ n + 1 2 2 + ζ n + 1 2 2 ( 1 + ( C + λ 4 ) τ ) ( ϵ 0 2 2 + ζ 0 2 2 ) + C 3 ( τ 2 + h 2 2 ) 2 + C 4 τ n = 0 m ϵ n 2 2 + ζ n 2 2 + τ ( C + λ + 1 2 ) ϵ m + 1 2 2 + ζ m + 1 2 2 .
Subtract the term τ ( C + λ + 1 2 ) ϵ m + 1 2 2 + ζ m + 1 2 2 from both ends of this series of inequalities. Notice now that Equation (45) is satisfied, for each m I N 1 , letting C 0 = C 4 , and
ω m = ϵ m 2 2 + ζ m 2 2 , m I N 1 ,
ρ m = ( 1 + ( C + λ 4 ) τ ) ( ϵ 0 2 2 + ζ 0 2 2 ) + C 3 ( τ 2 + h 2 2 ) 2 , m I ¯ N .
Lemma 4 guarantees now that
ϵ n 2 2 + ζ n 2 2 C 5 ( τ 2 + h 2 2 ) 2 , ( j , n ) J × I ¯ N ,
where C 5 = C 3 e C 0 T . From this, it follows that ϵ n 2 , ζ n 2 C 5 ( τ 2 + h 2 2 ) , for each ( j , n ) J × I ¯ N . The conclusion of this proposition readily follows now.    □

4. Illustrative Simulations

We provide illustrative simulations to show the main properties of Equation (16). In particular, we provide a numerical study of the convergence of the scheme, and we confirm the validity of the conclusion of Theorem 4. The experiments provided in this section make use of the Algorithm provided in Appendix A, which is a MATLAB implementation of Equation (16).
In this section, we consider only the spatially one-dimensional scenario. To that end, let p = 1 . We employ the set of parameters in Table 1. In particular, notice that Ω = ( 7 , 7 ) . The parameters α 1 and α 2 may change values from one simulation to the other, and the function V is defined by V ( x ) = 1 2 x 2 , for each x Ω . Finally, the initial data is defined by ϕ 1 ( x ) = ϕ 2 ( x ) = 1 π exp ( x 2 ) , for each x Ω .
Example 1.
Throughout this example, we let T = 10 and fix the computer parameters τ = 0.01 and h 1 = 0.1 . The results of our simulations are provided in Figure 1. In that figure, we provide the graphs of the approximate solutions of the problem Equation (4) as functions of ( x , t ) Ω ¯ × [ 0 , T ] using the scheme Equation (16). More precisely, graphs show the numerical solution of the real parts of (a) ψ 1 and (b) ψ 2 , the imaginary parts of (c) ψ 1 and ψ 2 , and the absolute values of (e) ψ 1 and (f) ψ 2 . The graphs show that the scheme is a stable technique. We repeated this experiment using all the parameters as before, and α 1 = α 2 = 1.5 . The results are provided in Figure 2, and they suggest that the scheme Equation (16) is stable, in agreement with Theorem 3. For illustration purposes, Figure 3 shows the simulations corresponding to α 1 = α 2 = 1 , and they also support the conclusions of this example.
Example 2.
We numerically studied the convergence of the scheme Equation (16), using the parameters in Table 1 along with the functions V, ϕ 1 and ϕ 2 defined in the paragraph before Example 1. In the following experiments, we fix T = 0.5 and α 1 = α 2 = 1.5 , though similar results have been obtained for other values of these parameters (we did not include them to avoid redundancy). The exact solution of the problem is calculated using τ = 5 × 10 5 and h = h 1 = 0.02 . Let ( u , v ) be the solution of the scheme Equation (16) at the time t N , obtained using the computer parameters τ and h. Under these circumstances, let us define the complex constants ϵ τ , h ( x j ) = ψ 1 ( x j ) u j N , for each j J . In turn, we define the temporal and spatial convergence rates, respectively, as
ρ τ , h = log 2 ϵ 2 τ , h 2 ϵ τ , h 2 ,
σ τ , h = log 2 ϵ τ , 2 h 2 ϵ τ , h 2 .
With these conventions, Table 2a,b provides numerical studies of convergence in time and in space, respectively, for the problem described in this example. Various values of τ and h were chosen to that end. However, in any case, the results seem to support that the scheme Equation (16) is quadratically convergent, in agreement with Theorem 4.

5. Conclusions

We proposed a three-level discrete model to solve a double-fractional coupled Bose–Einstein system considering Riesz spatial fractional operators along with two different orders of differentiation. The system is a fractional generalization of the two-component Gross–Pitaevskii system investigated in physics. The solutions of the continuous system are complex; therefore, the computational complexity for solving the system is a shortcoming to be considered.
To alleviate this problem, a three-level scheme was designed to solve our continuous model. This scheme is linear and explicit, which implies that the computational implementation is relatively easy to propose. Indeed, a numerical algorithm to solve this system is presented in the appendix of this work. The code is provided in MATLAB for the sake of convenience, and the nomenclature used to describe it follows that of the theoretical description.
In this manuscript, we proved the existence and uniqueness of the solutions using arguments of linear algebra. We also established that the algorithm is spatially and temporally quadratically consistent. We proved that the finite-difference scheme is spatially and temporally quadratically convergent to the solution of the differential problem. Computational simulations are presented in this work to illustrate this fact.

Author Contributions

Conceptualization, J.E.M.-D.; methodology, J.E.M.-D.; software, A.J.S.-R., J.E.M.-D. and N.R.; validation, A.J.S.-R., J.E.M.-D. and N.R.; formal analysis, A.J.S.-R. and J.E.M.-D.; investigation, A.J.S.-R., J.E.M.-D. and N.R.; resources, A.J.S.-R., J.E.M.-D. and N.R.; data curation, A.J.S.-R. and J.E.M.-D.; writing—original draft preparation, A.J.S.-R., J.E.M.-D. and N.R.; writing—review and editing, J.E.M.-D. and N.R.; visualization, J.E.M.-D.; supervision, J.E.M.-D.; project administration, J.E.M.-D.; funding acquisition, J.E.M.-D. and N.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Council for Science and Technology of Mexico (CONACYT) through grant A1-S-45928 and by Ministerio de Ciencia e Innovación and Regional Development European Funds through project PGC2018-101443-B-I00.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Acknowledgments

The authors would like to thank the anonymous reviewers and the Guest Editor in charge of handling this work for their criticisms. All of their suggestions and comments were followed in the revision of this work. As a consequence, this final version of our work resulted in a substantially improved manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. MATLAB Algorithm

In this appendix, we present a MATLAB implementation of the scheme in Equation (16). To that end, we will provide an alternative form of the linear equations. Equations (24) and (28) in order to avoid computational instabilities. Departing from Equation (21), we multiply both sides of that equation by 2 τ . Proceeding as in the proof of Theorem 1, it follows that expression may be expressed in vector form as
A u n + 1 = 2 τ a ( u n , v n ) + b ( u n 1 ) , n I N 1 .
Here,
a ( u j n , v j n ) = ( β 11 | u j n | 2 + β 12 | v j n | 2 ) u j n + λ v j n , ( j , n ) J × I N 1 ,
b ( u j n ) = i τ 2 δ x 1 ( α 1 ) + τ ( V j + D ) u j n , ( j , n ) J × I N 1 ,
and
A = i τ g 0 ( α 1 ) 2 h 1 α 1 τ ( V 0 + D ) τ g 1 ( α 1 ) 2 h 1 α 1 τ g 2 M 1 ( α 1 ) 2 h 1 α 1 τ g 1 ( α 1 ) 2 h 1 α 1 i τ g 0 ( α 1 ) 2 h 1 α 1 τ ( V 1 + D ) τ g 3 M 1 ( α 1 ) 2 h 1 α 1 τ g M 1 2 ( α 1 ) 2 h 1 α 1 τ g M 1 3 ( α 1 ) 2 h 1 α 1 i τ g 0 ( α 1 ) 2 h 1 α 1 τ ( V M 1 + D ) .
In similar fashion, it is possible to rewrite Equation (28) as
B v n + 1 = 2 τ c ( u n , v n ) + d ( u n 1 ) , n I N 1 .
In this case, we define
c ( u j n , v j n ) = ( β 12 | u j n | 2 + β 22 | v j n | 2 ) v j n + λ u j n , ( j , n ) J × I N 1 ,
d ( v j n ) = i τ 2 δ x 1 ( α 2 ) + τ V j v j n , ( j , n ) J × I N 1
and
B = i τ g 0 ( α 2 ) 2 h 1 α 2 τ V 0 τ g 1 ( α 2 ) 2 h 1 α 2 τ g 2 M 1 ( α 2 ) 2 h 1 α 2 τ g 1 ( α 2 ) 2 h 1 α 2 i τ g 0 ( α 2 ) 2 h 1 α 2 τ V 1 τ g 3 M 1 ( α 2 ) 2 h 1 α 2 τ g M 1 2 ( α 2 ) 2 h 1 α 2 τ g M 1 3 ( α 2 ) 2 h 1 α 2 i τ g 0 ( α 2 ) 2 h 1 α 2 τ V M 1 .
The description of these recursive formulas will be described next. The initial steps Equations (17) and (18) will be observed also in our implementation. Moreover, the computational implementation will make use of the fact that the coefficients ( g k ( α ) ) k = can be expressed in alternative form as
g 0 ( α ) = Γ ( α + 1 ) Γ ( α / 2 + 1 ) 2
and
g k + 1 ( α ) = 1 α + 1 α / 2 + k + 1 g k ( α ) , k N { 0 } .
The following conventions and definitions will be observed in our code:
  • Input:
    • a1 = a 1 ,
    • b1 = b 1 ,
    • T = T ,
    • betta11 = β 11 ,
    • betta12 = β 12 ,
    • betta22 = β 22 ,
    • lambda = λ ,
    • D = D ,
    • tau = τ ,
    • h1 = h 1 ,
    • P = P ,
    • u1 = ϕ 1 , and
    • v1 = ϕ 2 .
  • Output:
    • u3 = u N and
    • v3 = v N .
a1=−10;
b1=10;
T=0.5;
 
alpha1=1.5;
alpha2=1.5;
beta11=1.5;
beta12=0.5;
beta22=1.5;
lambda=−0.5;
D=2;
tau=0.01;
h1=0.1;
 
x=a1:h1:b1;
t=0:tau:T;
M=length(x);
N=length(t);
 
P=0.5.*x.2;
 
u1=exp(-x.2)./sqrt(pi);
v1=exp(-x.2)./sqrt(pi);
 
ga1=zeros(1,M);
ga2=zeros(1,M);
ga1(1)=gamma(alpha1+1)/gamma(0.5*alpha1+1)2/h1alpha1;
ga2(1)=gamma(alpha2+1)/gamma(0.5*alpha2+1)2/h1alpha2;
for k=1:M-1
ga1(k+1)=(1-(alpha1+1)/(0.5*alpha1+k))*ga1(k);
ga2(k+1)=(1-(alpha2+1)/(0.5*alpha2+k))*ga2(k);
end
 
Ha1=zeros(M,M);
Ha2=zeros(M,M);
for j=1:M
for k=1:M
Ha1(j,k)=ga1(abs(j-k)+1);
Ha2(j,k)=ga2(abs(j-k)+1);
end
end
 
A=diag(0.5*tau.*(P+D)-1i)+0.25.*tau.*Ha1;
C=diag(0.5*tau.*(P+D)+1i)+0.25.*tau*Ha1;
B=diag(0.5*tau.*P-1i)+0.25.*tau.*Ha2;
D=diag(0.5*tau.*P+1i)+0.25.*tau.*Ha2;
 
u2=u1+1i.*tau.*(-0.5.*u1*Ha1-P.*u1-u1*D-lambda*v1... -(beta11.*abs(u1).2+beta12.*abs(v1).2).*u1);
v2=u1+1i.*tau.*(-0.5.*v1*Ha2-P.*v1-lambda*u1... -(beta12.*abs(u1).2+beta22.*abs(v1).2).*v1);
 
for n=3:N
meanu=0.5.*(3.*u2-u1);
meanv=0.5.*(3.*v2-v1);
c=(beta11.*abs(meanu).2+beta12.*abs(meanv).2).*meanu+lambda.*meanv;
d=(beta12.*abs(meanu).2+beta22.*abs(meanv).2).*meanv+lambda.*meanu;
 
u3=linsolve(A,-C*u2’-tau.*c’)’;
v3=linsolve(B,-D*v2’-tau.*d’)’;
u3(1)=0;
u3(M)=0;
 
u1=u2;
u2=u3;
v1=v2;
v2=v3;
end
end

References

  1. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Elsevier: Amsterdam, The Netherlands, 1998; Volume 198. [Google Scholar]
  2. Bao, W.; Markowich, P.A.; Wang, H. Ground, symmetric and central vortex states in rotating Bose-Einstein condensates. Commun. Math. Sci. 2005, 3, 57–88. [Google Scholar] [CrossRef]
  3. 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. 2016, 325, 74–97. [Google Scholar] [CrossRef] [Green Version]
  4. Bao, W.; Cai, Y. Optimal error estimates of finite difference methods for the Gross-Pitaevskii equation with angular momentum rotation. Math. Comput. 2013, 82, 99–128. [Google Scholar] [CrossRef] [Green Version]
  5. Cerimele, M.M.; Chiofalo, M.L.; Pistella, F.; Succi, S.; Tosi, M.P. Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates. Phys. Rev. E 2000, 62, 1382. [Google Scholar] [CrossRef] [Green Version]
  6. Myatt, C.; Burt, E.; Ghrist, R.; Cornell, E.A.; Wieman, C. Production of two overlapping Bose-Einstein condensates by sympathetic cooling. Phys. Rev. Lett. 1997, 78, 586. [Google Scholar] [CrossRef]
  7. Hall, D.; Matthews, M.; Ensher, J.; Wieman, C.; Cornell, E.A. Dynamics of component separation in a binary mixture of Bose-Einstein condensates. Phys. Rev. Lett. 1998, 81, 1539. [Google Scholar] [CrossRef] [Green Version]
  8. Modugno, G.; Ferrari, G.; Roati, G.; Brecha, R.J.; Simoni, A.; Inguscio, M. Bose-Einstein condensation of potassium atoms by sympathetic cooling. Science 2001, 294, 1320–1322. [Google Scholar] [CrossRef]
  9. Mudrich, M.; Kraft, S.; Singer, K.; Grimm, R.; Mosk, A.; Weidemüller, M. Sympathetic cooling with two atomic species in an optical trap. Phys. Rev. Lett. 2002, 88, 253001. [Google Scholar] [CrossRef] [Green Version]
  10. Haas, M.; Leung, V.; Frese, D.; Haubrich, D.; John, S.; Weber, C.; Rauschenbeutel, A.; Meschede, D. Species-selective microwave cooling of a mixture of rubidium and caesium atoms. New J. Phys. 2007, 9, 147. [Google Scholar] [CrossRef]
  11. Ortigueira, M.D. Riesz potential operators and inverses via fractional centred derivatives. Int. J. Math. Math. Sci. 2006, 2006, 048391. [Google Scholar] [CrossRef]
  12. Wang, X.; Liu, F.; Chen, X. Novel second-order accurate implicit numerical methods for the Riesz space distributed-order advection-dispersion equations. Adv. Math. Phys. 2015, 2015, 590435. [Google Scholar] [CrossRef]
  13. Gorenflo, R.; Mainardi, F. Fractional calculus. In Fractals and Fractional Calculus in Continuum Mechanics; Springer: Berlin/Heidelberg, Germany, 1997; pp. 223–276. [Google Scholar]
  14. Nonnenmacher, T.F.; Metzler, R. On the Riemann–Liouville fractional calculus and some recent applications. Fractals 1995, 3, 557–566. [Google Scholar] [CrossRef]
  15. Dalir, M.; Bashour, M. Applications of fractional calculus. Appl. Math. Sci. 2010, 4, 1021–1032. [Google Scholar]
  16. Sun, H.; Zhang, Y.; Baleanu, D.; Chen, W.; Chen, Y. A new collection of real world applications of fractional calculus in science and engineering. Commun. Nonlinear Sci. Numer. Simul. 2018, 64, 213–231. [Google Scholar] [CrossRef]
  17. Bagley, R.L.; Torvik, P. A theoretical basis for the application of fractional calculus to viscoelasticity. J. Rheol. 1983, 27, 201–210. [Google Scholar] [CrossRef]
  18. Rossikhin, Y.A.; Shitikova, M.V. Application of fractional calculus for dynamic problems of solid mechanics: Novel trends and recent results. Appl. Mech. Rev. 2010, 63. [Google Scholar] [CrossRef]
  19. Fallahgoul, H.; Focardi, S.; Fabozzi, F. Fractional Calculus and Fractional Processes with Applications to Financial Economics: Theory and Application; Academic Press: Cambridge, MA, USA, 2016. [Google Scholar]
  20. Scalas, E.; Gorenflo, R.; Mainardi, F. Fractional calculus and continuous-time finance. Phys. Stat. Mech. Appl. 2000, 284, 376–384. [Google Scholar] [CrossRef] [Green Version]
  21. Zhang, Y.; Sun, H.; Stowell, H.H.; Zayernouri, M.; Hansen, S.E. A review of applications of fractional calculus in Earth system dynamics. Chaos Solitons Fractals 2017, 102, 29–46. [Google Scholar] [CrossRef]
  22. Ionescu, C.; Lopes, A.; Copot, D.; Machado, J.T.; Bates, J. The role of fractional calculus in modeling biological phenomena: A review. Commun. Nonlinear Sci. Numer. Simul. 2017, 51, 141–159. [Google Scholar] [CrossRef]
  23. Nowakowski, J.; Ostalczyk, P.; Sankowski, D. Application of fractional calculus for modelling of two-phase gas/liquid flow system. Inform. Autom. Pomiary Gospod. Ochr. ŚRodowiska 2017, 7, 42–45. [Google Scholar] [CrossRef]
  24. Tarasov, V.E. Continuous limit of discrete systems with long-range interaction. J. Phys. Math. Gen. 2006, 39, 14895. [Google Scholar] [CrossRef] [Green Version]
  25. Macías-Díaz, J.E.; Bountis, A. Supratransmission in β-Fermi–Pasta–Ulam chains with different ranges of interactions. Commun. Nonlinear Sci. Numer. Simul. 2018, 63, 307–321. [Google Scholar] [CrossRef]
  26. Christodoulidi, H.; Tsallis, C.; Bountis, T. Fermi–Pasta–Ulam model with long-range interactions: Dynamics and thermostatistics. EPL 2014, 108, 40006. [Google Scholar] [CrossRef] [Green Version]
  27. Bountis, T. From mechanical to biological oscillator networks: The role of long range interactions. Eur. Phys. J. Spec. Top. 2016, 225, 1017–1035. [Google Scholar] [CrossRef]
  28. Ortigueira, M.D. Fractional central differences and derivatives. Ifac Proc. Vol. 2006, 39, 58–63. [Google Scholar] [CrossRef]
  29. Hao, Z.P.; Sun, Z.Z.; Cao, W.R. A fourth-order approximation of fractional derivatives with its applications. J. Comput. Phys. 2015, 281, 787–805. [Google Scholar] [CrossRef]
  30. Lin, F.R.; Wang, Q.Y.; Jin, X.Q. Crank-Nicolson-weighted-shifted-Grünwald-difference schemes for space Riesz variable-order fractional diffusion equations. Numer. Algorithms 2020, 1–31. [Google Scholar] [CrossRef]
  31. Hendy, A.S.; Macías-Díaz, J.E. An efficient Hamiltonian numerical model for a fractional Klein–Gordon equation through weighted-shifted Grünwald differences. J. Math. Chem. 2019, 57, 1394–1412. [Google Scholar] [CrossRef]
  32. Macías-Díaz, J.E. A numerically efficient variational algorithm to solve a fractional nonlinear elastic string equation. Numer. Algorithms 2020, 1–28. [Google Scholar] [CrossRef]
  33. Martínez, R.; Macías-Díaz, J.E. An energy-preserving and efficient scheme for a double-fractional conservative Klein–Gordon–Zakharov system. Appl. Numer. Math. 2020, 158, 292–313. [Google Scholar] [CrossRef]
  34. Wang, Y.M.; Ren, L. A high-order L2-compact difference method for Caputo-type time-fractional sub-diffusion equations with variable coefficients. Appl. Math. Comput. 2019, 342, 71–93. [Google Scholar] [CrossRef]
  35. Diethelm, K. The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2010. [Google Scholar]
  36. Oldham, K.; Spanier, J. The Fractional Calculus Theory and Applications of Differentiation and Integration to Arbitrary Order; Elsevier: Amsterdam, The Netherlands, 1974. [Google Scholar]
  37. Gao, G.H.; Sun, Z.Z.; Zhang, H.W. A new fractional numerical differentiation formula to approximate the Caputo fractional derivative and its applications. J. Comput. Phys. 2014, 259, 33–50. [Google Scholar] [CrossRef]
  38. Yan, Y.; Sun, Z.Z.; Zhang, J. Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations: A second-order scheme. Commun. Comput. Phys. 2017, 22, 1028–1048. [Google Scholar] [CrossRef]
  39. Macías-Díaz, J.E. Nonlinear wave transmission in harmonically driven hamiltonian sine-Gordon regimes with memory effects. Chaos Solitons Fractals 2020, 110362. [Google Scholar] [CrossRef]
  40. Murillo, J.Q.; Yuste, S.B. An explicit difference method for solving fractional diffusion and diffusion-wave equations in the Caputo form. J. Comput. Nonlinear Dyn. 2011, 6. [Google Scholar] [CrossRef]
  41. Macías-Díaz, J.E. Numerical study of the process of nonlinear supratransmission in Riesz space-fractional sine-Gordon equations. Commun. Nonlinear Sci. Numer. Simul. 2017, 46, 89–102. [Google Scholar] [CrossRef]
  42. De Pablo, A.; Quirós, F.; Rodríguez, A.; Vázquez, J.L. A fractional porous medium equation. Adv. Math. 2011, 226, 1378–1409. [Google Scholar] [CrossRef] [Green Version]
  43. Vázquez, J.L.; de Pablo, A.; Quirós, F.; Rodríguez, A. Classical solutions and higher regularity for nonlinear fractional diffusion equations. J. Eur. Math. Soc. 2017, 19, 1949–1975. [Google Scholar] [CrossRef] [Green Version]
  44. Bonforte, M.; Sire, Y.; Vázquez, J.L. Optimal existence and uniqueness theory for the fractional heat equation. Nonlinear Anal. Theory Methods Appl. 2017, 153, 142–168. [Google Scholar] [CrossRef] [Green Version]
  45. Stan, D.; Vázquez, J.L. The Fisher–KPP equation with nonlinear fractional diffusion. Siam J. Math. Anal. 2014, 46, 3241–3276. [Google Scholar] [CrossRef] [Green Version]
  46. Segatti, A.; Vázquez, J.L. On a fractional thin film equation. Adv. Nonlinear Anal. 2020, 9, 1516–1558. [Google Scholar] [CrossRef]
  47. Díaz, J.I.; Gómez-Castro, D.; Vázquez, J.L. The fractional Schrödinger equation with general nonnegative potentials. The weighted space approach. Nonlinear Anal. 2018, 177, 325–360. [Google Scholar] [CrossRef] [Green Version]
  48. Adda, F.B. Geometric interpretation of the differentiability and gradient of real order. C. R. L’Academie Des. Sci. Ser. Math. 1998, 8, 931–934. [Google Scholar]
  49. Adda, F.B. The differentiability in the fractional calculus. Nonlinear Anal. 2001, 47, 5423–5428. [Google Scholar] [CrossRef]
  50. Meerschaert, M.M.; Mortensen, J.; Wheatcraft, S.W. Fractional vector calculus for fractional advection–dispersion. Phys. Stat. Mech. Appl. 2006, 367, 181–190. [Google Scholar] [CrossRef]
  51. Tarasov, V.E. Fractional vector calculus and fractional Maxwell’s equations. Ann. Phys. 2008, 323, 2756–2778. [Google Scholar] [CrossRef] [Green Version]
  52. Ortigueira, M.D.; Rivero, M.; Trujillo, J.J. From a generalised Helmholtz decomposition theorem to fractional Maxwell equations. Commun. Nonlinear Sci. Numer. Simul. 2015, 22, 1036–1049. [Google Scholar] [CrossRef]
  53. Ortigueira, M.; Machado, J. On fractional vectorial calculus. Bull. Pol. Acad. Sci. Tech. Sci. 2018, 66. [Google Scholar] [CrossRef]
  54. Desplanques, J. Théoreme d’algébre. J. Math. Spec. 1887, 9, 12–13. [Google Scholar]
  55. Pen-Yu, K. Numerical methods for incompressible viscous flow. Sci. Sin. 1977, 20, 287–304. [Google Scholar]
  56. Macías-Díaz, J.E. A structure-preserving method for a class of nonlinear dissipative wave equations with Riesz space-fractional derivatives. J. Comput. Phys. 2017, 351, 40–58. [Google Scholar] [CrossRef]
Figure 1. Numerical approximations to the solutions of Equation (4) versus x and t, using the constants of Table 1. The graphs were obtained using the finite-difference method of Equation (16) with T = 10 , α 1 = α 2 = 2 , τ = 0.01 , and h 1 = 0.1 . We also let V ( x ) = 1 2 x 2 and ϕ 1 ( x ) = ϕ 2 ( x ) = 1 π exp ( x 2 ) , for each x Ω . The graphs provide (a) Re ( ψ 1 ) , (b) Re ( ψ 2 ) , (c) Im ( ψ 1 ) , (d) Im ( ψ 2 ) , (e) | ψ 1 | , and (f) | ψ 2 | .
Figure 1. Numerical approximations to the solutions of Equation (4) versus x and t, using the constants of Table 1. The graphs were obtained using the finite-difference method of Equation (16) with T = 10 , α 1 = α 2 = 2 , τ = 0.01 , and h 1 = 0.1 . We also let V ( x ) = 1 2 x 2 and ϕ 1 ( x ) = ϕ 2 ( x ) = 1 π exp ( x 2 ) , for each x Ω . The graphs provide (a) Re ( ψ 1 ) , (b) Re ( ψ 2 ) , (c) Im ( ψ 1 ) , (d) Im ( ψ 2 ) , (e) | ψ 1 | , and (f) | ψ 2 | .
Mathematics 09 01412 g001
Figure 2. Numerical approximations to the solutions of Equation (4) versus x and t, using the constants of Table 1. The graphs were obtained using the finite-difference method of Equation (16) with T = 10 , α 1 = α 2 = 1.5 , τ = 0.01 and h 1 = 0.1 . We also let V ( x ) = 1 2 x 2 and ϕ 1 ( x ) = ϕ 2 ( x ) = 1 π exp ( x 2 ) , for each x Ω . The graphs provide (a) Re ( Ψ 1 ) , (b) Re ( Ψ 2 ) , (c) Im ( Ψ 1 ) , (d) Im ( Ψ 2 ) , (e) | Ψ 1 | , and (f) | Ψ 2 | .
Figure 2. Numerical approximations to the solutions of Equation (4) versus x and t, using the constants of Table 1. The graphs were obtained using the finite-difference method of Equation (16) with T = 10 , α 1 = α 2 = 1.5 , τ = 0.01 and h 1 = 0.1 . We also let V ( x ) = 1 2 x 2 and ϕ 1 ( x ) = ϕ 2 ( x ) = 1 π exp ( x 2 ) , for each x Ω . The graphs provide (a) Re ( Ψ 1 ) , (b) Re ( Ψ 2 ) , (c) Im ( Ψ 1 ) , (d) Im ( Ψ 2 ) , (e) | Ψ 1 | , and (f) | Ψ 2 | .
Mathematics 09 01412 g002
Figure 3. Numerical approximations to the solutions of Equation (4) versus x and t, using the constants of Table 1. The graphs were obtained using the finite-difference method of Equation (16) with T = 10 , α 1 = α 2 = 1 , τ = 0.01 and h 1 = 0.1 . We also let V ( x ) = 1 2 x 2 and ϕ 1 ( x ) = ϕ 2 ( x ) = 1 π exp ( x 2 ) for each x Ω . The graphs provide (a) Re ( Ψ 1 ) , (b) Re ( Ψ 2 ) , (c) Im ( Ψ 1 ) , (d) Im ( Ψ 2 ) , (e) | Ψ 1 | , and (f) | Ψ 2 | .
Figure 3. Numerical approximations to the solutions of Equation (4) versus x and t, using the constants of Table 1. The graphs were obtained using the finite-difference method of Equation (16) with T = 10 , α 1 = α 2 = 1 , τ = 0.01 and h 1 = 0.1 . We also let V ( x ) = 1 2 x 2 and ϕ 1 ( x ) = ϕ 2 ( x ) = 1 π exp ( x 2 ) for each x Ω . The graphs provide (a) Re ( Ψ 1 ) , (b) Re ( Ψ 2 ) , (c) Im ( Ψ 1 ) , (d) Im ( Ψ 2 ) , (e) | Ψ 1 | , and (f) | Ψ 2 | .
Mathematics 09 01412 g003
Table 1. Table of the parameter values employed in the numerical simulations of this manuscript.
Table 1. Table of the parameter values employed in the numerical simulations of this manuscript.
p a 1 b 1 β 11 β 12 β 22 λ D
1 7 7 1.5 0.5 1.5 0.5 2
Table 2. Absolute errors and the rates of convergence (a) in time and (b) in space. Example 2 describes the experimental setting for these results.
Table 2. Absolute errors and the rates of convergence (a) in time and (b) in space. Example 2 describes the experimental setting for these results.
(a) Temporal Study of Convergence.
τ 0.08 h = 0.04 h = 0.02
ϵ t , h 2 ρ τ , h ϵ t , h 2 ρ τ , h ϵ t , h 2 ρ τ , h
0.04 7.1074 × 10 2 1.9124 × 10 2 3.9163 × 10 3
0.02 1.2825 × 10 2 2.4703 2.9995 × 10 3 2.6726 8.6646 × 10 4 2.1763
0.01 2.6863 × 10 3 2.2553 6.0564 × 10 4 2.3082 2.3468 × 10 4 1.8844
(b) Spatial Study of Convergence.
h τ = 0.001 τ = 0.0005 τ = 0.00025
ϵ t , h 2 σ τ , h ϵ t , h 2 σ τ , h ϵ t , h 2 σ τ , h
0.08 1.6737 × 10 5 3.9923 × 10 6 9.3696 × 10 7
0.04 3.3586 × 10 6 2.3171 8.3193 × 10 7 2.2627 2.4757 × 10 7 1.9201
0.02 7.5983 × 10 7 2.1441 2.2384 × 10 7 1.8940 6.8871 × 10 8 1.8459
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Serna-Reyes, A.J.; Macías-Díaz, J.E.; Reguera, N. A Convergent Three-Step Numerical Method to Solve a Double-Fractional Two-Component Bose–Einstein Condensate. Mathematics 2021, 9, 1412. https://0-doi-org.brum.beds.ac.uk/10.3390/math9121412

AMA Style

Serna-Reyes AJ, Macías-Díaz JE, Reguera N. A Convergent Three-Step Numerical Method to Solve a Double-Fractional Two-Component Bose–Einstein Condensate. Mathematics. 2021; 9(12):1412. https://0-doi-org.brum.beds.ac.uk/10.3390/math9121412

Chicago/Turabian Style

Serna-Reyes, Adán J., Jorge E. Macías-Díaz, and Nuria Reguera. 2021. "A Convergent Three-Step Numerical Method to Solve a Double-Fractional Two-Component Bose–Einstein Condensate" Mathematics 9, no. 12: 1412. https://0-doi-org.brum.beds.ac.uk/10.3390/math9121412

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop