Finite Element Method for a Kind of Two-Dimensional Space-Fractional Diffusion Equation with Its Implementation

In this article, we consider a two-dimensional symmetric space-fractional diffusion equation in which the space fractional derivatives are defined in Riesz potential sense. The well-posed feature is guaranteed by energy inequality. To solve the diffusion equation, a fully discrete form is established by employing Crank-Nicolson technique in time and Galerkin finite element method in space. The stability and convergence are proved and the stiffness matrix is given analytically. Three numerical examples are given to confirm our theoretical analysis in which we find that even with the same initial condition, the classical and fractional diffusion equations perform differently but tend to be uniform diffusion at last.


Introduction
Fractional convection-diffusion equations are generalizations of classical convection-diffusion equations, which have come to be applied in Physics [1]- [4], hydrology [5] [6] and many other fields.As it is difficult to get the analytic solutions of these equations, numerical approaches to different type of fractional convection-diffusion equations are proposed in recent years.Tadjeran et al. [7] considered one-dimensional space-fractional diffusion equation with variable coefficient by fractional Crank-Nicholson method based on the shifted Grünwald formula, and obtained an unconditional stable second-order accurate numerical approximation by extrapolation.Later, Tadjeran and Meerschaert [8] utilized the classical alternating directions implicit (ADI) approach with a Crank-Nicholson discretization and a Richardson extrapolation to solve two-dimensional space-fractional diffusion equation, and proved it is unconditional stable second-order accurate.Sousa [9] derived an implicit second-order accurate numerical method which used a spline approximation for space-fractional diffusion equation and the consistency and stability were examined.A space-time spectral method for time fractional diffusion equation was developed by Li and Xu [10], in which the convergence was proven and priori error estimate was given.Xu [11] proposed a discontinuous Galerkin method for one-dimensional convection-subdiffusion equations with fractional Laplace operator and derived stability analysis and optimal convergence rate.Jin et al. [12] gave a full discretization scheme for multi-term time-fractional diffusion equation by using finite difference method in time and finite element method in space, and discussed its stability and error estimate.
The symmetric space-fractional convection-diffusion equation (including both left and right derivatives) was firstly proposed by Chaves [13] to investigate the mechanism of super-diffusion and was later generalized by Benson et al. [14] [15].It is a powerful approach for a description of transport dynamics in complex systems governed by anomalous diffusion.Zhang [16] et al. considered one-dimensional symmetric space-fractional partial differential equations with Galerkin finite element method in space and a backward difference technique in time, and the stability and convergency were proven.Sousa [17] derived a second order numerical method for one-dimensional symmetric space-fractional convection-diffusion equation and studied its convergence.
Recently, numerical methods for multi-dimensional problems of fractional differential equational are studied.For example, in [18], a semi-alternating direction method for a 2-D fractional reaction diffusion equation are proposed to solve FitzHugh-Nagumo model on an approximate irregular domain.In [19], Crank-Nicolson ADI spectral method is presented to approximate the two-dimensional Riesz space fractional nonlinear reactiondiffusion equation.In [20] [21], Wang and Du proposed fast finite difference methods to compute three-dimensional space-fractional diffusion equations, which reduce the computational cost a lot.
In this paper, we consider the following two-dimensional symmetric space-fractional diffusion equation (SSFDE) ( y t u x y t u x y t a f x y t t x y where 1 2 α < ≤ , 0 a > is a constant, ( ) , , , , , , , , Remark: In this paper, the default fractional derivative is Riemann-Liouville derivative.This article is organized as follows.In Section 2, we introduce some functional spaces.In Section 3 and Section 4, we prove existence and uniqueness of the variational solution.The full discretization of SSFDE is given in Section 5, where we apply Crank-Nicolson technique in time and Galerkin finite element method in space.Moreover, a detailed stability and convergence analysis is carried out.In section 6, we present the implementation of how to get the stiffness matrix.Finally, some numerical examples are given in Section 7 to confirm our theretical analysis and to compare the difference between fractional diffusion and integer order diffusion system.: where û denotes the Fourier tansform of u with variable . Also let In the following, a semi-norm is defined by integral ( ) Theorem 2.4 (Fractional Poincarà Friedrichs Inequality [23]).For ( ) The definitions and theorems above are basic frame of multi-dimensional fractional derivative spaces.In terms of Equation (1), we let M be atomic with atoms { } or and norm It is easy to derive that ( 12) is equivalent to (13) with using theorem 2.1 and Parseval equality.Lemma 2.1 (The relationship between R-L and Caputo fractional order derivatives [24]).Assume that the derivatives ( ) ( ) ( ) a t and where C p a t D denotes Caputo fractional order derivative, which is defined as , the two kinds of derivates is equivalent, i.e.

(
) Proof.Let m p σ = − , combining Lemma 2.2 and Lemma 2.3 we have , , , , , For convenience, we denote then Equation ( 1) can be written in the following form To derive the variational form of (23), we introduce two properties of ( ) where ( ) ˆ, e , d d .
Proof.In view of Theorem 2.1, we can derive the Fourier Transform Remark: Here, we use ( ) to make difference from the fractional Laplace operator ( ) , which defined as [25] [26] and its Fourier Transform is ( ) .

α
, the formula is the classical Green formula.Proof.Using Theorem 2.5 and taking notice that / 2 (0,1)

Variational Formulation
In order to derive the variational form of (23), we assume u is a sufficiently smooth solution of (23), and multiply by arbitrary ( ) The weak formulation of the equation is to find the ( ) With using property 2, the above formula could be written as

∫∫
Thus we define the associated bilinear form ( ) ( ) Proof.According to the definition of ( ) ( ) Associating the definition of the semi-norm of ( ) J α Ω and using Young's inequality it follows Combining the equivalence of ( ) According to the equivalence of ( ) θ Ω , and combining Theorem 2.4 we can obtain i.e., the form ( ) × Ω .Then, we can obtain the energy estimate i.e. the solution of ( 23) is well posed.
Proof.Multiply the first formula of ( 23) by u and integrate both sides of the equation in Ω , then we have As the coercivity of the form ( ) , B ⋅ ⋅ and Young's inequality, we obtain and integrating over ( ) to the above inequality we get

Crank-Nicolson-Galerkin Finite Element Fully Discrete System
Let h S denote a uniform of partition of Ω , with grid parameter h, and . h V denote the continuous functions on G.We define the finite dimensional subspace with the piecewise polynomials k P of order ( ) or less then k.Taking a uniform mesh for the time variable t and let ( ) where 0 τ > being the time step, and . Then by the Galerkin finite element method and Crank-Nicolson technique, ( 23) is transformed into the following problem: find ( ) , where 0,h Then the first formula of (32) can be written as , .
In view of Theorem 4.1, we have Therefore, the bilinear ( ) in (32), noticing the coercivity of the bilinear form ( ) , B ⋅ ⋅ and employing Hölder inequality, we have Then we can obtain So the result is valid.Lemma 4.1 (Approximation Property [27]) Let ( ),0 , : V which is defined as follows for each ( ) Looking back to the first formula of (32), we can derive , , , , and with using (37), we can obtain

∫∫
and combining Cauchy-Schwarz inequality, we can obtain In the following we will estimate the three parts of the above inequality respectively.The first part ( ) ( ) Hence we can obtain a recursive inequality Take 0 s = in (35), then we can obtain From ( 38) and (39), we can derive the following error estimate Finally, the formula (40) leads to (36).

Computational Implementation
Since the fractional derivative is a non-local operator, the implementation of finite element method for fractional differential equations is very complex.The main problem is how to obtain the stiffness matrix.In [28], Roop investigated the computational aspects of the Galerkin approximating using continuous piecewise polynomial basis functions on a regular triangulation of the domain.In this section we give the computational details, in which the bilinear functions are chosen as the basis functions.The computational domain is and the number of computational grid is 1 2 N N × .First of all, we consider the problem of finding the fractional derivative of each of the basis function  1).It is defined as where ( ) ( ) ( ) ( ) ,  ,  ,  ,  x y x y x y x y are the centers of the blocks 1 2 3 4 , , , K K K K , and .Assume the coordinate of the ith node (see Figure 1) is ( ) x y , then we can derive ( ) , , , we , , x x y y y , replacing 1 0 y by 4 0 y in the three cases above respectively, we can get 2 a x i D α ψ in the corresponding region.
Secondly, we consider the problem of calculating the inner product ( ) , is the number of inner points.Denote ( ) then the coordinate of the jth node is ( ) x y .It is easy to know when , 0 For the other cases, we present the results here.Please see the appendix for the expatiation.Case 1: ( ) ( ) ( ) Finally, we consider the problem of calculating the stiffness matrix A via the inner product obtained.Form Equation (29) we can see that A can be decomposed into four parts 1 2 3 4 , , , A A A A .With ignoring the coefficient ( ) , , A i j 3 In fact, for i ψ and j ψ , if we start numbering these nodes along the direction of y axis and rename the two basis functions to y i ψ and y j ψ , then we have where , , , m n p q are defined in (41) and (42).
( ) ( ) ( ) ( ) , in 0, , , 0 , , 0 500 0.25 0.25 , on 0.5, 0.5 0.5, 0.5    .Remark: The trial function in all of the numerical experiments is bilinear function.We can see that the results support our error estimate and ensure the numerical approximation is effective.In the following, we take fixed initial value and source term independent of α to try to describe the character of the solution with the change of α .We note that the initial condition in the fractional system affect wider area than integer order in a short period of time by comparing the first two contour maps.Moreover, the diffusion under the influence of initial condition last longer in the fractional system.So at 0.35 t = , the diffusion in classical system is almost uniform in every direction but this state needs more time to reach in the fractional system (see the left map at 1 t = ).

