Efficient Temporal Third/Fourth-Order Finite Element Method for a Time-Fractional Mobile/Immobile Transport Equation with Smooth and Nonsmooth Data

In recent years, the numerical theory of fractional models has received more and more attention from researchers, due to the broad and important applications in materials and mechanics, anomalous diffusion processes and other physical phenomena. In this paper, we propose two efficient finite element schemes based on convolution quadrature for solving the time-fractional mobile/immobile transport equation with the smooth and nonsmooth data. In order to deal with the weak singularity of solution near t=0, we choose suitable corrections for the derived schemes to restore the third/fourth-order accuracy in time. Error estimates of the two fully discrete schemes are presented with respect to data regularity. Numerical examples are given to illustrate the effectiveness of the schemes.


Introduction
In this paper, we are interested in the robust high-order numerical methods for the following time-fractional mobile/immobile transport equation [1]: which is subject to the boundary condition u = 0 on ∂Ω × (0, T] and the initial condition u(x, 0) = v(x) in Ω. Here, Ω ⊂ R 2 is a bounded convex polygonal domain and the parameters κ 1 , κ 2 , and µ are positive constants. The source term f and the initial condition v are both given functions. The notation RL D α 0,t with α ∈ (0, 1) denotes the Riemann-Liouville derivative defined by In recent years, Equation (1) has received extensive attention, and it has significant applications in diverse fields, especially in the modeling of groundwater solutes [2][3][4]. There are some numerical studies for Equation (1) or its variants; see [1,[5][6][7][8][9] and the references therein.
In our recent paper [1], we studied the solution regularity of Equation (1) and proposed two efficient finite element schemes by employing convolution quadrature based on the backward Euler and the second-order backward difference methods. It is proved that the proposed schemes are robust with respect to data regularity, including the nonsmooth initial data, i.e., v ∈ L 2 (Ω). One of the distinct features of the two schemes is that they do not impose any compatibility conditions on the source term f , so this avoids some unrealistic assumption for the solution regularity, i.e., the assumption of the sufficiently smooth solution. When dealing with the nonsmooth solution problem, we observe numerically that the convolution quadrature generated by the high-order backward difference formula (denoted by BDF) is only of first-order accuracy when no corrections are added. This motivates us to investigate how to choose some suitable corrections in order to restore the desired high-order accuracy.
Since the convolution quadrature has a nice stable discrete structure, for which it is easy to perform the analysis of the resulting numerical scheme [10,11], it seems that it is suitable for developing high-order numerical schemes of the fractional model with the Riemann-Liouville derivative. In order to restore the high-order accuracy, one usually needs to add some corrections to the original scheme. There are a few papers on this direction [12][13][14][15][16][17][18], just to name a few. One may refer to two review papers [19,20] for more details. In [14,16,17], the convolution quadratures generated by the backward Euler and the second-order backward difference methods were applied to obtain the robust scheme for the time-fractional Cattaneo equation, the two-term time-fractional diffusionwave equation, and the modified anomalous subdiffusion equation, respectively. In [12], the authors constructed and analyzed an efficient fractional Crank-Nicolson finite element scheme with two corrections at the starting time levels. They proved that the proposed scheme is of second-order accuracy in time for both smooth and nonsmooth data. Recently, Wang et al. improved the results by adding only one single correction at the starting time level and presented the optimal error estimates of their resulting scheme [18]. In [13], Jin et al. developed the BDFk (with the k denoting the order of BDF) convolution quadrature with suitable correction terms to deal with the temporal discretization of the fractional evolution equation. Later on, Shi and Chen extended the BDFk convolution quadrature to numerically solve the fractional Feynman-Kac equation with Lévy flight [15]. As far as we know, there is no convolution quadrature generated by high-order BDF to solve Equation (1) with smooth and nonsmooth problem data. The goal of this paper is to fill this gap.
The contributions of this paper are listed below. First, based on the convolution quadrature generated by BDF3/BDF4, we develop the robust temporal third/fourthorder finite element scheme for problem (1) by carefully choosing the suitable corrections at the starting two/three steps. Second, the error estimates of the resulting schemes are expressed in term of data regularity, and are proved rigorously, cf. Theorems 2-5. Third, the designed numerical tests show that the two schemes indeed yield the temporal third/fourth-order theoretical accuracy, so this further verifies numerically the correctness of the error estimates and illustrates the significant improvement in the temporal accuracy of our schemes compared to that of [1], cf. Tables 4-6 in Section 6.
The rest of this paper is organized as follows. In Section 2, we propose the temporal third-order finite element scheme based on the convolution quadrature generated by BDF3. In Section 3, the method to obtain the suitable correction terms is discussed in detail. In Section 4, we prove that the BDF3 is of third-order accuracy, with respect to smooth and nonsmooth data. The generalization of the BDF3 is given in Section 5, and numerical examples are provided in Section 6 to verify the effectiveness of the two proposed numerical schemes. Finally, we present the conclusion in the last section, i.e., Section 7. The symbol c in this paper denotes a positive constant that is independent of the temporal and spatial step sizes.

The Temporal Third-Order Numerical Scheme
In this section, we propose a temporal third-order fully discrete scheme for Equation (1) based on convolution quadrature in time and the standard Galerkin finite element method in space.

Spatial Semidiscrete Scheme
We first introduce the semidiscrete scheme by using the finite element method in space and recall the corresponding error estimates [1].
Denote the continuous piecewise linear finite element space V h as follows: where T h is the partition of the domain Ω with h = max K∈T h h K , and h K is the diameter of the triangles K [21]. The semidiscrete scheme for (1) is given by the following: Find u h (t) = u h (·, t) ∈ V h such that the following holds: Here, v h is the proper approximation to the function v. In this paper, we [21]. Here, P h and R h are the L 2 (Ω) orthogonal projection operator and Ritz projection operator, respectively. By using the discrete Laplacian −∆ h : where f h (t) = P h f (t).

Fully Discrete Schemes
We divide the time interval [0, T] into uniform grids {t n = nτ} N n=0 with τ = T/N. Here, N is a given positive integer. Based on the convolution quadrature generated by BDFk, the numerical approximation of RL D α 0,t u h (t) at t = t n is described by in which u n h = u h (t n ) and the weights w j are provided by Applying the convolution quadrature generated by BDF3/BDF4 for solving the semidiscrete scheme (2), one may obtain the following fully discrete finite element scheme: Find the numerical approximation U n h , n ≥ 1 of u(t n ), such that the following holds: where n = 1, 2, · · · , N, U 0 h = v h , and F n h = f h (t n ). It is known that the scheme (4) generally exhibits only temporal first-order accuracy, due to the low regularity of the solution near t = 0. Such a fact is also verified by the numerical examples in Section 6.
In what follows of this section, we consider the scheme (4) generated by BDF3 and use the idea presented in [13] to derive a modified scheme which can restore the temporal third-order accuracy. To this end, we define The semidiscrete scheme (2) can be rewritten as follows: In view of the construction of BDF2 in [1], we may naturally consider the following two corrections for (4) at the starting two steps: and where W n h = U n h − v, n ≥ 0. The four parameters λ ij (i, j = 1, 2) are unknown and need to be determined.
Proof. By using the Laplace transform technology, we can easily obtain the integral representation (9) from (5). One may refer to [1] for further details. It remains to derive the fully discrete solution in (4) based on discrete Laplace transform, i.e., generating function.
. The scheme (4) has the following equivalent form: with the two initial corrections (6) and (7). Multiplying ξ n on both sides of (11), and summing up for n = 1, 2, · · · , we obtain the following: From (12), we have Thus, with r ∈ (0, 1) and the notation Γ τ r := {z = − ln(r)/τ + iy : y ∈ R and |y| ≤ π/τ}. Here, we have used the change of variables ξ = e −zτ in the second equality of (14). Denote the complex region as Σ enclosed by Γ τ (14) is analytic with respect to z in Σ, so applying the Cauchy's theorem, we have This leads to the desired result (10) since there holds All this completes the proof of the lemma.
From the representations of (9) and (10), it is intuitive to select λ ij (i, j = 1, 2) in (10) such that and for any given z ∈ Γ τ θ,δ . For the left hand sides of (15) and (16), the triangle inequality leads to and We shall need the following lemma.

Lemma 2.
Let K(z) and δ τ (ξ) be defined as in (8) and (3) with k = 3, respectively. Then, for any In view of Lemma B.1. in [13], we derive that and Thus, we further deduce that for z ∈ Γ τ θ,δ , where the last inequality holds since |z| < 1. By Lemma 2, we only need to choose and It is derived in ( [13], cf. (16) and (17)) that the above error estimates can be valid by choosing λ 11 = 23/12, λ 12 = 13/12, λ 21 = 7/12, and λ 22 = 1/2. So, the suitable corrections in (6) and (7) can be written as follows: and The fully discrete scheme with the two corrections (23) and (24) as the starting two steps is referred to as the BDF3, and can restore third-order accuracy in time, which will be proved in the next section.

Error Estimates of the BDF3
In this part, we present the error estimates for the modified scheme (25). We first consider the error estimate for a homogenous problem with v ∈ D(∆) and f = 0. In the proof of error estimates, we may explicitly or implicitly employ the identity ∆ h R h = P h ∆ and the L 2 -stability of the orthogonal projection P h when needed. (1) and (25), respectively. Assuming v ∈ D(∆) and f = 0. Then, for v h := R h v, we have

Theorem 2. Let u and U n h be the solutions of
Proof. It suffices to consider the bound of u h (t n ) − U n h in view of Theorem 1. From (9) and (10) in Lemma 1, we can obtain Since K(δ τ (e −zτ )) ≤ c min{|z| −1 , |z| −α }, in view of (17) with Lemma 2 and letting δ ≤ 1/t n , we have For I 2 , we get e t n r cos θ r −α−1 dr.
Note that r ≥ π/(τ sin θ), one has r −1 ≤ r 2 τ 3 (sin θ/π) 3 . So, we have Putting the error estimates of the two parts, I 1 and I 2 , into (26), we complete the proof of the theorem.
Next, we state the error estimate for the homogenous problem with v ∈ L 2 (Ω) and f = 0. Theorem 3. Let u and U n h be the solutions of (1) and (25), respectively. Assuming v ∈ L 2 (Ω) and f = 0. Then, for v h := P h v we have In view of (26), we have ). By the triangle inequality, the termK(z) has the following error estimate: For the first term I 11 , using the result in (19) and the criterion (21), we have So, For the second term I 12 , we have Applying the error estimates in Lemma B.1. of [13] again, we derive that Since we obtain the error estimate for I 12 : I 12 ≤ cτ 3 |z| 2 .
Consequently, by choosing δ ≤ 1/t n , we obtain For the second term I 2 , similar to the derivation of I 2 in Theorem 2, we can easily derive the following inequality: where the first inequality is valid since This, together with the error estimate of I 1 and Theorem 1, ends the proof. Now, we consider the inhomogenous problem, i.e., f = 0 and v = 0. (1) and (25) with v = 0, respectively. Then, we have

