ERROR ESTIMATES FOR A ROBUST FINITE ELEMENT METHOD OF TWO-TERM TIME-FRACTIONAL DIFFUSION-WAVE EQUATION WITH NONSMOOTH DATA

. In this paper, we consider a two-term time-fractional diﬀusion-wave equation which in-volves the fractional orders α ∈ (1 , 2) and β ∈ (0 , 1), respectively. By using piecewise linear Galerkin ﬁnite element method in space and convolution quadrature based on second-order backward diﬀerence method in time, we obtain a robust fully discrete scheme. Error estimates for semidiscrete and fully discrete schemes are established with respect to nonsmooth data. Numerical experiments for two-dimensional problems are provided to illustrate the eﬃciency of the method and conform the theoretical results.


Introduction
During the last few decades, fractional calculus has received increasing attention in the modeling of practical problems, e.g., the anomalous diffusion phenomena in complex system [7], the spread of disease [10,11], the electrical RF circuit [24], and so on. Time-fractional diffusion equation is one of these important models, which can be rigoroully derived by the continuous time random walks (CTRWs) method [21] or formally obtained from the classical diffusion equation by replacing the first (or second) order time derivative by a Caputo derivative of order γ with 0 < γ ≤ 1 (or 1 < γ ≤ 2). In order to improve the modeling accuracy in depicting the anomalous diffusion process, some researchers proposed the multi-term time-fractional diffusion equations (mT-FDEs) and successfully applied in the modeling of some abnormal phenomenon, such as the anomalous behavior of tracer plumes in heterogeneous aquifers [23], the generalized Oldroyd-B fluid through a porous space [15], the unsteady flow of a generalized Maxwell fluid [9], etc.
In this paper, we focus on a two-term time-fractional diffusion-wave equation, which can be viewed as a special case of mT-FDEs [6]. Let Ω ⊂ R 2 be a bounded convex polygonal domain with smooth boundary ∂Ω, T > 0 be a fixed time. The equation is described as follows.
Here ∂ n ∂s n u(·, s) is the n-th partial derivative of u with respect to s and Γ(·) is the Gamma function. Equation (1) can be regarded as the fractional analogue of the classical telegraph equation (if α = 2 and β = 1, then it is just the classical telegraph equation). By choosing β = 1, model (1) recovers the time-fractional diffusion equation with damping [3]: When coefficient κ goes to zero, we obtain the classical time-fractional diffusion-wave equation: Since the term of C D α 0,t u(x, t) is nonzero and κ > 0, we refer to the equation (1) as the two-term time-fractional diffusion-wave equation, see [12] for more details.
Some efforts have been made on the equation (1) in theoretical analysis. Orsingher and Beghin shown that the fundamental solution to time-fractional telegraph equation (that is the case α = 2β in (1)) is the distribution of a telegraph process with Brownian time [22]. Later on, Stojanović and Gorenflo proved the existence and uniqueness of the solution to the nonlinear case of equation (1) by the linearization method [25]. Ferreira et al. studied the functional solution of the multidimensional time-fractional telegraph equation with Laplace and Fourier transforms [8].
For the numerical issue of the equation (1) or other related mult-term cases, it is known that there are some numerical schemes, for examples, see [5,6,18,19,26,30]. One can also refer to the book [17] or the review paper [16] and the references therein. We can see that the time fractional orders in most of the listed papers were restricted on the case (0, 1] or (1,2] except that in the three papers [5,6,19]. Based on the classical L1 method in time, Fan et al. proposed the unstructured mesh finite element method for the multi-term time-space fractional diffusion-wave equation with the time fractional orders belonging to the whole interval (0, 2) [5]. They also provided the stability and convergence analyses for the proposed numerical scheme in details. In [6], the authors employed the temporal finite difference method similar to the above to derive the fully discrete scheme. One can see that the derivation of error estimates in most of aforementioned works were highly relied on the unrealistic smooth assumptions for the exact solution, which is equivalent to rely on the smoothness of the initial data and the compatibility of the source term f at t = 0.
In order to deal with the nonsmooth solution problem in the fractional model, some efforts have been made. By assuming that the analytical solution of time-fractional diffusion-wave equation (i.e., 2T-FDEs) has a certain form, Zeng et al. proposed the second-order accurate in time numerical scheme based on a modified weighted shifted Grünwald-Letnikov formula [29]. That scheme can deal with the insufficiently smooth solution case by introducing some correction terms. This idea is applied to solve multi-term time-fractional mixed diffusion and diffusion-wave equation with fractional orders range from 0 to 2 in [19]. That is, based on an equivalent form of the original equation, Liu et al. derived an efficient numerical scheme which has the second-order convergence in time by adding some corrections. By splitting the right-hand side function f , approximating f (0) and the initial data by the equivalent forms, Jin et al. developed the fully discrete schemes for numerically solving two fractional diffusion equations [13]. Their schemes are second-order accurate in time for both smooth and nonsmooth data if the second-order backward difference method (SBD) is chosen in convolution quadrature. This method was recently systematically developed in [14] for higher order schemes. Extension of this idea to mT-FDEs, especially to (1), is possible but to the best of our knowledge, it is still unavailable in the literature so far.
The goal of this work is to develop a Galerkin finite element method for problem (1) and derive the corresponding error estimates for the nonsmooth initial data and the incompatible source term f . Stimulated by the idea in [13], we propose an efficient fully discrete numerical scheme based on the piecewise linear finite element approximation in space, and convolution quadrature with SBD in time. Specifically, the essential contributions of this study are as follows. First, we develop a robust numerical scheme based on convolution quadrature for the model (1), and derive optimal error estimates without any requirement of regularity of exact solution, only few simple assumptions on source term f , for nonsmooth initial data u 0 , u 1 ∈ L 2 (Ω); see Theorem 4.2. Second, we generalize some error estimate results in [13] to the case of the two-term time-fractional diffusion-wave equation (1); see Theorem 3.3. Especially, we improve an error estimate result in [13]; see Remark 3.4. Third, we theoretically clarify the reason why we need to use the equivalent forms of the initial data u 0 and the source f , and numerically verify this clarification; see Remark 4.4 and Case (c) in Example 5.2 of Section 5.
The rest of the paper is organized as follows. In Section 2 we introduce the preliminary. In Section 3, we present the semidiscrete scheme based on Galerkin finite element method in space and derive the corresponding error estimate. In Section 4, we develop the fully discrete scheme based on convolution quadrature with SBD in time, and establish the error estimate. Two numerical examples are presented in Section 5 to confirm the theoretical results. Finally, some concluding remarks are given in Section 6.
Throughout this paper, the notation c is denoted as a constant which may vary at different occurrences, but is always independent of the mesh size h and the time step size τ .

