A Mixed Finite Volume Element Method for Time-Fractional Reaction-Diffusion Equations on Triangular Grids

: In this article, the time-fractional reaction-diffusion equations are solved by using a mixed ﬁnite volume element (MFVE) method and the L 1-formula of approximating the Caputo fractional derivative. The existence, uniqueness and unconditional stability analysis for the fully discrete MFVE scheme are given. A priori error estimates for the scalar unknown variable (in L 2 ( Ω ) -norm) and the vector-valued auxiliary variable (in ( L 2 ( Ω )) 2 -norm and H ( div, Ω ) -norm) are derived. Finally, two numerical examples in one-dimensional and two-dimensional spatial regions are given to examine the feasibility and effectiveness.

In recent years, finite volume (FV) and finite volume element (FVE) methods have been used to solve FPDEs by many scholars. Wang and Du [23] constructed a fast locally conservative FV method to solve steady-state space-fractional diffusion equations with Riesz fractional derivatives. Hejazi et al. [24,25] provided an FV method for the time-space two-sided fractional advection-dispersion equation and space fractional advection-dispersion equation, respectively. Zhuang et al. [26] proposed FV and FE methods for the space-fractional Boussinesq equation with Riemann-Liouville derivatives on a one-dimensional domain. Yang et al. [27] constructed an FV scheme with pretreatment Lanczos method to solve two-dimensional space fractional reaction-diffusion equations. Cheng et al. [28] provided an Eulerian-Lagrangian control volume method to solve the space-fractional advection-diffusion equations with Caputo fractional derivatives. Feng et al. [29] proposed an FV scheme to treat the two-sided space-fractional diffusion equation with Riemann-Liouville fractional derivatives. Jia and Wang [30,31] constructed fast FV schemes for two kinds of space-fractional differential equations, respectively. Simmons et al. [32] proposed an FV method to solve two-sided fractional diffusion equations with Riemann-Liouville derivatives. Jiang and Xu [33] constructed a monotone FV scheme to solve the time-fractional Fokker-Planck equation. Karaa et al. [34] proposed an FVE method to solve two-dimensional fractional subdiffusion problems with the Riemann-Liouville derivative. Karaa and Pani [35] applied an FVE method to solve the time-fractional diffusion equations with nonsmooth initial data.
In this article, we will develop a mixed finite volume element (MFVE) method to solve the time-fractional reaction-diffusion equations as follows ∂t α − ε∆u(X, t) + p(X)u(X, t) = f (X, t), (X, t) ∈ Ω × J, u(X, t) = 0, (X, t) ∈ ∂Ω ×J, u(X, 0) = u 0 (X), where Ω ⊂ R 2 is a bounded convex polygonal domain with boundary ∂Ω, J = (0, T] with 0 < T < ∞. The functions u 0 (X), p(X) and f (X, t) are smooth enough, the diffusion coefficient ε > 0. And there exist constants p 0 > 0 and p 1 > 0 such that 0 Our aim is to construct an MFVE scheme [36][37][38][39][40][41] by combining the MFE method [42][43][44][45] with the FVE method [46][47][48][49][50][51] to treat the time Caputo fractional reaction-diffusion equation. The MFVE method can not only simultaneously compute several different quantities but also maintain the local conservativity, which is very important in computational fluid dynamics. In this article, we apply the L1-formula [52,53] to approximate the Caputo fractional derivative, introduce a vector-valued auxiliary variable to rewrite the model equation as the first-order coupled system, and approximate this coupled system by the MFVE method. We use the lowest order Raviart-Thomas (RT 0 ) space and the piecewise constant function space as the trial function spaces of the vector-valued auxiliary variable and scalar unknown variable, respectively. In terms of theoretical analysis, we give the existence, uniqueness and stability analysis, and obtain a priori error estimates for the scalar unknown variable (in L 2 (Ω)-norm) and the vector-valued auxiliary variable (in (L 2 (Ω)) 2 -norm and H(div, Ω)-norm). Further, we give two numerical examples in one-dimensional and two-dimensional spatial regions to verify the feasibility and effectiveness of the proposed scheme.
The remaining structure of the article is as follows. In Section 2, we present a fully discrete MFVE scheme for the time-fractional reaction-diffusion equation by introducing transfer operator γ h and the L1-formula of approximating the Caputo fractional derivative. In Section 3, the properties of transfer operator and truncation errors of L1-formula are given, and the generalized MFVE projection and its estimates are introduced. The existence, uniqueness, stability and convergence analysis for the MFVE scheme are analyzed in Sections 4 and 5, respectively. Two numerical examples are given to verify the feasibility and effectiveness in Section 6. Throughout the article, we use the standard notations and definitions of the Sobolev spaces as in Reference [54]. The mark C will denote a generic positive constant which is independent of the spatial mesh parameter h and time discretization parameter τ.

