Efficient Galerkin finite element methods for a time-fractional Cattaneo equation

In this paper, we develop two efficient fully discrete schemes for solving the time-fractional Cattaneo equation, where the fractional derivative is in the Caputo sense with order in (1, 2]. The schemes are based on the Galerkin finite element method in space and convolution quadrature in time generated by the backward Euler and the second-order backward difference methods. Error estimates are established with respect to data regularity. We further compare our schemes with the L2-1σ scheme. Numerical examples are provided to show the efficiency of the schemes.


Introduction
Let ⊂ R 2 be a bounded convex domain with a boundary ∂ , and T > 0 be a fixed time. The time-fractional Cattaneo equation considered in this paper is described as with a homogeneous Dirichlet boundary condition u(x, t) = 0 on ∂ × (0, T], and initial conditions u(x, 0) = v(x), ∂ t u(x, 0) = b(x) in . Here, the parameter κ is some fixed positive constant related to the relaxation time and μ is a diffusion constant. The source term f and the initial conditions v and b are given functions. The fractional order α belongs to (1,2] and the fractional derivative C D α 0,t is in the Caputo sense defined by where the Riemann-Liouville (R-L) integral RL D -ν 0,t with ν > 0 is given by The time-fractional problem (1.1) has attracted considerable attention in recent decades due to its capacity of modeling the anomalous dynamics of physical diffusion processes. The Caputo derivative C D α 0,t in (1.1) may capture the anomalous diffusion processes. One can refer to [1,2] for derivation details. Equation (1.1) becomes a special case of the timefractional partial differential equation with damping by choosing suitable κ [3]. It is well known that the classical Cattaneo model is introduced to address the insufficiency of the Fickian law in modeling the necessary finite velocity of propagation in heat and signal transport. When α = 2, the time-fractional Cattaneo equation (1.1) recovers the classical one. So the model (1.1) provides a flexible way to model some more general dynamic crossover behaviors [4].
The theoretical studies for the fractional Cattaneo equation have appeared in some literature; see [1,2,4,5], just to name a few. In most of this work, the solution of the fractional Cattaneo equation is obtained in the series form in terms of the H-functions by employing the Laplace and Fourier transform methods. It is very complicated and may not be applicable for numerical simulation. So, we need to resort to the numerical solution.
The numerical schemes for various fractional models have received a lot of attention in recent years; see [6][7][8][9][10][11] for finite difference methods, [12][13][14] for finite element methods, and [15][16][17][18][19][20] for spectral methods [21,22]. However, it seems that the numerical studies for the fractional model (1.1) are still limited. In [23], Ren and Gao presented a compact alternating direction implicit (ADI) difference scheme for solving a two-dimensional Cattaneo equation with the time-fractional derivative. Zhao and Sun proposed a compact difference scheme by applying the Crank-Nicolson discretization in time [24]. Wei developed a fully discrete local discontinuous Galerkin scheme and studied the corresponding stability and convergence [25]. In [3], Chen and Li proposed an ADI Galerkin method for a time-fractional partial differential equation with damping based on the L2-1 σ method in time. We remark that the discretization of Caputo derivative in the above-mentioned numerical schemes are based on the difference method. Recently, Li et al. presented a space-time spectral method for the one-dimensional time-fractional Cattaneo equation [26]. One can see that the error estimates in this work were derived based on the assumption that the solution is sufficiently smooth. However, such assumptions are somewhat restrictive, and in some cases only the regularity of the source term f and the initial data v and b is available [27]. This motivates us to further develop some more reliable numerical schemes for solving (1.1).
In this paper, we aim to develop the robust and efficient Galerkin finite element method for the fractional model (1.1) by employing the convolution quadrature method in time and derive error estimates that are expressed by problem data. With these considerations in mind, we present the semidiscrete Galerkin scheme for (1.1) based on the Galerkin finite element method in space. However, the numerical schemes are only first-order time accurate if one uses convolution quadrature method directly even high-order backward differentiation formulas are employed, since the required compatibility of the problem data is usually not satisfied. In order to restore the desired convergent rate for any t n > 0, we add a few of corrections to the first time step when constructing the second-order backward difference method, which is based on the prevailing treatment [28]. Therefore, stimulated by the idea in [29], we reformulate the semidiscrete Galerkin scheme using the relationship between Caputo derivative and Riemann-Liouville derivatives, and we obtain two fully discrete schemes by employing the convolution quadrature in time generated by backward Euler and second-order backward difference methods. The error estimates are derived with respect to data regularity which involves both smooth and nonsmooth data. We further compare our schemes with the L2-1 σ scheme to show the robustness of the schemes.
The paper is organized as follows. In Sect. 2, we introduce the preliminaries about the notations and lemmas used in this paper. In Sect. 3, we present two fully discrete schemes based on the Galerkin finite element method in space and convolution quadrature in time. The error estimates are established in Sect. 4. In Sect. 5, we compare our schemes with the L2-1 σ scheme. Extensive numerical examples with two-dimensional cases are demonstrated in Sect. 6. Finally, we draw our conclusions in Sect. 7.