Preliminary
In this section, we introduce a useful lemma and some notations about the operational calculus in [4], which are fundamental in the development and error estimates of fully discrete scheme.
Let K(s) be a complex valued or operator valued function that is analytic and bounded in a sector Σ θ where θ ∈ ( π 2 , π). Then K(s) is the Laplace transform of a distribution k on the real line, which is an analytic function for t > 0. One can refer to [4,20] for more details. We define K(∂ t ) as the operator of convolution with the kernel k: K(∂ t )g = k * g where ∂ t is the time differentiation. Let the time step size τ = T /N with a positive integer N . The convolution quadrature K(∂ τ )g of K(∂ t )g is given by where the quadrature weights ω j are obtained by the generating function: Here ϑ is the quotient of the generating polynomials of the stable and consistent linear multistep method [20]. In this paper, we focus on ϑ(ζ) = 3 2 − 2ζ + 1 2 ζ 2 for SBD scheme. The backward Euler or other higher order schemes can be considered similarly.
In the following, we derive the representation of the solution in equation (1). Let u(x, z) be the Laplace transform of u(x, t). To simplify the notation, we always denote u(t) = u(x, t) and u(z) = u(x, z) when no ambiguity arises. Similar notations are also used for f (x, t) and f (x, z). An application of the Laplace transform to (1) yields So, we formally have with the kernel function K(z) given by K(z) = (g(z)I + A) −1 . Combining the resolvent estimate (4) and Lemma 2.1, we readily have . By the inversion of Laplace transform to (11), the solution of problem (1) can be written by where the operators E, E, and E are respectively defined by and Here, the contour Γ θ,δ is defined by Γ θ,δ = ρe ±iθ : ρ ≥ δ ∪ δe iψ : |ψ| ≤ θ with θ ∈ (π/2, π) and δ > 0. The representation (12) serves to define the mild solution of the equation (1). We do not intend to go further for solution theory since it will be off the focus of the current paper. We leave the rigorous analysis of the mild solution of (1) as one of the future study.