Fully Discrete MFVE Scheme
In order to formulate the MFVE approximate scheme, we introduce an auxiliary variable σ(X, t) = −∇u(X, t), and rewrite the problem (1) as The associated mixed variational formulation of (3) is to find {u, σ}: where H(div, Ω) = {w ∈ (L 2 (Ω)) 2 : divw ∈ L 2 (Ω)}. Now, we choose a quasi-uniform triangulation mesh T h = {K B } of the region Ω, where K B represents a triangular element with B as the centre of gravity, referring to Figure 1. Let h K B be the diameter of triangle K B , and h = max{h K B }. The nodes are defined as the midpoint of the three sides of the triangle. P 1 , P 2 , · · · , P M 0 represent the inner nodes and P M 0 +1 , · · · , P M S represent the boundary nodes.
We select the RT 0 space H h as the trial function space for variable σ, where and select L h as the trial function space for variable u, where Next, we construct the dual partition T * h , which is a union of quadrilaterals and triangles. For the interior node such as P 3 , referring to Figure 1, the dual element K * P 3 is the quadrilateral A 1 B 1 A 3 B 2 , and contains triangles K L (with ∆A 1 B 1 A 3 ) and K R (with ∆A 1 A 3 B 2 ). For the boundary node such as P 6 , the dual element K * P 6 is the triangle K I (with ∆A 4 A 5 B 3 ). Integrate (3) on the primal and dual partitions to obtain Similar to Reference [37], the transfer operator γ h : H h → (L 2 (Ω)) 2 is defined by where χ K represents the characteristic function of a set K. We select the range Y h of the transfer operator γ h as the test function space. Then, we apply the transfer operator γ h and rewrite (7) to obtain the following form Making use of the Green theorem, we have where n represents the unit out-normal direction. We can easily obtain that Then, we rewrite (9) to the following form Now, let 0 = t 0 < t 1 < · · · < t N = T be a equidistant grid of the time interval [0, T] with step length τ = T/N, for a positive integer N, and denote t n = nτ, n = 0, 1, · · · , N. For a given function ϕ . Following References [52,53], we can approximate the time-fractional derivative ∂ α u(X, t) ∂t α at t = t n as follows Denote Let U n and Z n be the approximate solutions of u and σ at t = t n , respectively. Then we get the fully-discrete MFVE scheme to find {U n , Z n } ∈ L h × H h , n = 1, 2, · · · , N, such that Remark 1.
(1) In practical numerical calculations, we need to apply the definition of D α t ϕ n to rewrite (13)(b) into the following or other equations The above equation is also used when the existence and uniqueness of discrete solutions are proved in Section 4.
(2) Compared with the standard FE methods [17][18][19][20] and FVE methods [34,35], the MFVE methods can simultaneously compute the different physical quantities about the scalar unknown variable u (such as displacement) and the vector-valued auxiliary variable σ (such as flux).
(3) In this article, we consider the Caputo fractional reaction-diffusion equations. Actually, we can apply the conversion relationship between the Caputo derivative and the Riemann-Liouville derivative to approximate the Riemann-Liouville derivative.

Some Lemmas and Notations
From (11) and (12), we can know that the truncation errors R n 1 (x), R n 2 (x) and R n t (x) will be estimated by the following lemma.

Lemma 1 ([52] (Inequalities (3.2) and (3.3))). The truncation errors R n
1 (x), R n 2 (x) and R n t (x) defined by (12) are bounded by where C > 0 is a generic constant which is independent of space-time grid step h and τ.
with ζ > 0, then there exists a constant C > 0 independent of space-time grid step h and τ such that Next, we give the properties of the transfer operator γ h . For ∀w h = (w 1 h , w 2 h ) ∈ H h , the discrete seminorm and the norm are defined as follows