Preliminary
We first list some useful notations in this part.
Next, we introduce the notations from the operational calculus [31]. Let K(z) be a complex valued or operator valued function that is analytic in a sector θ where θ ∈ ( π 2 , π). Besides, the function K(z) is bounded by for some real numbers λ and M. Then K(z) is the Laplace transform of a distribution k on the real line, which vanishes for t < 0, has its singular support empty or concentrated at t = 0, and which is an analytic function for t > 0 (see [31,32] for more details). By the inversion Laplace transform for K(z) with t > 0, we get Here the contour lies in θ , and parallel to its boundary and oriented with increasing imaginary part. We denote that θ, We define K(∂ t ) as the operator of convolution with the kernel k: K(∂ t )g = k * g. Here ∂ t is the time differentiation and g is suitably smooth. We divide the time interval [0, T] into a uniform grid with a grid point t k = kτ . Here the time step size τ = T/N with a positive integer N . The convolution quadrature K(∂ τ )g of K(∂ t )g at t = t n is given by where the quadrature weights ω k (τ ) are determined by the generating function Here ϑ is the quotient of the generating polynomials of linear multistep method [32].
In this paper, we focus on the two cases: the backward Euler (BE) and the second-order backward difference methods (SBD). The convolution quadrature has the following error estimates [33]. In the following, we use c to denote a constant that may vary at different occurrences but is always independent of the grid size h and the time step size τ .
Lemma 2.1 Let K(z) be analytic in θ and (2.1) hold. Then, for g(t) = ct σ -1 , the convolution quadrature based on BE (p = 1) or SBD (p = 2) satisfies Finally, we state some fundamental properties for the functions with z ∈ θ by the following lemma.

Lemma 2.2
Let θ ∈ (π/2, π/α) be fixed. Then, for any z ∈ θ , we have g(z) ∈ θ , wherē θ = αθ < π , Proof The proof is trivial and one can complete the proof by a similar idea to that presented in Lemma 2.1 of Ref. [34]. So we omit the details here.

Fully discrete schemes
In this section, we develop two fully discrete schemes by using the standard Galerkin finite element method in space and convolution quadrature in time.

Semidiscrete Galerkin scheme in space
A partition of the domain is denoted by T h in which h is the maximal length of the sides of the triangulation T h . We denote the continuous piecewise linear finite element space V h as The L 2 ( ) orthogonal projection P h : From the definitions of P h and R h , one can see that they are stable in L 2 and H 1 0 , respectively. These properties would be used frequently and may not be mentioned explicitly in the error estimates. We also need the following approximation properties of these two operators P h and R h [30].

Lemma 3.1 The following approximation properties of the operators P h and R h hold:
The semidiscrete Galerkin scheme for (1.1) reads: Here v h and b h are proper approximations to the functions v and b, respectively. The bilinear form a(u, v) is given by (∇u h , ∇χ). By introducing the discrete Laplacian h : V h → V h with the definition: and letting A h =h , we can rewrite the semidiscrete scheme (3.1) as follows:

