Finite volume element method with the WSGD formula for nonlinear fractional mobile/immobile transport equations

In this article, a finite volume element method with the second-order weighted and shifted Grünwald difference (WSGD) formula is proposed and studied for nonlinear time fractional mobile/immobile transport equations on triangular grids. By using the WSGD formula of approximating the Riemann–Liouville fractional derivative and an interpolation operator I∗ h , a second-order fully discrete finite volume element (FVE) scheme is formulated. The existence, uniqueness, and unconditional stability for the fully discrete FVE scheme are derived, the optimal a priori error estimates in L∞(L2(Ω )) and L2(H1(Ω )) norms are obtained, in which the convergence orders are independent of the fractional parameters. At the end of this article, two numerical examples with different nonlinear terms are given to verify the feasibility and effectiveness.


Introduction
The theories and numerical methods of fractional differential equations (FDEs) have become hot topics, which attract more and more scholars in science and engineering. This is because the FDEs are widely used in many fields, such as physics, chemistry, biology, and ecology [1][2][3][4][5][6][7]. Many practical problems with some properties such as memory, heterogeneity, or heredity can be described well by the corresponding FDEs, including fractional cable equations, fractional diffusion equations, fractional Allen-Cahn equations, and so on. However, due to the existence and complexity of fractional derivative, it is difficult to obtain the analytical solutions for most of the FDEs. Thus, many numerical methods [8][9][10][11][12][13][14][15][16][17][18][19] have been proposed and studied to solve different types of the FDEs.
In this article, we focus on the following nonlinear time fractional mobile/immobile transport equations for (x, t) ∈ Ω × J, with boundary and initial conditions u(x, t) = 0, (x, t) ∈ ∂Ω ×J, u(x, 0) = u 0 (x), x ∈Ω, (2) where Ω ⊂ R 2 is a bounded convex polygonal domain with boundary ∂Ω, J = (0, T] with 0 < T < ∞, β 1 > 0 and β 2 ≥ 0 are two given nonnegative constants. The source function f (x, t) and initial data u 0 (x) are smooth enough. Moreover, we assume that the coefficient A(x) = {a i,j (x)} 2×2 is a sufficiently smooth matrix function, which is symmetric and uniformly positive definite, that is, there exists a constant β 0 > 0 such that The nonlinear term g(u) satisfies |g (u)| ≤ C, where C > 0 is a constant. Also ∂ α u(x,t) ∂t α is the Riemann-Liouville time fractional derivative defined by The mobile/immobile transport equations have been widely used in unsaturated transport through homogeneous media [20][21][22][23][24]. Due to some partitioning between the phases [25], mobile/immobile formulations equate the divergence of advective and dispersive flux of a mobile phase to the change in concentration in both the mobile and immobile zones. Fractional mobile/immobile transport equations, as described in [26] by Schumer et al., are equivalent to the mobile/immobile model with power law memory function, and are considered to be the limiting equation which governs continuous time random walks with heavy tailed random waiting times. Recently, some numerical methods have been designed to solve the fractional mobile/immobile models. Liu et al. [27] proposed an implicit finite difference (FD) method for the fractional mobile/immobile advectiondispersion equation, in which the Caputo time fractional derivative is discretized by L1 formula [12,13] and the Riemann-Liouville space fractional derivative is discretized by shifted Grünwald-Letnikov formula [1,28]. Zhang et al. [29] provided a stable implicit Euler approximation scheme to solve the fractional mobile/immobile advection-dispersion equation with the Coimbra variable-order derivative. Liu et al. [30] proposed a meshless method to treat the two-dimensional fractional mobile/immobile transport equation based on radial basis functions for the spatial discretization. Wang [31] constructed a highorder compact FD scheme to solve the fractional mobile/immobile convection-diffusion equations, and gave a Richardson extrapolation algorithm to improve the temporal convergence accuracy. Yin et al. [32] constructed a generalized BDF2-θ scheme for the fractional mobile/immobile transport equations with the initial singularity of the time fractional derivative. However, we find that there is no report about the finite volume element method with the weighted and shifted Grünwald-Letnikov difference (WSGD) formula to solve the nonlinear time fractional mobile/immobile transport equations.
The finite volume element (FVE) method, as an important numerical method, has been widely used to solve various differential equations [33][34][35][36][37][38][39] in the field of science and engineering. This method can preserve the local conservation laws for some physical quantities, which is very important in scientific computing. Recently, the FVE method has been used to solve some FDEs by some scholars. Sayevand and Arjang [40] designed an FVE scheme to solve the time fractional subdiffusion equation on a rectangular partition, where the Caputo fractional derivative was approximated by using L1 formula. Karaa et al. [41] constructed an FVE scheme for the time fractional subdiffusion equation with the Riemann-Liouville fractional derivative, and applied a piecewise linear discontinuous Galerkin method in time, where the convergence rate in time was k 1+α (0 < α < 1). In order to improve the results in [41], Karaa and Pani [42] considered smooth and nonsmooth initial data, and gave two fully discrete numerical schemes by using convolution quadrature in time generated by the backward Euler difference method and the second-order backward difference method. Recently, numerical methods based on the WSGD formula have attracted more and more scholars, and have been studied to solve many FDEs [43][44][45][46][47][48][49]. Compared with the L1 formula, the WSGD formula can obtain the second-order convergence rate, which is independent of the fractional parameters. This motivates us to find a way to combine the FVE methods with the WSGD formula so that we can use their advantages to solve more FDEs.
In this article, our purpose is to construct an FVE scheme for the nonlinear time fractional mobile/immobile transport equations on triangular girds by using the WSGD formula. In spatial discretization, we construct the primal and dual partitions, select the piecewise linear polynomial space and the piecewise constant function space as the trial and test function spaces, respectively, then construct the FVE scheme by using the interpolation operator. In temporal discretization, we adopt the second-order WSGD formula to approximate the Riemann-Liouville fractional derivative ∂ α u/∂t α , apply a second-order three-level difference scheme to approximate the time derivative ∂u/∂t, and give a secondorder approximation formula for the nonlinear term g(u). We give the existence, uniqueness, and unconditional stability analyses for the FVE scheme in detail, and obtain the optimal a priori error estimates in L ∞ (L 2 (Ω)) and L 2 (H 1 (Ω)) norms. Compared with the discrete schemes in [27,29,30], our scheme can achieve second-order temporal convergence rate.
This article is organized as follows. In Sect. 2, a fully discrete FVE scheme for the nonlinear time fractional mobile/immobile transport equation (1)-(2) is proposed. In Sect. 3, the existence, uniqueness, and unconditional stability analyses are derived. In Sect. 4, the optimal a prioir error estimates are obtained. Finally, in Sect. 5, two examples with different nonlinear terms are given to illustrate the feasibility and effectiveness. Furthermore, we use general definitions and notations of the Sobolev spaces as in [50], and adopt the symbol C to represent a generic positive constant, which is independent of temporal and spatial mesh.

