Error analysis of a finite volume element method for fractional order evolution equations with nonsmooth initial data

# Error analysis of a finite volume element method for fractional order evolution equations with nonsmooth initial data

Samir Karaa and Amiya K. Pani111 Department of Mathematics, Industrial Mathematics Group, Indian Institute of Technology Bombay, Powai, Mumbai-400076. Department of Mathematics and Statistics, Sultan Qaboos University, P. O. Box 36, Al-Khod 123, Muscat, Oman. Email: skaraa@squ.edu.om. The research of this author is supported by Sultan Qaboos University under Grant IG/SCI/DOMS/16/01.
###### Abstract

In this paper, a finite volume element (FVE) method is considered for spatial approximations of time-fractional diffusion equations involving a Riemann-Liouville fractional derivative of order in time. Improving upon earlier results (Karaa et al., IMA J. Numer. Anal. 2016), optimal error estimates in - and -norms for the semidiscrete problem with smooth and middly smooth initial data, i.e., and are established. For nonsmooth data, that is, , the optimal -error estimate is shown to hold only under an additional assumption on the triangulation, which is known to be satisfied for symmetric triangulations. Superconvergence result is also proved and as a consequence, a quasi-optimal error estimate is established in the -norm. Further, two fully discrete schemes using convolution quadrature in time generated by the backward Euler and the second-order backward difference methods are analyzed, and error estimates are derived for both smooth and nonsmooth initial data. Based on a comparison of the standard Galerkin finite element solution with the FVE solution and exploiting tools for Laplace transforms with semigroup type properties of the FVE solution operator, our analysis is then extended in a unified manner to several time-fractional order evolution problems. Finally, several numerical experiments are conducted to confirm our theoretical findings.

Key words. fractional order evolution equation, subdiffusion, finite volume element method, Laplace transform, backward Euler and second-order backward difference methods, convolution quadrature, optimal error estimate, smooth and nonsmooth data.

AMS subject classifications. 65M60, 65M12, 65M15

## 1 Introduction

Let be a bounded, convex polygonal domain in with boundary , and let be a given function (initial data) defined on . We now consider the following time-fractional diffusion problem: find in such that

 u′(x,t)+∂1−αtAu(x,t)=0 in Ω×(0,T], (1.1a) u(x,t)=0 on ∂Ω×(0,T], (1.1b) u(x,0)=v(x) in Ω, (1.1c)

where , is the partial derivative of with respect to time, and is the Riemann-Liouville fractional derivative in time defined for by:

 ∂1−αtφ(t):=ddtIαφ(t):=ddt∫t0ωα(t−s)φ(s)dswithωα(t):=tα−1Γ(α). (1.2)

Here, denotes the temporal Riemann-Liouville fractional integral operator of order . This class of problems describes the model of an anomalous subdiffusion, see [9], [10] and [25].

Over the last two decades, considerable attention from both practical and theoretical point of views has been given to fractional diffusion models due to their various applications. Several numerical techniques for the problem (1.1) have been proposed with different types of spatial discretizations. The finite element (FE) method has, in particular, been given a special attention in approximating the solution of the problem (1.1), see [24, 22, 23, 26, 12, 13, 11, 2] and references, there in. Most recently, a FVE method is analyzed in [14] and a prior error estimates with respect to data regularity have been derived.

Although the numerical study of (1.1) has been discussed in a large number of papers, optimal error estimates with respect to the smoothness of the solution expressed through initial data have been established only in few papers recently. This is due to the presence of time-fractional derivative, and hence, deriving sharp error bounds under reasonable regularity assumptions on the exact solution has become a challenging task.

To motivate our results, we begin by recalling some facts on the spatially semidiscrete standard Galerkin FE method for the problem (1.1) in the piecewise FE element space

 Vh={χ∈C0(¯¯¯¯Ω):χ|Kis linear for % all K∈Thandχ|∂Ω=0},

where is a family of regular triangulations of the domain into triangles with denoting the maximum diameter of the triangles . With denoting the bilinear form associated with the operator , and the inner product in , the semidiscrete Galerkin FE method is to seek satisfying

 (u′h,χ)+a(∂1−αtuh,χ)=0∀χ∈Vh,t∈(0,T],uh(0)=vh, (1.3)