Semidiscrete Galerkin FEM in space
In this section, we propose the semidiscrete finite element method and derive its error estimate.

Semidiscrete Galerkin scheme
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 . Then the continuous piecewise linear finite element space V h is given by If u 0 and u 1 both belong to L 2 (Ω), then the semidiscrete Galerkin FEM for (1) reads: where the initial value conditions u h (0) = P h u 0 ∈ V h and ∂ t u h (0) = P h u 1 ∈ V h . Here, the operator P h : and letting A h = −∆ h , we can rewrite the semidiscrete scheme (13) as follows: where

Error estimate for semidiscrete scheme
Similar to the derivation of (12), the solution to (14) can be represented by where respectively. We have the error estimate between K(z) and K h (z) [2].
we present the error estimates for these operators F h (t), F h (t), and F h (t). Unless otherwise specified, we always set δ = 1/t for the contour Γ θ,δ in the following proofs.
Similarly, we can derive that employing Lemma 3.1 again, we obtain The error estimates of ∇F h (t)φ L 2 (Ω) and ∇ F h (t)φ L 2 (Ω) can be proved in a similar way, so the proof is completed.
Subtracting (15) from (12), we obtain the error equation: Now we state the error estimate of the semidiscrete scheme (14) for the nonsmooth initial data u 0 , u 1 ∈ L 2 (Ω).
. Let u and u h be the solutions of problems (1) and (14) with u 0 , u 1 ∈ L 2 (Ω), respectively. Then Proof. We consider the following two cases.
(a) Case u 0 , u 1 = 0, but f = 0. The error equation (16) becomes By Lemma 3.2, we readily get (b) Case u 0 = u 1 = 0, but f = 0. In this situation, we consider the two cases: t ≤ h 2/β and t > h 2/β for the error estimate of e h (t) by following the idea in [1].
For the first case: Since the term I in (17) yields Since 2α/β − 2 > 0, we have h 2α/β−2 < 1 for any given mesh size h < 1. It follows that Similarly, we can obtain a bound for II which is the same with I. Combining the error estimates for I and II, we have e h (t) ≤ ch 2 min(1, κ −1 ) f L ∞ (0,T ;L 2 (Ω)) .
Now we consider the second case: t > h 2/β . We have Hence, for t ∈ (0, T ], we get Putting these error bounds together, we arrive at the desired result.
Remark 3.4. From Theorem 3.3, when κ goes to zero, we can readily obtain the estimate: for the time-fractional diffusion-wave equation (3). Comparison with the error estimate with Theorem 3.3 in [13], we improve that result by involving only one | ln h| instead of | ln h| 2 .

Fully discrete scheme
In this section we apply convolution quadrature based on SBD formula in time to obtain a fully discrete scheme and derive the corresponding error estimate.
Then using the relationship between Caputo and R-L derivatives: (0)) , and combining with the semidiscrete scheme (14), we derive that Applying K(∂ t ) = ( RL D α 0,t + κ RL D β 0,t + A h ) −1 on both sides of the above semidiscrete scheme, we obtain Combining the convolution quadrature with the associativity of convolution, we have the approximation of u h (t n ) at t n = nτ with u n h by where n = 1, 2, · · · , N, and f n h = P h f (t n ). Here K(∂ τ ), ∂ α τ and ∂ β τ are the convolution quadratures generated by the SBD formula.
The above numerical scheme (19) has a more familiar form: Next, we present a more robust numerical scheme which can maintain the second-order accuracy for the original SBD scheme (20) (the reason why we need to do the correction, please refer to Remark 4.4, or [13]).
Since (18) can be rewritten as Using the splitting It follows that wheref n h = f h (t n ) − f h (0). The exact contribution RL D −1 0,t is kept in order to preserve the second-order time accuracy.

