Numerical Approximation of SpaceTime Fractional Parabolic Equations
Abstract.
In this paper, we develop a numerical scheme for the spacetime fractional parabolic equation, i.e., an equation involving a fractional time derivative and a fractional spatial operator. Both the initial value problem and the nonhomogeneous forcing problem (with zero initial data) are considered. The solution operator for the initial value problem can be written as a DunfordTaylor integral involving the MittagLeffler function and the resolvent of the underlying (nonfractional) spatial operator over an appropriate integration path in the complex plane. Here denotes the order of the fractional time derivative. The solution for the nonhomogeneous problem can be written as a convolution involving an operator and the forcing function .
We develop and analyze semidiscrete methods based on finite element approximation to the underlying (nonfractional) spatial operator in terms of analogous DunfordTaylor 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 semidiscrete approximation by applying an exponentially convergent sinc quadrature technique to approximate the DunfordTaylor integral of the discrete operator and is free of any time stepping.
To approximate the convolution appearing in the semidiscrete approximation to the nonhomogeneous problem, we apply a pseudo midpoint quadrature. This involves the average of , (the semidiscrete approximation to ) over the quadrature interval. This average can also be written as a DunfordTaylor integral. We first analyze the error between this quadrature and the semidiscrete approximation. To develop a fully discrete method, we then introduce sinc quadrature approximations to the DunfordTaylor integrals for computing the averages.
1991 Mathematics Subject Classification:
65M12, 65M15, 65M60, 35S11, 65R201. 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
(1) 
Here the fractional derivative in time with is defined by the leftsided Caputo fractional derivative of order ,
(2) 
Note that (2) holds for smooth and extends by continuity to a bounded operator on satisfying
where denotes the RiemannLiouville 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
(3) 
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
(4) 
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
(5) 
Here, for ,
(6) 
and
(7) 
with denoting the MittagLeffler 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 factionalinspace operator was approximated as a DirichlettoNeumann mapping via a CaffarelliSilvestre extension problem [8] on . In [17], Nochetto et. al. analyze an L1 time stepping scheme for (4) in the context of the CaffarelliSilvestre 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 DunfordTaylor 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 semidiscrete schemes. In Section 4, we study the semidiscrete 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 DunfordTaylor integral representation of the semidiscrete 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 semidiscrete 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 nonhomogeneous forcing problem. The approximation in time is based on a pseudomidpoint quadrature applied to the convolution in (5), i.e., given a partition on ,
(8) 
where is the semidiscrete 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 DunfordTaylor 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 seminorm 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 LaxMilgram lemma) of
(9) 
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 nonincreasing 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
(10) 
In addition, we define the associated sesquilinear form by
(11) 
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
(12) 
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

is a bounded map of into ;

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 Lshaped domain and smooth , Assumption 2.1 only holds for .
2.6. The MittagLeffler Function
The MittagLeffler 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 twoparameter MittagLeffler function is defined by
(13) 
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
(14) 
and
(15) 
Recall that always denotes the leftsided 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
(16) 
2.7. Solution via superposition
The solution of (4) is the superposition of two solutions: the homogeneous solution and the nonhomogeneous solution ,
(17) 
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 DunfordTaylor integral representation of ,
(18) 
Here and with the logarithm defined with branch cut along the negative real axis. Given , the contour consists of three segments (see Figure 1):
(19)  
3. Finite Element Approximations
3.1. Subdivisions and Finite Element Spaces
Let be a sequence of globally shape regular and quasiuniform 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
(21)  
(22) 
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
(23) 
In addition, for any , there exists a constant such that
(24) 
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
(25) 
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
(26) 
for all .
For any , the dotted spaces described in Section 2.3 also have discrete counterparts , which are characterized by their norms
(27) 
On , the two dotted norms are equivalent: For , there exists a constant independent of such that for all ,
(28) 
(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
(29) 
3.3. The Semidiscrete Scheme in Space
We propose a Galerkin finite element method for the space discretization of (5). This is to find satisfying
(30) 
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
(31) 
where
(32) 
and is as in (19).
3.4. A semidiscrete estimate
The purpose of this section (Theorem 3.3) is to obtain estimates for
(33) 
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
(34) 
We are now in a position to prove the error estimate for the semidiscrete approximation in space. Before doing so, for and , we set
(35) 
We assume that
(36) 
The assumption (36) is sufficient to guarantee that the solution is in and we have the following theorem.
Theorem 3.3 (Space Discretization of ).
Proof.
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:
(38) 
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
(39) 
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