where and is an approximation of the initial data . Upon introducing the discrete operator defined by

 (Ahψ,χ)=(∇ψ,∇χ)∀ψ,χ∈Vh,

the semidiscrete FE scheme (1.3) is rewritten in an operator form as

 u′h(t)+∂1−αtAhuh(t)=0,t>0,uh(0)=vh. (1.4)

In [22], McLean and Thomée have established the following estimate for the Galerkin FE approximation to (1.1): with , there holds for

 ∥uh(t)−u(t)∥≤Ch2t−α(2−q)/2|v|q,0≤q≤2, (1.5)

where is the -norm of and is a weighted norm defined on the space to be described in Section 2. Here, is the -projection given by : for all . For a smooth initial data, that is, , the estimate (1.5) is still valid for the initial approximation where is the standard Ritz projection defined by the relation: for all . The estimate (1.5) extends results obtained for the standard parabolic problem, i.e, , which has been thoroughly studied, see [27]. In the recent work [2], an approach based on Laplace transform and semigroup type theory has been exploited to derive a priori error estimates of the type (1.5), and most recently, a delicate energy analysis has been developed in [15] to obtain similar estimates.

Regarding the optimal estimate in the gradient norm, the following result

 ∥∇(uh(t)−u(t))∥≤Cht−α(2−q)/2|v|q,0≤q≤2, (1.6)

holds with on quasi-uniform meshes. For the cases , one can also choose . However, without the quasi-uniformity assumption on the mesh, the estimate (1.6) remains valid only for , see [15].

Optimal convergence rate up to a logarithmic factor in the stronger -norm has been derived in [23, 15]. While in [23], Laplace transform technique combined with semigroup type theoretic approach is used to derive maximum norm estimates, in [15] a novel energy argument combined with Sobolev inequality for D-problems is employed to establish, under quasi-uniformity assumption on the mesh, the following -error estimate for and

 ∥u(t)−uh(t)∥L∞(Ω)≤C|lnh|52h2t−α(3−q)/2(|v|q+∥v∥L∞(Ω)),1≤q≤2. (1.7)

In this article, we discuss the error analysis of the approximate solution satisfying the following FVE method:

 (¯u′h,χ)h+a(∂1−αt¯uh,χ)=0∀χ∈Vh,t∈(0,T],¯uh(0)=vh, (1.8)

where is a discrete inner product on to be defined in Section 3. Here, one of our objective is to establish the analogous of estimates (1.5) and (1.6) for the solution of the FVE semidiscrete problem (1.8), namely; with the appropriate choices of ,

 ∥¯uh(t)−u(t)∥+h∥∇(¯uh(t)−u(t))∥≤Ch2t−α(2−q)/2|v|q,0≤q≤2. (1.9)

We shall derive this estimate for in Section 4.1 and for in Section 4.2. For the latter case, we are only able to prove the a priori estimate under an additional hypothesis on , which is known to be satisfied for symmetric triangulations. Without any such condition, only sub-optimal order convergence is obtained, which is similar to the result proved in [5] for linear parabolic problems. For the stronger -norm, a quasi-optimal error estimate analogous to (1.7) is established for .

Our analysis provides improvements of earlier results in [14], where the initial data is required to be in with . Unlike the classical FE error analysis in which an intermediate projection, usually, a Ritz projection, is introduced to derive optimal error estimates, our approach, here, shall combine the error estimates for the standard Galerkin FE solution stated above with new bounds for the difference . A similar idea has been used in [4] and [5] for the approximation of the standard parabolic problem by the lumped mass FE method and the FVE method, respectively, leading to an improvement of their earlier results in [3].

Our second objective is to analyze two fully discrete schemes for the semidiscrete problem (1.8) based on convolution quadrature in time generated by the backward Euler and the second-order backward difference methods. Error estimates with respect to the data regularity are provided in Theorems 5.1 and 5.2. For instance, it is shown that the discrete solution obtained by the backward Euler method with a time step size satisfies the following a priori error estimate

 ∥Unh−¯uh(tn)∥≤C(τ−1+αq/2+h2t−α(1−q/2)n)|v|q,q=0,1,2.

When an additional restriction on the triangulation is imposed. A similar type of error bound is shown to hold for the second-order backward difference scheme in Subsection 5.2.