Theorem 4. Let u and U n h be the solutions of
Proof. It follows from (9) and (10) in Lemma 1 that where From the derivation in Theorem 2, we can easily obtain the following error estimate of I 1 : Next, we consider the two terms I 2 and I 3 , respectively. Through the triangle inequality, we can obtain Since (κ 1 δ τ (e −zτ ) + κ 2 δ τ (e −zτ ) α − µ∆ h ) −1 ≤ c min{|z| −1 , |z| −α }, using the criterion (16), we have I 21 ≤ cτ 3 |z| 1−α . For I 22 , the following inequality holds: Analogous to the derivation of the error estimates of I 121 and I 122 for I 12 in Theorem 3, we can obtain I 22 ≤ cτ 3 |z| 1−α . So, we choose δ = 1/t n and get

For the last term I 3 , by using the Taylor expansion of
. It suffices to analyze the error estimates for source terms of the form f 1 (t) and f 2 (t), respectively. Assume that f h (t) ≡ f 1 (t). Then the term I 3 has the following form: 3 . By repeating the similar argument in the above error estimate for I 2 , one can derive that I 3 ≤ cτ 3 t α−1 n f (0) . When f h (t) ≡ f 2 (t), the corresponding semidiscrete Galerkin solution in (9) can be written as follows: where the operator E(t) is given by Furthermore, by (13), we have from which we can derive the representation of the fully discrete solution U n h below: with E(δ τ (ξ)) = ∑ ∞ n=0 E n τ ξ n . In view of (14), we get Denote E τ (t) = ∑ ∞ n=0 E n τ δ t n (t) with the Dirac delta function δ t n (t) at t n , we can rewrite (27) as with the following three corrections: The corresponding error estimates for BDF4 is described by the following theorem.