Lemma 3 ([37] (Lemma 2.4)). The transfer operator γ h is bounded
Moreover, there exists a constant C > 0 independent of h such that

Lemma 4 ([39] (Lemmas 3.3 and 3.4)
). The following symmetry relation holds and there exists a constant µ 0 > 0 independent of h such that Proof. Noting that − ∑ n−1 k=0 b n k = 1, and applying Lemma 4, we obtain Applying Lemma 4, we can have the following equation Substituting the above equation into (15), we obtain the desired result.

Existence, Uniqueness and Stability Analysis for the MFVE Scheme
We first give the detailed proof of the existence and uniqueness for the MFVE scheme (13).
There exists a unique solution for the MFVE scheme (13).
be the basis functions of the space H h and L h , then Z n ∈ H h and U n ∈ L h are expressed as follows Substituting the above expressions into (13), and taking w h = φ i (i = 1, 2, · · · , M S ) and v h = ψ j (j = 1, 2, · · · , M B ), then (13) can be expressed as a matrix form: find {z n , u n } such that where z n = (z n 1 , z n 2 , · · · , z n M S ) T , u n = (u n 1 , u n 2 , · · · , u n M B ) T , It should be noted that A 1 , A 2 and A 3 are symmetric positive definite matrices. Noting that where E is the identity matrix, G = 1 Γ(2−α) A 2 + τ α A 3 + ετ α C T A −1 1 C, we have that G is a symmetric positive definite matrix and the coefficient matrix of linear equations (21) is invertible. Thus, linear Equations (21) have a unique solution, that is to say the MFVE scheme (13) has a unique solution. The proof of Theorem 1 has been completed.
Next, we consider the stability for the MFVE scheme (13). Theorem 2. Let {U n , Z n } be the solution of the MFVE scheme (13), then there exists a constant C > 0 independent of h and τ such that

Proof.
Choosing v h = U n and w h = Z n in (13), we can obtain By virtue of the definition of D α t U n , we derive Substituting (24) into (23), and applying Lemma 4, we have The above inequality leads to the following result Applying Lemma 2, we obtain Now, making use of (13)(a) and (14) Choosing w h = Z n in (28) and Apply Lemma 5 in (29) and note that Applying the Cauchy-Schwarz inequality and the Young inequality in (30), we have Then we have Multiply (31) by Applying Lemma 2 and (27), we derive Thus, we conclude the desired result.

Convergence Analysis for the MFVE Scheme
This section mainly studies the problem of convergence analysis for the MFVE scheme (13). First, let {ũ h ,σ h } be the generalized MFVE projection of {u, σ} defined by (17), thus the errors are expressed as Making use of (4), (13) and (14), and applying the generalized MFVE projection, we obtain the error equations Theorem 3. Let {u, σ} and {U n , Z n } be the solutions of the system (4) and (13), respectively, then there exists a constant C > 0 independent of h and τ such that

Proof.
Choosing v h = ξ n and w h = θ n in (33), we have Noting the fact that we can rewrite (35) as follows Apply Lemma 4 in (37) to obtain Now, we estimate the four terms F i (i = 1, 2, 3, 4) on the right-hand side of (38). Noting that b n k < 0(0 ≤ k < n), and applying the Cauchy-Schwarz inequality, we can obtain Applying D α t η n = τ −α Γ(2−α) ∑ n−1 k=0b k ∂ t η n−k , ∑ n−1 k=0b k = n 1−α and Lemma 6, we obtain Thus, the F 3 can be estimated by the following result Then, combining the estimates of F i (i = 1, 2, 3, 4), we have Applying Lemma 2, we obtain To estimate σ n − Z n , we apply (33)(a) and (34) Choose v h = D α t ξ n in (33)(b) and w h = θ n in (47) to obtain Applying the Cauchy-Schwarz inequality and the Young inequality in (48), we can get With the help of Lemma 5 in (49), we have Multiply the above inequality by 2 Applying Lemma 2 and Lemma 4, we obtain Next, we estimate σ n − Z n H(div,Ω) . Choose v h = divθ n in (33)(b) and w h = θ n in (47) to obtain Then, (53) leads to the following result Apply the Cauchy-Schwarz inequality and the Young inequality in (53) to obtain Noting that (θ n , γ h θ n ) ≥ 0 and b n k < 0(0 ≤ k < n), we have Replace −(D α t θ n , γ h θ n ) with the above result in (54) to obtain Finally, apply Lemma 1 and Lemma 6 with (46), (52) and (57) to complete the proof.

