A Crank–Nicolson Finite Volume Element Method for Time Fractional Sobolev Equations on Triangular Grids

In this paper, a finite volume element (FVE) method is proposed for the time fractional Sobolev equations with the Caputo time fractional derivative. Based on the L1-formula and the Crank–Nicolson scheme, a fully discrete Crank–Nicolson FVE scheme is established by using an interpolation operator I∗ h . The unconditional stability result and the optimal a priori error estimate in the L2(Ω)-norm for the Crank–Nicolson FVE scheme are obtained by using the direct recursive method. Finally, some numerical results are given to verify the time and space convergence accuracy, and to examine the feasibility and effectiveness for the proposed scheme.


Introduction
Fractional partial differential equations (FPDEs) have received extensive attentions by more and more scholars, and have been widely applied in many fields of science and engineering [1,2]. Many practical problems can be portrayed very well by the some FPDEs, such as fractional (reaction) diffusion equations [3][4][5][6][7][8][9][10][11][12], fractional Allen-Cahn equations [13][14][15], fractional Cable equations [16][17][18], and fractional mobile/immobile transport equations [19][20][21]. In the past few decades, a large number of numerical methods [22][23][24] have been proposed and used to solve the FPDEs, which have achieved excellent theoretical and numerical results. Recently, some scholars began to propose some numerical methods to solve fractional Sobolev equations. Liu et al. [25] constructed a modified reduced-order finite element scheme to solve the time fractional Sobolev equations by using the proper orthogonal decomposition technique [26,27], and obtained the stability and convergence results. Haq and Hussain [28] proposed a meshfree spectral method for the time fractional Sobolev equations by using radial basis functions and point interpolation method. Beshtokov [29] proposed a finite difference method to solve the time fractional Sobolev-type equations with memory effect.
In this paper, we study a Crank-Nicolson finite volume element method to solve the following time fractional Sobolev equations with the initial and boundary conditions u(x, t) = 0, (x, t) ∈ ∂Ω ×J, u(x, 0) = u 0 (x), where J = (0, T] with a positive constant T. Here, we assume that Ω ⊂ R 2 is a bounded convex polygonal region with boundary ∂Ω, the coefficients κ 1 and κ 2 are two positive constants. Moreover, we assume that the initial data u 0 (x) and source function f (x, t) are smooth enough. In (1), ∂ α /∂t α is the Caputo time fractional derivative with order α (0 < α < 1) defined by ∂φ(x, s) ∂s ds (t − s) α , for a given function φ(x, t). ( In recent years, the finite volume element (FVE) methods [30][31][32][33][34] have been applied by more and more scholars to solve some FPDEs numerically. Sayevand and Arjang [35] proposed an FVE method for the sub-diffusion equation with the Caputo fractional derivative. Karaa et al. [36] adopted a piecewise linear discontinuous Galerkin method in time, and constructed an FVE scheme for the sub-diffusion equation with the Riemann-Liouville fractional derivative. Karaa and Pani [37] constructed two fully discrete FVE schemes to solve the time fractional sub-diffusion equation with smooth and nonsmooth initial datas, in which the Riemann-Liouville fractional derivative was approximated by using convolution quadrature generated by two different difference schemes. Badr et al. [38] proposed a B-spline FVE method for the time fractional advection-diffusion problem with the Caputo fractional derivative, and gave the stability analysis. Zhao et al. [39] designed an fully discrete FVE scheme to solve the nonlinear time fractional mobile/immobile transport equations based on the weighted and shifted Grünwald difference (WSGD) formula, and obtained the unconditional stability and optimal error estimates.
In this paper, the main aim was to establish a fully discrete FVE scheme to solve the time fractional Sobolev equation by using the Crank-Nicolson scheme. In spatial discretization, we use the classical FVE method, which can maintain the local conservation of some physical quantities, and this is very important in numerical computing. In temporal discretization, we apply the Crank-Nicolson scheme to discretize the equation at the time node t n+1/2 on the whole, and combine the L1-formula [40,41] to discretize the Caputo time fractional derivative ∂ α ∆u/∂t α , so that the time convergence accuracy is O(τ 2−α ). In our theoretical analysis, we specially adopt the recursive analysis method, and obtain the unconditional stability result and the optimal a priori error estimate in the L 2 (Ω)-norm. Moreover, we provide some numerical results to verify the theoretical result, and to examine the feasibility and effectiveness of the fully discrete FVE scheme.
The remainder part of this paper is organized as follows. In Section 2, a Crank-Nicolson FVE scheme for the time fractional Sobolev Equation (1) is proposed. In Section 3, the unconditional stability analysis for the Crank-Nicolson FVE scheme is derived in detail. The a prioir error estimate is given in Section 4. Finally, in Section 5, some numerical results are given to illustrate the feasibility and effectiveness. Furthermore, throughout the article, we adopt the mark C to denote a generic positive constant, which is independent of the mesh parameters.
where E α,ϕ (t n+ 1 2 ) is the truncation error. Now, we apply (3) and (7) to rewrite the original equation in (1) at t = t n+ 1 2 as the following equivalent equation∂ where the truncation error E α (t n+ 1 2 ) is as follows: Following References [23,41], and making use of Taylor's expansion, we can obtain the estimate of the truncation error E α (t n+ 1 2 ), that is, if u ∈ C 2 (J, H 2 (Ω)) and u ttt ∈ L ∞ (J, L 2 (Ω)), then there exist a positive constant C independent of h and τ such that E α (t n+ 1 2 ) ≤ C(τ 2 + τ 2−α ). Making use of (8), we can obtain the variational formulation at t = t n+ 1 2 as follows where a(u, v) = Ω ∇u · ∇vdx, ∀u, v ∈ H 1 0 (Ω). Now, as shown in Figure 1, let T h = {K} be a quasi-uniform triangular partition of the domain Ω with h = max{h K }, where h K is the diameter of the triangular element K ∈ T h . Then Ω = ∪ K∈T h K and Z h denotes all vertices, that is  Next, based on the primal partition T h , we construct a dual partition T * h . With z 0 ∈ Z 0 h as an interior node, let z j (j = 1, 2, · · · , s) be its adjacent nodes (as shown in Figure 1, s = 6). Let M j (j = 1, 2, · · · , s) denote the midpoints of z 0 z j and let Q j (j = 1, 2, · · · , s) be the barycenters of the triangle z 0 z j z j+1 with z n+1 = z 1 . We construct the control volume K * z 0 by joining successively M 1 , Q 1 , · · · , M s , Q s , M 1 . With Q j , j = 1, 2, · · · , s as the nodes of control volume K * z j , let Z * h be the set of all dual nodes. Then, we construct the dual partition T * h through the union of all control volumes K * z j . Then, we define the trial function space V h and test function space V * h as follows: With a node z, let Φ z be the standard nodal linear basis function, and Ψ z be the characteristic function of the control volume K * z . It is obvious that We define the piecewise linear interpolation operator I h : C(Ω) → V h and the piecewise constant interpolation operator I * h : C(Ω) → V * h as follows: From Reference [30], we can see that I h and I * h have the following approximation property Now, integrating (8) over each control volume K * z , and applying the Green formula, we can get where n means the outer-normal direction on ∂K * z . Let u h ∈ V h be the discrete solution of u, make use of the operator I * h to rewrite (13) as the following variational formulation: where the bilinear form a(·, ·), following References [30,31], can be rewritten as follows: Let u n h be the discrete solution of u at t = t n . We establish a fully discrete Crank-Nicolson FVE scheme to seek u n h ∈ V h , (n = 1, · · · , N) such that And we can rewrite the Crank-Nicolson FVE scheme (16) as the following equivalent formulation: Case n = 0: Case n ≥ 1: Remark 1. In the construction of the dual partition T * h , we can also choose Q j as the circumcenter of the triangular element, as shown in Figure 3.2.2 in Reference [30]. Based on this circumcenter dual partition, we can also construct the Crank-Nicolson FVE scheme, and obtain the same theoretical analysis.

Remark 2.
When the barycenter dual partition is selected, there is no essential difference between the structured and unstructured primal partition in the proposed Crank-Nicolson FVE scheme. However, when the circumcenter dual partition is selected and the space domain Ω is rectangular, based on the structured primal partition, the proposed scheme will become more simple, and is similar to the finite difference scheme.

Remark 3.
Making use of Lemmas 1 and 2 in Section 3, we can easily have that the coefficient matrices of linear equations generated by (17) and (18) are invertible. Then, there exists a unique discrete solution for the Crank-Nicolson FVE scheme (16).

Stability Analysis
In order to give the stability analysis, we first introduce some properties of the bilinear forms (·, I * h ·) and a(·, I * h ·). Lemma 1. [30] The following symmetry relation for the bilinear form (·, I * h ·) holds: Moreover, there exist two positive constants µ 1 and µ 2 independent of h such that Lemma 2. [30] The following symmetry relation for the bilinear form a(·, I * h ·) holds: Moreover, there exist two positive constants µ 3 and µ 4 independent of h such that Next, we give the unconditional stability result for the Crank-Nicolson FVE scheme (16).
be the solutions of the Crank-Nicolson FVE scheme (16), then there exists a positive constant C independent of h and τ such that Proof. For the case of n = 0, choosing v h = 2u (17), we obtain Applying Lemma 2, we have Then, making use of (26) in (25), we obtain Applying Lemmas 1 and 2, we obtain For the case of n ≥ 1, choosing v h = 2u Noting that and substituting (30) into (29), we obtain Applying Lemma 2, we have and (32) and (33) into (31), we have Apply Lemma 2 and the Young inequality to obtain Now, let Λ n = (u n h , We can rewrite (35) as follows: By using the direct recursive method in (36), we have Making use of (28), we estimate Λ 1 as follows: Substituting (38) into (37), we have Thus, apply the definition of Λ n and Lemma 1 to complete the proof.

Convergence Analysis
In order to obtain the error estimates for the Crank-Nicolson FVE scheme (16), we need to introduce a projection operator P h : H 1 0 (Ω) ∩ H 2 (Ω) → V h , which is defined by the following formulation: Following Reference [30], we give the following estimates for the projection operator P h .
Next, we give the optimal a priori error estimate for the Crank-Nicolson FVE scheme (16).

Numerical Examples
In this section, we will give some numerical results to examine the feasibility and effectiveness of the proposed Crank-Nicolson FVE scheme. Example 1. We consider the following time fractional Sobolev equation: where κ 1 = κ 2 = 1, Ω = (0, 1) × (0, 1) and J = (0, 1]. We choose the source function f (x, t) as follows: Thus, we can obtain the corresponding exact solution In the actual numerical calculation, in order to reduce the influence of numerical integration on the calculation accuracy as much as possible, we use the composite fifth-order Gauss quadrature formula in triangular domain to calculate some definite integral of space variable, including the calculation of the error u(t n ) − u n h . In Table 1, for the fixed space step h = √ 2/160, we choose different fractional parameters α = 0.1, 0.3, 0.5, 0.7, 0.9 and time step τ = 1/10, 1/20, 1/40, and get the corresponding error results. From these results, it is easy to see that the time convergence rates in L 2 (Ω)-norm are close to 2 − α on the whole, which is a little higher when α = 0.1, 0.3 and a little lower when α = 0.5, 0.7, 0.9. In Table 2, we select the same fractional parameters α and different meshes with h = √ 2/10, √ 2/20, √ 2/40 and time step τ = 1/1000, and get the space convergence rates, which are close to 2. These numerical results are consistent with the theoretical analysis.    Table 2. Error max n u(t n ) − u n h with τ = 1 1000 for Example 1. Furthermore, we give the error results at different time points t = 0, 0.1, 0.2, · · · , 1. In Tables 3 and 4, we take the fractional parameter α = 0.1, 0.9, respectively, fix the time step τ = 1/1000, and list the error results in L 2 (Ω)-norm at each time point with different space step h = √ 2/10, √ 2/20, √ 2/40. From these error results, it is easy to see that the errors in L 2 (Ω)-norm are gradually increasing as time goes on. For the fractional parameter α = 0.3, 0.5, 0.7, we also calculate corresponding error results, which are similar to the numerical behaviors in the cases α = 0.1, 0.9, so we will not repeat them.   Table 4. Errors in L 2 -norm with τ = 1 1000 and α = 0.9 for Example 1. Example 2. We consider the time domainJ, space domainΩ, initial function u(x, 0), coefficients κ 1 and κ 2 as in Example 1. In this example, we choose the source function f (x, t) as follows: Then, the exact solution in this example is as follows: u(x, t) = t 2 sin(2πx 1 ) sin(2πx 2 ), ∀x = (x 1 , x 2 ) ∈Ω.
In this example, we also choose some different fractional parameters α and mesh sizes to carry out numerical experiments, and obtain the corresponding error results in Tables 5-8. In Table 5, we fix the space step h = √ 2/200, take different time steps and fractional parameters as in Example 1, obtain the corresponding error results, and the time convergence rates in L 2 (Ω)-norm are close to 2 − α on the whole, which is a little higher when α = 0.3, 0.5 and a little lower when α = 0.1, 0.7, 0.9. In Table 6, we take the mesh sizes and fractional parameters as in Example 1, and find that the space convergence rates in L 2 (Ω)-norm are also close to 2, which is consistent with the theoretical analysis. In Tables 7 and 8, we give the error results at different points t = 0, 0.1, 0.2, · · · , 1 with fractional parameters α = 0.3, 0.7, which have the same error behaviors as Example 1. Based on the above error results and analysis, we can easily see that the constructed Crank-Nicolson FVE scheme for the time fractional Sobolev equations is feasible and effective.    3.69351278 × 10 −2 9.22110708 × 10 −3 2.29129165 × 10 −3

Remark 4.
Based on the L1-formula, we can also construct a backward Euler FVE scheme, which is to find u n h ∈ V h , (n = 1, · · · , N) such that Next, we compare the numerical results of the Crank-Nicolson FVE scheme (16) and backward Euler FVE scheme (63), and also consider Examples 1 and 2. In Table 9, we select the same mesh sizes as in Table 1 of Example 1, and obtain the error results, in which the time convergence rates of the backward Euler FVE scheme are significantly lower than 2 − α when α = 0.1, 0.3. As can be seen from Tables 5 and 10, there are similar error behaviors in the calculation of the two schemes in Example 2. It can be seen from the comparison results that the Crank-Nicolson FVE scheme is better than the backward FVE scheme.

Conclusions
In this paper, we study the FVE methods based on the Crank-Nicolson scheme and L1-formula to solve the time fractional Sobolev equations with the Caputo fractional derivative. We derive the unconditional stability result which is only dependent on the source term function f (x, t) and the initial data u 0 (x). We also obtain the optimal a priori error estimates in L 2 (Ω)-norm by using the properties of the operator I * h and the direct recursive method. Moreover, we give two numerical examples with some error results to examine the accuracy and efficiency of the proposed Crank-Nicolson FVE scheme. In this paper, we consider the Dirichlet boundary condition. In the future, we will try to use the FVE method to solve more fractional models with different boundary conditions.