Our third objective is to generalize our results on FVE method for both smooth and nonsmooth initial data to other classes of fractional order evolution equations in Section 6. Say for example, we can extend our FVE analysis to the following class of time-fractional problems:

 u′(x,t)+JαAu(x,t)=0 in Ω×(0,T], (1.10)

with homogeneous Dirichlet boundary conditions and initial condition for When this class of problems is known as fractional diffusion-wave equation or evolution equation with positive memory, see [20, 22] and references, therein. The case corresponds to the PIDE with singular kernel, refer to [21]. Now if then this class of problems is known as the Rayleigh-Stokes problems for generalized second grade fluid, see [2]. Even our FVE analysis can be directly applied to the following time-fractional order diffusion problem:

 C∂αtu(x,t)+Au(x,t)=0, (1.11)

where is the fractional Caputo derivative of order For the semidiscrete FE analysis of (1.11), we refer to Jin et al. [12]. The unifying analysis of all these classes of evolution problems is based on comparing the FVE solution with the corresponding FE solution and exploiting the Laplace transform technique along with semigroup type properties of the FVE solution operator.

The rest of the paper is organized as follows. In the next section, we introduce notation, recall the solution representation for the continuous problem (1.1) and some smoothing properties of the solution operator, which play an important role in our subsequent error analysis. Section 3 deals with a brief description of the spatially semidiscrete FVE scheme and their properties. In Section 4, we derive error estimates for the semidiscrete FVE scheme for smooth and nonsmooth initial data , in Subsections 4.1 and 4.2. For , i.e., , we show an optimal error bound under an additional assumption on the triangulation. Superconveregence result is proved in Subsection 4.3 and as a consequence, a quasi-optimal error estimate is established in the -norm. In Section 5, two fully discrete schemes based on convolution quadrature approximation of the fractional derivative are presented and error estimates are established. Section 6 focuses on possible generalization of the present FVE error analysis to various types of time-fractional evolution problems. Finally, in Section 7, we present numerical results to confirm our theoretical findings.

Throughout the paper, denotes a generic positive constant that may depend on and , but is independent of the spatial mesh element size .

## 2 Representation of exact solution and properties

We first introduce some notations. Let be the Dirichlet eigenpairs of the selfadjoint and positive definite operator , with being an orthonormal basis in . For , we denote by the Hilbert space induced by the norm

 |v|2r=∥Ar/2v∥2=∞∑j=1λrj(v,ϕj)2,

with being the inner product on . Then, it follows that , see [27, Lemma 3.1]. In particular, is the norm on , is also the norm on and is the equivalent norm in . Note that , , form a Hilbert scale of interpolation spaces. Motivated by this, we denote by the norm on the interpolation scale between and for in the interval . Then, the and norms are equivalent for any for by interpolation.

For and , we introduce the contour defined by

 Γθ,δ={ρe±iθ:ρ≥δ}∪{δeiψ:|ψ|≤θ},

oriented with an increasing imaginary part. Further, we denote by the sector

 Σθ={z∈C,z≠0,|argz|<θ}.

For , it is clear that as . Since the operator is selfadjoint and positive definite, its resolvent satisfies the bound

 ∥(zαI+A)−1∥≤Mθ|z|−α∀z∈Σθ, (2.1)

where . We now make use of the Laplace transform of the solution defined by

 ^u(z,x)=∫∞0e−ztu(t,x)dt.

The boundary condition on transforms into on . Taking Laplace transforms in (1.1a), we, then, arrive at

 (zI+z1−αA)^u(z)=v, (2.2)

and hence,

 ^u(z)=^E(z)v,^E(z)=zα−1(zαI+A)−1. (2.3)

In view of (2.1) and (2.3), satisfies the following bound

 ∥^E(z)∥≤Mθ|z|−1∀z∈Σθ. (2.4)

From (2.3), the Laplace inversion formula yields an integral representation for the solution of (1.1) as

 u(t)=12πi∫Cezt^E(z)vdz,t>0, (2.5)

where the contour of integration , known as Bromwich contour, is any line in the right-half plane parallel to the imaginary axis and with Im increasing. Since is analytic in and satisfies the bound (2.4), the path of integration may, therefore, be deformed into the curve so that the integrand has an exponential decay property.

