Numerical Approximation of Space-Time Fractional Parabolic Equations

Andrea Bonito Department of Mathematics, Texas A&M University, College Station, TX 77843-3368. Wenyu Lei Department of Mathematics, Texas A&M University, College Station, TX 77843-3368.  and  Joseph E. Pasciak Department of Mathematics, Texas A&M University, College Station, TX 77843-3368.
July 6, 2019

In this paper, we develop a numerical scheme for the space-time fractional parabolic equation, i.e., an equation involving a fractional time derivative and a fractional spatial operator. Both the initial value problem and the non-homogeneous forcing problem (with zero initial data) are considered. The solution operator for the initial value problem can be written as a Dunford-Taylor integral involving the Mittag-Leffler function and the resolvent of the underlying (non-fractional) spatial operator over an appropriate integration path in the complex plane. Here denotes the order of the fractional time derivative. The solution for the non-homogeneous problem can be written as a convolution involving an operator and the forcing function .

We develop and analyze semi-discrete methods based on finite element approximation to the underlying (non-fractional) spatial operator in terms of analogous Dunford-Taylor integrals applied to the discrete operator. The space error is of optimal order up to a logarithm of . The fully discrete method for the initial value problem is developed from the semi-discrete approximation by applying an exponentially convergent sinc quadrature technique to approximate the Dunford-Taylor integral of the discrete operator and is free of any time stepping.

To approximate the convolution appearing in the semi-discrete approximation to the non-homogeneous problem, we apply a pseudo midpoint quadrature. This involves the average of , (the semi-discrete approximation to ) over the quadrature interval. This average can also be written as a Dunford-Taylor integral. We first analyze the error between this quadrature and the semi-discrete approximation. To develop a fully discrete method, we then introduce sinc quadrature approximations to the Dunford-Taylor integrals for computing the averages.

1991 Mathematics Subject Classification:
65M12, 65M15, 65M60, 35S11, 65R20

1. Introduction

In this paper, we investigate the numerical approximation to the following time dependent problem: given a bounded Lipschitz polygonal domain , a final time , an initial value (a complex valued Sobolev space) and a forcing function , we seek satisfying


Here the fractional derivative in time with is defined by the left-sided Caputo fractional derivative of order ,


Note that (2) holds for smooth and extends by continuity to a bounded operator on satisfying

where denotes the Riemann-Liouville fractional derivative. The differential operator appearing in (1) is an unbounded operator associated with a Hermitian, coercive and sesquilinear form on . For , the fractional differential operator is defined by the following eigenfunction expansion


where denotes the inner product and is an -orthonormal basis of eigenfunctions of with eigenvalues . The above definition is valid for , where denotes the functions such that . A weak formulation of (1) reads: find and satisfying


Here the bilinear form and denotes the duality pairing between and . As a consequence of [17, Theorem 6], the above problem has a unique solution, which can be explicitly written as


Here, for ,




with denoting the Mittag-Leffler function (see the defintion (13)). We also refer to Theorem 2.1 and 2.2 of [19] for a detailed proof of the above formula when , noting that the argument is similar for any .

A major difficulty in approximation solutions of (4) involves time stepping in the presence of the fractional time derivative. The L1 time stepping method was developed in [14] and applied for the case . Letting be the time step, it was shown in [14] that the L1 scheme gives the rate of convergence provided that the solution is twice continuously differentiable in time. For the homogeneous problem (), the L1 scheme is guaranteed to yield first order convergence assuming the initial data is in (see [10]). See also [11] and the reference therein for other time discretization methods and error analyses. We also refer to [13] for the backward time stepping scheme for the case .

The numerical approximation to the solution (5) has been studied recently in [17]. The main difficulty is to discretize the fractional differential operators and simultaneously. In [16], the factional-in-space operator was approximated as a Dirichlet-to-Neumann mapping via a Caffarelli-Silvestre extension problem [8] on . In [17], Nochetto et. al. analyze an L1 time stepping scheme for (4) in the context of the Caffarelli-Silvestre extension problem and obtain a rate of convergence in time of with (see Theorem 3.11 in [17]).

The goal of the paper is to approximate the solution of (4) directly based on the solution formula (5). Our approximation technique and its numerical analysis relies on the Dunford-Taylor integral representation of the solution formula (5). Such a numerical method has been developed for the classical parabolic problem [3, 13] (i.e. the case ) and the stationary problem [4]; see also [5] when the differential operator is regularly accretive [12].

