A Mixed Finite Element Formulation for the Conservative Fractional Diffusion Equations

We consider a boundary-value problem of one-side conservative elliptic equation involving Riemann-Liouville fractional integral. The appearance of the singular term in the solution leads to lower regularity of the solution of the equation, so to the lower order convergence rate for the numerical solution. In this paper, by the dividing of equation, we drop the lower regularity term in the solution successfully and get a new fractional elliptic equation which has full regularity. We present a theoretical framework of mixed finite element approximation to the new fractional elliptic equation and derive the error estimates for unknown function, its derivative, and fractional-order flux. Some numerical results are illustrated to confirm the optimal error estimates.

Equations (1a) and (1b), delicately derived from random walk model [1], fractional conservation of mass [2], and the impact of boundary [3], describe subsurface fluid flow and solute transport processes taking place in a multiscale heterogeneous medium/aquifer (see [4] and the references cited therein).These processes exhibit anomalous diffusion or non-Fickian diffusion that cannot be modeled adequately and properly by classical second-order diffusion equation and thus have contracted considerable attention in practical applications.
Considering the conservation property of the diffusion processes and the application of the fractional-order flux  = − 0     in engineering, an ideal numerical procedure should recognize both the unknown function and its flux and obey the conservation of mass to reflect the physical character of the diffusion model ((1a) and (1b)).Hence, one should involve the state variable  and the fractional-order flux  = − 0     to form a saddle-point formulation then design a locally conservative numerical procedure.For this purpose, by introducing two intermediate variables  =  and the flux  = − 0    , Chen and Wang [19] proposed a saddle-point framework and developed a locally conservative expanded mixed finite element method recently.Although a rigorous numerical analysis theory was established in [19], the 2 Advances in Mathematical Physics convergence rates for the unknown  and the flux  are not optimal with respect to regularity requirement for sufficiently smooth solution .In particular, for the solution  only in  1−+ (Ω) with  ∈ [/2, 1/2), the convergence rates are heavily destroyed to  − /2, due to the appearance of the singular term  1− .
The aim of this paper is to express the solution of (1a) and (1b) by a full-regular solution satisfying a general fractional diffusion equation with -dependent right-hand side term and a  1− -type solution satisfying an analytic-solved fractional equation.Then, the expanded mixed finite element method proposed in [19] can be applied to solve the dependent fractional equation and obtain good error estimates whatever the regularity of the solution  is.
The rest of the paper is organized as follows.In Section 2, we revisit the fractional differential operators and their properties.In Section 3, we split the fractional diffusion equations (1a) and (1b) into two equations: one is an -dependent fractional diffusion equation that permits a fully regular solution, and the other is an analytic-solved fractional equation from which the  1− -type analytic solution is easily solved.By doing so, the solution  of (1a) and (1b) is expressed as the sum of the -dependent solution and the  1− -type solution.We are devoted to discretizing the -dependent fractional equation by the expanded mixed finite element procedure in Section 4 and deriving the related error estimates in Section 5. Obviously, the error estimates for the original fractional diffusion equations (1a) and (1b) obey those for the -dependent fractional equation since the  1− -type solution is analytically solved.In Section 6, numerical experiments are performed to confirm our theoretical findings.

Fractional Operators and Fractional Sobolev Spaces
In this section, we will recall the definitions and properties of Riemann-Liouville fractional differential operators and fractional Sobolev spaces.
The left and right Riemann-Liouville fractional derivatives of order  are We note that fractional integral operator [20][21][22] satisfies semigroup property; that is, for ,  > 0, and adjoint property We also observe that the left (right) fractional derivative operator is a left inverse of the left (right) fractional integral operator [20][21][22]: We introduce the fractional derivative spaces and present their well-established equivalence to fractional-order Sobolev spaces [9,19,23].
Definition 3 (see [9]).For  > 0, define the norm with seminorm Let    (R) be the left fractional derivative space defined as the closures of  ∞ 0 (R) with respect to norm ‖ ⋅ ‖    (R) .Definition 4 (see [23]).For  > 0, define the norm and let The right fractional derivative spaces    (R) and  −  (R) and their respective norms are defined analogously.
Definition 5 (see [9]).For  > 0, define the norm with seminorm where û() denotes the Fourier transform of Let   (R) be the left fractional derivative space defined as the closures of  ∞ 0 (R) with respect to norm ‖ ⋅ ‖   (R) .
We then reiterate the equivalence theories established in [9,23] for the fractional derivative spaces.
We now restrict the fractional derivative spaces to Ω. Definition 9 (see [9]).Let  > 0. Define the spaces  For simplicity of presentation, we use ‖⋅‖  and |⋅|  to stand for the norms and seminorms of fractional-order Sobolev spaces   (Ω) and ‖ ⋅ ‖ − to stand for the norms of negative fractional-order Sobolev spaces  − (Ω), respectively.
Remark 15.From Theorem 13, the regularity of original problem (1a) and (1b) is only  (Ω), due to the presence of the singular term  1− , while the regularity to ( 18) is In particular, for the right term () ∈   (Ω),  ≥ 0, the regularity can be improved to  +2− 0 (Ω) by choosing a proper decomposition.For example, we can replace  2 by   with  ≥  + 2 −  in the decomposition of the solution when () ∈   (Ω).This is the advantage of the splitting.

Saddle-Point Formulation and Mixed Finite Element Method
In this section, we shall formulate a saddle-point variational formulation and establish a mixed finite element formulation to the -dependent fractional diffusion equations (16a) and (16b).

Convergence Analysis
In this section, we are devoted to the convergence analysis.
In error analysis below, we shall make use of three projection operators.The first operator is Π ℎ :   →  ℎ defined by The other two operators are the  2 -projections  ℎ :   →  ℎ and  ℎ :   →  ℎ : Lemmas 18-20 hold in the integer-order Sobolev spaces (if  = 1 see [19] and if  = 2, 3 see [24][25][26][27]).We can prove that they also hold in the fractional-order spaces by applying standard real interpolation methods [26].

Numerical Experiments
In this part, we present some numerical experiments to verify the theoretically proven convergence results.In all the following examples, we let  = 1.For each example, we consider three different  values, that is, 1/4, 1/2, and 2/3, and present the estimates and the convergence orders of ‖ −  ℎ ‖, ‖− ℎ ‖, and ‖− ℎ ‖  −/2 (Ω) , where the number in the bracket under the column order is the theoretical convergence rate derived in Theorem 26.
Remark 5. From the above four examples, we can see that the error rates for ‖ −  ℎ ‖  are identical to the predicted rates of convergence in Theorem 26.However, numerical convergence rates for ‖ −  ℎ ‖ and ‖ −  ℎ ‖ are much higher than the theoretical prediction of Theorem 26.

Table 1 :
Numerical result of Example 1 on a uniform mesh of size ℎ for  = 1.

Table 2 :
Numerical result of Example 1 on a uniform mesh of size ℎ for  = 2.

Table 3 :
Numerical result of Example 2 on a uniform mesh of size ℎ for  = 1.

Table 4 :
Numerical result of Example 2 on a uniform mesh of size ℎ for  = 2.

Table 5 :
Numerical result of Example 3 on a uniform mesh of size ℎ for  = 1.

Table 6 :
Numerical result of Example 3 on a uniform mesh of size ℎ for  = 2.

Table 7 :
Numerical result of Example 4 on a uniform mesh of size ℎ for  = 1.

Table 8 :
Numerical result of Example 4 on a uniform mesh of size ℎ for  = 2.