In the next lemma, we present some smoothing properties of the operator which play a key role in our error analysis. The estimates are proved for instance in [7, Lemma 2.2]. Note that the first estimate (2.6) given below is obtained by interpolation technique.

###### Lemma 2.1.

The following estimates hold:

 ∥A^E(z)χ∥≤Cθ|z|α(1−p/2)−1|χ|p∀z∈Σθ,0≤p≤2, (2.6)
 ∥∇^E(z)χ∥≤Cθ|z|α/2−1∥χ∥∀z∈Σθ, (2.7)

where depends only on .

In the next section, we introduce the semidiscrete finite volume element scheme.

## 3 Semidiscrete FVE scheme and its properties

To describe the finite volume element formulation, we first introduce the dual mesh on the domain . Let be the set of nodes or vertices, that is,

 Nh:={Pi:Pi   is a vertex of the element  K∈Th and Pi∈¯¯¯¯Ω}

and let be the set of interior nodes in Further, let be the dual mesh associated with the primary mesh which is defined as follows. With as an interior node of the triangulation let be its adjacent nodes (see, Figure  1 with ). Let denote the midpoints of and let be the barycenters of the triangle with . The control volume is constructed by joining successively . With as the nodes of let be the set of all dual nodes . For a boundary node , the control volume is shown in Figure  1. Note that the union of the control volumes forms a partition of .

The dual volume element space on the dual mesh is defined as

 V∗h={χ∈L2(Ω):χ|K∗P0is constant% for allK∗P0∈T∗handχ|∂Ω=0}.

The semidiscrete FVE formulation for (1.1) is to seek such that

 (¯u′h,χ)+ah(∂1−αt¯uh,χ)=0∀χ∈V∗h,t>0,¯uh(0)=vh, (3.8)

where the bilinear form is defined by

 ah(ψ,χ)=−∑Pi∈N0hχ(Pi)∫∂K∗Pi∇ψ⋅nds∀ψ∈Vh,χ∈V∗h (3.9)

with denoting the outward unit normal to the boundary of the control volume . For and , a use of Green’s formula yields

 (Aw,χ)=ah(w,χ).

To rewrite the Petrov-Galerkin method (3.8) as a Galerkin method in , we introduce the interpolation operator by

 Π∗hχ=∑Pi∈N0hχ(Pi)ηi(x),

where is the characteristic function of the control volume . The operator is selfadjoint and positive definite, see [6], and hence, the following relation

 (ψ,χ)h=(ψ,Π∗hχ)∀ψ,χ∈Vh

defines an inner product on . Also, the corresponding norm is equivalent to the -norm on , uniformly in , see [16]. Furthermore, from the following identity [1, 8]

 ah(χ,Π∗hv)=(∇χ,∇v)∀χ,v∈Vh,

the bilinear form is symmetric and for .

We now introduce the discrete operator corresponding to the inner product by

 (¯Ahψ,χ)h=(∇ψ,∇χ)∀ψ,χ∈Vh.

Then, the FVE method (1.8) is written in an operator form as

 ¯u′h(t)+∂1−αt¯Ah¯uh(t)=0,t>0,¯uh(0)=vh. (3.10)

An appropriate modification of arguments in [5, 12] yields the following discrete analogous of Lemma 2.1 and therefore, we skip the proof.

###### Lemma 3.1.

Let . With , the following estimates hold:

 ∥¯Ah^Eh(z)χ∥≤Cθ|z|α(1−p/2)−1∥¯Ap/2hχ∥∀z∈Σθ,0≤p≤2, (3.11)
 |^Eh(z)χ|1≤Cθ|z|α/2−1∥χ∥∀z∈Σθ, (3.12)

where is independent of the mesh size .

Moreover, an analogous of Lemma 3.1 holds for , when we replace in Lemma 3.1 by

## 4 Error analysis

This section deals with a priori optimal error estimates for the semidiscrete FVE scheme (1.8) with initial data , . To do so, we first introduce the quadrature error defined by

 (∇Qhχ,∇ψ)=ϵh(χ,ψ):=(χ,ψ)h−(χ,ψ)∀ψ∈Vh. (4.1)