Error estimate for fully discrete scheme
In this subsection, we only derive the error estimate for the fully discrete scheme (26). The error estimates for (20) and (23) will be obtained as the byproducts, see Remarks 4.4 and 4.5.
We shall need the following important lemma [20].
Lemma 4.1. Let K(z) be analytic in Σ θ and K(z) ≤ c|z| −λ , ∀z ∈ Σ θ for some real λ. Then for g(t) = ct γ−1 , the convolution quadrature based on the SBD method satisfies We obtain the error estimate for the fully discrete scheme (26) as follows.
(a)Case u 0 , u 1 = 0, and f = 0. Let K(·) = K(·)(·), and K(·) = K(·)(·) α , with notation · denoted by ∂ t or ∂ τ . Subtracting (24) from (25) at grid-point t n , we get Letting g(z) = z α + κz β , we derive that Then Lemma 4.1(with λ = −1, γ = 2 and λ = 0(or β − α), γ = 2) and the L 2 (Ω) stability of P h yield , again applying Lemma 4.1 (with λ = α − 1, γ = 2 and λ = β − 1, γ = 2), we get in which the singular integral term included (t n − s) β−2 is dropped . Taking these error bounds together, we complete the proof by triangle inequality.  (26) is second-order accurate both in time and in space for any fixed t n away from 0. It is worth noting that, the terms with the power of t n appearing in the error estimate reflect the singularity behavior of the solution near t = 0 for nonsmooth data. Furthermore, when κ goes to zero, one can have the following error estimate of the time-fractional diffusion-wave equation (3): which improves the one presented in [13] since there's only one | ln h| here.
Remark 4.4. In this remark, we discuss the error estimate for original numerical scheme (20). Suppose that the source term f has the expression: f (x, t) = g(x)( j c j t γj ), where g ∈ L 2 (Ω), c j and γ j are some nonnengative constants, then from the idea of proof in Theorem 4.2, we readily get So, the temporal convergent order of the original numerical scheme (20) at most first-order accuracy if u 0 = 0 or f h (0) = 0. This is the reason why we need to do some modification for the initial condition u 0 and the source term f in order to obtain a higher order in time.
Remark 4.5. Along the line of proof in Theorem 4.2, we can derive the error estimate for the fully discrete scheme (23), that is Thus the numerical scheme (23) is less competitive than (26) since (23) requires f (0) exists.
Remark 4.6. It is interesting to see that the schemes (20) and (23) coincide with each other when u 0 = f (0) = 0. Furthermore, when u 0 = f = 0, the three numerical schemes (that is, (20), (23), and (26)) have the same form and their corresponding error estimates also coincide with each other.