Numerical Examples
For examining the feasibility and effectiveness of the MFVE scheme, we consider two numerical examples with one-dimensional and two-dimensional spatial regions. Example 1. We consider the following time-fractional reaction-diffusion equation in one-dimensional spatial regions where Ω = (a, b) and J = (0, T]. We also introduce an auxiliary variable σ(x, t) = −u x (x, t), and rewrite the Equation (58) as the first-order coupled system. As in Reference [41], we establish the primal partition T h and dual partition T * h , and choose the space H * h × L * h as the trial function space, where As in Reference [41], we also define transfer operator γ * h : where A * i ∈ T * h and M T is the number of spatial nodes.
The numerical simulation results with parameter α = 0.1, 0.3, 0.5, 0.7, 0.9 are given in Tables 1-6. We fix the spatial step length h = 1/10,000, select the time step length τ = 1/10, 1/20, 1/40, 1/80, and give the error results of u (in L ∞ (L 2 (Ω))-norm) and σ (in L ∞ (L 2 (Ω))-norm and L ∞ (H 1 (Ω))-norm) in Tables 1-5. It is easy to see that the orders of time convergence are approximate to 2 − α which are consistent with the theoretical results in Theorem 3. Moreover, taking different parameter α, fixing the time step length τ = 1/1000, and selecting the spatial step length h = 1/10, 1/20, 1/40, 1/80, we also do some numerical experiments to examine the orders of spatial convergence. The numerical results show that the orders of spatial convergence are approximate to 2 which are larger than the theoretical results. The numerical simulation results with parameter α = 0.1 are given in Table 6, and the other results with different parameter α are similar, and we will not repeat the description.    Table 3. Error results with α = 0.5 and h = 1 × 10 −4 for Example 1.  Table 4. Error results with α = 0.7 and h = 1 × 10 −4 for Example 1.  Table 5. Error results with α = 0.9 and h = 1 × 10 −4 for Example 1.  Table 6. Error results with α = 0.1 and τ = 1 × 10 −3 for Example 1.
In this example, we do some numerical experiments with different parameter α and mesh size h = √ 2τ, and apply the fifth-order Gauss integral formula in triangular regions to calculate the errors of u (in L ∞ (L 2 (Ω))-norm) and σ (in L ∞ (L 2 (Ω)-norm and L ∞ (H(div, Ω))-norm, where L 2 (Ω) = (L 2 (Ω)) 2 ). The numerical results show that the orders of convergence are approximate to 1. The numerical simulation results are given in Table 7  are given in Figures 2 and 3, respectively. Because the trial function spaces L h (defined in (6)) and H h (defined in (5)) are not continuous function spaces, the surfaces in Figures 2 and 3 agree with the properties of the trial function spaces. The figures and numerical results show that the proposed MFVE method for the time-fractional reaction-diffusion equations in two-dimensional spatial regions is feasible and effective.       Table 9. Error results with α = 0.3 and h = 1×10 −4 for Example 1 (ε = 1/100).

Conclusions
We apply the MFVE method to treat the time-fractional reaction-diffusion equations. By introducing the auxiliary variable σ and applying the L1-formula, we construct the MFVE scheme, prove the existence and uniqueness of the proposed scheme, and obtain the stability results which only depend on the source term function f (X, t). We also obtain a priori error estimates for the variable u (in L 2 (Ω)-norm) and the auxiliary variable σ (in (L 2 (Ω)) 2 -norm and H(div, Ω)-norm) by using generalized MFVE projection and the properties of the transfer operator γ h . Moreover, we briefly describe that the MFVE method can also be constructed in one-dimensional spatial regions, and give two numerical examples to examine the feasibility and effectiveness. In this article, we consider a linear model. In the future, we will try to study some nonlinear fractional models by using our method, which will be a challenge for us from the point of view of numerical analysis.
Author Contributions: All authors contributed to the draft of the manuscript, all authors read and approved the final manuscript.