The operator , introduced in [4] for the lumped mass FE element, represents the quadrature error in a special way. It satisfies the following error estimates, see [4, 5].

###### Lemma 4.1.

Let be defined by (4.1). Then, there holds

 ∥∇Qhχ∥+h∥¯AhQhχ∥≤Chp+1∥∇pχ∥∀χ∈Vh,p=0,1. (4.2)

Note that, by Lemma 4.1, and without additional assumptions on the mesh, the following estimate holds:

 ∥Qhχ∥≤C∥∇Qhχ∥≤Ch∥χ∥∀χ∈Vh.

This estimate cannot be improved in general, see [4, 5] for some counter examples. However, on some special meshes, one can derive a better approximation. For instance, if the mesh is symmetric (see [4, 5] for the definition and examples), the operator is shown to satisfy

 ∥Qhχ∥≤Ch2∥χ∥∀χ∈Vh. (4.3)

To derive optimal error estimates for the FVE solution , we split the error into , where and being the standard Galerkin FE solution. Then, from the definitions of , and , satisfies

 ξt(t)+∂1−αt¯Ahξ(t)=−¯AhQhuht(t),t>0,ξ(0)=0. (4.4)

### 4.1 Error estimates for smooth initial data

In the following theorem, optimal error estimates are derived for smooth initial data with

###### Theorem 4.1.

Let and be the solutions of and , respectively, with for and . Then, there is a positive constant , independent of such that

 ∥¯uh(t)−u(t)∥+h∥∇(¯uh(t)−u(t))∥≤Ct−α(2−q)/2h2|v|q,t>0. (4.5)
###### Proof.

Since the estimates for are given in (1.5) and (1.6), it is sufficient to show

 ∥ξ(t)∥+h∥∇ξ(t)∥≤Ct−α(2−q)/2h2|v|q,q∈[1,2]. (4.6)

By taking Laplace transforms in (4.4) and following the analysis in Section 2, we represent by

 ξ(t)=−12πi∫Γezt^Eh(z)¯AhQhˆuht(z)dz. (4.7)

Here and also throughout this article, is the particular contour chosen as with . From (4.7), it follows that

 ∥ξ(t)∥+h∥∇ξ(t)∥≤12π∫Γ|ezt|(∥^Eh(z)¯AhQhˆuht(z)∥+h∥∇^Eh(z)¯AhQhˆuht(z)∥)|dz|. (4.8)

To complete the proof of the estimate, we need to compute the terms under the integral sign on the right of side of (4.8). Now, we discuss two cases for and separately.

When , that is, apply (3.11) with and (3.12) in Lemma 3.1 to obtain

 ∥^Eh(z)¯AhQhˆuht(z)∥≤C|z|α/2−1∥∇Qhˆuht(z)∥, (4.9)

and

 ∥∇^Eh(z)¯AhQhˆuht(z)∥≤C|z|α/2−1∥¯AhQhˆuht(z)∥. (4.10)

Then, by (4.2), it follows that

 ∥^Eh(z)¯AhQhˆuht(z)∥+h∥∇^Eh(z)¯AhQhˆuht(z)∥≤Ch2|z|α/2−1∥∇ˆuht(z)∥. (4.11)

Since

 ˆuht(z)=−z1−αAh^uh(z)=−z1−αAh^Fh(z)vh,

an estimate analogous to (3.12) yields

 ∥∇ˆuht(z)∥=|z|1−α∥∇^Fh(z)Ahvh∥≤C|z|1−α||z|α/2−1∥Ahvh∥≤C|z|−α/2∥Ahvh∥. (4.12)

On substitution of (4.11) and (4.12) in (4.8), we use (4.7) to obtain

 ∥ξ(t)∥+h∥∇ξ(t)∥ ≤ Ch2(∫Γ|ezt||z|−1|dz|)∥Ahvh∥ (4.13) ≤ Ch2(∫∞1/teρtcosθρ−1dρ+∫θ−θecosψdψ)∥Ahvh∥ ≤ Ch2∥Ahvh∥.

Now, by the identity , we have

 ∥AhRhv∥=∥PhAv∥≤∥Av∥=|v|2,

which shows the estimate (4.6) for

For the case , that is, consider (4.11) and the identity

 ˆuht(z)=z^uh(z)−vh

