Numerical Approximation of Space-Time Fractional Parabolic Equations
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
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  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  and applied for the case . Letting be the time step, it was shown in  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 ). See also  and the reference therein for other time discretization methods and error analyses. We also refer to  for the backward time stepping scheme for the case .
The numerical approximation to the solution (5) has been studied recently in . The main difficulty is to discretize the fractional differential operators and simultaneously. In , the factional-in-space operator was approximated as a Dirichlet-to-Neumann mapping via a Caffarelli-Silvestre extension problem  on . In , 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 ).
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 ; see also  when the differential operator is regularly accretive .
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
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
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 )
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
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 ).
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  for more details.
For and , the two-parameter Mittag-Leffler function is defined by
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 ,
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):
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  guarantees the existence of a constant independent of such that
In addition, for any , there exists a constant such that
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 ).
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 ,
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
and is as in (19).
3.4. A semi-discrete estimate
The purpose of this section (Theorem 3.3) is to obtain estimates for
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 ).
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