Fully discrete schemes with convolution quadrature in time
In this part we apply convolution quadrature based on BE and SBD methods in time to obtain two fully discrete schemes.
Using the relationship between Caputo and R-L derivatives: we derive from the semidiscrete scheme (3.2) that with t > 0. Combining the convolution quadrature with the associativity of convolution, we have the approximation of u h (t) at t = nτ with U n h by Here, K(∂ τ ) and ∂ α τ are the convolution quadratures with ϑ(ζ ) = 1ζ for BE scheme, or ϑ(ζ ) = 3 2 -2ζ + 1 2 ζ 2 for SBD scheme in (2.2).
Thus, the BE scheme for (3.5) is stated thus: Find U n h for n ≥ 1 such that Next, we consider a robust numerical scheme which can maintain the second-order accuracy for the scheme (3.5). Since the semidiscrete scheme (3.4) can be recast as Set 1 τ = (0, 3/2, 1, 1, . . .). Using the identity 1 τ = ∂ τ RL D -1 0,t , we obtain the SBD scheme as follows: Find U n h for n ≥ 1 such that That is,

Error estimates
In this section, we establish the error estimates for the proposed numerical schemes by using the operator theoretic approach presented in [29] and originating from [31].

Error estimates for the semidiscrete scheme
We first consider the integral representation of the solution for the homogeneous problem with f = 0. Employing the Laplace transform to (1.1), one has and the functions h(z) and g(z) are given by (2.3). It follows that the solution u(t) is represented by Similarly, the solution u h (t) to (3.2) can be represented by the following: where From the structures of E(z) and E h (z), we let F h (z) = (g(z)I + A) -1 -(g(z)I + A h ) -1 P h for notation convenience. The operator F h (z) has the following property, which plays a key role in the error estimate [34].
Now we are ready to state the error estimates for the semidiscrete problem. Unless otherwise noted, we will always choose δ = 1/t for contour θ,δ when deriving the error bound in what follows. Theorem 4.1 Let u be the solution of problem (1.1) with f = 0, and set u h be the solution of (3.2). Then the following estimates hold: Proof We first consider the second case (b). Subtracting (4.3) from (4.2), we have By Lemma 4.1 and choosing δ = 1/t, we derive that which yields the L 2 -error estimates for the case (b). Next for the first case (a), by (4.3) and (4.2), we have By Lemmas 4.1, 3.1, and 2.2, respectively, we derive that Putting the bounds of I and II together yield the L 2 -error estimates of case (a). A similar argument yields the H 1 -error estimates. So the proof is completed.
Now we consider the inhomogeneous problem with f = 0. By applying the Laplace transform, the solution u(t) can be represented by where By a similar argument, the solution u h (t) to (3.7) is represented by Subtracting (4.4) from (4.5), we get The operator G h has the following error estimate.
Proof By Lemma 4.1, we derive that A similar argument also yields the H 1 -estimates. The proof is thus completed.
We are in a position to state the error estimate for the inhomogeneous problem.
Proof Using the idea of the proof in Theorem 4.4 in [35], we consider the two cases: t ≤ h 2/α and t > h 2/α for the error estimate of e h (t).
For the first case: t ≤ h 2/α , we have So, Next we consider the second case: t > h 2/α . By Lemma 4.2, we have Hence, for t ∈ (0, T], we get the desired results with L 2 -estimate for e h (t) by taking the absolute value of ln(t α /h 2 ). A similar argument yields the same estimate for h ∇e h (t) .
Remark 4.1 For t ∈ (0, T], we may simplify the error bound in Theorem 4.2 as e h (t) ≤ ch 2 | ln h| f L ∞ (0,T;L 2 ( )) , which is an improved version of that in [29] by involving | ln h| instead of | ln h| 2 . Note that the factor | ln h| reflects the limited smoothing property in space.

Error estimates for the fully discrete schemes
In this part, we derive the L 2 -error estimates for the fully discrete schemes (3.6) and (3.9).

Error estimates for the BE method
We first consider the homogeneous problem with f = 0. Using the pervading strategy for the error analysis, we split the error u(t n ) -U n h of (3.6) as a sum of two terms u(t n )u h (t n ) and u h (t n ) -U n h . Since the spatial error u(t)u h (t) has been discussed in the previous part, we focus on the temporal error u h (t n ) -U n h . We state the error estimates as below.

