A Priori Estimates of a Finite Element Method for Fractional Diffusion Problems by Energy Arguments

Samir Karaa222Department of Mathematics and Statistics, Sultan Qaboos University, P. O. Box 36, Al-Khod 123, Muscat, Oman. Email: skaraa@squ.edu.om    Kassem Mustapha222Department of Mathematics and Statistics, Sultan Qaboos University, P. O. Box 36, Al-Khod 123, Muscat, Oman. Email: skaraa@squ.edu.om    Amiya K. Pani333Department of Mathematics, Industrial Mathematics Group, Indian Institute of Technology Bombay, Powai, Mumbai-400076. Email: akp@math.iitb.ac.in

In this article, the Galerkin piecewise-linear finite element (FE) method is applied to approximate the solution of time-fractional 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 quasi-optimal 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.

22footnotetext: Department of Mathematics and Statistics, King Fahd University of Petroleum and Minerals, Dhahran, 31261, Saudi Arabia. Email:kassem@kfupm.edu.sa

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 time-fractional diffusion problems of the form: find such that


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 time-fractional derivative defined by: for ,


( is the Riemann–Liouville time-fractional 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


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 time-stepping 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 quasi-uniform 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 time-stepping 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 time-fractional 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 quasi-optimal 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 quasi-uniform FE mesh, an error estimate in -norm is obtained, see Remark LABEL:rem:_1. By dropping the quasi-uniformity 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 quasi-uniform, a quasi-optimal error estimate of order in the stronger -norm is proved, see Theorem LABEL:thm:_smooth_and_nonsmooth-5.

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 quasi-uniform 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_nonsmooth-5 to the time-fractional diffusion equation: for ,


where is the time-fractional 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 quasi-uniform, 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 time-fractional Rayleigh-Stokes problem (described by the time-fractional 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:Semi-discrete_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:sup-conv-1, and as a consequence, an optimal gradient FE error estimate in the -norm is derived in Theorem LABEL:thm:_H1_bound, and a quasi-optimal FE error bound in the -norm is achieved for , see Theorem LABEL:thm:_smooth_and_nonsmooth-5. Particularly relevant to this a priori error analysis is the appropriate use of several properties of the time-fractional 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


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


where denotes the -norm.

By [13, Lemma A.1] and again since the bilinear form is symmetric positive definite, the following holds: for ,


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


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 Semi-discrete 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


with given Thus, the standard semidiscrete FE formulation for (LABEL:a) is to seek such that


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


and hence, by using the regularity property in (LABEL:eq:_regularity_property), we observe


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


Then, the following result holds, whose proof can be found in [15, 18].

Theorem 3.1

For we have


Noting that, the solution of problem (LABEL:a) has singularity near . For instance, if and , then by (LABEL:rho-estimate) and the regularity property for

This leads to a quasi-optimal 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),


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:sup-2) over the time interval and obtain


But because Therefore,


Multiply by and use by Lemma LABEL:Ia (b) to find that

However, from (LABEL:sup-3_ph), we get


and thus,


Consequently, an application of Lemma LABEL:lem:_reg_use (with ) yields


To complete our proof, we rewrite (LABEL:eq:_I_alpha+1) as

Again, an application of Lemma LABEL:lem:_reg_use (with ) shows


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


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 ,

Lemma 4.3

Let . Then, the following estimate holds

Proof. We multiply (LABEL:sup-2) by so that


where From the fractional Leibniz formula, we have

Hence, we rearrange (LABEL:k-3) as


and then, by equations (LABEL:sup-3_ph) and (LABEL:eq:_I_alpha+1),


Hence, by Lemma LABEL:lem:_reg_use (with ), we obtain

and thus, an application of the Cauchy-Schwarz 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:rho-estimate) (with ) and the following regularity property (which follows from [12, Theorems 4.2 and 5.6])


where we arrive to

However, by (LABEL:eq:_rho) and the inequality for , we find that


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



Remark 4.6

Under the quasi-uniformity 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:rho-estimate) (with and ) and the regularity property (LABEL:eq:_regularity_property), we obtain the following optimal error estimate in the -norm:

By removing the quasi-uniformity 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 line-by-line, 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 quasi-optimal 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:sup-2) by and then using the identity (Lemma LABEL:Ia (a)), we obtain


Then, a use of (LABEL:sup-3_ph) yields after simplifying


Now, set in (LABEL:k3-1), 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