to obtain using (2.4)

 ∥∇^uht(z)∥=∥∇(z^Fh(z)vh−vh)∥≤(M+1)∥∇vh∥. (4.14)

From the estimate (4.8), using (4.11) and (4.14) with we deduce that

 ∥ξ(t)∥+h∥∇ξ(t)∥ ≤ Ch2(∫Γ|ezt|z|α/2−1|dz|)|v|1 ≤ Ch2(∫∞1/teρtcosθρα/2−1dρ+∫θ−θecosψt−α/2dψ)|v|1 ≤ Ct−α/2h2|v|1.

This completes the proof for the case

Since estimates for and are known, then interpolation technique provides result for This concludes the rest of the proof. ∎

###### Remark 4.1.

Note that the estimate (4.5) in Theorem 4.1 remains valid when . Indeed, for , let denote the solution of (1.8) with Then satisfies

 ζt+∂1−αt¯Ahζ=0,t>0,ζ(0)=Phv−Rhv.

Since

 ζ(t)=−12πi∫Γezt^Eh(z)(Phv−Rhv)dz,

we deduce

 ∥ζ(t)∥≤C∥Phv−Rhv∥∫Γ|ezt||z|−1|dz|≤Ch2|v|2.

Thus, the estimate (4.5) with follows by the triangle inequality. If the inverse inequality holds, which is the case if the mesh is quasi-uniform, then the estimate in the gradient norm follows directly for .

If the -projection operator is stable in , i.e., , then the estimate (4.5) holds for the case and the choice . A sufficient condition for such stability of is the quasi-uniformity of the mesh. Now, by interpolation the estimate (4.5) holds for and .

### 4.2 Error estimates for nonsmooth initial data

In this subsection, we establish optimal error estimates for the semidiscrete FVE scheme (1.8) for nonsmooth initial data .

###### Theorem 4.2.

Let and be the solution of and , respectively, with and . Then, there exists a positive constant , independent of such that

 ∥¯uh(t)−u(t)∥+∥∇(¯uh(t)−u(t))∥≤Cht−α∥v∥,t>0. (4.15)

Furthermore, if the quadrature error operator satisfies (4.3), then the following optimal error estimate holds:

 ∥¯uh(t)−u(t)∥≤Ch2t−α∥v∥,t>0. (4.16)
###### Proof.

As before, it is sufficient to prove estimates for . We first apply (3.11) with to arrive at

 ∥^Eh(z)¯AhQhˆuht∥≤C|z|α−1∥Qhˆuht∥.

Then, the following bound follows from the integral representation (4.7):

 ∥ξ(t)∥≤C∫Γ|ezt||z|α−1∥Qhˆuht(z)∥|dz|. (4.17)

To estimate the gradient of , we note that

 ∥∇^Eh(z)¯AhQhˆuht∥≤C|z|α−1∥∇Qhˆuht∥,

and hence,

 ∥∇ξ(t)∥≤C∫Γ|ezt||z|α−1∥∇Qhˆuht(z)∥|dz|. (4.18)

Note that holds on a general mesh, and by (4.2). Since by (2.4), a substitution into (4.17) and (4.18) yields the first estimate (4.15). Finally, if (4.3) holds, then (4.16) follows immediately from (4.17), which completes the proof. ∎

### 4.3 L∞(Ω)-error estimates

In the following, we obtain a superconvergence result for the gradient of in the -norm. As a consequence, assuming and the quasi-uniformity on the mesh, a quasi-optimal error estimate in the stronger -norm is derived for the semidiscrete FVE solution . We first prove the following Lemma by refining some of the estimates derived in the proofs of Theorem 4.1.

###### Lemma 4.2.

For , and with , there is a positive constant independent of such that

 ∥∇ξ(t)∥≤Ch2t−α(3−q)/2|v|q,t>0.

The estimate is still valid for but with quasi-uniform assumption on the mesh.

###### Proof.

By using bounds (3.11) and (4.2), we obtain instead of (4.10) the following estimate

 ∥∇^Eh(z)¯AhQhˆuht(z)∥≤C|z|α−1∥∇Qhˆuht(z)∥≤Ch2|z|α−1∥∇ˆuht(z)∥.

Since