A Priori Estimates of a Finite Element Method for Fractional Diffusion Problems by Energy Arguments
Abstract
In this article, the Galerkin piecewiselinear finite element (FE) method is applied to approximate the solution of timefractional diffusion equations with variable diffusivity on bounded convex domains. Standard energy arguments do not provide satisfactory results for such a problem due to the low regularity of its exact solution. Using a delicate energy analysis, a priori optimal error bounds in , , and quasioptimal in norms are derived for the semidiscrete FE scheme for cases with smooth and nonsmooth initial data. The main tool of our analysis is based on a repeated use of an integral operator and use of a type of weights to take care of the singular behavior at The generalized Leibniz formula for fractional derivatives is found to play a key role in our analysis. Numerical experiments are presented to illustrate some of the theoretical results.
Key words. Fractional diffusion equation, finite element approximation, energy argument, error analysis, nonsmooth data, generalized Leibniz formula.
AMS subject classifications. 65M60, 65M12, 65M15
1 Introduction
In this paper, we aim to investigate the error analysis via energy arguments of a semidiscrete Galerkin finite element (FE) method for timefractional diffusion problems of the form: find such that
\hb@xt@.01(1.1a)  
\hb@xt@.01(1.1b)  
\hb@xt@.01(1.1c) 
where , is a bounded, convex polygonal domain in with boundary , , and are given functions defined on their respective domains. Here, is the partial derivative of with respect to time and is the Riemann–Liouville timefractional derivative defined by: for ,
\hb@xt@.01(1.2) 
( is the Riemann–Liouville timefractional integral). As , converges to , and thus, problem (LABEL:a) reduces to the classical diffusion equation.
Throughout the paper, we assume that the source term and the diffusivity coefficient function are sufficiently regular and
\hb@xt@.01(1.3) 
Several numerical techniques for the problem (LABEL:a) (with constant diffusivity coefficient) in one and several space variables have been proposed with various types of spatial discretizations including finite difference, FE or spectral element methods, see [2, 10, 9, 18]. For the time discretization, different timestepping schemes (implicit and explicit) have been investigated including finite difference, convolution quadrature, and discontinuous Galerkin methods, see [3, 4, 5, 17, 18, 20, 22]. The error analyses in most studies in the existing literature typically assume that the solution of (LABEL:a) is sufficiently regular including at , which is not practically the case. Indeed, assuming high regularity on imposes additional compatibility conditions on the given data, which are not reasonable in many cases.
Though the numerical approximation of the solution of (LABEL:a) was considered by many authors over the last decade, the optimality of the estimates with respect to the solution smoothness expressed through the problem data, and , was considered in a few papers for the case of constant diffusivity and quasiuniform FE meshes. Obtaining sharp error bounds under reasonable regularity assumptions on has proved challenging. McLean and Thomée [15] established the first optimal error estimates for the Galerkin FE solution of (LABEL:a) with respect to the regularity of initial data. More precisely, for , convergence rates of order ( denoting the maximum diameter of the spatial mesh elements) were proved assuming that the initial data for (see, Section LABEL:sec:WRT for the definition of these spaces). The proof was based on some refined estimates of the Laplace transform in time for the error. In [16] and by using a similar approach, the same authors derived convergence rates in the stronger norm. For , assumed to be in , while for . Recently, McLean and Mustapha [14] studied the error analysis of a first order semidiscrete timestepping scheme for problem (LABEL:a) with allowing nonsmooth using discrete Laplace transform technique. Since standard energy arguments are used heavily in the error analysis of Galerkin FEs for classical diffusion equations, it is more pertinent to extend the analysis to these timefractional order diffusion problems with a variable diffusivity. Since and do not commute, extending these arguments to problem (LABEL:a) is not a straightforward task, especially in the case of nonsmooth .
The main motivation of this work is to derive optimal error estimates of the semidiscrete Galerkin FE method for the problem (LABEL:a) for both smooth and nonsmooth initial data using energy arguments. Since the solution has limited smoothing property [12], a repeated use of the integral operator like (see [6, 7]) along with type weights is an essential tool to provide optimal error estimates. Earlier, for smooth initial data (), error analysis of the Galerkin FE method for the fractional diffusion equation (LABEL:a) was considered in [18] using energy arguments and they have derived quasioptimal error estimates of order in norm. Below, we briefly summarize our main results obtained in this article: for and when with ,