Lemma 4.3 Let u h and U n
h be the solutions of (3. 2) with f = 0 and (3.6) with f h = 0, respectively. Then the following error estimate holds: Proof By (3.4) and (3.5), we have it follows from Lemma 2.1 (let σ = 1, σ = 2 with p = 1 and the suitable chosen parameters λ) that which leads to the desired result. The case (a) can be derived similarly, so the proof is completed.
Remark 4.2 When proving the error estimate in terms of smooth data inḢ 2 ( ), we only can obtain the terms with v h and b h . One may wonder whether this result can be improved to A h v h and A h b h where the identity A h R h = P h A and the L 2 -stability of P h can be applied. However, it seems that this situation may not happen unless the term ∂ t u in (1.1) vanishes. In such situation, the fractional model (1.1) reduces to the fractional diffusion-wave equation which has been discussed in [29]. So in this sense, we further extend the results presented in [29] to some extent.
We now summarize the results for the error estimates of the fully discrete BE scheme (3.6) in the next theorem.

Theorem 4.3 Let u and U n
h be the solutions of (1.1) with f = 0 and (3.6) with f n h = 0, respectively. Then the following error estimates hold: Next we state the error estimates for the inhomogeneous problem with f = 0 but v = b = 0.
Proof It suffices to bound U n hu h (t n ) . From (3.4) and (3.5), we have For the first term I, note that F(z) ≤ c|h(z)|, we apply Lemma 2.2 and Lemma 2.1 (with σ = 1, p = 1, and some suitable parameters λ) to derive that The second term II has the following estimate: Combining the above estimates I and II with Theorem 4.2, we complete the proof.

Error estimate for the SBD method
We now present the error estimates for the homogeneous problem in the next lemma.

Lemma 4.4 Let u h and U n h be the solutions of (3.
2) with f = 0 and (3.9), respectively. Then the following error estimates hold: We first consider the second case: v, b ∈ L 2 ( ). By the identity (g(z)I + A h ) -1 A h = Ig(z)(g(z)I + A h ) -1 and Lemma 2.2, we, respectively, get F 1 (z)A h ≤ c|z| and So employing Lemma 2.1 (with σ = 2, p = 2, and the suitable chosen parameters λ), we derive that Now we turn back to the first case: v, b ∈Ḣ 2 ( ). By noting that use of Lemma 2.1 again yields The error estimates for the term II when b ∈Ḣ 2 ( ) and b h = R h v is similar to that in (b), so the proof is completed.
Remark 4.3 It seems that the error estimates for the smooth data inḢ 2 ( ) in the proof of Lemma 4.4 may not be sharp due to the appearance of the factor t -1 n in b Ḣ2 ( ) . One can refer to Remark 4.2 for similar discussions. Now we summarize the error estimate of the SBD scheme (3.9) for the homogeneous problem.

Theorem 4.5 Let u and U n
h be the solutions of (1.1) with f = 0 and (3.9), respectively. Then the following error estimates holds: Finally, we present the following error estimate for the inhomogeneous problem with v = b = 0.
Proof It suffices to bound e n h = u h (t n ) -U n h . By the expansion f h (t) = f h (0) + tf h (0) + t * f h , we rewrite the expressions of u h (t n ) and U n h in (3.7) and (3.8), respectively, as For ∀z ∈ θ , one has F(z) = E h (z) ≤ c min{|z| -1 , |z| -α /κ}. By Lemma 2.1 (with σ = 2, p = 2, and the suitable chosen parameters λ), we obtain For the third term III, we similarly have Putting the bounds I, II, and III together, we thus finish the proof.
By setting φ h = ∂ t u h , we may rewrite the semidiscrete equation (3.2) as the following first-order system: Applying the L2-1 σ method for the above equation on time level t n+σ for n ≥ 0, we obtain the fully discrete scheme (denoted as the L2-1 σ scheme) for solving (1.1): Find φ n h , n ≥ 1, such that from which we obtain the numerical solution U n h to (1.1): Similar to the strategy of the proof for Theorem 4.1 and 4.2 in [3], one can readily see that the above numerical scheme (5.1) is stable and accurate with order O(τ 2 + h 2 ), under the assumption that u ∈ C 3 ([0, T];Ḣ 2 ( )). The assumption may be too strong which would lead to the limitation of the application for this scheme (5.1). Note that from Theorems 4.3, 4.4, 4.5, and 4.6, one can see that the error estimates of BE and SBD schemes depend only on the regularity of the problem data rather than the regularity of the solution. In addition, the problem data can be incompatible or nonsmooth. From this point of view, the proposed numerical schemes (BE and SBD) are more competitive than the L2-1 σ scheme (5.1).

