Positivity and Boundedness Preserving Schemes for Space-Time Fractional Predator-Prey Reaction-Diffusion Model
aa r X i v : . [ m a t h . NA ] N ov doi:Volume X , Number , XX pp. X–XX
POSITIVITY AND BOUNDEDNESS PRESERVING SCHEMESFOR SPACE-TIME FRACTIONAL PREDATOR-PREYREACTION-DIFFUSION MODEL
YANYAN YU, WEIHUA DENG, AND YUJIANG WU
School of Mathematics and Statistics, Lanzhou University,Lanzhou, 730000, People’s Republic of China
Abstract.
The semi-implicit schemes for the nonlinear predator-prey reaction-diffusion model with the space-time fractional derivatives are discussed, wherethe space fractional derivative is discretized by the fractional centered differ-ence and WSGD scheme. The stability and convergence of the semi-implicitschemes are analyzed in the L ∞ norm. We theoretically prove that the numer-ical schemes are stable and convergent without the restriction on the ratio ofspace and time stepsizes and numerically further confirm that the schemes havefirst order convergence in time and second order convergence in space. Thenwe discuss the positivity and boundedness properties of the analytical solu-tions of the discussed model, and show that the numerical solutions preservethe positivity and boundedness. The numerical example is also presented. Introduction.
The predator-prey model, also known as the Lotka-Volterraequation, is a pair of first-order nonlinear ordinary differential equations frequentlyused to describe the dynamics of biological system in which two species interact, oneas the predator and the other as prey [1]; and the unknown variables usually denotecertain measure of total population. To interpret the unknown variables as spacialdensities so as to allow the population size to vary throughout the considered region,Conway and Smoller introduce diffusion terms to the predator-prey model whichallows diffusing as well as interacting each other [2]; and then the predator-preyreaction-diffusion model is obtained. With the development of the predator-preymodel, it seems nature to introduce diffusion terms to the corresponding model,e.g., the Michaelis-Menten-Holling predator-prey reaction-diffusion model [3, 4]; andmore general other models [5, 6, 7, 8].Anomalous diffusion, including subdiffusion and superdiffuion, is also a diffusionprocess, but its mean squared displacement (MSD) is nonlinear with respect totime t , in contract to the classical diffusion process, in which the MSD is linear [9];nowadays, it is widely recognized that the anomalous diffusion is ubiquitous, e.g.,diffusion through porous media, protein diffusion within cells, and also being foundin many other biological systems. So it seems reasonable/nature to introduce theanomalous diffusion to the predator-prey model. The subdiffusion is introduced tothe Michaelis-Menten-Holling predator-prey model in [10] to get a new model; it is Mathematics Subject Classification.
Primary: 65M06, 26A33; Secondary: 45M20.
Key words and phrases. semi-implicit scheme, fractional predator-prey model, fractional cen-tered difference, WSGD scheme, convergence. proved that the solution of the model is positive and bounded; and the numericalschemes preserving the positivity and boundedness are detailedly discussed. Here weintroduce both the subdiffusion and superdiffusion to the Michaelis-Menten-Hollingpredator-prey model, and the system can be written as ∂ α N∂t α = D ∂ β N∂ | x | β + N (cid:18) − N − ̺PP + N (cid:19) , x ∈ ( l, r ) , t > ,∂ α P∂t α = D ∂ β P∂ | x | β + σP (cid:18) − γ + κδP κP + NP + N (cid:19) , x ∈ ( l, r ) , t > , (1)with the Caputo derivative in time and Riesz space fractional derivetive, where ̺ , σ and κ are positive real numbers and N and P denote the population densities of preyand predator respectively. Based on the practical applications, we are interested inthe solutions of (1) with the nonnegative initial conditions N ( x,
0) = g ( x ) ≥ , P ( x,
0) = g ( x ) ≥ , x ∈ ( l, r ) , (2)and the homogeneous Neumann boundary conditions( ∂N ( x, t ) /∂x ) | x = l and r , respectively = ( ∂P ( x, t ) /∂x ) | x = l and r , respectively = 0 . (3)The positive constants γ and δ in the coupled equations denote the minimal mortal-ity and the limiting mortality of the predator, respectively. Throughout the paper,we assume that γ satisfies the natural condition 0 < γ ≤ δ and consider the case ofthe diffusion constants D i > , i = 1 , RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 3 two schemes preserve the positivity and boundedness. The numerical experimentsto confirm the convergent orders and the positivity and boundedness preserving areperformed in Section 5. And we conclude the paper with some discussions in thelast section.2.
Numerical schemes for the space-time fractional predator-prey reaction-diffusion model.
As mentioned in the introduction section, the model we discussis as follows: ∂ α N∂t α = D ∂ β N∂ | x | β + N (cid:18) − N − ̺PP + N (cid:19) , x ∈ ( l, r ) , < t ≤ T,∂ α P∂t α = D ∂ β P∂ | x | β + σP (cid:18) − γ + κδP κP + NP + N (cid:19) , x ∈ ( l, r ) , < t ≤ T, (4)with 0 < α <
1, 1 < β <
2, and use the following initial conditions N ( x,
0) = g ( x ) , P ( x,
0) = g ( x ) , x ∈ [ l, r ] (5)and the homogeneous boundary conditions ∂N ( x, t ) ∂x | x = l = ∂N ( x, t ) ∂x | x = r = ∂P ( x, t ) ∂x | x = l = ∂P ( x, t ) ∂x | x = r = 0 , < t ≤ T. (6)In (4), the time fractional operator ∂ α u ( x,t ) ∂t α denotes the Caputo fractional derivative ∂ α u ( x, t ) ∂x α = 1Γ(1 − α ) Z t ∂u ( x, s ) ∂s t − s ) α ds, and ∂ β u∂ | x | β denotes the Riesz fractional derivative ∂ β u ( x, t ) ∂ | x | β = − cos ( βπ/ − β ) d dx Z rl | x − ξ | − β u ( ξ, t ) dξ. (7)For ease of presentation, we uniformly divide the spacial domain [ l, r ] into M subintervals with stepsize h (= ( r − l ) /M ) and the time domain [0 , T ] into N subintervals with steplength τ (= T /N ). Let x i = l + ih ( i = 0 , , · · · , M ), t k = kτ ( k = 0 , , · · · , N ), and denote the grid function as u ki = u ( x i , t k ). We use thediscrete scheme of the Caputo derivative in the following form [16] D ατ u ki = τ − α Γ(2 − α ) (cid:2) u ki − k − X n =1 ( b k − n − − b k − n ) u ni − b k − u i (cid:3) , (8)where b n = ( n + 1) − α − n − α >
0, and 1 = b > b > · · · > b k − > (1 − α ) k − α .Note that the above discrete scheme has 2 − α order accuracy in time, i.e., ∂ α u ( x i , t ) ∂t α (cid:12)(cid:12) t = t k = D ατ u ki + O ( τ − α ) . (9)For the Riesz space fractional derivative, we introduce two second order finite dif-ference approximation schemes (fractional centered difference scheme and secondorder WSGD scheme) in the next parts. YANYAN YU, WEIHUA DENG, AND YUJIANG WU
Fractional centered difference scheme.
In this part, we use the fractionalcentered difference [19] to discretize the spatial factional derivative. The Rieszfractional derivative can be redefined as ∂ β u ( x, t ) ∂ | x | β = lim h → − h β [ x/h ] X j =[( l − r + x ) /h ] ( − j Γ( β + 1)Γ( β − j + 1)Γ( β + j + 1) u ( x − jh, t ) . (10)Denoting the weights as g j = ( − j Γ( β + 1)Γ( β − j + 1)Γ( β + j + 1) , we get the fractional centered difference approximation for the Riesz fractionalderivative δ βx u ki = − h β i X j = − M + i g j u ki − j , (11)where the weights g j satisfy the following properties:1. g ≥ g − j = g j ≤ g j +1 = (cid:16) − β +1 β/ j +1 (cid:17) g j , and | g j +1 | < | g j | ,4. ∞ P j = −∞ g j = 0, and i P j = − M + ij =0 | g j | < g .From [20], we know that the accuracy of fractional centered difference approximationis as, Lemma 2.1. If u ( x, t i ) ∈ C [ l, r ] , then we have ∂ β u ( x, t i ) ∂ | x | β | x = x i = δ βx u ki + O ( h ) , with < β < . If the fractional centered difference approximation is substituted into the systemof reaction-diffusion equation (4), then the resulting semi-implicit finite differencescheme is D ατ N k +1 i = D δ βx N k +1 i + N ki (cid:18) − N ki − ̺P ki P ki + N ki (cid:19) ,D ατ P k +1 i = D δ βx P k +1 i + σP ki (cid:18) − γ + κδP ki κP ki + N ki P ki + N ki (cid:19) , (12)with initial conditions N i = g ( x i ) , P i = g ( x i ) . (13)For guaranteeing the second-order accuracy in space, we give the three point inter-polation scheme for the Neumann boundary conditions: N k − N k + 3 N k = 0 , P k − P k + 3 P k = 0 , ≤ k ≤ N,N kM − − N kM − + 3 N kM = 0 , P kM − − P kM − + 3 P kM = 0 , ≤ k ≤ N. (14) RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 5
Eq. (12)-(14) may be rearranged and written as matrix form and solved to get thenumerical solution at the time step t k +1 :( I + A ) N k +1 = k X j =1 ( b j − b j +1 ) N k − j + b k N + µN k (cid:18) − N k − ̺P k P k + N k (cid:19) , ( I + A ) P k +1 = k X j =1 ( b j − b j +1 ) P k − j + b k P + µσP k (cid:18) − γ + κδP k κP k + N k P k + N k (cid:19) , where N k = ( N k , N k , · · · , N kM − ), P k = ( P k , P k , · · · , P kM − ) and A i is a ( M − × ( M −
1) square matrix A i = µD i h β g + g g − − g g − · · · g − M +2 + g − M +1 g + g g − g g − g − M +3 + g − M +2 ... . . . g M − + g M − g M − − g M − g M − g + g − . (15)2.2. Second order WSGD scheme.
The other equivalent definition of the Rieszfractional derivative ∂ β u ( x,t ) ∂ | x | β is ∂ β u ( x, t ) ∂ | x | β = − cos ( βπ/
2) ( l D βx u ( x, t ) + x D βr u ( x, t )) , (16)where l D βx u ( x, t ) and x D βr u ( x, t ) denote the left and right Riemann-Liouville frac-tional derivatives of the function u ( x, t ), respectively. According to the discus-sions of [22], we present the discrete approximations of the left and right Riemann-Liouville fractional derivatives as follows: l D βx u ( x i ) = 1 h β i +1 X j =0 w j u ( x i − j +1 ) + O ( h ) , x D βr u ( x i ) = 1 h β M − i +1 X j =0 w j u ( x i + j − ) + O ( h ) , (17)where w = α/ w j = Γ( j − β − / (Γ( − β )Γ( j )) · (1 − ( β/ · ( β + 1) /j ). Af-ter substituting (17) into (16) and merging the same terms, we get the discreteapproximation of the Riesz fractional derivative δ x u ki = − h β i X j = − M + i θ j u ki − j , where the weights θ j satisfy the following properties1. θ = w cos ( βπ/ = − β − β cos ( βπ/ > θ = θ − = w + w cos ( βπ/ = β ( β +2)( β − cos ( βπ/ < | j | > θ j = θ − j = w j +1 cos ( βπ/ = Γ( j − β )2 cos ( βπ/ − β )Γ( j +1) (cid:16) − β ( β +1)2( j +1) (cid:17) ≤ YANYAN YU, WEIHUA DENG, AND YUJIANG WU θ j +1 = ( − β + j )( β + β − j ))(2+ j )( β + β − j )) θ j , for | j | ≥
2, and | θ | > | θ | ≥ | θ | ≥ · · · ,4. ∞ P j = −∞ θ j = 0, and i P j = − M + ij =0 | θ j | < θ .In view of the approximations (17), we know that the above discrete scheme of theRiesz fractional derivative also has second-order accuracy ∂ β u ( x, t ) ∂ | x | β (cid:12)(cid:12) x = x k ,t = t j = − h β i X j = − M + i θ j u ki − j + O ( h ) . (18)If the second-order WSGD discretization is substituted into the reaction-diffusionequations (4), we obtain the new semi-implicit WSGD finite difference scheme. Andthe same initial and boundary discretizations (13) and (14) are used. Then we canget the same formulations for the WSGD scheme as the ones for the fractionalcentered difference scheme (12)-(15).The semi-implicit fractional centered difference and the semi-implicit WSGDscheme have the same structure, and in particular, their discretized weights havethe similar properties. So in the following, the two schemes have the same analyses,and we just focus on the first scheme.3. Stability and convergence of the numerical schemes.
Now we first denotethe two nonlinear terms by f ( N, P ) = N (cid:18) − N − ̺PP + N (cid:19) (19)and f ( N, P ) = σP (cid:18) − γ + κδP κP + NP + N (cid:19) (20)for convenience. And assume that they satisfy local Lipschitz condition, i.e., thereexists a positive constant L such that | f ( N , P ) − f ( N , P ) | ≤ L ( | N − N | + | P − P | )and | f ( N , P ) − f ( N , P ) | ≤ L ( | N − N | + | P − P | ) , RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 7 when | N − N | ≤ ǫ and | P − P | ≤ ǫ for a given positive constant ǫ . Then thediscrete scheme (12) can be rewritten as the following form (cid:18) µ + D g h β (cid:19) N k +1 i = 1 µ k X n =1 ( b k − n − b k − n +1 ) N ni + 1 µ b k N i − D h β i X j = − M + ij =0 g j N k +1 i − j + f ( N ki , P ki ) , (cid:18) µ + D g h β (cid:19) P k +1 i = 1 µ k X n =1 ( b k − n − b k − n +1 ) P ni + 1 µ b k P i − D h β i X j = − M + ij =0 g j P k +1 i − j + f ( N ki , P ki ) , (21)where µ = Γ(2 − α ) τ α . In this paper, we use L ∞ norm of { u ki } M − i =1 which is definedas k u k k = max ≤ i ≤ M − | u ki | . Let ( ¯ N ki , ¯ P ki ) be the approximate solution to the numerical scheme and denote ǫ ki = N ki − ¯ N ki and ε ki = P ki − ¯ P ki . Then under the small perturbations, there existsthe following numerical stability result. Theorem 3.1.
The numerical schemes (12)-(14) are stable and there exist k ǫ k k ≤ − α ) − T α Γ(2 − α ) L ( k ǫ k + k ε k ) , k ε k k ≤ − α ) − T α Γ(2 − α ) L ( k ǫ k + k ε k ) , (22) when T < (1 / (Γ(1 − α ) L )) /α .Proof. We prove this theorem by mathematical induction. Using the first equationof (21), we have (cid:18) µ + D g h β (cid:19) ǫ k +1 i = 1 µ k X n =1 ( b k − n − b k − n +1 ) ǫ ni + 1 µ b k ǫ i − D h β i X j = − M + ij =0 g j ǫ k +1 i − j + ( f ( N ki , P ki ) − f ( ¯ N ki , ¯ P ki )) . YANYAN YU, WEIHUA DENG, AND YUJIANG WU
According to the numerical approximations (14), when i = 2 , M −
2, the aboveequation can be rewritten as (cid:18) µ + D g h β (cid:19) ǫ k +1 i = 1 µ k X n =1 ( b k − n − b k − n +1 ) ǫ ni + 1 µ b k ǫ i − D h β i − X j = − M + i +1 j =0 g j ǫ k +1 i − j + D h β ( g i ǫ k +12 + g − M + i ǫ k +1 M − ) − D h β ( g i ǫ k +11 + g − M + i ǫ k +1 M − )+ ( f ( N ki , P ki ) − f ( ¯ N ki , ¯ P ki ))= 1 µ k X n =1 ( b k − n − b k − n +1 ) ǫ ni + 1 µ b k ǫ i − D h β i − X j = − M + i +1 j =0 ,i − , − M + i +2 g j ǫ k +1 i − j + D h β (( g i / − g i − ) ǫ k +12 + ( g − M + i / − g − M + i +2 ) ǫ k +1 M − ) − D h β ( g i ǫ k +11 + g − M + i ǫ k +1 M − ) + ( f ( N ki , P ki ) − f ( ¯ N ki , ¯ P ki )) . Since g i / − g i − > g − M + i / − g − M + i +2 >
0, we get (cid:18) µ + D g h β (cid:19) k ǫ k +1 i k≤ µ k X n =1 ( b k − n − b k − n +1 ) k ǫ ni k + 1 µ b k k ǫ i k − D h β i − X j = − M + i +1 j =0 ,i − , − M + i +2 g j k ǫ k +1 i − j k + D h β (( g i / − g i − ) k ǫ k +12 k + ( g − M + i / − g − M + i +2 ) k ǫ k +1 M − k ) − D h β ( g i k ǫ k +11 k + g − M + i k ǫ k +1 M − k ) + L ( k N ki − ¯ N ki k + k P ki − ¯ P ki k ) . From the above inequality we know that (cid:18) µ + D g h β (cid:19) k ǫ k +1 k ≤ µ k X n =1 ( b k − n − b k − n +1 ) k ǫ n k + 1 µ b k k ǫ k + L ( k ǫ k k + k ε k k ) − D h β i X j = − M + ij =0 g j k ǫ k +1 k . Together with the fourth property of coefficient g i , we can obtain the followinginequality (cid:18) µ + D g h β (cid:19) k ǫ k +1 k ≤ µ k X n =1 ( b k − n − b k − n +1 ) k ǫ n k + 1 µ b k k ǫ k + D g h β k ǫ k +1 k + L ( k ǫ k k + k ε k k ) . RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 9
Therefore, k ǫ k +1 k ≤ k X n =1 ( b k − n − b k − n +1 ) k ǫ n k + b k k ǫ k + Lµ ( k ǫ k k + k ε k k ) . (23)When i = 2, we have (cid:18) µ + D g h β (cid:19) ǫ k +12 − D h β g / ǫ k +12 = 1 µ k X n =1 ( b k − n − b k − n +1 ) ǫ n + 1 µ b k ǫ − D h β X j = − M +3 j =0 , − M +4 g j ǫ k +12 − j + D h β ( g − M +2 / − g − M +4 ) ǫ k +1 M − − D h β ( g ǫ k +11 + g − M +2 ǫ k +1 M − ) + ( f ( N k , P k ) − f ( ¯ N k , ¯ P k )) . Then in the similar way as above, we can also get the estimate (23) for i = 2. If i = M −
2, there exists a similar argument.By almost the same deduction as k ǫ k +1 k , we can get k ε k +1 k ≤ k X n =1 ( b k − n − b k − n +1 ) k ε n k + b k k ε k + Lµ ( k ǫ k k + k ε k k ) . (24)For k = 1, it’s easy to check that (22) holds. When k >
1, suppose (22) holds for n = 1 , · · · , k . Then there exists k ǫ k +1 k ≤ (cid:18) − b k + 2 Lµ (1 − α ) − T α Γ(2 − α ) L + b k (cid:19) ( k ǫ k + k ε k ) . We just need to prove1 − b k + 2 Lµ + b k ((1 − α ) − T α Γ(2 − α ) L ) ≤ , i.e., 2 Lµ · b − k ≤ α + T α Γ(2 − α ) L. Using µ = Γ(2 − α ) τ α = Γ(2 − α ) T α N α , b − k < ( k +1) α − α and Γ(1 − α ) T α L <
1, we canobtain 2 Lµ · b − k < L · Γ(2 − α ) T α N α · ( k + 1) α − α ≤ L · Γ(1 − α ) T α , and α + T α Γ(2 − α ) L > L · Γ(1 − α ) T α . This implies that 1 − b k + 2 Lµ + b k ((1 − α ) − T α Γ(2 − α ) L ) ≤ . So k ǫ k +1 k ≤ − α ) − T α Γ(2 − α ) L ( k ǫ k + k ε k ) , and k ε k +1 k ≤ − α ) − T α Γ(2 − α ) L ( k ǫ k + k ε k ) . Theorem 3.2.
Let { ( N ( x i , t k ) , P ( x i , t k )) } M − i =1 and { ( N ki , P ki ) } M − i =1 be the exactsolutions of the subdiffusive reaction-diffusion equation (4)-(6) and of the numericalschemes (12)-(14) respectively and define ρ ki = N ( x i , t k ) − N ki and η ki = P ( x i , t k ) − P ki . Then ρ k and η k satisfy the following error estimates: k ρ n k ≤ C ( τ + h ) and k η n k ≤ C ( τ + h ) , (25) when T < (1 / (Γ(1 − α ) L )) /α .Proof. According to the first equation of (21), when i = 2 , M −
2, we can get thefollowing inequality (cid:18) µ + D g h β (cid:19) | ρ k +1 i |≤ µ k X n =1 ( b k − n − b k − n +1 ) | ρ ni | + 1 µ b k | ρ i | − D h β i − X j = − M + i +1 j =0 ,i − , − M + i +2 g j | ρ k +1 i − j | + D h β (( g i / − g i − ) | ρ k +12 | + ( g − M + i / − g − M + i +2 ) | ρ k +1 M − | ) − D h β ( g i | ρ k +11 | + g − M + i | ρ k +1 M − | ) + L ( | ρ ki | + | η ki | ) + r k +11 i , where r k +11 = O ( τ + h ). From the above inequality we know that (cid:18) µ + D g h β (cid:19) k ρ k +1 k≤ µ k X n =1 ( b k − n − b k − n +1 ) k ρ n k + 1 µ b k k ρ k − D h β i X j = − M + ij =0 g j k ρ k +1 k + L ( k ρ k k + k η k k ) + k r k +11 k . In the view of the fourth property of g i , we obtain k ρ k +1 k ≤ k X n =1 ( b k − n − b k − n +1 ) k ρ n k + b k k ǫ k + Lµ ( k ρ k k + k η k k ) + µ k r k +11 k , and k η k +1 k ≤ k X n =1 ( b k − n − b k − n +1 ) k η n k + b k k η k + Lµ ( k ρ k k + k η k k ) + µ k r k +12 k , where r k +11 = O ( τ + h ) and r k +12 = O ( τ + h ). Similar to the discussions of thenumerical stability for the special cases i = 2 , M −
2, we know that the above
RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 11 inequalities hold for i = 1 , , · · · , M −
1. Now suppose that k ρ m k ≤ b − m − (1 − α ) − T α Γ(2 − α ) L ! ( k ρ k + k η k + max ≤ n ≤ mi =1 , µ || r ni || ) (26)and k ρ m k ≤ b − m − (1 − α ) − T α Γ(2 − α ) L ! ( k ρ k + k η k + max ≤ n ≤ mi =1 , µ || r ni || ) (27)hold for m = 1 , , · · · , k . Next we prove that the estimates (26) and (27) hold for m = k + 1.When m = 1, it’s easy to check that (26) holds. When m = k + 1, together with b − k − n − < b − n , we have || ρ k +1 || ≤ k X n =0 ( b n − b n +1 ) b − k − n − (1 − α ) − T α Γ(2 − α ) L ! ( k ρ k + k η k + max ≤ n ≤ ki =1 , µ || r ni || ) + b n || ρ || + 2 Lµb − n − (1 − α ) − T α Γ(2 − α ) L ( k ρ k + k η k + max ≤ n ≤ ki =1 , µ || r ni || ) + µ || r k +1 ||≤ (cid:18) b − n (1 − α ) − T α Γ(2 − α ) L (cid:19) ( k ρ k + k η k + max ≤ n ≤ ki =1 , µ || r ni || ) + µ || r k +1 || . Notice that ρ i = 0 , η i = 0 , ≤ i ≤ M. Together with b − n ≤ ( n +1) α − α and ( n + 1) τ ≤ T , we obtain || ρ k +1 || ≤ C ( τ + h ) . In a similar way, we can get || η k +1 || ≤ C ( τ + h ) . Positivity and boundedness of the analytical and numerical solutionsof the space-time fractional predator-prey reaction-diffusion model.
Inthis section, we prove that the analytical solutions of (4) have positivity and bound-edness. And then we demonstrate that the numerical solutions of the given schemes(12)-(14) can preserve the properties: positivity and boundedness.4.1.
Positivity and boundedness of the analytical solutions.
Firstly we in-troduce the maximum principle which is necessary for getting the positivity andboundedness of the analytical solutions. Here we consider the following one dimen-sional space-time fractional equation Lu = ∂ α u∂t α − D ∂ β u∂ | x | β + c ( x, t ) u = f ( x, t ) , ( x, t ) ∈ Ω T = [ l, r ] × (0 , T ] , (28)where Ω T is a bounded domain with Lipschitz continuous boundary. Theorem 4.1.
Let the coefficient c ( x, t ) ≥ and the non-homogeneous term f ( x, t ) ≤ (resp. f ( x, t ) ≥ ) in Ω T . If u ∈ C , (Ω T ) T C (Ω T ) is the solution of (28), thenthe non-negative maximum (resp. non-positive minimum) of u ( x, t ) in Ω T (if exists)must reach at the parabolic boundary Γ T , i.e., max Ω T u ( x, t ) ≤ max Γ T { u ( x, t ) , } (resp . min Ω T u ( x, t ) ≥ min Γ T { u ( x, t ) , } ) . (29)In fact, if the non-negative maximum value of u ( x, t ) is not at the boundary Γ T ,then there exists an interior point ( x ∗ , t ∗ ) ∈ Ω T such that u ( x ∗ , t ∗ ) > max Γ T { u ( x, t ) , } and u ( x ∗ , t ∗ ) ≥ max Ω T u ( x, t ) . (30)We introduce the auxiliary function v ( x, t ) in the following form v ( x, t ) = u ( x, t ) − εE α, ( bt α )with ε > b >
0. Choosing sufficient small ε , v can be positive at the point( x ∗ , t ∗ ). We know that, for any ( x, t ) ∈ Ω T , v ( x, t ) satisfies ∂ α v∂t α − D ∂ β v∂ | x | β + c ( x, t ) v = f ( x, t ) − ε (cid:0) b + D C ( x ) + c ( x, t ) (cid:1) E α, ( bt α ) , where C ( x ) = (( r − x ) − β +( x − l ) − β ) sec( πβ/ − β )2Γ(2 − β ) . Note that when 1 < β <
2, sec( πβ/ − β ) >
0. And because of f ( x, t ) ≤
0, we deduce that ∂ α v∂t α − D ∂ β v∂ | x | β + c ( x, t ) v < , (31)at all interior points. While at the point ( x ∗ , t ∗ ), we already have the result [10] ∂ α v ( x ∗ , t ) ∂t α | t = t ∗ = ∂ α u ( x ∗ , t ∗ ) ∂t α − εbE α, ( b ( t ∗ ) α ) ≥ , (32)when ε ≤ (1 − α ) T − α m ∗ Γ(2 − α ) bE α, ( b ( t ∗ ) α ) . For the Riesz fractional derivative, we know that u satisfies ∂ β u∂ | x | β | x = x ∗ = lim h → h β (cid:0) − [ x ∗ /h ] X j =[( l − r + x ∗ ) /h ] j =0 g j u ( x ∗ − jh, t ∗ ) − g u ( x ∗ , t ∗ ) (cid:1) ≤ lim h → h β (cid:0) [ x ∗ /h ] X j =[( l − r + x ∗ ) /h ] j =0 g j ( u ( x ∗ , t ∗ ) − u ( x ∗ − jh, t ∗ ) (cid:1) < . Then for sufficient small ε , ∂ β v∂ | x | β | x = x ∗ = ∂ β u ( x, t ∗ ) ∂ | x | β | x = x ∗ + εD E α, ( b ( t ∗ ) α ) C ( x ) | x = x ∗ ≤ . (33)Together with (32) and (33), we obtain ∂ α v∂t α − D ∂ β v∂ | x | β + c ( x, t ) v ≥ at ( x ∗ , t ∗ ) , which is contradictory with (31).Similar analysis can be done for the case f ( x, t ) ≥
0. So from the above analysiswe arrive at Theorem 4.1.
RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 13
Next we introduce the definitions of the upper and lower solutions, from whichwe can get positivity and boundedness of the analytical solutions.
Definition 4.2.
Assume u i ( i = 1 ,
2) solve the following equations ∂ α u i ∂t α − D i ∂ β u i ∂ | x | β = f i ( u , u ) , x ∈ Ω , t ∈ (0 , T ] ,Bu i = ∂u i ∂x = g i ( x, t ) , x ∈ ∂ Ω , t ∈ (0 , T ] , (34) u i ( x,
0) = ϕ i ( x ) , x ∈ Ω , and the nonlinear functions f ( · , · ) is quasi-monotone decreasing and f ( · , · ) is quasi-monotone increasing. If there exist two functions ˜ u i ( x, t ) and u e i ( x, t ) which satisfy B ˜ u i − g i ( x, t ) ≥ ≥ Bu e i − g i ( x, t ) , x ∈ ∂ Ω , t ∈ (0 , T ] , (35)˜ u i ( x, − ϕ i ( x ) ≥ ≥ u e i ( x, − ϕ i ( x ) , x ∈ ¯Ω , (36)and ∂ α ˜ u ∂t α − D ∂ β ˜ u ∂ | x | β − f (˜ u , u e ) ≥ ≥ ∂ α u e ∂t α − D ∂ β u e ∂ | x | β − f ( u e , ˜ u ) , (37) ∂ α ˜ u ∂t α − D ∂ β ˜ u ∂ | x | β − f (˜ u , ˜ u ) ≥ ≥ ∂ α u e ∂t α − D ∂ β u e ∂ | x | β − f ( u e , u e ) . (38)Then we call U ( x, t ) = (˜ u ( x, t ) , ˜ u ( x, t )) and V ( x, t ) = ( u e ( x, t ) , u e ( x, t )) upperand lower solutions of the system (34). Theorem 4.3.
Suppose { f , f } is Lipschitz continuous and mixed quasi-monotonous.If the upper and lower solutions, U ( x, t ) and V ( x, t ) , satisfy the inequality V ( x, t ) ≤ U ( x, t ) , then (34) has a unique solution in [ V ( x, t ) , U ( x, t )] .Proof. Let us define the following iteration ∂ α ¯ u ( k )1 ∂t α − D ∂ β ¯ u ( k )1 ∂ | x | β + L · ¯ u ( k )1 = L · ¯ u ( k − + f (¯ u ( k − , u ( k − ) ,∂ α ¯ u ( k )2 ∂t α − D ∂ β ¯ u ( k )2 ∂ | x | β + L · ¯ u ( k )2 = L · ¯ u ( k − + f (¯ u ( k − , ¯ u ( k − ) ,∂ α u ( k )1 ∂t α − D ∂ β u ( k )1 ∂ | x | β + L · u ( k )1 = L · u ( k − + f (u ( k − , ¯ u ( k − ) ,∂ α u ( k )2 ∂t α − D ∂ β u ( k )2 ∂ | x | β + L · u ( k )2 = L · u ( k − + f (u ( k − , u ( k − ) ,B ¯ u ( k ) i | ∂ Ω × (0 ,T ] = B u ( k ) i | ∂ Ω × (0 ,T ] = g i ( x, t ) | ∂ Ω × (0 ,T ] , ( i = 1 , u ( k ) i ( x,
0) = u ( k ) i ( x,
0) = ϕ i ( x ) , ( x ∈ ¯Ω , i = 1 , , (39) where L is the maximum of Lipschiz constants of f and f , and denote the initialiteration function as (¯ u (0)1 , ¯ u (0)2 ) = (˜ u , ˜ u ) , (u (0)1 , u (0)2 ) = ( u e , u e ) , where u e ≤ ˜ u and u e ≤ ˜ u . According to the maximum principle Theorem 4.1 andusing the iteration (39), similar to the analysis [10], we can get u e i ≤ u (1) i ≤ · · · ≤ u ( k ) i ≤ ¯ u ( k ) i ≤ · · · ≤ ¯ u (1) i ≤ ˜ u i , ( i = 1 , . Note that f is quasi-monotone decreasing and f is quasi-monotone increasing,then lim k → + ∞ ¯ u ( k ) i = ¯ u i ( x, t ) , lim k → + ∞ u ( k ) i = u i ( x, t ) , ( i = 1 , , which satisfy ¯ u i ≥ u i . Next we will check that¯ u i = u i = u i , ( i = 1 , . Denote w = ¯ u − u ( ≥ w = ¯ u − u ( ≥ w i ≥ , i = 1 ,
2. Then weget the following inequalities ∂ α w ∂t α − D ∂ β w ∂ | x | β − L · ( w + w ) ≤ ,∂ α w ∂t α − D ∂ β w ∂ | x | β − L · ( w + w ) ≤ ,Bw = 0 , Bw = 0 ,w ( x,
0) = 0 , w ( x,
0) = 0 . (40)The below arguments show that w = w = 0 in Ω T . Suppose that there ex-ist (ˆ x i , ˆ t i ) ∈ Γ T ( i = 1 ,
2) and w and w can obtain their maximum values at(ˆ x i , ˆ t i ), ( i = 1 , w (ˆ x , ˆ t ) ≤ w (ˆ x , ˆ t ). Since w i ( x,
0) = 0 and w i ( x, t ) ≥
0, we can get t ∗ i ( < ˆ t i ) such that w i ( x, t ) < w i (ˆ x i , ˆ t i )for any ( x, t ) ∈ Ω × (0 , t ∗ i ). Suppose the function u i solve the equations ∂u i ( t ) ∂t = − w i (ˆ x i , ˆ t i ) t ∗ i in (0 , t ∗ i ) ,u i (0) = 12 w i (ˆ x i , ˆ t i ) ,u i ( t ) = 0 , t ∈ [ t ∗ i , T ] , where i = 1 ,
2, and let¯ w i ( x, t ) = w i ( x, t ) + u i ( t ) in Ω × [0 , T ] . Then ∂ α ¯ w ∂t α − ∂ β ¯ w ∂ | x | β ≤ L ( w + w ) + (1 + D C ( x )) ∂ α u ( t ) ∂t α , where C ( x ) ≥
0. So at (ˆ x , ˆ t ), ¯ w satisfies ∂ α ¯ w (ˆ x , ˆ t ) ∂t α − ∂ β ¯ w (ˆ x , ˆ t ) ∂ | x | β ≤ Lw (ˆ x , ˆ t ) + (1 + D C ( x )) ∂ α u (ˆ t ) ∂t α . RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 15
While ∂ α ¯ w (ˆ x ,t ) ∂t α | t =ˆ t ≥ − ∂ β ¯ w ( x, ˆ t ) ∂ | x | β | x =ˆ x ≥
0, and since ∂ α u ( t ) ∂t α | t =ˆ t < − w (ˆ x , ˆ t )2Γ(1 − α ) ˆ t − α , we get 2 Lw (ˆ x , ˆ t ) + ∂ α u (ˆ t ) ∂t α < , for T ≤ ( L Γ(1 − α ) ) α . This gives a contradiction, so w ≡
0. Now going back to ¯ w ,we know that ∂ α ¯ w ∂t α − ∂ ¯ w ∂x ≤ L ( w + w ) + (1 + D C ( x )) ∂ α u ( t ) ∂t α = Lw + (1 + D C ( x )) ∂ α u ( t ) ∂t α . Then w ≡
0, when T ≤ ( L Γ(1 − α ) ) α . Repeating the same process as done in [10],we have w = w in Ω T for any T . So ¯ u = u and ¯ u = u . Set u i ( x, t ) = ¯ u i = u i , i = 1 , . Then u ( x, t ) = ( u , u ) solves (34). And the uniqueness of the solution can besimilarly proved as [10].From [10], we know that the monotone properties of the nonlinear terms f ( · , · )and f ( · , · ) of (4) defined in (19) and (20), i.e., f ( · , · ) is quasi-monotone decreasing, f ( · , · ) is quasi-monotone increasing. Now we take U ( x, t ) = (˜ u ( x, t ) , ˜ u ( x, t )) = (1 , L ) , ( L > /γ )and V ( x, t ) = ( u e ( x, t ) , u e ( x, t )) = (0 , . By calculating, we know U ( x, t ) and V ( x, t ) satisfy the inequalities (37) and (38).If we specify the initial condition of (4) such that the given initial N ( x, ∈ [0 , P ( x, ∈ [0 , L ] for any x ∈ ( l, r ), then the initial conditions satisfy (36). Thuswe know that V ( x, t ) = (0 ,
0) and U ( x, t ) = (1 , L ) are lower and upper solutionsof (4)-(6); and then from Theorem 4.3, that the analytical solutions of (4)-(6) arepositive and bounded is obtained.4.2. Positivity and boundedness of the numerical solutions.
In this sub-section, we show that the numerical schemes (12) can preserve the positivity andboundedness of the corresponding analytical solutions, i.e., for any i and k , 0 1, and(1 − b ) N ki ≤ k − X n =0 ( b n − b n +1 ) N k − ni + b k N i ≤ − b N ki + 1 + b , (1 − b ) P ki ≤ k − X n =0 ( b n − b n +1 ) P k − ni + b k P i ≤ − b P ki + 1 + b L . Together with the above estimates and the numerical scheme (12), we obtain0 < N k +1 i + D µh β i X j = − M + i g j N k +1 i − j ≤ µ < − b ̺ , and 0 < P k +1 i + D µh β i X j = − M + i g j P k +1 i − j ≤ L (42)for µ < − b σδ . The inequalities above can yield 0 < N k +1 i ≤ < P k +1 i ≤ L .In fact, if 0 < N k +1 i ≤ i such that solution N k +1 i satisfies N k +1 i ≤ or N k +1 i > . Case 1. If N k +1 i ≤ 0, then we choose the minimum in { N k +1 i } M − i =1 . Denote it as N k +1 s which is non-positive. Along with the left inequality of (41), there exists N k +1 s + D µh β s − X j = − M + s +1 g j N k +1 s − j + D µg s h β (4 N k +11 − N k +12 )+ D µg − M + s h β (4 N k +1 M − − N k +1 M − ) > . Then we have (cid:18) D µg h β (cid:19) N k +1 s > D µh β s − X j = − M + s +1 j =0 − g j N k +1 s − j − D µg s h β N k +11 − D µg − M + s h β N k +1 M − + D µg s h β N k +12 + D µg − M + s h β N k +1 M − . (43)When s = 2 , M − 2, we rewrite the inequality in the following form (cid:18) D µg h β (cid:19) N k +1 s > D µh β s − X j = − M + s +1 j =0 ,s − , − M + s +2 − g j N k +1 s − j − D µg s h β N k +11 − D µg − M + s h β N k +1 M − + D µh β ( g s − g s − ) N k +12 + D µh β ( g − M + s − g − M + s +2 ) N k +1 M − . RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 17 Since − g s ≤ − g s − and − g − M + s ≤ − g − M + s +2 , then g s − g s − ≥ g − M + s − g − M + s +2 ≥ 0. In view of − g j ≥ j = 0, we get (cid:18) D µg h β (cid:19) N k +1 s > D µh β s − X j = − M + s +1 j =0 ,s − , − M + s +2 − g j N k +1 s − D µg s h β N k +1 s − D µg − M + s h β N k +1 s + D µh β ( g s − g s − ) N k +1 s + D µh β ( g − M + s − g − M + s +2 ) N k +1 s . This implies D µh β s X j = − M + s g j N k +1 s > . Combining with the properties of g i , we know 1+ D µh β s P j = − M + s g j > 0. So N k +1 s > s = 2, then along with the inequality (43) we have (cid:18) D µg h β (cid:19) N k +12 − D µg h β N k +12 > D µh β X j = − M +1 j =0 − g j N k +12 − j − D µg h β N k +11 − D µg − M +2 h β N k +1 M − + D µg − M +2 h β N k +1 M − . After calculating, we get D µh β X j = − M +2 − g j N k +12 > . Repeating the arguments above, we know that there exists a contradiction. If s = M − 2, the similar argument works.Case 2. If N k +1 i > 1, then there exists a maximum N k +1 s in { N k +1 i } M − i =1 . Sincethe right inequality of (41) implies N k +1 s + D µh β s X j = − M + s g j N k +1 s − j ≤ , we have (cid:18) D µh β (cid:19) N k +1 s ≤ D µh β s X j = − M + sj =0 − g j N k +1 s . Rewriting the inequality, it follows that D µh β s X j = − M + s g j N k +1 s ≤ . Noting that 1 + D µh β s P j = − M + s g j > N k +1 s needs to satisfy N n +1 s < 1, which iscontradictory with the assumption.The argument for P k +1 i is similar to the discussion of N k +1 i .5. Numerical experiments. To verify the above theoretical results, we presentsome numerical results of the two numerical schemes (fractional centered differencescheme and WSGD scheme) with Neumann boundary. And the initial conditionsare taken as [8] N ( x, 0) = 0 . . cos ( πx ) ,P ( x, 0) = 0 . . cos ( πx ) . In this section, we fix the parameters σ = 1 , ̺ = 1 . , γ = 0 . , κ = 1, and δ = 0 . 5, use the diffusion coefficients D = 0 . 005 and D = 0 . , × (0 , h = τ = 0 . 01. We observe that the orders α and β affect the shape of the solutionsobviously. Fig. 1 shows the solutions of classical predator-prey reaction-diffusionmodel with α = 1 and β = 2. When α tends to 1 and β to 2, the numerical solutionsof the fractional predator-prey reaction-diffusion equations are also convergent tothe solutions of the classical ones. N P Figure 1. Numerical solution for the case α = 1 and β = 2To get the convergent order in time and space, we consider the nonhomogeneousequations: ∂ α N∂t α = D ∂ β N∂ | x | β + N (cid:18) − N − ̺PP + N (cid:19) + f ( x, t ) ,∂ α P∂t α = D ∂ β P∂ | x | β + σP (cid:18) − γ + κδP κP + NP + N (cid:19) + g ( x, t ) . RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 19 N P Figure 2. Numerical solution for the case α = 1 and β = 1 . N P Figure 3. Numerical solution for the case α = 1 and β = 1 . N ( x, t ) = P ( x, t ) = ( t + 1) x (1 − x ) , N P Figure 4. Numerical solution for the case α = 0 . β = 2 N P Figure 5. Numerical solution for the case α = 0 . β = 1 . f ( x, t ) =2( − x ) x ( t − α Γ(2 − α ) + t − α Γ(3 − α ) )+ D (1 + t ) Γ(5 − β ) x − β ( ( − x ) x β (1 − x ) β (12 x − xβ + ( − β ) β )+ x (12( − x ) + ( − x ) β + β )) sec( πβ − ( t + 1) x (1 − x ) (1 − ( t + 1) x (1 − x ) − ̺/ , RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 21 g ( x, t ) =2( − x ) x ( t − α Γ(2 − α ) + t − α Γ(3 − α ) )+ D (1 + t ) Γ(5 − β ) x − β ( ( − x ) x β (1 − x ) β (12 x − xβ + ( − β ) β )+ x (12( − x ) + ( − x ) β + β )) sec( πβ − σ ( t + 1) x (1 − x ) ( − γ + κδ ( t + 1) x (1 − x ) κ ( t + 1) x (1 − x ) + 1 / . Table 1 shows that the method has first-order accuracy for the time discretiza-tions when α = 0 . β = 1 . 5. In the computations of the example, we takethe spacial steplength h = 0 . τ = 0 . Table 1. The numerical errors and convergent orders for the timediscretization with α = 0 . β = 1 . τ ( h = 0 . e N ( h, τ ) rate e P ( h, τ ) rate0 . . . . 05 0 . . . . . 025 8 . e − 04 0 . . . . . e − 04 0 . . e − 04 0 . Table 2. The numerical errors and convergent orders of the frac-tional centered difference scheme with α = 0 . β = 1 . .h ( τ = 0 . e N ( h, τ ) order e P ( h, τ ) order0 . . . . 05 4 . e − 04 2 . . . . 025 4 . e − 05 3 . . e − 04 3 . . . e − 06 2 . . e − 05 2 . Conclusion. More than three decades ago, the classical diffusion is introducedinto the predator-prey model, which leads to the predator-prey reaction-diffusionmodel. With the deep research on the anomalous diffusion, it is noticed that theanomalous diffusion is ubiquitous. It seems nature to introduce the anomalous Table 3. The numerical errors and convergent orders of the frac-tional centered difference scheme with α = 0 . β = 1 . .h ( τ = 0 . e N ( h, τ ) order e P ( h, τ ) order0 . . . . 05 1 . e − 04 3 . . e − 04 3 . . 025 1 . e − 05 2 . . e − 05 2 . Table 4. The numerical errors and convergent orders of the frac-tional centered difference scheme with α = 0 . β = 1 . .h ( τ = 0 . e N ( h, τ ) order e P ( h, τ ) order0 . . . . 05 0 . . . . . 025 4 . e − 04 2 . . . . . e − 05 2 . . e − 04 2 . Table 5. The numerical errors and convergent orders of theWSGD scheme with α = 0 . β = 1 . .h ( τ = 0 . e N ( h, τ ) order e P ( h, τ ) order0 . . e − 04 0 . . 05 4 . e − 04 0 . . . . 025 1 . e − 04 1 . . e − 04 1 . . . e − 05 2 . . e − 04 1 . Table 6. The numerical errors and convergent orders of theWSGD scheme with α = 0 . β = 1 . h ( τ = 0 . e N ( h, τ ) order e P ( h, τ ) order0 . . . . 05 0 . . . . . 025 0 . . 222 0 . . . . . . . Acknowledgements. This work was supported by the National Natural ScienceFoundation of China under Grant No. 11271173. RACTIONAL PREDATOR-PREY REACTION-DIFFUSION MODEL 23 REFERENCES [1] F. Brauer and C. Castillo-Chavez, “Mathematical Models in Population Biology and Epi-demiology”, Springer-Verlag, 2000.[2] E. D. Conway and J. A. Smoller, Diffusion and the predator-prey interaction , SIAM J. Appl.Math., (1977), 673-686.[3] M. Cavani and M. Farkas, Bifurcation in a predator-prey model with memory and diffusionII: Turing bifurcation , Acta Math. Hungar., (1994), 375-393.[4] S. Aly, I. Kim and D. Sheen, Turing instability for a ratio-dependent predator-prey modelwith diffusion , Appl. Math. Comp., (2011), 7265-7281.[5] F. Bartumeus, D. Alonso, and J. Catalan, Self organized spatial structures in a ratio dependentpredator-prey model , Phys. A, (2001), 53-57.[6] P. Y. H. Pang and M. X. Wang, Qualitative analysis of a ratio-dependent predator prey systemwith diffusion , Proc. R. Soc. Edinb., (2003), 919-942.[7] M. Wang, Stationary patterns for a prey predator model with prey-dependent and ratio-dependent functional responses and diffusion , Physica D, (2004), 172-192.[8] A. Yun, D. Jeong, and J. Kim, An Efficient and Accurate Numerical Scheme for TuringInstability on a Predator-Prey Model , Int. J. Bifur. Chaos, (2012), 1250139.[9] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractionaldynamics approach , Phys. Rep., (2000), 1-77.[10] Y. Y. Yu, W. H. Deng, and Y. J. Wu, Positivity and boundedness preserving schemes for thefractional reaction-diffusion equation , Sci. China Math., (2013), 2161-2178.[11] M. M. Meerschaert, H.P. Scheffler, and C. Tadjeran, Finite difference methods for two-dimensional fractional dispersion equation , J. Comput. Phys., (2006), 249-261.[12] Q. Yang, F. Liu, and I. Turner, Numerical methods for fractional partial differential equationswith Riesz space fractional derivatives , Appl. Math. Model., (2010), 200-218.[13] C. Chen, F. Liu, I. Turner, and V. Anh, A Fourier method for the fractional diffusion equationdescribing sub-diffusion , J. Comput. Phys., (2007), 886-897.[14] M. R. Cui, Compact finite difference method for the fractional diffusion equation , J. Comput.Phys., (2009), 7792-7804.[15] G. Gao and Z. Sun, A compact difference scheme for the fractional sub-diffusion equations ,J. Comput. Phys., (2011), 586-595.[16] Y. Lin and C. Xu, Finite difference/spectral approximations for the time-fractional diffusionequation , J. Comput. Phys., (2007), 1533-1552.[17] S. B. Yuste and L. Acedo, An explicit finite difference method and a new Von-NeumannType stability analysis for fractional diffusion equations , SIAM J. Numer. Anal. (2005),1862-1874.[18] W. H. Deng, Numerical algorithm for the time fractional Fokker-Planck equation , J. Comput.Phys., (2007), 1510-1522.[19] M. D. Ortigueira, Riesz potential operators and inverses via fractianal centred derivatives ,International J. Math. Math. Sci., (2006), 48391.[20] C. Celik and M. Duman, Crank-Nicolson method for the fractional diffusion equation withthe Riesz fractional derivative , J. Comput. Phys., (2012), 1743-1750.[21] D. Wang, A. Xiao, and W. Yang, Crank-Nicolson difference scheme for the coupled nonlinearSchr¨odinger equations with the Riesz space fractional derivative , J. Comput. Phys., (2013), 670-681.[22] W. Y. Tian, H. Zhou, and W. H. Deng, A class of second order difference approximation forsolving space fractional diffusion equations , arXiv:1201.5949 [math.NA] .[23] H. Zhou, W. Y. Tian, and W. H. Deng, Quasi-compact finite difference schemes for spacefractional diffusion equations , J. Sci. Comput., (2013), 45-66. E-mail address : [email protected] E-mail address : [email protected] E-mail address ::