Conclusion
Many different numerical methods for fractional convection-diffusion equation have been discussed by researchers in recent 10 years.In this paper, we discussed one kind of space-fractional diffusion equation which could be derived through replacing the second order derivative of x and y by corresponding Riesz fractional derivative in the classical diffusion equation.A numerical approximation for the equation was presented by using C-N technique in time direction and Galerkin finite method in space.Furthermore, a detailed stability and convergence analysis was carried out for the fully discrete system.Then, some numerical examples were given and the differences between fractional and classical diffusion were presented.It is known that the stiffness matrix of fractional differential equation is rather complex, so to make the approach applicatory.We give the implementation of computational aspect.However, because of the non-local property of fractional derivative, the stiffness matrix is not sparse (almost dense) which challenges the computational resources.1, p n q m =− = , i.e.
With noticing that i ψ and j ψ are both symmetrical about the straight lines m y y = , we have ( ) are equivalent with equivalent semi-norms and norms.
The solution of variational formulation(28) exists and is unique.Proof.The existence can be derived directly from Theorem 3.6 with Lax-Milgram theorem and Theorem 3.7 ensure the uniqueness.

D
α ψ .The support set of the ith basis function i

Figure 1 .
Figure 1.Sketch for the element and node number.
the numerical solution and error visually, we present the surfaces of n
the right side of (33) is continuous.According to Lax-Milgram theorem, the fully discrete approximating system (32) has unique solution n 36) d

Table 1 .
Experimental error and convergence rate in 2