Numerical experiments
In this part, numerical examples are presented to verify the convergence results derived in previous sections. In all tests, we fix the coefficient κ = 1, and examine the temporal and spatial convergence rates at a fixed final time T. The L 2 -norm errors on Ω = (0, 1) 2 are measured by U n h − u h (t n ) L 2 (Ω) at t n = T . Example 5.1 (The solution is known). Consider problem (1) with u 0 = u 1 = 0 and the source term: f = sin(πx) sin(πy) t γ + g(γ, α)t γ−α + κg(γ, β)t γ−β , where γ ≥ α and g(x, y) = Γ(x + 1)/Γ(x + y − 1), The corresponding solution is provided by u(x, y, t) = sin(πx) sin(πy)t γ . Here, we let γ = 3.1. To examine the temporal convergence rates of the proposed SBD scheme (26), we choose a sufficiently small mesh size h = 1/512, so that the error generated by spatial discretization is negligible. Similarly, in order to ignore the errors arising from temporal discretization, we take the time step size τ = T /4096 when testing spatial convergence. The numerical results are shown in Tables 1 and 2. We can observe that the convergence order in both time and space is two when the finial time T = 1. Below we are also interested in the performance of the SBD scheme (26) at time T close to initial moment. One can see that the SBD scheme still has 2nd-order accuracy in the temporal and spatial directions at T = 0.1. Thus these numerical results validate our theoretical analysis in Theorem 4.2.  (1)  Here, χ S is the characteristic function of the set S, u 0 ∈Ḣ (Ω) for 0 ≤ < 1/2 in (a), and u 1 ∈Ḣ (Ω) in (b).
In the following, we focus on the temporal convergence rate at T = 0.1. The analytic solution of (1) with the above data can be expressed in terms of H-function with two variables [8], which however is difficult to compute. In order to observe the convergence order in the time direction, we only need to compute a semidiscrete reference solution here on very refined mesh, i.e., the semidiscrete solution u h (t n ) with h = 1/10 and τ = T /4096. The numerical results for cases (a) and (b) are demonstrated on Table 3. In both cases (a) and (b), the data for the initial conditions are not smooth, we still observe a second-order convergence accuracy in the time direction, which fits well with our convergence theoretical analysis in Theorem 4.2 again. Next, we test the effect of the regularity of the right-hand side function f on the numerical solution. To this end, we compute and compare the numerical results using three different numerical schemes (i.e., (20), (23), and (26)) for case (c). See Tables 4 and 5. The relevant discussions are as follows.
(I) Case c = 1 and 0 < γ < 2. From the numerical results in Table 4, we observe that the accuracy in time for the schemes (20) and (23) (23), when γ = 0.1, we still can observe some convergent order which is greater than one, although f (0) does not exist. This scheme exhibits second-order accuracy when γ = 1.1 which is in agreement with Remark 4.5. (II) Case c = 0 and 0 < γ < 2. In this case, u 0 = f (0) = 0, so the schemes (20) and (23) have the same form. Thus here we only present the numerical results for (20) and (26), see Table 5. The scheme (20) has better performance than that in the above case (I), and it can achieve second-order convergent order when γ = 1.1, which is in good agreement with our analysis in Remark 4.4. In the above two cases (I) and (II), the scheme (26) has second-order convergence order in time, so it is more robust than the other two numerical schemes (20) and (23), which verifies the correctness of our error estimate in Theorem 4.2. Finally, we perform numerical simulation of the fractional model (1) using data from case (a) in Example 5.2. Since numerical simulations of the fractional-time diffusion-wave equation (3) have been given in other articles, e.g., see [13], we focus here on the effect of parameters κ and fractional order β on the solution behavior of the model. The numerical results are presented in Figures 1 and 2. We can see that whether it is fixed κ, change β, or fixed β, change κ, the corresponding numerical profiles of the model (1) are not much different when the time T = 0.01 or T = 0.1, see the first two lines in Figures 1 and 2. However, as time progresses, the smaller the fractional order β, the slower the numerical solution decay, and the behavior of the solution at this time behaves more like wave propagation, which suggests that the term RL D α 0,t u in (1) is playing a dominant role, while the term RL D β 0,t u is playing a damping role. Similar results are also observed in the case of fixed fractional order β, varying coefficient k. See the last two lines in Figures 1 and 2, respectively. Based on this, we conclude that the equation (1) possess the great potential flexible capability in modeling abnormal diffusion than the other two equations (2) and (3). For the corresponding theoretical analysis, one may refer to [8].

Conclusions
In this paper, a robust fully discrete scheme for two-term time-fractional diffusion-wave equation with second order convergence both in time and space are developed. The scheme is based on the piecewise linear Galerkin finite element method in space and convolution quadrature with second-order backward difference in time. Without any regularity requirement of exact solution, the error estimate of the fully discrete scheme for nonsmooth data is established rigorously. In the process of deriving the numerical scheme, two other numerical schemes (20) and (23) were also obtained, and the corresponding error estimates were discussed and presented as the byproducts. Numerical tests for two-dimensional problems fully confirm our error estimates, further illustrating that the proposed scheme (26) outperforms the other two schemes (20) and (23) in terms of effectiveness In subsequent studies, we shall complete the solution theory of the equation (1), as well as construct effective finite element schemes based on convolution quadrature of other generating functions, and give error estimates with respect to data regularity (including smooth data). The nonlinear problems are also one of our interests. In [28], based on the convolution quadrature generated by k-step backward differentiation formula, Wang and Zhou developed and analyzed high-order time stepping finite element schemes for the semilinear subdiffusion equations. Whether their approach can be extended to the corresponding nonlinear problems of (1) is worthy of further study.