Fully discrete finite volume element scheme
In order to construct the FVE scheme, we first give primal and dual partitions. Let T h = {K} be a set of quasiuniform triangulation mesh of the domain Ω with h = max{h K }, referring to Fig. 1, where h K denote the diameter of the triangle K ∈ T h . Then we have Ω = K∈T h K . Moreover, let Z h = {Z : Z is a vertex of element K, K ∈ T h } represent all vertices, and Z 0 h ⊂ Z h represent the set of interior vertices in T h .
Next, let T * h be the dual mesh based on the primary mesh T h . With Z 0 ∈ Z 0 h as an interior node, let Z i (i = 1, 2, . . . , m) be the corresponding adjacent nodes (as shown in Fig. 1, m = 6). We denote the midpoints of Z 0 Z i by M i (i = 1, 2, . . . , m), and denote the barycenters of the and define the piecewise constant function space V * h as the test function space, that is, Let Φ Z be the general nodal linear basis function associated with the node Z ∈ Z 0 h , and Ψ z be the characteristic function of the control volume K * be the classical piecewise linear interpolation operator and I * h : C(Ω) → V * h be the piece constant interpolation operator, that is, where n denotes the outer-normal direction on ∂K * Z . We apply the operator I * h to rewrite (5) as the following formulation: where the bilinear form a(·, ·), following [33,34], can be taken as follows: In our theoretical analysis, we also need to give the variational formulation of the problem (1)-(2) to find u(t) ∈ H 1 0 (Ω) such that where a(u, v) = Ω A∇u · ∇v dx, ∀u, v ∈ H 1 0 (Ω). Now, we introduce the mesh of the temporal intervalJ = [0, T] given by 0 = t 0 < t 1 < · · · < t N = T for some positive integer N , where t n = nτ and τ = T/N , n = 0, 1, . . . , N . We denote ϕ n = ϕ(t n ) for a function ϕ and Next, following [43,44], we apply the WSGD formula to approximate the Riemann-Liouville fractional derivative ∂ α u(x,t) ∂t α at time t = t n+1 as follows: where the truncation error E n+1 u,α = O(τ 2 ), and For approximating the nonlinear term g(u) at t = t n+1 , as in [51], we use the linearized formulation denoted by G[u n+1 ] as follows: Then, we can obtain the equivalent formulation of variational formulation (8) at t = t n+1 which is to find u n+1 = u(t n+1 ) ∈ H 1 0 (Ω) such that where and Now, we denote by U n the approximate solution of u at t = t n , make use of (9), (10), and (13), and obtain the fully discrete FVE scheme to find U n+1 ∈ V h (n = 0, 1, . . . , N -1) such that We can also split the FVE scheme (17) into the following equivalent iterative formulation: Case n = 0: Case n ≥ 1: Then we will prove the existence and uniqueness of the fully discrete solutions for the FVE scheme (17) (or (18)- (19)) in the next section.
Remark 2.1 For the fully discrete FVE scheme (17) or equivalent formulations (18)- (19), in the practical calculations, we can obtain U 1 by using U 0 and solving (18), where U 0 = P h u 0 defined in Sect. 4. Thus, when U n and U n-1 (n ≥ 1) have been obtained, we can solve (19) to obtain U n+1 .