The outline of the remainder of the paper is as follows. Section 2 provides some notation and preliminaries related to (1). In Section 3, we review some classical results from the finite element discretization and provide a key result (Theorem 3.3) instrumental to derive error estimates for semi-discrete schemes. In Section 4, we study the semi-discrete approximation to . Here is the Galerkin finite element approximation of in the continuous piecewise linear finite element space and denote the projection onto . We subsequently apply a sinc quadrature scheme to the Dunford-Taylor integral representation of the semi-discrete solution. For the sinc approximation, we choose the hyperbolic contour for , with . Here denotes the smallest eigenvalue of . Theorem 3.3 directly gives an error estimate for the semi-discrete approximation in fractional Sobolev spaces of order , with . As expected, the rate of convergence depends on the smoothness of the solution which, in term, depends on the smoothness of the initial data and the regularity pickup associated with the spatial exponent . Theorem 4.3 proves that for a quadrature of points with quadrature spacing and depending on , the sinc quadrature error is bounded by , where the constant is independent of and . In Section 5, we focus on the approximation scheme for the non-homogeneous forcing problem. The approximation in time is based on a pseudo-midpoint quadrature applied to the convolution in (5), i.e., given a partition on ,


where is the semi-discrete approximation to . Assuming that the forcing function is in , We show in Theorem 5.3 that the error in the approximation (8) in time is under a geometric partition refined towards (with subintervals). We then apply an exponentially convergent sinc quadrature scheme to approximate the Dunford-Taylor integral representation of the discrete operator . Theorem 5.5 shows that the sinc quadrature leads to an additional error which is . Some technical proofs are given in Appendices A and B.

Throughout this paper, and denote generic constants. We shall sometimes explicity indicate their dependence when appropriate.

2. Notation and Preliminaries

2.1. Notation

Let be a bounded polygonal domain with Lipschitz boundary. Denote by and (or in short and ) the standard Sobolev spaces of complex valued functions equipped with the norms

The scalar product is denoted :

We also denote by the closed subspace of consisting of functions with vanishing traces. Thanks to the Poincaré inequality, we will use the semi-norm as the norm on . The dual space of is denoted and is equipped with the dual norm:

where stands for the duality pairing between and .

The norm of an operator between two Banach spaces and is given by

and in short when .

2.2. The Unbounded Operator

Let us assume that is a Hermitian, coercive and sesquilinear form on . We denote by and the two positive constants such that

Furthermore, we let be the solution operator, i.e. for , , where is the unique solution (thanks to Lax-Milgram lemma) of


Following Section 2 of [12], see also Section 2.3 in [5], we denote to be the inverse of and define .

2.3. The Dotted Spaces

The operator is compact and symmetric on . Fredholm theory guarantees the existence of an -orthonormal basis of eigenfunctions with non-increasing real eigenvalues . For every positive integer , is also an eigenfunction of with corresponding eigenvalue . The decay of the coefficients in the representation

characterizes the dotted spaces . Indeed, for , we set

On , we consider the natural norm

We also denote by the dual space of for . It is known that (see for instance [5])

Note that, we identify functions with functionals by for .

2.4. Fractional Powers of Elliptic Operators

Let be defined from a Hermitian, coercive and sesquilinear form on as described in Section 2.2. For , the fractional power of is given by


In addition, we define the associated sesquilinear form by


which satisfies .

2.5. Intermediate Spaces and the Regularity Assumption

As we saw above, the dotted spaces relies on the eigenfunction decomposition of a compact operator. These are natural spaces to consider fractional powers of operators but are less adequate to describe standard smoothness properties. The latter are better characterized by the intermediate spaces defined for by real interpolation


In order to link the two set of functional spaces introduced above, we assume the following elliptic regularity condition:

Assumption 2.1 (Elliptic Regularity).

There exists so that

  1. is a bounded map of into ;

  2. is a bounded operator from to .

Under the above assumption we have the following equivalence property:

Proposition 2.1 (Equivalence, Proposition 4.1 in [4]).

Suppose that Assumption 2.1 holds for . Then the spaces and coincide for with equivalent norms.

Notice that Assumption 2.1 is quite standard and holds for a large class of sesquilinear forms . An important example is the diffusion process given by

defined on , where satisfies

The in Assumption  2.1 is related to the domain and the smoothness of the coefficients. For example, if is convex and is smooth, Assumption  2.1 holds for any in . In contrast, for the two dimensional L-shaped domain and smooth , Assumption  2.1 only holds for .

2.6. The Mittag-Leffler Function

The Mittag-Leffler functions are instrumental to represent the solution of fractional time evolution, see (6) and (7). We briefly introduce them together with their properties used in our argumentation. We refer to Section 1.8 in [20] for more details.

