Numerical solution of fractional order diffusion problems with Neumann boundary conditions
aa r X i v : . [ m a t h . NA ] N ov Numerical solution of fractional order diffusionproblems with Neumann boundary conditions
Béla J. Szekeres and Ferenc Izsák , ∗ September 11, 2018 Department of Applied Analysis and Computational Mathematics,Eötvös Loránd University H-1117, Budapest, Pázmány P. s. 1/C, Hungary.
Abstract
A finite difference numerical method is investigated for fractional order diffu-sion problems in one space dimension. For this, a mathematical model is devel-oped to incorporate homogeneous Dirichlet and Neumann type boundary conditions.The models are based on an appropriate extension of the initial values. The well-posedness of the obtained initial value problems is proved and it is pointed out thatthe extensions are compatible with the above boundary conditions. Accordingly, afinite difference scheme is constructed for the Neumann problem using the shiftedGrünwald–Letnikov approximation of the fractional order derivatives, which is basedon infinite many basis points. The corresponding matrix is expressed in a closed formand the convergence of an appropriate implicit Euler scheme is proved.
Keywords: fractional order diffusion; Grünwald–Letnikov formula; non-local deriva-tive; Neumann boundary conditions; implicit Euler scheme.
A widely accepted constitutive relation, the first Fick’s law leads to the standard diffu-sion model. At the same time, more observations confirmed the presence of super- and ∗ Corresponding author. E-mail address: [email protected] , tel.: +36 13722500-8428, fax: +3613812158. u is observed on a physical volume then we have information onthe density only on the closure of this. On the other hand, in the mathematical model,nonlocal operators are used, which require the value of u also outside of the domain. Foran accurate finite difference approximation we also need values outside of the domain. Inthis way, it is natural to look for an extension of the density u .The majority of the authors consider homogeneous Dirichlet boundary conditions andassume zero values outside of the domain. A similar approach is applied in the non-localcalculus [5].Regarding Neumann type boundary conditions, at our best knowledge, there is onlyone attempt [12], dealing the boundary condition at the operator level and proposing amatrix transformation technique.The aim of this paper is to • develop a meaningful mathematical approach to model homogeneous Neumann (andDirichlet) boundary conditions • analyze the well-posedness of the corresponding PDE’s • construct a corresponding finite difference scheme with a full error analysis. Fractional calculus
We summarize some basic notions and properties of the fractionalcalculus. For more details and examples we refer to the monographs [13],[15],[18].2o define an appropriate function space for the fractional order derivatives on the realaxis we introduce for arbitrary a, b ∈ R the function spaces ¯ C ( a, b ) = C [ a, b ] | ( a,b ) and ¯ C ( a, b ) / R = { f ∈ ¯ C ( a, b ) : Z ba f = 0 } . With these we define C I ( R ) = { b − a - periodic extension of f : f ∈ ¯ C ( a, b ) / R } . Remarks:
1. Functions in C I ( R ) are bounded and they are continuous except of the possible dis-continuity points { a + k ( b − a ) : k ∈ Z } .2. In the points { a + k ( b − a ) : k ∈ Z } we define f to be f ( a )+ f ( b )2 . (cid:3) Definition 1
For the exponent β ∈ (0 , the fractional order integral operators −∞ I βx and x I β ∞ on the space C I ( R ) are defined with −∞ I βx f ( x ) = 1Γ( β ) Z x −∞ f ( s )( x − s ) − β d s and x I β ∞ f ( x ) = 1Γ( β ) Z ∞ x f ( s )( s − x ) − β d s. With this, the left and right-sided
Riemann–Liouville derivatives of order α ∈ R + \ Z aregiven by RL −∞ ∂ αx f ( x ) = ∂ nx −∞ I n − αx f ( x ) and RLx ∂ α ∞ f ( x ) = ( − n ∂ nx x I n − α ∞ f ( x ) , where n is the integer with n − < α < n .Accordingly, for a function f = f + C f , where f ∈ C I ( R ) and C f is a constant functionwe define the Riesz derivative with ∂ α | x | f ( x ) = C σ ( RL −∞ ∂ αx f ( x ) + RLx ∂ α ∞ f ( x )) , where C σ = − σ α π with a given positive constant σ . emarks:
1. For the simplicity, the constant σ – which corresponds to the intensity of the superdif-fusive process in the real-life phenomena – does not appear in the notation ∂ α | x | .2. In the original definition, the Riemann–Liouville derivatives are given on a boundedinterval ( a, b ) ⊂ R such that in the above definition −∞ and ∞ are substituted with a and b , respectively. In this case a I βx and x I βb can be defined for the exponent β ∈ R + oreven for β ∈ C with Re β > , in case of complex valued functions, see Section 2.1 in[13]. Moreover, the definition can be extended to be a bounded operator on the functionspace L ( a, b ) , see Theorem 2.6 in [18]. For an overview of the alternating notations anddefinitions, we refer to the review paper [10].3. An advantage of the above approach is that alternative definitions on real axis coincide.For instance, fractional derivatives can be interpreted using the fractional power of thenegative Laplacian − ∆ , see Lemma 1 in [26]. This can be introduced via Fourier trans-form, see Section 2.6 in [18]. For more information and multidimensional extension of theRiesz derivative see also Section 2.10 in [13]. Note that finite difference discretizations forRiesz fractional derivatives has been studied also in [19] highlighting its connection withprobabilistic models.4. We have left open the question for which functions does Definition 1 make sense. Thegeneral answer requires the discussion of the Triebel–Lizorkin spaces [1], [18], which isbeyond the scope of this paper. Some sufficient conditions on a bounded interval ( a, b ) are also discussed in Lemma 2.2 in [13]. At the same time, since we will approximatethe Riesz derivative with finite differences of the fractional integrals and verify that thefractional integrals make sense on the function space C I ( R ) . Lemma 1
For each function f ∈ C I ( R ) and any exponent α ∈ (1 , the fractional orderintegral operators −∞ I − αx f and x I − α ∞ f make sense.Proof: We prove the statement for the right-sided approximation, the left sided can behandled in a similar way. In concrete terms, we prove that Z ∞ x f ( s )( s − x ) α − d s < ∞ . (1)Obviously, there is a k ∈ Z such that a + ( b − a ) k > x and accordingly, Z ∞ x f ( s )( s − x ) α − d s = Z a +( b − a ) kx f ( s )( s − x ) α − d s + Z ∞ a +( b − a ) k f ( s )( s − x ) α − d s. (2)4ere, using the condition < α − < we have that the first term is finite. To estimatethe second one, we introduce the function F : ( a + ( b − a ) k, ∞ ) → R with F ( s ) = Z sa +( b − a ) k f ( s ∗ ) d s ∗ such that F ′ ( s ) = f ( s ) on ( a + ( b − a ) k, ∞ ) . Also, since f ∈ C I ( R ) , we have that F ( a + ( b − a ) k ) = F ( a + ( b − a )( k + 1)) = . . . such that F is bounded.With this the second term on the right hand side of (2) can be rewritten as Z ∞ a +( b − a ) k F ′ ( s )( s − x ) α − d s = (cid:20) F ( s )( s − x ) α − (cid:21) ∞ a +( b − a ) k − Z ∞ a +( b − a ) k F ( s ) · (1 − α )( s − x ) α d s = − Z ∞ a +( b − a ) k F ( s ) · (1 − α )( s − x ) α d s, which is also finite since F is bounded and − β > . Therefore, (1) is finite, which provesthe lemma. (cid:3) Fundamental solutions
For the forthcoming analysis we analyze the Cauchy problem ( ∂ t u ( t, x ) = ∂ α | x | u ( t, x ) t ∈ R + , x ∈ R u (0 , x ) = u ( x ) x ∈ R (3)with a given initial function u ∈ C I ( R ) and α ∈ (1 , . Lemma 2
The Cauchy problem in (3) has a unique solution and can be given by u ( t, x ) = (Φ α,t ∗ u )( x ) , (4) where Φ α,t denotes the fundamental solution corresponding to the Riesz fractional differ-ential operator ∂ α | x | , furthermore, u ( t, · ) ∈ C ∞ ( R ) for all t ∈ R + .Proof: Applying the spatial Fourier transform F to the equations in (3) we have ( ∂ t F u ( t, s ) = −| s | α F u ( t, s ) t ∈ R + , s ∈ R F u (0 , s ) = F u ( s ) s ∈ R , F h ∂ α | x | u ( t, x ) i ( s ) = −| s | α F u ( t, s ) , see [10], p. 38. There-fore, F u ( t, s ) = e − t | s | α F u ( s ) such that an inverse Fourier transform F − implies that u ( t, x ) = F − ( e − t | s | α F u ( s ))( x ) = F − ( e − t | s | α ) ∗ u ( x ) . In this way, the fundamental solution of (3) can be given as Φ α,t ( s ) = F − ( e − t | x | α )( s ) = F ( e − t | x | α )( s ) . Using the fact that the function to transform is even and the Fourier transform F exp {− a |·|} ( s ) can be given (see, e.g. , [9], p. 1111) we have Φ α,t ( s ) = F ( e − t | x | α )( s ) = F ( e − t | x | e − t | x | α − )( s ) = F ( e − t | x | )( s ) ∗ F ( e − t | x | α − )( s )= q π tt + s ∗ F ( e − t | x | α − )( s ) . (5)Here for any fixed t the real function given by s → √ π tt + s is in C ∞ ( R ) and also all ofits derivatives are in L ( R ) . On the other hand, for α > the real function given with e − t | x | α − is in L ( R ) , therefore F ( e − t | x | α − ) is bounded and continuous. Consequently,using (5) the right hand side of the equality ∂ k Φ α,t ( s ) = ∂ k q π tt + s ∗ F ( e − t | x | α − )( s ) makes sense, which gives statement in the lemma. (cid:3) Discretization
The finite difference approximation of the fractional order derivatives isnot straightforward. It turns out that an obvious one-sided finite difference approximationof the one-sided Riemann–Liouville derivatives results in an unstable method even if animplicit Euler method is applied for the time marching scheme [14]. To stabilize theseschemes, we need to use the translated Grünwald–Letnikov formula , which is given for f ∈ C ( R ) with D α,p,h −∞ ,GL f ( x ) = 1Γ( − α ) lim M →∞ h α M X k =0 Γ( k − α )Γ( k + 1) f ( x − ( k − p ) h )= 1 h α ∞ X k =0 g k f ( x − ( k − p ) h ) (6)6 * * * * * * * ♦ ♦ ♦ ♦ ♦ ♦ ♦ ♦ · · · x − h x − h x − h x − h x x + h x + 2 h x + 3 h x + 4 h · · · Figure 1: Basis points for the the left-sided ( ∗ ) and the right-sided ( ♦ ) translatedGrünwald–Letnikov formula (given in (6) and (7)) applied in x with the translation pa-rameter p = 2 .and D α,p,h ∞ ,GL f ( x ) = 1Γ( − α ) lim M →∞ h α M X k =0 Γ( k − α )Γ( k + 1) f ( x + ( k − p ) h )= 1 h α ∞ X k =0 g k f ( x + ( k − p ) h ) (7)depending on the translation parameter p ∈ N , the order of the differentiation α ∈ (1 , and the discretization parameter h ∈ R + , where we used the coefficients g k = Γ( k − α )Γ( − α )Γ( k + 1) = ( − k (cid:18) αk (cid:19) . The principle of the two-sided translated discretizations is depicted in Figure 1. Thesecoefficients satisfy the following: ∞ X k =0 g k = 0 ∀ α ∈ (1 , g = − α, g j ≥ for j = 1 . (8)We use the same notation for the discrete differential (or difference) operators, i.e. foreach v = ( . . . , v − , v , v , . . . ) ∈ R Z we write h D α,p,h −∞ ,GL v i j = 1 h α ∞ X k =0 ( − k g k v j + p − k and h D α,p,h ∞ ,GL v i j = 1 h α ∞ X k =0 ( − k g k v j + k − p , (9)where the superscript j denotes the j th component. Remarks:
7. One can prove [14] that the integrals in (6) and (7) approximate the correspondingRiemann–Liouville derivatives in the following sense: RL −∞ ∂ αx f ( x ) = D α,p,h −∞ ,GL f ( x ) + O ( h ) (10)and similarly, RLx ∂ α ∞ f ( x ) = D α,p,h ∞ ,GL f ( x ) + O ( h ) , (11)provided that both of the Fourier transform of f and that of RL −∞ ∂ αx f ( x ) and RLx ∂ α ∞ f ( x ) are in L ( R ) . We will point out that in the present framework no such assumption isnecessary.2. Higher-order finite difference approximations can be obtained as a linear combinationof first order ones using different translation parameters. For example, the sum definedby GL −∞ D α,p,qx,h u ( x ) = 2 q − α p − q ) D α,p,h −∞ ,GL u ( x ) + 2 p − α p − q ) D α,p,h −∞ ,GL u ( x ) (12)provides a second order accurate approximation of the Riemann–Liouville derivative. Sim-ilar statement holds if −∞ is switched to ∞ . For the details, see [27].3. The nonlocal effect of the differential operators result in full matrices. At the sametime, one can save some computing efforts with an appropriate decomposition of the abovematrix [25].4. An alternative approximation of fractional order elliptic operators (which can be ap-plied in multidimensional cases) was proposed in [3], which can be a basis also for finiteelement discretizations. Extensions are not only necessary to have well-posed problems involving nonlocal diffusionoperators, but also essential at the discrete level. In order to have sufficient accuracy in thefinite difference approximation near to the boundary and at the boundary of a nonlocaldifferential operator, it is necessary to have (virtual) gridpoints outside of the originalcomputational domain
Ω = ( a, b ) . This is clearly shown in (6), (7) and (12).Summarized, the sketch of our approach is the following: • we extend the problem to R to get rid of the boundary conditions • we solve the corresponding Cauchy problem (3) (we will approximate this with finitedifferences) 8 we verify that the desired homogeneous Neumann (or homogeneous Dirichlet) bound-ary conditions are satisfied for the restriction of the solution. Definition 2
We say that the extension ˜ · : L (Ω) → C I ( R ) is compatible with the homogeneous Neumann (no-flux) or Dirichlet boundary conditionsand the operator ∂ α | x | if the function ˜ u is the unique solution of the following problem ( ∂ t ˜ u ( t, x ) = ∂ α | x | ˜ u ( t, x ) t ∈ R + , x ∈ R ˜ u (0 , x ) = ˜ u ( x ) x ∈ R and ∂ x ˜ u ( t, x ) = 0 or ˜ u ( t, x ) = 0 for x ∈ ∂ Ω and t > . In rough terms, one can say that a correct extension of the solution from Ω is the functionwhich solves the extended problem on the whole real axis such that the boundary conditionon ∂ Ω is still satisfied. Extension corresponding to homogeneous Dirichlet boundary condition
As amotivation, we use the idea in [20] which is generalized to the case of two absorbing walls.For the simplicity, the definition is given for functions u : (0 , → R . Definition 3
We call the 2-periodic extension of the function ˆ u D , ( x ) = ( u ( x ) for x ∈ (0 , − u ( − x ) for x ∈ ( − , the Dirichlet type extension of u and we denote it with ˆ u D .Remarks:
1. This is sometimes called the odd extension of u and can be obtained first with reflectingthe graph of u to (0 , and (1 , ∈ R to extend it to ( − , and reflecting this graphfurther to ( − , and (2 , ∈ R to extend it to ( − , and continuing this process.2. A natural physical interpretation of this extension is the following: To ensure zero-concentration at the end points in the consecutive time steps, we have to force a skew-symmetric concentration profile around and .3. We may define u D ( k ) = 0 for k ∈ Z but as we will see this is not essential.9 simple calculation shows that we have ˆ u D ( y ) = − ˆ u D ( − y ) and ˆ u D (1 + y ) = − ˆ u D (1 − y ) y ∈ R . (13)We define similarly the extension ˆ v D ∈ R Z of the vector v = ( v , v , . . . , v N +1 ) ∈ R N +2 with v = v N +1 = 0 : [ˆ v D ] j = v j j = 0 , , . . . , N + 1 − v N +1) − j j = N + 2 , N + 3 , . . . , N + 1 v j +2 N j
6∈ { , , . . . , N + 1 } . (14) Extension corresponding to homogeneous Neumann boundary condition
Sim-ilarly to the previous case, the definition is given for functions u : (0 , → R . Definition 4
We call the 2-periodic extension of the function ˆ u N , ( x ) = ( u ( x ) for x ∈ (0 , u ( − x ) for x ∈ ( − , the Neumann type extension of u and we denote it with ˆ u N .Remarks:
1. This is sometimes called the even extension of u and can be obtained first with reflectingthe graph of u to the vertical line given with x = 0 to extend it to ( − , then to thevertical line given with x = 1 to extend it to (1 , and continuing this process.2. A natural physical interpretation of this extension is the following: To ensure zero-flux,we force a symmetric concentration profile around and .A simple calculation shows that ∂ x ˆ f N (1 + y ) = ∂ x ˆ f N (1 − y ) and ∂ x ˆ f N ( y ) = ∂ x ˆ f N ( − y ) y ∈ R . (15)Similar notations are used for the “extended” vector ˆ v N ∈ R Z of v = ( v , v , . . . , v N ) ∈ R N +1 which is defined as follows: [ˆ v N ] j = v j j = 0 , , . . . , Nv − j − j = − N − , − N, . . . , − v j +2 N +2 j
6∈ {− N − , . . . , , . . . , N } . (16)A natural physical interpretation of this extension is that particles are reflected at theboundaries to ensure zero flux. In this way, we also reflect the concentration profile inthe model. The principle of the Neumann extension for a vector is visualized in Figure 2.We verify that the above extensions meet the requirements in Definition 2.10 ⋄ ⋄ ⋄ • • • • ⋄ ⋄ ⋄ ⋄ v − v − v − v − v v v v v v v v Figure 2: Entries of the Neumann extension ˆ v N of the vector v = ( v , v , v , v ) . Lemma 3
The extensions ˆ u N and ˆ u D are compatible with the Dirichlet and the Neumannboundary conditions, respectively and with the differential operator ∂ α | x | .Proof: According to Lemma (2) we can express the solution of ( ∂ t u = ∂ α | x | u on R + × R u (0 , · ) = ˆ u N with the convolution u ( t, · ) = Φ α,t ∗ ˆ u N = ˆ u N ∗ Φ α,t , such that u ( t, x ) = Z R ˆ u N ( x − y )Φ α,t ( y ) d y . Since Φ α,t ∈ C ∞ ( R ) , the same holds for u ( t, · ) . Accordingly, the right and left limit of ∂ x u ( t, · ) in coincide. Using this, (15) and the fact that Φ α,t is even we obtain ∂ x u ( t,
1) = lim ǫ n → − ∂ x u ( t, − ǫ n ) = lim ǫ n → − ( ∂ x ˆ u N ∗ Φ α,t )(1 − ǫ n )= lim ǫ n → − Z R ∂ x ˆ u N (1 − ǫ n − y )Φ α,t ( y ) d y = − lim ǫ n → − Z R ∂ x ˆ u N (1 + ǫ n + y )Φ α,t ( y ) d y = − lim ǫ n → − Z R ∂ x ˆ u N (1 + ǫ n + y )Φ α,t ( − y ) d y = − lim ǫ n → − Z R ∂ x ˆ u N (1 + ǫ n − y )Φ α,t ( y ) d y = − lim ǫ n → − ( ∂ x ˆ u N ∗ Φ α,t )(1 + ǫ n )= − lim ǫ n → − ∂ x u ( t, ǫ n ) = − ∂ x u ( t, , which gives that the homogeneous Neumann boundary condition is satisfied in . With anobvious modification, using (15), we can also verify the homogeneous Neumann boundarycondition ∂ x u ( t,
0) = 0 . 11imilarly, we can express the solution of ( ∂ t u = ∂ α | x | u on R + × R u (0 , · ) = ˆ u D with u ( t, x ) = Z R ˆ u D ( x − y )Φ α,t ( y ) d y . Since Φ α,t ∈ C ∞ ( R ) , the same holds for u ( t, · ) . Accordingly, the right and left limit of u ( t, · ) in coincide. Using this and (13) we obtain u ( t,
1) = lim ǫ n → − u ( t, − ǫ n ) = lim ǫ n → − ˆ u D ∗ Φ α,t (1 − ǫ n )= lim ǫ n → − Z R ˆ u D (1 − ǫ n − y )Φ α,t ( y ) d y = − lim ǫ n → − Z R ˆ u D (1 + ǫ n + y )Φ α,t ( y ) d y = − lim ǫ n → − Z R ˆ u D (1 + ǫ n + y )Φ α,t ( − y ) d y = − lim ǫ n → − Z R ˆ u D (1 + ǫ n − y )Φ α,t ( y ) d y = − lim ǫ n → − ˆ u D ∗ Φ α,t (1 + ǫ n ) = − lim ǫ n → − u ( t, ǫ n ) = − u ( t, , which gives that the homogeneous Dirichlet boundary condition is satisfied in . With anobvious modification, using (13), we can verify homogeneous Dirichlet boundary conditionalso in .In both cases, the equality ∂ t u = ∂ α | x | u is satisfied in (0 , which gives the statementin the lemma. (cid:3) Using the above extension, we will investigate the numerical solution of the problem ( ∂ t u ( t, x ) = ∂ α | x | u N ( t, · )( x ) t ∈ (0 , T ) , x ∈ (0 , u (0 , x ) = u x ∈ (0 , . (17)The following theorem is the basis of our numerical method, which again confirms thefavor of the Neumann type extension. Theorem 1
For any u ∈ C [0 , the problem in (17) is well-posed, and its unique solu-tion u ∈ C ∞ [0 , satisfies the boundary conditions ∂ x u ( t,
0) = ∂ x u ( t,
1) = 0 .Proof:
We first note that the problem ( ∂ t u ( t, x ) = ∂ α | x | u ( t, x ) t ∈ (0 , T ) , x ∈ R u (0 , x ) = ˆ u N ( x ) x ∈ R . (18)12s well-posed and corresponding to Lemma 3 its solution can be given by u ( t, x ) = Z R ˆ u N ( x − y )Φ α,t ( y ) d y . This implies that u ( t, − x ) = Z R ˆ u N ( − x − y )Φ α,t ( y ) d y = Z R ˆ u N ( x + y )Φ α,t ( − y ) d y = Z R ˆ u N ( x − y )Φ α,t ( y ) d y = u ( t, x ) . and u ( t, x + 2) = Z R ˆ u N ( x + 2 − y )Φ α,t ( y ) d y = Z R ˆ u N ( x − y )Φ α,t ( y ) d y = u ( t, x ) . Therefore, the restriction of the solution u of (18) solves the problem in (17). Here wehave also used throughout that for the Neumann extension: u N ∈ C I ( R ) such that theRiesz derivative ∂ α | x | u N makes sense.To prove the uniqueness, we assume that v solves (17). According to (15) we obviouslyhave that the Neumann extension satisfies ∂ t ˆ v N ( t, − x ) = ∂ t ˆ v N ( t, x ) x ∈ R . (19)To compute the fractional order integral we introduce the function C I ( R ) ∋ ˆ v N =ˆ v ( t, · ) N − R v ( t, · ) such that we have Γ(2 − α ) ∞ I − αx ˆ v N ( − x ) = Z − x −∞ v N ( y ∗ )( − x − y ∗ ) α d y ∗ + Z ∞− x v N ( y ∗ )( y ∗ + x ) α d y ∗ = Z ∞ x v N ( − y )( − x + y ) α d y + Z x −∞ v N ( − y )( − y + x ) α d y = Z ∞ x v N ( y )( − x + y ) α d y + Z x −∞ v N ( y )( − y + x ) α d y = Γ(2 − α ) ∞ I − αx ˆ v N ( x ) . Therefore, we also have ∂ α | x | v N ( t, − x ) = ∂ α | x | v N ( t, x ) . (20)Consequently, (19) and (20) imply ∂ t ˆ v N ( t, − x ) = ∂ t ˆ v N ( t, x ) = ∂ α | x | v N ( t, x ) = ∂ α | x | v N ( t, − x ) and the periodicity obviously gives ∂ t ˆ v N ( t, x ) = ∂ t ˆ v N ( t, x ) = ∂ α | x | v N ( t, x ) = ∂ α | x | v N ( t, x ) such that the Neumann extension u N solves (18). This, however, has a unique solution,which gives the uniqueness of the solution of (17). (cid:3) .2 Numerical methods Following the classical method of lines technique we first discretize the spatial variablesin the extended problem and choose a time stepping scheme for the full discretization.To streamline the forthcoming computations, the interval [0 , will be transformed to [0 , π ] , where we use the following grid points: x j := πh (cid:18) j + 12 (cid:19) = π (cid:18) N + 1) + jN + 1 (cid:19) . (21) u nj denotes the numerical approximation at time nτ in the grid point x j and u ( nδ, · ) the values of the analytic solution at time nδ in the grid points. For the spatial discretization we use the Grünwald–Letnikov approximations in (9) intro-ducing A α,h ∈ R ( N +1) × ( N +1) with A α,h u = D α, ,h −∞ , GL u N + D α, ,h ∞ , GL u N . (22)This is combined with an implicit Euler time stepping to obtain u n +1 j − u nj τ = (cid:2) A α,h u n +1 (cid:3) j . (23)To make the consecutive formulas more accessible, we expand (22) in a concrete example. Example
We give the first component of A α,h v for v = ( v , v , v , v ) . [ A α,h v ] = g v + g v + g v + g v + g v + g v + g v + g v + g v + . . . + g v + g v + g v + g v + g v + g v + g v + g v + g v + . . . . In the general case, using (22), (9) and (16) we obtain that [ A α,h v ] j = C σ h α j +1 X k =0 g k v j − k +1 + N − j +1 X k =0 g k v k + j − + ∞ X l =0 j + N +2 X k = j +2 g k +2 l ( N +1) v k − j − + N − j +2 X k = N − j +2 g k +2 l ( N +1) v N − k − j +2 ! + ∞ X l =0 j +2 N +3 X k = j + N +3 g k +2 l ( N +1) v N + j − k +3 + N − j +3 X k =2 N − j +3 g k +2 l ( N +1) v k − N + j − !! j = 1 , , . . . , N − (24)14 A v ] = 2 C σ h α ( g v + g v + ∞ X l =0 N +2 X k =2 g k +2 l ( N +1) v N +2 − k + N +3 X k = N +3 g k +2 l ( N +1) v N − k +3 !! , (25) [ A v ] N = 2 C σ h α ( g v N − + g v N + ∞ X l =0 N +2 X k =2 g k +2 l ( N +1) v N − k +2 + N +3 X k = N +3 g k +2 l ( N +1) v k − N − !! . (26)For K = 2 m +1 the corresponding matrix can be given with a slight modification. Observethat all of the coefficients g i arises once on the right hand side of (24), (25) and (26). Proposition 1
The matrix A α,h has negative diagonal and non-negative off-diagonal el-ements.Proof: Observe that in (24), (25) and (26) the coefficient of v j , v and v N , respectivelyis g and g appears only here. Therefore, using also (8) we obtain that A has negativediagonals and positive off-diagonals. (cid:3) We analyze the properties of the matrix A α,h and the corresponding differential oper-ator. Lemma 4
The eigenvectors of the matrix A α,h ∈ R ( N +1) × ( N +1) are given as v kh = (cos kx , cos kx , . . . , cos kx N − , cos kx N ) T for each k = 0 , , . . . , N with the corresponding eigenvalues − σ cos (cid:0) α π (cid:1) (cid:18) h (cid:19) α sin α kπh · cos (cid:16) kπh + α π − kπh ) (cid:17) . Lemma 5
The eigenfunctions of the operator ∂ α | x | on (0 , with homogeneous Neumannboundary conditions are given by { cos kπx } ∞ k =1 with the corresponding eigenvalues { σ · ( kπ ) α } ∞ k =1 . The technical proofs of these statements are postponed to the Appendix.
Proposition 2
The numerical approximation defined in the scheme (23) is consistentwith the problem in (17) in the maximum norm and the order of the consistency is O ( τ ) + O ( h ) . roof: It is sufficient to prove that the right hand side of (23) provides a first orderapproximation for the Riesz derivative of order | α | . According to the proof of Theorem1 the analytic solution u N of (18) is periodic, smooth and it satisfies the homogeneousNeumann boundary conditions in and . Therefore, its cosine Fourier series is pointwiseconvergent: u ( t, x ) = ∞ X k =0 F k cos kπx x ∈ (0 , , t ∈ R + , (27)where we do not denote the time dependence of the cosine Fourier coefficients F k . Usingagain that u ( t, · ) ∈ C ∞ ( R ) we also have - using the regularity theory of Fourier series -that for any r ∈ N there exists a constant C r such that | F k | ≤ C r ( kπ ) r ∀ k ∈ N + . (28)Using (27) componentwise for t = nδ we have that u ( hδ, · ) = [ u ( nδ, x π ) u ( nδ, x π ) . . . u ( nδ, x N π )] T = [ ∞ X k =0 F k cos kx ∞ X k =0 F k cos kx . . . ∞ X k =0 F k cos kx N ] T = ∞ X k =0 F k v kh . (29)Using Lemma 5 for the expansion in (27) and the matrix A α,h for (29) according to (41)we obtain the following equality: ∂ α | x | u ( nδ, x j π ) − [ A α,h u ( nδ, · )] j = ∞ X k =0 − σ ( kπ ) α F k cos kπx j − " A α,h ∞ X k =0 F k v kh j = ∞ X k =0 F k cos kπx j (cid:18) − σ ( kπ ) α − C σ h α α sin α kπh · cos (cid:16) kπh + α π − kπh ) (cid:17)(cid:19) = − ∞ X k =0 σF k cos kπx j ( kπ ) α − sin kπh kπh ! α cos (cid:0) kπh + α ( π − kπh ) (cid:1) cos kπh ! . (30)To prove the proposition we first verify that for a mesh-independent constant C thefollowing inequality is valid: (cid:12)(cid:12)(cid:12)(cid:12) − (cid:18) sin ss (cid:19) α (cid:16) cos s ( α − − sin s ( α −
2) tan απ (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) ≤ Cs, (31)16hich will be applied with s = kπh . We first verify that h ( s ) := 1 − (cid:18) sin ss (cid:19) α (cid:16) cos s ( α − − sin s ( α −
2) tan απ (cid:17) ≤ Cs, (32)where lim s → h ( s ) = 0 . Therefore, it is sufficient to prove that h ′ is bounded on [0 , π ] .Obviously, h ′ ( s ) = α (cid:18) sin ss (cid:19) α − sin s − s · cos ss · ( α − h cos s ( α − − sin s ( α −
2) tan απ i + (cid:18) sin ss (cid:19) α (cid:16) ( α − (cid:16) sin s ( α − − cos s ( α −
2) tan απ (cid:17)(cid:17) , where lim s · cos s − sin ss = lim s sin s s = 0 . Hence, all components in the expansion of h ′ ( s ) are bounded, which really verifies (32).Using (32) in (30) and applying (28) we obtain the following estimation: (cid:12)(cid:12)(cid:12) ∂ α | x | u ( nδ, x j π ) − [ A α,h u ( nδ, · )] j (cid:12)(cid:12)(cid:12) ≤ C h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ∞ X k =0 σF k cos kπx j ( kπ ) α kπ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ C · σ h ∞ X k =0 | F k | ( kπ ) α +1 ≤ C C · σ h ∞ X k =1 kπ ) ( kπ ) α +1 ≤ C C · σ h ∞ X k =1 kπ ) , (33)which completes the proof. (cid:3) Theorem 2
The numerical approximation defined in the scheme (23) converges to thesolution of (17) in the maximum norm for α ∈ (1 , and the order of the convergence is O ( τ ) + O ( h ) .Proof: Using Proposition 2 we only have to verify the stability of (23). For this, werewrite it into a linear system ( I − τ A α,h ) u n +1 = u n . Using Proposition 1 we obtain that the diagonal of I − τ A α,h is strictly positive and hasnon-positive off-diagonals. Moreover, using (24), a simple calculation shows that for the17ndices j = 1 , , . . . , N − we have (cid:2) ( I − τ A α,h ) · (1 , , . . . , T (cid:3) j = 1 + τh α − j +1 X k =0 k =1 g k − ∞ X l =0 j +2+ N X k = j +2 g k +2 l ( N +1) − ∞ X l =0 j +2 N +3 X k = j + N +3 g k +2 l ( N +1) + τh α − N − j +1 X k =0 k =1 g k − ∞ X l =0 2 N − j +2 X k = N − j +2 g k +2 l ( N +1) − ∞ X l =0 3 N − j +3 X k =2 N − j +3 g k +2 l ( N +1) . (34)Observe that in the brackets in (34) each coefficient g i , i = 1 appears once (see theExample and the remark after (26)). According to (8), we obtain α − ∞ X k =0 k =1 g k = − ∞ X k =0 g k = 0 and therefore, the sum in both brackets on the right hand side of (34) is zero such thatthe entire right hand side is one . Using (25) and (26), an obvious modification of (34)gives its positivity both for the indices j = 0 and j = N . Therefore, (cid:2) ( I − τ A α,h ) · (1 , , . . . , T (cid:3) j = 1 is valid for all j = 0 , , . . . , N . In this way, ( I − τ A α,h ) − elementwise positive such that k ( I − τ A α,h ) − k ∞ = 1 , and consequently, the scheme in (23) is unconditionally stable. (cid:3) Construction of the matrix A α,h Whenever the coefficient in the matrix A α,h arebased on an infinite number of grid points, it can be computed in concrete terms.For this we introduce B α,h ∈ R ( N +1) × ( N +1) with B α,h = (cid:0) v h , v h , . . . , v Nh (cid:1) , which consists of the eigenvectors of A α,h , see Lemma 4. Then ( I − τ A α,h ) B α,h = (cid:0) (1 − τ λ ) v h , (1 − τ λ ) v h , . . . , (1 − τ λ ) v Nh (cid:1) , and therefore, τ A α,h = I − (cid:0) (1 − τ λ ) v h , (1 − τ λ ) v h , . . . , (1 − τ λ ) v Nh (cid:1) B − α,h , (35)where on the right hand side all terms can be computed.18 omplexity and extension to higher-order methods Since we have explored theeigenvectors of A α,h we do not have to compute its components in the practice as a series.We simply obtain the matrix using (35) such that the extension does not result in extracomputational costs.Higher order methods can be obtained in the same fashion. In the semidiscretization,we should then choose a higher order spatial approximation, e.g. , the one in (12) and theaccuracy of the time stepping can also be increased, e.g. , using a Crank–Nicolson scheme.For homogeneous Dirichlet boundary conditions, such a study is performed (even for themultidimensional case) in [4]. A homogeneous model problem
We first investigate the model problem ∂ t u ( t, x ) = 0 . ∂ . | x | u ( t, x ) x ∈ (0 , , t ∈ (0 , u (0 , x ) = x − x x ∈ (0 , ∂ x u ( t,
0) = ∂ x u ( t,
1) = 0 t ∈ (0 , , which is converted to the well-posed extended problem ∂ t u ( t, x ) = 0 . ∂ . | x | u ( t, x ) x ∈ (0 , , t ∈ (0 , u (0 , x ) = \ x − x N x ∈ R . (36)The analytic solution of (36) on (0 , × (0 , ) is u ( t, x ) = − ∞ X k =1 ( − k +1 kπ ) e − t ( kπ ) . cos( kπx ) x ∈ (0 , , t ∈ (0 , , which has been computed in the grid points with a high accuracy to verify the convergenceof the implicit Euler method. The results of the computations are summarized in Table 1.We computed the error k e τ,h k ∞ of the approximation in maximum-norm for various timesteps and discretization parameters at the final time t = 1 . One can clearly see the firstorder convergence which was predicted by the theory, see Theorem 2. The convergencerate was estimated in the consecutive refinement steps using the formula log (cid:16) k e τ, h k ∞ k e τ,h k ∞ (cid:17) .19able 1: Convergence results for the implicit Euler method in (23) applied to (36)Grid parameter ( h ) Time step ( τ ) Error in k · k ∞ -norm Convergence rate1/2 1/2 1.4 · − ∅ · − · − · − · − · − · − · − · − l norm should be constant, which is an easy consequenceof the fact, that the sum of the elements in the columns of A α,h, ∞ is zero.Therefore, we examine the boundary conditions in course of the simulations. For this,we use the second order accurate approximation ∂ x u ( t,
0) = 12 h (3 u ( t, − u ( t, h ) + u ( t, h )) (37)and the results are shown in Figure 3. For the simplicity, we applied the same number ofgrid points in the spatial and the time coordinates. This accurate approximation can berecognized as the numerical equivalent of Lemma 3. An inhomogeneous model problem
We second investigate the model problem ∂ t u ( t, x ) = 0 . · ∂ . | x | u ( t, x ) + e − t cos πx x ∈ (0 , , t ∈ (0 , u (0 , x ) = 2 x − x x ∈ (0 , ∂ x u ( t,
0) = ∂ x u ( t,
1) = 0 t ∈ (0 , , (38)which is converted to the well-posed extended problem ( ∂ t u ( t, x ) = 0 . · ∂ . | x | u ( t, x ) + e − t cos πx x ∈ (0 , , t ∈ (0 , u (0 , x ) = \ x − x N x ∈ R . (39)20
10 15 2000.0010.0020.0030.004 N u ’ h ( ) Figure 3: The approximation (37) of the derivative ∂ x u (1 , in the numerical simulationsvs. the number of gridpoints.We split the original equation in (38) into two ones as ∂ t u ( t, x ) = 0 . · ∂ . | x | u ( t, x ) x ∈ (0 , , t ∈ (0 , u (0 , x ) = 2 x − x x ∈ (0 , ∂ x u ( t,
0) = ∂ x u ( t,
1) = 0 t ∈ (0 , and ∂ t u ,t ∗ ( t, x ) = 0 . · ∂ . | x | u ,t ∗ ( t, x ) x ∈ (0 , , t ∈ (0 , u ,t ∗ (0 , x ) = e − t ∗ cos πx x ∈ (0 , ∂ x u ,t ∗ ( t,
0) = ∂ x u ,t ∗ ( t,
1) = 0 t ∈ (0 , , the analytic solution of which be given as u ( t, x ) = 13 + ∞ X k =1 kπ ) e − t ( kπ ) . (( − k −
1) cos( kπx ) x ∈ (0 , , t ∈ (0 , , and u ,t ∗ ( t, x ) = e − t . − tπ x ∈ (0 , , t ∈ (0 , . These have been computed in the grid points with a high accuracy to verify the conver-gence of the implicit Euler method. 21hen, according to the Duhamel’s principle (see [8], p. 49) we have u ( t, x ) = u ( t, x ) + Z t u ,t ∗ ( t − t ∗ , x ) d t ∗ . Accordingly, for the numerical solution at t = 1 we compute first the approximationof u ( t, · ) at the grid points. Then with the same time steps and spatial accuracy weapproximate u , (1 , · ) , u ,τ (1 − τ, · ) , u , τ (1 − τ, x ) , . . . , u , (0 , x ) , where indeed, the last term is already given. Then a composite trapezoidal numericalintegration u (1 , x ) ≈ τ u , (1 , · ) + τ ( u ,τ (1 − τ, · ) + u , τ (1 − τ, x ) + · · · + u , − τ ( τ, x )) + τ u , (0 , x ) gives the desired approximation in x .The results of the computations are summarized in Table 2, where in all cases τ =1 / such that the numerical integration does not harm the predicted order of conver-gence.Table 2: Convergence results for the implicit Euler method in (23) applied to (38)Grid parameter ( h ) Time step ( τ ) Error in k · k ∞ -norm Convergence rate1/2 1/2 6.2 · − ∅ · − · − · − · − · − · − · − · − Appendix
Proof of Lemma 4:
We first observe that the Neumann extension is the natural one for v kh in the sense that (cid:20)c v kh N (cid:21) j = cos kx j j ∈ Z . (40)22hen according to (22), (40), (21), (6) and (7) we obtain (cid:2) A α,h v kh (cid:3) j = ∞ X l =0 C σ h α g l (cid:20)c v kh N (cid:21) j + l − + g l (cid:20)c v kh N (cid:21) j − l +1 ! = C σ h α ∞ X l =0 g l (cos kx j + l − + cos kx j − l +1 )= C σ h α ∞ X l =0 g l (cid:18) cos kπh (cid:18) j + 12 + l − (cid:19) + cos kπh (cid:18) j + 12 − ( l − (cid:19)(cid:19) = 2 C σ cos kx j h α ∞ X l =0 g l cos kπh ( l − C σ cos kx j h α R ∞ X l =0 exp {− ikπh } g l exp { ikπhl } ! = 2 C σ cos kx j h α R exp {− ikπh } ∞ X l =0 ( − l (cid:18) αl (cid:19) exp { ikπhl } ! = 2 C σ cos kx j h α R (exp {− ikπh } (1 − exp { ikπh } ) α )= 2 C σ cos kx j h α R (cid:18) exp {− ikπh } · α sin α kπh · (cid:16) cos α π − kπh ) − i sin α π − kπh ) (cid:17)(cid:19) = 2 C σ cos kx j h α α sin α kπh · (cid:16) cos kπh cos α π − kπh ) − sin kπh sin α π − kπh ) (cid:17) = 2 C σ cos kx j h α α sin α kπh · cos (cid:16) kπh + α π − kπh ) (cid:17) (41)where the have used the identity (1 − exp { ikπh } ) α = (cid:18) kπh · (cid:18) sin kπh − i cos kπh (cid:19)(cid:19) α = 2 α sin α k πh · (cid:18) cos (cid:18) π − kπh (cid:19) − i sin (cid:18) π − kπh (cid:19)(cid:19) α = 2 α sin α kπh · (cid:16) cos α π − kπh ) − i sin α π − kπh ) (cid:17) . The definition of C σ gives then the statement in the lemma. (cid:3) roof of Lemma 5: In the proof, we use the identities Z ∞ x n − cos bx d x = Γ( n ) b n cos nπ , Z ∞ x n − sin bx d x = Γ( n ) b n sin nπ , (42)which can be found in [9], 3.761/9.Observe that the even extension of the cos( kπ · ) | (0 , function to the real axis is the cos( kπ · ) function itself. On the other hand, as it was pointed out in [26], we can differen-tiate the integrals in the Riemann–Liouville formula to obtain ∂ α | x | cos( kπx ) = − σ · ( kπ ) (cid:0) απ (cid:1) (cid:18) − Z x ∞ cos kπs ( x − s ) α − d s − Z x ∞ cos kπs ( s − x ) α − d s (cid:19) . (43)Using (42), the first term can be rewritten as Z x −∞ cos kπs ( x − s ) α − d s = Z ∞ cos kπ ( x − y ) y α − d y = cos kπx Z ∞ cos kπyy α − d y + sin kπx Z ∞ sin kπyy α − d y = Γ(2 − α )( kπ ) − α (cid:18) cos kπx cos π (2 − α )2 + sin kπx sin π (2 − α )2 (cid:19) . A similar computations gives that Z ∞ x cos kπs ( s − x ) α − d s = Z ∞ cos kπ ( x + y ) y α − d y = Γ(2 − α )( kπ ) − α (cid:18) cos kπx cos π (2 − α )2 − sin kπx sin π (2 − α )2 (cid:19) . Therefore, the equality in (43) can be rewritten as ∂ α | x | cos( kπx ) = σ · ( kπ ) (cid:0) απ (cid:1) Γ(2 − α ) Γ(2 − α )( kπ ) − α kπx cos π (2 − α )2= σ · ( kπ ) α (cid:0) απ (cid:1) kπx cos π (2 − α )2 = − σ · ( kπ ) α cos kπx. On the other hand, the system { cos kπx } ∞ k =1 is complete in L (0 , such that no furthereigenfunctions can exist. (cid:3) cknowledgements The authors acknowledge the financial support of the Hungarian National Research FundOTKA (grants K104666 and K112157).
References [1] R. A. Adams and J. J. Fournier.
Sobolev spaces . Academic Press, Amsterdam, secondedition, 2003. Pure and Applied Mathematics, Vol. 140.[2] D. A. Benson, S. W. Wheatcraft, and M. M. Meerschaert. Application of a frac-tional advection-dispersion equation.
Water Resources Research , 36(6):1403–1412(electronic), 2000.[3] A. Bonito and J. Pasciak. Numerical approximation of fractional powers of ellipticoperators. Submitted, arXiv:1307.0888.[4] W. H. Deng and M. Chen. Efficient numerical algorithms for three-dimensionalfractional partial diffusion equations.
J. Comp. Math. , 32(4):371–391, 2014.[5] Q. Du, M. Gunzburger, R. Lehoucq, and K. Zhou. Analysis and approximationof nonlocal diffusion problems with volume constraints.
SIAM Rev. , 54(4):667–696,2012.[6] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou. A nonlocal vector calculus,nonlocal volume-constrained problems, and nonlocal balance laws.
Math. Mod. Meth.Appl. Sci. , 23(03):493–540, 2013.[7] A. M. Edwards, R. A. Phillips, N. W. Watkins, M. P. Freeman, E. J. Murphy,V. Afanasyev, S. V. Buldyrev, M. G. E. da Luz, E. P. Raposo, H. E. Stanley, andG. M. Viswanathan. Revisiting Lévy flight search patterns of wandering albatrosses,bumblebees and deer.
Nature , 449:1044–1048, 2007.[8] L. C. Evans.
Partial differential equations , volume 19 of
Graduate Studies in Math-ematics . American Mathematical Society, Providence, RI, 1998.[9] I. S. Gradshteyn and I. M. Ryzhik.
Table of integrals, series, and products . Aca-demic Press Inc., San Diego, CA, sixth edition, 2000. Translated from the Russian,Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger.2510] R. Hilfer. Threefold introduction to fractional derivatives. In R. Klages, G. Radons,and I. Sokolov, editors,
Anomalous Transport: Foundations and Applications , pages17–73. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, 2008.[11] J.-F. Huang, N. Nie, and Y.-F. Tang. A second order finite difference-spectral methodfor space fractional diffusion equations.
Sci. Chin. Math. , 57(6):1303–1317, 2014.[12] M. Ilic, F. Liu, I. Turner, and V. Anh. Numerical approximation of a fractional-in-space diffusion equation. I.
Fract. Calc. Appl. Anal. , 8(3):323–341, 2005.[13] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo.
Theory and applications of frac-tional differential equations , volume 204 of
North-Holland Mathematics Studies . El-sevier Science B.V., Amsterdam, 2006.[14] M. M. Meerschaert and C. Tadjeran. Finite difference approximations for fractionaladvection-dispersion flow equations.
J. Comput. Appl. Math. , 172(1):65–77, 2004.[15] K. S. Miller and B. Ross.
An introduction to the fractional calculus and fractionaldifferential equations . A Wiley-Interscience Publication. John Wiley & Sons Inc.,New York, 1993.[16] R. Nochetto, E. Otárola, and A. Salgado. A pde approach to fractional diffusion ingeneral domains: A priori error analysis.
Foundations of Computational Mathematics ,2014.[17] I. Podlubny.
Fractional differential equations , volume 198 of
Mathematics in Scienceand Engineering . Academic Press Inc., San Diego, CA, 1999. An introduction tofractional derivatives, fractional differential equations, to methods of their solutionand some of their applications.[18] S. G. Samko, A. A. Kilbas, and O. I. Marichev.
Fractional integrals and derivatives .Gordon and Breach Science Publishers, Yverdon, 1993. Theory and applications,Edited and with a foreword by S. M. Nikol ′ ski˘ı, Translated from the 1987 Russianoriginal, Revised by the authors.[19] S. Shen, F. Liu, V. Anh, and I. Turner. The fundamental solution and numericalsolution of the Riesz fractional advection-dispersion equation. IMA J. Appl. Math. ,73(6):850–872, 2008.[20] P. Szymczak and A. J. C. Ladd. Boundary conditions for stochastic solutions of theconvection-diffusion equation.
Phys. Rev. E , 68:036704, Sep 2003.2621] C. Tadjeran and M. M. Meerschaert. A second-order accurate numerical method forthe two-dimensional fractional diffusion equation.
J. Comput. Phys. , 220(2):813–823,2007.[22] C. Tadjeran, M. M. Meerschaert, and H.-P. Scheffler. A second-order accurate numeri-cal approximation for the fractional diffusion equation.
J. Comput. Phys. , 213(1):205–213, 2006.[23] W. Tian, H. Zhou, and W. Deng. A class of second order difference approximationfor solving space fractional diffusion equations.
Math. Comp. , 2014. to appear, arXiv:1201.5949[math.NA] .[24] R. A. Treumann. Theory of super-diffusion for the magnetopause.
Geophys. Res.Lett. , 24(14):1727–1730, 1997.[25] H. Wang and T. S. Basu. A fast finite difference method for two-dimensional space-fractional diffusion equations.
SIAM J. Sci. Comput. , 34(5):A2444–A2458, 2012.[26] Q. Yang, F. Liu, and I. Turner. Numerical methods for fractional partial differentialequations with Riesz space fractional derivatives.
Appl. Math. Model. , 34(1):200–218,2010.[27] H. Zhou, W. Tian, and W. Deng. Quasi-compact finite difference schemes for spacefractional diffusion equations.