Theorem 5. Let u and U n
h be the solutions of (1) and (28), respectively. Then, we have where β = 0 and v = v if v ∈ L 2 (Ω) and β = α and v = ∆v if v ∈ D(∆).
Proof. Similar to the idea discussed in the error estimates of Theorems 2-4, one can easily deduce the desired error estimate results for BDF4, and thus we omit the details here.

Numerical Examples
Now, we perform the numerical examples to test the error estimates of BDF3 (25) and BDF4 (28). Let Ω = (0, 1) 2 and T = 0.1. We consider the following three types of problem data: Here, χ S is denoted as the characteristic function of the set S. We remark that the three types of data considered above cover all cases of error theories, and therefore this is sufficient for our numerical tests. Moreover, we focus on the temporal convergence rate since the spatial one is well understood. Since the analytical solution of equation with the above data is hard to obtain, we use reference solution instead and the reference solution is computed by BDF4 with fixing h = 1/10 and τ = T/4096. The L 2 (Ω) norm errors are measured by u(t n ) − U n h . The numerical results are demonstrated in Tables 1-6. From these tables, we can see that when no correction terms are added to BDF3 and BDF4, the temporal convergence accuracy is only first-order accuracy for both (see Tables 1-3), while after adding suitable correction terms (cf. Tables 4-6), we observe the theoretical temporal third-/fourth-order accuracy, which numerically verifies the importance of adding correction terms and the correctness of our constructed correction terms in BDF3/BDF4.
In addition, we also numerically compare the BDF3 and BDF4 with the BDF2 proposed in [1]; see Tables 4-6. The numerical results further indicate the correctness of our error estimates presented in Theorems 2-5 and reveal that our schemes greatly improve the temporal convergence accuracy of the BDF2.

Conclusions
In this paper, we developed an efficient temporal third/fourth-order finite element scheme for model (1) with smooth and nonsmooth data. The schemes were obtained by using the Galerkin finite element method in space and the convolution quadrature generated by the BDF3 and BDF4 in time, so this is why we refer to the resulting schemes as BDF3 and BDF4, respectively. In order to restore the desired accuracy for both smooth and nonsmooth data, we carefully chose the suitable corrections for the BDF3 and BDF4 at the starting two/three steps. Error estimates of the two modified schemes were established, with respect to the data regularity. Finally, the numerical examples confirmed the robustness and accuracy of the proposed schemes. With the BDF3/BDF4 modified at the starting two/three steps, we greatly improved the convergence results presented in [1].

Informed Consent Statement: Not applicable.
Data Availability Statement: The data of numerical simulation used to support the findings of this study are included within the article.

Conflicts of Interest:
The authors declare no conflict of interest.