Existence, uniqueness, and stability analyses
In this section, we will derive the results of the existence, uniqueness, and stability for the fully discrete FVE scheme (17). First, we give some properties of the coefficients in WSGD formula and the interpolation operator I * h in Sect. 3.1.

Some lemmas
The interpolation operators I h and I * h satisfy the following properties: The bilinear form (·, I * h ·) satisfies the following property: and there exist two positive constants μ 1 and μ 2 independent of h such that Making use of Lemma 3.3, similar to the proof in [43,44], we can also obtain the following property of the sequence {q α (k)} ∞ k=1 .
Proof When n = 0, making use of Lemma 3.3 and the definition of ∂ 2 t φ n , we obtain the following result: For the case of n ≥ 1, we have Noting that we substitute (33) into (32), and obtain the desired result.

Existence and uniqueness
Theorem 3.1 There exists a unique discrete solution for the fully discrete FVE scheme (17).
Proof Let M 0 Z be the number of the vertices in Z 0 h , and {Φ i : i = 1, 2, . . . , M 0 Z } be the basis functions of the space V h , then U n ∈ V h can be expressed as follows: We substitute (34) into the formulations (18) and (19) equivalent to the FVE scheme (17), take v h = Φ j (j = 1, 2, . . . , M 0 Z ), then rewrite (18) and (19) in the matrix form, and search for u n such that whereũ n = ũ n 1 ,ũ n 2 , . . . ,ũ n It is easy to see that A 1 is a symmetric positive definite matrix. Let Next, we will prove B 1 and B 2 are invertible. Applying Lemma 3.4, for ∀Y = (y 1 , y 2 , . . . , is a positive definite quadratic form generated by the asymmetric matrix A 2 . Therefore, Y T B 1 Y and Y T B 2 Y (for Y ∈ R M 0 Z ) are positive definite quadratic forms generated by asymmetric matrices B 1 and B 2 , respectively. Then we have that B 1 and B 2 are invertible. Hence, the linear equations (35) have a unique solution, and so the FVE scheme (17) has a unique solution. Thus, we complete the proof of Theorem 3.1.

Stability
The fully discrete FVE scheme (17) satisfies the following unconditional stability results. (17), then there exists a constant C > 0 independent of h and τ such that

Theorem 3.2 Let {U n } N n=1 be the solutions of the fully discrete FVE scheme
Proof Taking v h = U n+1 in (17) yields the following result: Making use of Lemma 3.4, and applying the Cauchy-Schwarz and Young inequalities, we have For the term G[U n+1 ] 2 in (37), when n ≥ 1, we have G U n+1 2 ≤ 2g U ng U n-1 2 ≤ C U n 2 + U n-1 2 .
When n = 0, we have Next, for the case of n ≥ 1 in (37), we apply Lemma 3.6 to obtain Multiplying (40) by 4τ , and summing over n from 1 to m, we have For the case of n = 0 in (37), taking v h = U 1 in the FVE system (17), we have Multiplying (42) by 2τ , and applying Lemma 3.4, we obtain Substitute (39) into (43) to obtain We rewrite (44) as follows: Use Lemma 3.5 to obtain Now, making use of (41) and (45), we have Applying Lemma 3.6, we obtain Substituting (48) into (47), and making use of (46), we rewrite (47) as follows: Applying Lemma 3.5 and the discrete Gronwall lemma yields Thus, we apply Lemma 3.3 to complete the proof.
Remark 3.1 From the results (50), we can also see that the fully discrete solution is unconditional stable in the discrete L 2 (H 1 (Ω)) norm, that is, In the next section, we also give the optimal a priori error estimate in this discrete L 2 (H 1 (Ω)) norm.

A priori error analysis
In order to give the error analysis for the fully discrete FVE scheme (17), we need to introduce an elliptic projection operator P h : Following [33], the above projection operator P h satisfies the following estimates.

Lemma 4.1 There exists a positive constant C such that
Next, we give the main results in this paper about the error estimates. Moreover, we can also obtain the following error estimate: Proof We split the error u(t n ) -U n = (u(t n ) -P h u(t n )) + (P h u(t n ) -U n ) = ξ n + η n . Applying Lemma 4.1, we only need to estimate η n . It is easy to see that η n satisfies the following error equation: Choose v h = η n+1 in (54) to obtain Making use of Lemma 3.4, we can obtain Next, we estimate some terms on the right-hand side of the inequality (56). Applying Lemma 4.1 and the triangle inequality, for the case of n ≥ 1, we have For the case of n = 0, we have We discuss the term G[u n+1 ] -G[U n+1 ] 2 separately. For the case of n ≥ 1, we have and for the case of n = 0, we have From the definitions of E n+1 u,t , E n+1 u,α , and E n+1 g , we easily get E n+1 Now, for the case of n ≥ 1, making use of (57)-(61) in (56), and applying Lemma 3.6, we obtain Multiplying (62) by 4τ , and summing over n from 1 to m, we have For the case of n = 0, taking v h = η 1 in (54), and applying Lemma 3.6, we obtain Multiplying (64) by 2τ , we have Making use of (58), (60), and (61) in (65), we can obtain Thus, we rewrite (66) as follows: Noting that η 0 = 0, and applying Lemma 3.5, we have Next, adding (63) and (67), we have Making use of (30) and (68), we have Substituting (70) into (69), we rewrite (69) as follows: Noting that η 0 = 0, and applying Lemma 3.5 and the discrete Gronwall lemma, we obtain Finally, we apply Lemma 4.1 and the triangle inequality to complete the proof.

Numerical examples
In order to illustrate the feasibility and effectiveness for the proposed FVE scheme, we give two numerical examples with different nonlinear terms g(u) = sin(u) and g(u) = u 3u and different exact solutions.
Then we can get the corresponding source function f (x, t).
In this example, we choose some different mesh sizes and parameters α to conduct numerical experiments, and give the error results in L ∞ (L 2 (Ω)) and L 2 (H 1 (Ω)) norms for u(x, t), in which we use the Gauss integral formula with fifth-order algebraic accuracy to calculate the space norms of the errors. In Table 1, we give the corresponding error results with different parameters α = 0.01, 0.5, 0.99 and mesh sizes (h, τ ) = ( . We point out here that error behaviors with other different parameters α such as α = 0.1, 0.3, 0.7, 0.9 are similar, so we will not repeat them. From the error results, we can easily see that the convergence order for u in L ∞ (L 2 (Ω)) norm is approximately equal to 2, and the convergence order for u(x, t) in L 2 (H 1 (Ω)) norm is approximately equal to 1. In order to observe the approximation effect intuitively, we choose the fractional parameter α = 0.5 in (73), and give the graphs of the exact solution and the numerical solution at time t = 1 in Fig. 2 20 , 1 20 )), respectively. It is easy to see that the graph of the numerical solution is also consistent with that of the exact solution.
Example 2 We consider the space domainΩ, time domainJ, initial data u 0 (x) and coefficient A(x) as in Example 1. In this example, we choose the nonlinear term g(u) = u 3u and the exact solution is then u(x, t) = t 1+α sin(2πx 1 ) sin(2πx 2 ), ∀x = (x 1 , x 2 ) ∈Ω.  In this example, we also take some different mesh sizes and parameters α to conduct numerical experiments, and give the corresponding error results in L ∞ (L 2 (Ω)) and L 2 (H 1 (Ω)) norms for u(x, t) in Table 2 with parameters α = 0.01, 0.5, 0.99 and mesh sizes . We can see that the error behaviors are consistent with those in Example 1. In Figs. 3(a) and 3(b), we describe the graphs of the exact and numerical solutions with parameter α = 0.5 for u(x, t) at time t = 1, respectively, where the mesh sizes are same as in Example 1. We also find that the exact solution u is approximated well by the fully discrete FVE solution U. The numerical behaviors and figures show that the constructed FVE scheme with second order WSGD formula for the nonlinear time fractional mobile/immobile transport equations in two-dimensional spatial regions is feasible and effective.

Conclusions
We apply the FVE methods based on the second-order WSGD formula to treat the nonlinear time fractional mobile/immobile equations with the Riemann-Liouville time fractional derivative. We construct the second-order fully discrete FVE scheme, give the existence, uniqueness, and unconditional stability results, derive the optimal a priori error estimates in L ∞ (L 2 (Ω)) and L 2 (H 1 (Ω)) norms, and give two numerical examples with different nonlinear terms to verify the theoretical results. The proposed method by combin-ing the FVE method with the WSGD formula can not only make use of the advantages of the FVE method, but also obtain the second-order convergence accuracy in time direction independent of the fractional parameters. In the future, we will extend and apply the proposed method to solve more fractional differential equations.