For and , the two-parameter Mittag-Leffler function is defined by


These functions are entire functions (analytic in ). We note that [20, equation (3.1.42)] (see also [9]) implies that for satisfies

i.e., is a solution of the scalar homogeneous version of the first equation of (1). For this reason, the function will play a major role in our analysis. We also note that




Recall that always denotes the left-sided Caputo fractional derivative (2).

Another critical property for our study is their decay when in a positive sector: For , and , there is a constant only depending on so that


2.7. Solution via superposition

The solution of (4) is the superposition of two solutions: the homogeneous solution and the non-homogeneous solution ,


where is defined by (6) and by (7). Following [19], we have that and in particular .

We discuss the approximation of each term in the decomposition separately. For the homogeneous problem (), we use the Dunford-Taylor integral representation of ,


Here and with the logarithm defined with branch cut along the negative real axis. Given , the contour consists of three segments (see Figure 1):

Figure 1. The contour given by (19).

We use an analogous representation for , namely,


The justification of (18) and (20) are a consequence of (16) and standard Dunford-Taylor integral techniques, see, [21, 2] for additional details.

3. Finite Element Approximations

3.1. Subdivisions and Finite Element Spaces

Let be a sequence of globally shape regular and quasi-uniform conforming subdivisions of made of simplexes, i.e. there are positive constants and independent of such that if for , denotes the diameter of and denotes the radius of the largest ball which can be inscribed in , then


Fix and denote by the space of continuous piecewise linear finite element functions with respect to and by the dimension of .

The projection onto is denoted by and satisfies

For and satisfying , Lemma 5.1 in [5] guarantees the existence of a constant independent of such that


In addition, for any , there exists a constant such that


The case follows from the definition of the -projection, the case is treated in [1, 6] and the general case follows by interpolation.

3.2. Discrete Operators

The finite element analogues of the operators and given in Section 2.2 are defined as follows: For , is defined by

and is given by

so that .

We now recall the following finite element error estimates.

Proposition 3.1 (Lemma 6.1 in [5]).

Let Assumption 2.1 (a) holds for some . Let and set . There is a constant independent of such that


Similar to the operator , has positive eigenvalues with corresponding -orthonormal eigenfunctions . The eigenvalues of are denoted as for . Then the discrete fractional operator is then given by

Its associated sesquilinear form reads


for all .

For any , the dotted spaces described in Section 2.3 also have discrete counterparts , which are characterized by their norms


On , the two dotted norms are equivalent: For , there exists a constant independent of such that for all ,


(see Appendix A.2 in [7]). From the property (cf. [7, equation (2.8)]), we also deduce an inverse inequality in discrete dotted spaces: For , we have


3.3. The Semi-discrete Scheme in Space

We propose a Galerkin finite element method for the space discretization of (5). This is to find satisfying


where the bilinear form is defined by (26) and is the -projection onto . Similarly to the continuous case (see discussion in Section 2.7), the solution of the above discrete problem is given by




and is as in (19).

3.4. A semi-discrete estimate

The purpose of this section (Theorem 3.3) is to obtain estimates for


which, in view of representations (17) and (31), will be instrumental to derive error estimates for the space discretization.

The following lemma assesses the discrepancy between the resolvant and its finite element approximation. Its somewhat technical proof is postponed to Appendix A.

Lemma 3.2 (Space Discretization of Resolvant).

Assume that Assumption 2.1 holds for some . Let and . Then, there exists a positive constant independent of such that for all with , and


We are now in a position to prove the error estimate for the semi-discrete approximation in space. Before doing so, for and , we set


We assume that


The assumption (36) is sufficient to guarantee that the solution is in and we have the following theorem.

Theorem 3.3 (Space Discretization of ).

Let , , and be as in (35). Assume that Assumption 2.1 holds for , and that satisfies (36). Then there exists a constant such that



Without loss of generality we assume that as the case follows from the continuous embedding

Also, we use the notation , and decompose the error in two terms:


For the first term on the right hand side above, we note that the assumptions on the parameters imply that and so the approximation property (23) of yields


We estimate by expanding in Fourier series with respect to the eigenfunctions of (see Section 2.3) and denote by the Fourier coefficient of so that

Two cases need to be considered:

Case 1: . Here, the regularity of the initial condition is large enough to directly use the bound deduced from (16) to get

Case 2: . In this case, we need to rely on the parabolic regularity for . We apply (16) again and obtain

Noting that , a Young’s inequality implies


Returning to (39) after gathering the estimates obtained for the two different cases, we obtain


We return to (38) and estimate now . This time we use the integral representations and the resolvant approximation (Lemma 3.2) to get