Numerical experiments
The numerical tests are demonstrated to verify the convergence theory in this part. Two sets of the numerical examples with known solution and unknown solution are considered. The first example with known solution is to confirm whether the error estimates of the three numerical schemes (BE, SBD, and L2-1 σ schemes) are correct when the problem data satisfy certain compatibility condition, i.e., when the solution is sufficiently smooth; while the second example with unknown solution is to test whether the convergences of the two numerical schemes (BE and SBD) are robust and superior to the L2-1 σ scheme for different regularity of the source term f and the initial condition v and b, which can be incompatible and nonsmooth. We test the temporal and spatial convergence rates at a fixed final time T with the L 2 norm errors defined by e n = U n hu(t n ) . The domain = (0, 1) 2 is divided by means of the linear Lagrange triangle finite element. All the tests are done with fixed κ = μ = 1.
In this example, we let γ = 3.5 and test the temporal and spatial convergence rates by employing the proposed schemes (3.6), (3.9), and (5.1), respectively. The numerical results are shown in Tables 1-2. Since the regularity of the solution meet the requirement of all three numerical schemes (3.6), (3.9), and (5.1), one can see that temporal and spatial convergence rates are agree well with the theoretical results.
Example 6.2 (The solution is unknown) Consider problem (1.1) with the following data:   In this example, we examine the temporal convergence rates for the proposed BE and SBD schemes with the final time T = 0.1. The spatial mesh size is chosen as h = 1/10 which is sufficient to observe the convergence rate in time for the semidiscrete scheme (3.2). Since it is difficult to obtain the exact solution, we compute the reference solution instead. The reference solution u h (t n ) is obtained by using the SBD scheme (3.9) with h = 1/10 and τ = 4096. The numerical results are demonstrated in Tables 3-7. From these tables, we can see that the BE and SBD schemes exhibit a steady convergence for both smooth and nonsmooth data, which coincide with the theoretical results. In addition, one can observe from the numerical results that the L2-1 σ scheme achieves less than the expected convergence rate O(τ 2 ) even for smooth data; see Tables 3 and 5. The possible reason is that the solution in such situation does not satisfy the requirement for a sufficiently smooth assumption in (5.1). This shows our schemes are more robust than the L2-1 σ scheme.   Finally, we perform numerical simulation on the fractional model (1.1) using certain data. We consider the following data: f = 0, b = 0, and v = 1 2πσ 2 expx 2 + y 2 2σ 2 ,   Fig. 1, we observe that the numerical solutions with κ = 1 decay slower than these with κ = 0.01 at the time t = 0.01, and behave more like waves as time evolves. So the parameter κ may reflect the propagation velocity of the particle transport [1]. Note that one recovers Fick's second law with an infinite velocity when κ → 0. For the fixed parameter κ = 1 in Fig. 2, it seems that wave feature with α = 1.9 are more obvious than the case with α = 1.1. So the fractional order α essentially indicates the anomalous nature of diffusive transport processes.

Conclusions
In this paper, two efficient and robust fully discrete schemes for solving the time-fractional Cattaneo equation (1.1) are proposed. The schemes are based on the Galerkin finite element method in space and convolution quadrature generated by the backward Euler and second-order backward difference methods in time. Error estimates are established with respect to data regularity, which includes both smooth and nonsmooth initial data. We further compare the proposed BE and SBD schemes with the L2-1 σ scheme and numerically verify the robustness of our schemes. The effects of parameter κ and fractional order α are numerically studied which may provide us a deep insight on the dynamics behaviors of the time-fractional Cattaneo model (1.1) [1,4].