a priori optimal error estimate in norm of order is established, see Theorem LABEL:thm:_smooth_and_nonsmooth. Consequently, for a quasiuniform FE mesh, an error estimate in norm is obtained, see Remark LABEL:rem:_1. By dropping the quasiuniformity mesh assumption, we showed that this error bound remains valid for see Theorem LABEL:thm:_H1_bound. However, for we derived an error estimate, which is reduced to for .

For , and assuming that and the FE mesh is quasiuniform, a quasioptimal error estimate of order in the stronger norm is proved, see Theorem LABEL:thm:_smooth_and_nonsmooth5.
The proposed technique has several attractive features. Some of these are: (1) allowing variable coefficients, and smooth and nonsmooth initial data in the error analysis, (2) the quasiuniform FE mesh assumption is not required to show the convergence results in norm for , and (3) the proposed technique can be applied to other fractional model problems with smooth and nonsmooth initial data. For instance, we discussed in Section LABEL:sec:_Caputo the extension of the achieved error bounds in Theorems LABEL:thm:_smooth_and_nonsmooth, LABEL:thm:_H1_bound and LABEL:thm:_smooth_and_nonsmooth5 to the timefractional diffusion equation: for ,
\hb@xt@.01(1.4) 
where is the timefractional Caputo derivative. The error estimate in Theorem LABEL:thm:_smooth_and_nonsmooth provides an improvement of the result obtained by Jin et al. [8, Theorem 3.7], where the error analysis of the lumped mass FE method was also considered. For constant diffusivity and under the assumption that the mesh is quasiuniform, the FE error bound in [8] involves a logarithmic factor which was derived using a semigroup approach. The error analysis approach in this work can also be used to investigate the error estimates for the FE method applied to the timefractional RayleighStokes problem (described by the timefractional differential equation) presented in the recent work [1], which is close to the one in [15].
Outline of the paper. In Section LABEL:sec:WRT, we recall some smoothness properties of the solution , we also state and derive some technical results. In Section LABEL:sec:Semidiscrete_FE, we introduce our semidiscrete FE scheme and recall some FE error results. We claim that a direct application of energy arguments to problem (LABEL:a) does not lead to optimal convergence rates even when the initial data . In Section LABEL:sec:_LinftyL2, for both smooth and nonsmooth initial data, we derive error estimates for the FE problem in the norm, see Theorem LABEL:thm:_smooth_and_nonsmooth. The generalized Leibniz formula is an essential ingredient in our error analysis. In Section LABEL:sec:_LinftyH1, a superconvergence result in the norm is obtained, see Theorem LABEL:supconv1, and as a consequence, an optimal gradient FE error estimate in the norm is derived in Theorem LABEL:thm:_H1_bound, and a quasioptimal FE error bound in the norm is achieved for , see Theorem LABEL:thm:_smooth_and_nonsmooth5. Particularly relevant to this a priori error analysis is the appropriate use of several properties of the timefractional integral and derivative operators. In Section LABEL:sec:_Caputo, we show that the achieved error estimates for problem (LABEL:a) are valid for the FE discretization of the fractional diffusion model (LABEL:D). Numerical tests are presented in Section LABEL:sec:_Numerical to confirm some of our theoretical findings. Throughout the paper, is a generic positive constant that may depend on and , but is independent of the spatial mesh element size .
2 Regularity and technical results
Smoothness properties of the solution of the fractional diffusion problem (LABEL:a) play a key role in the error analysis of the Galerkin FE method, particularly, since has singularity near , even for smooth given data. Below, we state the required regularity results for problem (LABEL:a) in terms of the initial data and the source term . Over the convex domain , by combining the results of Theorems 4.1, 4.2 and 5.6 in [12], for , we obtain
\hb@xt@.01(2.1)  
for , where Here, denotes the norm on the Hilbert space defined by
where (with ) are the eigenvalues of the operator elliptic operator (subject homogeneous Dirichlet boundary conditions) and are the associated orthonormal eigenfunctions. Noting that, for and for convex polygonal domains, for where (with ) is the standard Sobolev space.
Next, we state the positivity properties of the fractional operators and , and derive some technical results that will be used in the subsequent sections. By [19, Lemma 3.1(ii)] and since the bilinear form associated with the operator (that is, ) is symmetric positive definite on the Sobolev space , it follows that for piecewise time continuous functions
\hb@xt@.01(2.2) 
where denotes the norm.
By [13, Lemma A.1] and again since the bilinear form is symmetric positive definite, the following holds: for ,
\hb@xt@.01(2.3) 
The next lemma will be used frequently in our convergence analysis. In the proof, we use the following integral inequality: if for any , then
Lemma 2.1
Let and let or Assume that
\hb@xt@.01(2.4) 
for . Then
Proof. Choose in (LABEL:eq:_reg_1), and then, integrate over the interval to obtain
Now, use the positivity properties of in (LABEL:eq:_positive_of_Ia) and of in (LABEL:eq:_positive_of_Ba) to find that
For , an application of the integral inequality (stated above) yields the desired inequality. However, for , we use the inequality and the desired inequality follows after simplifying.
3 Semidiscrete FE method
This section focuses on a semidiscrete Galerkin FE scheme for problem (LABEL:a). To define the scheme, let be a family of regular triangulations (made of simplexes ) of the domain and let where denotes the diameter of the element The FE space on is given by
The weak formulation for problem (LABEL:a) is to find such that
\hb@xt@.01(3.1) 
with given Thus, the standard semidiscrete FE formulation for (LABEL:a) is to seek such that
\hb@xt@.01(3.2) 
with given to be defined later.
To derive a priori error estimates for the FE scheme (LABEL:semi), we split the error
where the Ritz projection is defined by the following relation: for all For , the projection error satisfies the following estimates [21]: for
\hb@xt@.01(3.3) 
and hence, by using the regularity property in (LABEL:eq:_regularity_property), we observe
\hb@xt@.01(3.4) 
Next, we show that a direct application of energy arguments to problem (LABEL:a) does not yield satisfactory results due to the low regularity of the continuous solution. From (LABEL:weak) and (LABEL:semi), the decomposition , and the property of the elliptic projection, we obtain the equation in as
\hb@xt@.01(3.5) 
Then, the following result holds, whose proof can be found in [15, 18].
Theorem 3.1
For we have
\hb@xt@.01(3.6) 
Noting that, the solution of problem (LABEL:a) has singularity near . For instance, if and , then by (LABEL:rhoestimate) and the regularity property for
This leads to a quasioptimal convergence for the spatial discretization by linear FEs. To achieve an optimal convergence, a stronger regularity assumption on is required and that, in turn imposes severe restrictions on initial data. Thus, the upper bound in Theorem LABEL:th1 is not sharp even for the case of smooth initial data, that is, regularity on the initial data is not sufficient to get an optimal convergence rate. Furthermore, it is clear that this upper bound is not suitable for the case of nonsmooth initial data. Therefore, we propose in the next section an approach via delicate energy arguments that provides valid error bounds for the problem with both smooth and nonsmooth initial data.
4 error estimates
For convenience, we introduce the notations:
In the next lemma, based on the generalized Leibniz formula for fractional derivatives, we state and show some identities for our subsequent use.
Lemma 4.1
For , the followings hold:
(a) ,
(b) .
Proof. The first identity follows from the fractional Leibniz formula. To show the second identity, noting first that Hence, by (a),
\hb@xt@.01(4.1) 
Now, we replace by in (LABEL:k3a) to obtain the second identity in the lemma.
Next, we derive an upper bound of . To do so, we let where denotes the projection defined by for all
Lemma 4.2
Let Then, we have
Proof. We integrate (LABEL:sup2) over the time interval and obtain
\hb@xt@.01(4.2) 
But because Therefore,
\hb@xt@.01(4.3) 
Multiply by and use by Lemma LABEL:Ia (b) to find that
However, from (LABEL:sup3_ph), we get
\hb@xt@.01(4.4) 
and thus,
\hb@xt@.01(4.5) 
Consequently, an application of Lemma LABEL:lem:_reg_use (with ) yields
\hb@xt@.01(4.6) 
To complete our proof, we rewrite (LABEL:eq:_I_alpha+1) as
Again, an application of Lemma LABEL:lem:_reg_use (with ) shows
\hb@xt@.01(4.7) 
Substitute (LABEL:eq:_theta_less_rho) in (LABEL:eq:_second_last_step) yields the desired bound.
An upper bound of the term will be derived in the next lemma. Again, for convenience, we introduce the following notation
\hb@xt@.01(4.8) 
For later use, by using the projection error estimates in (LABEL:eq:_rho) (with and ) for upper bounds of and , and then integrating, we find that for ,
\hb@xt@.01(4.9) 
Lemma 4.3
Let . Then, the following estimate holds
Proof. We multiply (LABEL:sup2) by so that
\hb@xt@.01(4.10) 
where From the fractional Leibniz formula, we have
Hence, we rearrange (LABEL:k3) as
\hb@xt@.01(4.11) 
and then, by equations (LABEL:sup3_ph) and (LABEL:eq:_I_alpha+1),
\hb@xt@.01(4.12) 
Hence, by Lemma LABEL:lem:_reg_use (with ), we obtain
and thus, an application of the CauchySchwarz inequality yields
Therefore, by using the identity , the inequality in (LABEL:eq:_theta_less_rho) and Lemma LABEL:lem:_estimate_of_Theta_in_l2_norm will complete the rest of the proof.
In the next theorem, we derive optimal convergence results of the FE scheme (LABEL:semi) in the norm for both smooth and nonsmooth initial data . For with we show that the error is bounded by for each . Recall that, for , while for . Noting that, in the limiting case , we recover the convergence rates for the parabolic equation .
Theorem 4.4
Let and be the solutions of and , respectively, with . Then, for
Proof. The desired result follows from the decomposition , the estimate of in Lemma LABEL:H1, the bound in (LABEL:eq:_bound_of_B1), and the estimate of in (LABEL:eq:_rho).
Remark 4.5
In the proof of the above theorem, we used (LABEL:eq:_bound_of_B1) which follows from the projection estimate in (LABEL:eq:_rho) for For , we follow similar steps where will be replaced with , to obtain for Now, for , we notice first that for with ,
Substitute this in the definition of defined in (LABEL:def:_B1), we observe
To estimate the first two terms, we use (LABEL:rhoestimate) (with ) and the following regularity property (which follows from [12, Theorems 4.2 and 5.6])
\hb@xt@.01(4.13) 
where we arrive to
However, by (LABEL:eq:_rho) and the inequality for , we find that
Therefore,
Consequently, by using the above bound of in Theorem LABEL:thm:_smooth_and_nonsmooth, we get the error estimate below that will be used in the forthcoming section to show the convergence of the gradient FE solution. For
\hb@xt@.01(4.14) 
where
Remark 4.6
Under the quasiuniformity condition on , for , from the decomposition , the inverse inequality, the estimate of in Lemma LABEL:H1, and the estimate (follows from the Ritz projection bound in (LABEL:rhoestimate) (with and ) and the regularity property (LABEL:eq:_regularity_property), we obtain the following optimal error estimate in the norm:
By removing the quasiuniformity mesh assumption, this error bound remains valid for see Theorem LABEL:thm:_H1_bound.
Remark 4.7
For smooth initial data , one may choose . An optimal convergence rate can be shown by following the proof of Theorem LABEL:thm:_smooth_and_nonsmooth linebyline, where the term in Lemma LABEL:lem:_estimate_of_Theta_in_l2_norm should be replaced with .
5  and error estimates
In this section, for each , we show optimal convergence results for the gradient FE error in the norm, and quasioptimal error bounds for the FE error in the norm, for both smooth and nonsmooth initial data We start our analysis by deriving an upper bound of .
Lemma 5.1
For and for , we have
Proof. Multiplying (LABEL:sup2) by and then using the identity (Lemma LABEL:Ia (a)), we obtain
\hb@xt@.01(5.1) 
Then, a use of (LABEL:sup3_ph) yields after simplifying
\hb@xt@.01(5.2) 
Now, set in (LABEL:k31), integrate the resulting equation over , and use the positivity property of in (LABEL:eq:_positive_of_Ba), to find that
By the integral inequality (stated before Lemma LABEL:lem:_reg_use), we observe