Optimal spatial H1 -norm analysis of a finite element method for a time-fractional diffusion equation

Journal of Computational and Applied Mathematics 367 (2020) 112435

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Optimal spatial H 1 -norm analysis of a finite element method for a time-fractional diffusion equation Chaobao Huang a,b , Martin Stynes b , a b

∗

School of Mathematics and Quantitative Economics, Shandong University of Finance and Economics, Jinan 250014, PR China Applied and Computational Mathematics Division, Beijing Computational Science Research Center, Beijing 100193, PR China

article

info

Article history: Received 3 January 2019 Received in revised form 12 August 2019 Keywords: Fractional diffusion equation Finite element method Optimal error in H 1 (Ω ) Graded mesh

a b s t r a c t A time-fractional initial–boundary value problem Dαt u − ∆u = f , where Dαt is a Caputo fractional derivative of order α ∈ (0, 1), is considered on the space–time domain Ω × [0, T ], where Ω ⊂ Rd (d ≥ 1) is a bounded Lipschitz domain. Typical solutions u(x, t) of such problems have components that behave like a multiple of t α as t → 0+ , so the integer-order temporal derivatives of u blow up at t = 0. The numerical method of the paper uses a standard finite element method in space on a quasiuniform mesh and considers both the L1 discretisation and Alikhanov’s L2-1σ discretisation of the Caputo derivative on suitably graded temporal meshes. Optimal error bounds in H 1 (Ω ) are proved; no previous analysis of a discretisation of this problem using finite elements in space has established such a bound. Furthermore, the optimal grading of the temporal mesh can be deduced from our analysis. Numerical experiments show that our error bounds are sharp. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Numerical methods for time-fractional initial–boundary value problems (IBVPs) have received a huge amount of attention in the research literature. In this paper we shall consider those fractional-derivative problems that generalise classical parabolic IBVPs such as the heat equation; that is, where the classical time derivative ut is replaced by a Caputo fractional derivative Dαt u of order α with 0 < α < 1. For this class of problems, several numerical methods have been proposed and analysed—see, e.g., [1–3] and their references. Many of these papers focus on the discretisation and analysis of the Caputo time derivative; the error in the spatial discretisation has received perhaps less attention, but it is our main topic here. As a fairly typical example of what has been done, consider [4], where the authors use the well-known L1 discretisation in time, combined with a direct discontinuous Galerkin finite element method in the spatial domain Ω . Optimal L2 (Ω ) errors in space are proved, but [4, Remark 4.1] the error estimate in H 1 (Ω ) is suboptimal unless one is sufficiently far from the initial time t = 0. A similar loss of optimality has appeared in all numerical analyses that investigated H 1 (Ω ) errors for this class of fractional-derivative problems, on uniform and graded temporal meshes; see, e.g., [5, Remark 5.6]. It is a disappointing weakness of these analyses, since for classical parabolic IBVPs the derivation of optimal error bounds in H 1 (Ω ) for all time t > 0 is standard—see for instance [6, Theorem 1.3]. Very recently, a way around this difficulty was found in [7], where the authors consider two different discretisations of the Caputo time derivative, while for the spatial discretisation [7, p. 2] ‘‘for simplicity, we only consider the finite difference ∗ Corresponding author. E-mail address: [email protected] (M. Stynes). https://doi.org/10.1016/j.cam.2019.112435 0377-0427/© 2019 Elsevier B.V. All rights reserved.

2

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

method’’ and ‘‘the theoretical results . . . are extendable for some other spatial discretisation such as the spectral method’’. It is surprising that in [7] no mention is made of finite element discretisations in space, even though such discretisations form the most natural setting for H 1 (Ω ) error analysis. The purpose of our paper is to fill this gap by showing how the optimal H 1 (Ω )-norm analysis of [7] can be extended to standard finite element discretisations in space. We consider two discretisations of the Caputo fractional derivative on a graded temporal mesh (with arbitrary userchosen mesh grading r ≥ 1): the well-known L1 discretisation, and Alikhanov’s higher-order L2-1σ discretisation. Our analysis proves optimal H 1 (Ω )-norm errors for the solutions computed by both discretisations, and also obtains the values of r that will deliver the optimal temporal error in each discretisation: r ≥ (2 − α )/α gives order 2 − α convergence for the L1 scheme, and r ≥ 2/α gives order 2 for the L2-1σ scheme. The paper is structured as follows. Section 2 presents the initial–boundary value problem and the properties of its solution. In Section 3 we describe the finite element discretisation of the spatial component of the differential operator, which will be used throughout the paper, and we define the graded temporal mesh on which we specify the L1 discretisation of the Caputo temporal derivative. The H 1 (Ω )-norm stability of this numerical method is shown in Section 4; then in Section 5 we derive an optimal H 1 (Ω ) convergence bound for the computed solution. Next, in Section 6, we examine a higher-order discretisation of the Caputo derivative on our graded mesh: Alikhanov’s L2-1σ scheme. Here we again show the H 1 (Ω )-norm stability of the method, but our approach differs from Section 4. Then we prove an optimal H 1 (Ω ) convergence bound. Finally, numerical results showing the sharpness of our theoretical results are presented in Section 7. Notation. C denotes a generic constant that depends on the data of the problem but is independent of the mesh parameters N and h; it can take different values in different places in the paper. The constant K , defined in (10) is a fixed constant. Set N = {1, 2, 3, . . .} and N0 = {0, 1, 2, . . .}. We write ∥ · ∥ for the norm in L2 (Ω ). For each q ∈ N, the notation H q (Ω ) is used for the standard Sobolev space with its associated norm ∥ · ∥q and seminorm |·|q , and H01 (Ω ) is the subspace of H 1 (Ω ) comprising those functions whose traces vanish on the boundary ∂ Ω . Furthermore, we write ∥ · ∥L∞ (H q ) = sup0 0, consider the time-fractional initial–boundary value problem Dαt u − ∆u = f on Q := Ω × (0, T ],

(1a)

u(x, 0) = u0 (x) for x ∈ Ω ,

(1b)

u|∂ Ω = 0 for 0 < t ≤ T .

(1c)

Here we assume that the spatial domain Ω ⊂ R (where d ∈ {1, 2, 3}) is bounded, with a Lipschitz-continuous boundary ∂ Ω . In (1a), Dαt is the standard Caputo fractional derivative operator, which is defined for all functions v (x, ·) that are absolutely continuous on [0, T ] by d

t

∫

1

Dαt v (x, t) =

Γ (1 − α )

(t − s)−α 0

∂v (x, s) ds. ∂s

(2)

¯ × [0, T ]. We assume that u0 ∈ C (Ω ¯ ) and f ∈ C (Q¯ ). Further regularity ¯ denote the closure of Ω in Rd . Set Q¯ := Ω Let Ω assumptions will be made later for these functions. Let {(λi , ψi ) : i = 1, 2, . . . } be the eigenvalues and eigenfunctions for the boundary value problem ∆ψi = λi ψi on Ω ,

ψi = 0 on ∂ Ω ,

where the eigenfunctions are normalised by requiring ∥ψi ∥ = 1 for all i. From the theory of sectorial operators [8,9], the fractional power ∆γ of the operator ∆ is defined for each γ ∈ R with domain

{ D(∆γ ) :=

g ∈ L2 ( Ω ) :

∞ ∑

2γ i

}

λ |(g , ψi )|2 < ∞ , and we set ∥g ∥∆γ :=

i=1

(

∞ ∑

2γ i

λ |(g , ψi )|2

)1/2 .

i=1

By imitating [5, Section 6.1], one can obtain the following a priori bounds on the solution of problem (1). Lemma 2.1. Let q ∈ N0 . Assume that u0 ∈ D(∆q+2 ) and ∂tl f (·, t) ∈ D(∆q ) for all t ∈ (0, T ], and that ∥u0 ∥∆q+2 + ∥∂tl f (·, t)∥∆q+1 ≤ C1 for l = 1, 2, where C1 is a constant independent of t. Then (1) has a unique solution u that satisfies (1a)–(1c) pointwise, and there exists a constant C such that

  u(·, t) ≤ C , q  l  ∂ u(·, t) ≤ C (1 + t α−l ) for l = 0, 1, 2, 3, t q

(3a) (3b)

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

3

and

  α D u(·, t) ≤ C , t q

(3c)

for all t ∈ (0, T ]. 2.1. Semidiscretisation in space To begin, we discretise (1a) only in space using a standard finite element method. Let M be a positive integer. Partition Ω by a quasiuniform mesh of M elements {Km : m = 1, . . . , M }. Set hm = diam(Km ) for each m and h = max1≤m≤M {hm }. Then the weak formulation of (1) is: Find u(·, t) ∈ H01 (Ω ) for each t ∈ (0, T ], such that

{

(Dαt u, v ) + (∇ u, ∇v ) = (f , v ) ∀ v ∈ H01 (Ω ), (u(0, ·), v (·)) = (u0 , v ) ∀ v ∈ H01 (Ω ),

(4)

where we have used the L2 (Ω ) inner products

∫ ∑ d ∂v ∂w (v, w ) := vw dx and (∇v, ∇w) := dx. ∂ xj ∂ xj Ω Ω ∫

j=1

For each k ∈ N and m = 1, 2, . . . , M, let P (Km ) denote the space of polynomials of total degree at most k defined on Km . Define the finite element spaces k

{

Vh := vh ∈ H01 (Ω ) : vh ⏐K ∈ P k (Km ), m = 1, 2, . . . , M

⏐

}

and V0h = {vh ∈ Vh : vh |∂ Ω = 0} .

m

Based on (4), our semi-discrete finite element method (FEM) is: Find uh (·, t) ∈ V0h for each t ∈ (0, T ], such that

{

(Dαt uh , vh ) + (∇ uh , ∇vh ) = (f , vh ) ∀ vh ∈ V0h , (uh (0, ·), vh (·)) = (u0 , vh ) ∀ vh ∈ V0h .

(5)

3. Temporal graded meshes; the L1 FEM To obtain a final fully-discrete numerical method, one must discretise (5) in time. We shall investigate two temporal discretisations, but both use graded meshes of the class that we now describe. Let N be a positive integer. Set tn = T (n/N)r for n = 0, 1, . . . , N, where the temporal mesh grading constant r ≥ 1 is chosen by the user. If r = 1, then the mesh is uniform. Similar graded meshes have been used in several papers that solve (1) numerically, e.g., [3]. The theoretical convergence results of Sections 5 and 6 will show the influence of r and will show how to choose it in an optimal way for each of the temporal discretisations that we consider. Set τn = tn − tn−1 for n = 0, 1, . . . , N. Note that

τn+1 = T

(

n+1 N

)r −T

( n )r N

≤ CTN −r nr −1 for n = 1, 2, . . . , N − 1.

(6)

3.1. The L1 FEM For n ≥ 1, we approximate the Caputo fractional derivative Dαt u(x, tn ) of (2) by the well-known L1 formula: Dαt u(x, tn ) ≈ DαN un :=

∫ n−1 i+1 ∑ u − ui

1

Γ (1 − α )

i=0

τi+1

ti+1

(tn − s)−α ds

ti

n−1 i+1 ∑ ] u − ui [ = (tn − ti )1−α − (tn − ti+1 )1−α Γ (2 − α ) τi+1

1

i=0

=

[

dn,1

where dn,i := (tn − tn−i )1−α shows easily that

dn,n

1

n−1 ∑

un−i (dn,i+1 − dn,i ), (7) Γ (2 − α ) Γ (2 − α ) i=1 ] − (tn − tn−i+1 )1−α /τn−i+1 for i = 1, . . . , n. Note that dn,1 = τn−α . The mean value theorem

Γ (2 − α )

un −

dn,i+1 < dn,i for 0 ≤ i ≤ n − 1 ≤ N − 1.

u0 +

(8)

4

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

Now use the L1 scheme (7) to discretise (5) in time. The fully discrete method (which we call the L1 FEM) is: find unh ∈ V0h for n = 0, 1, . . . , N such that

{

(DαN unh , vh ) + (∇ unh , ∇vh ) = (f n , vh ) for n = 1, . . . , N and all vh ∈ V0h ,

(u0h , vh ) = (u0 , vh ) ∀vh ∈ V0h ,

(9)

where f n (·) := f (·, tn ). To facilitate our analysis, we now introduce three operators that are often used in finite element analyses of time-dependent problems [6]. First, define the L2 projector Ph : L2 (Ω ) → V0h by (Ph w, vh ) = (w, vh ) ∀ vh ∈ V0h . Then [10] there exists a constant K , which is independent of g and of the mesh, such that

∥∇ Ph v∥ ≤ K ∥∇v∥ for all v ∈ H01 (Ω ). Next, define the Ritz projector Rh :

H01 (Ω )

(10)

→ V0h by

(∇ Rh w, ∇vh ) = (∇w, ∇vh ) ∀ vh ∈ V0h . Since V0h ⊂ H01 (Ω ) is the space of piecewise polynomials of degree at most k, it is well known [6, Lemma 1.1] that

∥w − Rh w∥ + h∥w − Rh w∥1 ≤ Chk+1 |w|k+1 ∀ w ∈ H k+1 (Ω ) ∩ H01 (Ω ).

(11)

Finally, define the discrete Laplacian ∆h : V0h → V0h by [6, (1.33)] (∆h v, w ) = −(∇v, ∇w ) ∀ v, w ∈ V0h .

(12)

These three operators are related [6, (1.34)] by

∆ h R h v = P h ∆ v ∀ v ∈ H 2 (Ω ) .

(13)

Using (12), the L1 FEM (9) takes the form: find

unh

∈ V0h for n = 0, 1, . . . , N such that

{

(DαN unh , vh ) − (∆h unh , vh ) = (Ph f n , vh ) for n = 1, . . . , N and all vh ∈ V0h , (u0h , vh ) = (Ph u0 , vh ) ∀ vh ∈ V0h ,

for n = 1, . . . , N. But since DαN unh , ∆h unh and Ph f n all lie in V0h , this formulation of our L1 FEM can be written more succinctly as: find unh ∈ V0h for n = 0, 1, . . . , N such that

{

DαN unh − ∆h unh = Ph f n for n = 1, . . . , N ,

(14)

u0h = Ph u0 . 4. H 1 Norm stability of the L1 FEM Our stability result for the L1 FEM will be presented in a general framework resembling that of [3, Section 4]. Suppose that for n = 0, 1, . . . , N, the function µn ∈ V0h satisfies

{

DαN µn − ∆h µn = Ph g n for n = 1, . . . , N ,

(15)

µ0 = Ph u0 .

By imitating part of the proof of [11, Lemma 4.2] (see also the proof of [5, Lemma 2.1]), we obtain the following coercivity property of the L1 scheme, which is important in our later analysis. Lemma 4.1.

Let the functions v j = v (·, t j ) be in L2 (Ω ) for j = 0, 1, . . . , N. Then the discrete L1 scheme satisfies

(DαN v n , v n ) ≥ DαN ∥v n ∥ ∥v n ∥ for n = 1, 2, . . . , N .

(

)

Proof. Let n ∈ {1, 2, . . . , N }. The definition of DαN v n gives (DαN v n , v n ) =

≥

dn,1

Γ (2 − α ) dn,1

(v n , v n ) −

∥v n ∥2 −

Γ (2 − α ) ( α n ) n = DN ∥v ∥ ∥v ∥,

dn,n

Γ (2 − α ) dn,n

Γ (2 − α )

(v 0 , v n ) −

∥v 0 ∥ ∥v n ∥ −

1

n−1 ∑

Γ (2 − α )

i=1

(dn,i − dn,i+1 )(v n−i , v n )

1

n−1 ∑

Γ (2 − α )

i=1

where we used Cauchy–Schwarz inequalities and dn,i > dn,i+1 > 0. □

(dn,i − dn,i+1 )∥v n−i ∥ ∥v n ∥

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

5

Note that the proof of Lemma 4.1 depends only the signs of the dn,i for i = 1, 2, . . . , n, so it could also be used for some other discretisations of Dαt u. The proof of the next lemma uses the key idea of [7]: multiply the discrete differential equation DαN µn − ∆h µn = Ph g n by −∆h µn . This contrasts with the related classical technique for H 1 (Ω )-norm analysis of the semidiscrete problem in the proof of [6, Theorem 1.3], where the discrete differential equation is multiplied by (µn )t . Lemma 4.2.

The solution µn of (15) satisfies

⎡ ∥∇µn ∥ ≤ τnα ⎣Γ (2 − α )K ∥∇ g n ∥ + dn,n ∥∇µ0 ∥ +

⎤

n−1 ∑

(dn,i − dn,i+1 )∥∇µj ∥⎦

(16)

j=1

for n = 1, 2, . . . , N. Proof. Fix n ∈ {1, 2, . . . , N }. Multiplying (15) by −∆h µn and integrating over Ω , one has

−(DαN µn , ∆h µn ) + ∥∆h µn ∥2 = −(Ph g n , ∆h µn ). Recalling the definition (12) of ∆h and discarding the non-negative term ∥∆h µn ∥2 , we get DαN (∇µn ), ∇µn ≤ (∇ Ph g n , ∇µn ).

(

)

Lemma 4.1 and a Cauchy–Schwarz inequality now yield DαN ∥∇µn ∥ ∥∇µn ∥ ≤ ∥∇ Ph g n ∥ ∥∇µn ∥.

(

)

Cancelling ∥∇µn ∥ then invoking (10), we get DαN ∥∇µn ∥ ≤ K ∥∇ g n ∥. Now using (7) to write out DαN ∥∇µn ∥ in full gives dn,1

n

Γ (2 − α )

∥∇µ ∥ ≤ K ∥∇ g ∥ +

which is equivalent to (16).

[

1

n

0

Γ (2 − α )

dn,n ∥∇µ ∥ +

n−1 ∑

]

(dn,i − dn,i+1 )∥∇µ ∥ , i

i=1

□

As in [3, (4.6)], define the positive real numbers θn,j , for n = 1, 2, . . . , N and j = 1, 2, . . . , n − 1, by

θn,n = 1,

θn,j =

n−j ∑

τnα−i (dn,i − dn,i+1 )θn−i,j .

i=1

Observe that (8) implies θn,j > 0 for all n, j. Furthermore, as in [3, Lemma 4.3], for n = 1, 2, . . . , N one has

τnα

n ∑

j−η θn,j ≤

j=1

T α N −η 1−α

, provided that η ≤ r α.

(17)

By imitating the proof of [3, Lemma 4.2], one can deduce the following stability bound for (15) from Lemma 4.2:

∥∇µn ∥ ≤ ∥∇µ0 ∥ + K τnα Γ (2 − α )

n ∑

θn,j ∥∇ g j ∥ for n = 1, 2, . . . , N .

j=1

Now we get our main stability result. Theorem 4.3 (H 1 (Ω )-stability of the L1 FEM). Let µn be the solution of (15). Then

∥∇µn ∥ ≤ ∥∇µ0 ∥ +

KT α Γ (2 − α )

Proof. From (17) one has τnα

1−α

∑n

j=1

max ∥∇ g j ∥ for n = 1, 2, . . . , N .

1≤j≤n

θn,j ≤ T α /(1 − α ). The result now follows immediately from (18). □

5. Error analysis of the L1 FEM From (2) and (9), the temporal truncation error in our method is

ϕ n (x) := Dαt u(x, tn ) − DαN u(x, tn ) [ ] ∫ ti+1 n−1 ∑ 1 ∂u u(x, ti+1 ) − u(x, tk ) = (tn − s)−α (x, s) − ds. Γ (1 − α ) ti ∂s τi+1 i=0

(18)

6

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

For all (x, tn ) ∈ Q one has  n ϕ  ≤ Cn− min{2−α,r α} . 1

Lemma 5.1.

Proof. One can obtain this bound similarly to [3, Lemma 5.2], using the bound (3b) of Lemma 2.1. □ Let un and unh be the solutions of (4) and (9) respectively at time t = tn for n = 0, 1, . . . , N. To facilitate the error analysis, we write un − unh = (Rh un − unh ) − (Rh un − un ) = ζ n − ρ n ,

(19)

where ζ := Rh u − and ρ := Rh u − u . (This splitting of u − into the error of the Ritz projection and a remainder is routine in the analysis of classical IBVPs; see [6, p.8].) We can estimate ρ n immediately using (11), but the analysis of ζ n is more difficult and we consider it now. From (1a), (13) and (14), one has n

unh

n

n

n

n

n

unh

DαN ζ n − ∆h ζ n = (Rh DαN un − ∆h Rh un ) − (DαN unh − ∆h unh )

= (Rh − Ph )DαN un + Ph (DαN un − ∆un ) − Ph f n = Ph (Rh − I)DαN un + Ph (f n + ϕ n ) − Ph f n = Ph DαN ρ n + Ph ϕ n ,

where the definitions of φ n and ρ n were used. Thus DαN ζ n − ∆h ζ n = Ph (DαN ρ n + ϕ n ) for n = 1, 2, . . . , N .

(20) ∞

1

We can now prove the optimal-rate convergence of our method in L (H ). Theorem 5.2 (Error Estimate for the L1 FEM). Let un and unh be the solutions of (4) and (9), respectively. Then for n = 1, 2, . . . , N, there exists a constant C such that

( ) ∥∇ un − ∇ unh ∥ ≤ C hk + N − min{2−α,r α} .

(21)

Proof. Observe that (20) is a particular case of (15). Thus we can invoke the stability bound (18) to get

∥∇ζ n ∥ ≤ ∥∇ζ 0 ∥ + K τnα Γ (2 − α )

n ∑

θn,j ∥∇ (DαN ρ j + ϕ j )∥.

(22)

j=1

The inequality ∥∇ Rh w∥ ≤ ∥∇w∥ ∀w ∈ H01 (Ω ) follows easily from the definition of Rh . Hence

∥∇ DαN ρ j ∥ = ∥∇ DαN ρ j − ∇ Dαt ρ j + ∇ Dαt ρ j ∥ ( ) ( ) ≤ ∥∇ Dαt uj − DαN uj − ∇ Rh Dαt uj − DαN uj ∥ + ∥∇ Dαt ρ j ∥ ≤ ∥∇ϕ j ∥ + ∥∇ Rh ϕ j ∥ + Chk |Dαt uj |k+1 ≤ 2∥∇ϕ j ∥ + Chk |Dαt uj |k+1 , where we used (11). Hence ∥∇ (DαN ρ j + ϕ j )∥ ≤ 3∥∇ϕ j ∥ + Chk |Dαt uj |k+1 . Substituting this inequality into (22) and recalling Lemma 5.1 yields

∥∇ζ n ∥ ≤ ∥∇ζ 0 ∥ + K τnα Γ (2 − α )

n ∑

( ) θn,j Chk |Dαt uj |k+1 + 3∥∇ϕ j ∥

j=1

≤ ∥∇ζ 0 ∥ + CK τnα Γ (2 − α )

n ∑

[ ] θn,j hk |Dαt uj |k+1 + j− min{2−α,r α}

j=1

≤ Chk + C τnα

(

k

≤C h +N

n ∑

θn,j hk + j− min{2−α,r α}

[

j=1 ) − min{2−α,r α}

]

,

0

where we used (3c) and ∥∇ζ ∥ = ∥∇ (Rh u0 − Ph u0 )∥ = ∥∇ Ph (Rh u0 − u0 )∥ ≤ K ∥∇ρ 0 ∥ ≤ Chk by (10) and (11), then invoked (17) with η = 0 for the hk term and η = min{2 − α, r α} for the j− min{2−α,r α} term. Combining this bound and (11) with (19), we get (21). □ Corollary 5.3.

Suppose that r ≥ (2 − α )/α . Then the L1 FEM solution unh satisfies

( ) ∥un − unh ∥H 1 (Ω ) ≤ C hk + N −(2−α) for n = 0, 1, . . . , N .

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

Proof. Theorem 5.2 and r ≥ (2 − α )/α imply that ∥∇ un − ∇ unh ∥ ≤ C hk + N yields the same bound in H 1 (Ω ). □

(

) −(2−α )

7

for each n. Then a Poincaré inequality

From [3] one sees that the rate N −(2−α ) is best possible for temporal errors when the L1 scheme is used, and for H (Ω ) approximation with piecewise polynomials of degree k it is well known that hk is the optimal rate; consequently Corollary 5.3 is an optimal error estimate. 1

6. The L2-1σ FEM In this section, we continue to use our FEM for the spatial discretisation, but for the temporal discretisation we replace the L1 scheme by Alikhanov’s L2-1σ scheme [12] on the same class of graded meshes, with the aim of obtaining a higher order of convergence in the temporal discretisation variable N. We continue to use the graded mesh of Section 3. For n = 0, . . . , N − 1 and 0 ≤ σ ≤ 1, set tn+σ = tn + σ τn+1 . Now approximate the Caputo fractional derivative Dαt v in (1a) at the point tn+σ by the L2−1σ formula of Alikhanov [12]: n ∑

Dαt v (tn+σ ) ≈ δtαn+σ v := gn,n v n+1 −

(gn,j − gn,j−1 )v j for n = 0, . . . , N − 1.

(23)

j=0

Here g0,0 = τ1−1 a0,0 , gn,−1 = 0, and for n ≥ 1 one has

gn,j

⎧ −1 ⎪ ⎨τj+1 (an,0 − bn,0 ) = τj−+11 (an,j + bn,j−1 − bn,j ) ⎪ ⎩ −1 τj+1 (an,n + bn,n−1 )

if j = 0, if 1 ≤ j ≤ n − 1, if j = n,

where an,n = an , j = bn,j =

tn+σ

∫

1

Γ (1 − α )

∫

1

Γ (1 − α )

tn tj+1

(tn+σ − η)−α dη =

(tn+σ − η)−α dη

for n ≥ 1

for n ≥ 0,

and 0 ≤ j ≤ n − 1,

tj

2

1

σ 1−α τ 1−α Γ (2 − α ) n+1

∫

Γ (1 − α ) tj+2 − tj

tj+1

(tn+σ − η)−α (η − tj+1/2 ) dη

for n ≥ 1

and 0 ≤ j ≤ n − 1.

tj

Next we state four technical lemmas that will be used in our analysis. Lemma 6.1 ([12, Corollary 1], [13, Corollary 1]). Suppose 0 ≤ σ ≤ 1. For any function V (·, t) defined on the graded mesh {tj }Nj=0 , one has

(

) 1 σ V n+1 + (1 − σ )V n , δtαn+σ V ≥ δtαn+σ ∥V ∥2 for n = 0, . . . , N − 1. 2

Set

ψvσ = τ1α sup

η∈(0,t1 )

( 1−α ) η |δt v (t0 ) − v ′ (η)|

α ψvj+σ = τj3+−α 1 tj+σ

ψvj,1 = τ1α sup ψvj,s =

sup

η∈(tj ,tj+1 )

|v ′′′ (η)|

for j = 1, . . . , N − 1,

( 1−α ) η |(I2,1 v (η))′ − v ′ (η)|

η∈(0,t1 ) 2 α τj−α +1 τs (τs + τs+1 )ts

sup

η∈(ts−1 ,ts+1 )

for j = 1, . . . , N − 1,

|v ′′′ (η)|

for 2 ≤ s ≤ j ≤ N − 1.

Lemma 6.2 ([13, Lemma 7]). Let v ∈ C [0, T ] ∩ C 3 (0, T ]. Assume that |v (l) (t)| ≲ 1 + t α−l for l = 0, 1, 2, 3 and t ∈ (0, T ]. Then

ψvj+σ ≤ CN − min{r α,3−α} for j = 0, . . . , N − 1,

ψvj,s ≤ CN − min{r α,3−α} for s = 1, . . . , j, when j ≥ 1.

Lemma 6.3 ([13, Lemma 5]). Suppose 1 − α/2 ≤ σ ≤ 1. Then for any function {V j }Nj=0 , one has

|V s+1 | ≤ |V 0 | + Γ (1 − α ) max

n=0,...,s

{

tnα+σ δtαn+σ |V |

}

for s = 0, . . . , N − 1.

8

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

Lemma 6.4 ([13, Lemma 1]). Suppose σ = 1 − α/2. Then for any function v (t) ∈ C 3 (0, T ], one has α

α

|δtj+σ v − Dt v (tj+σ )| ≤

Ctj−α +σ

(

ψv

j+σ

) + max {ψv } for j = 0, . . . , N − 1, j,s

s=1,...,j

We now define the L2-1σ FEM, imitating [13]. Let unh ∈ V0h for n = 0, 1, . . . , N denote the computed solution. Set f = f (·, tn+σ ) and un,σ = σ unh+1 + (1 − σ )unh for n = 0, 1, . . . , N − 1. Use the L2-1σ scheme (23) to discretise Dαt u(tn+σ ) in (5). Then the L2-1σ FEM is: for n = 0, . . . , N − 1, solve for unh+1 from n+σ

n,σ

(δtαn+σ uh , vh ) + (∇ uh , ∇vh ) = (f n+σ , vh ) ∀ vh ∈ V0h ,

(24)

where we set u0h = Rh u0 . Similarly to (13), the L2-1σ FEM can be expressed as

δtαn+σ uh − ∆h unh,σ = Ph f n+σ for n = 0, . . . , N − 1, with

u0h

(25)

= Rh u0 .

6.1. H 1 (Ω )-stability and error bound for the L2-1σ FEM We begin by proving a stability bound (cf. Theorem 4.3). Theorem 6.5 (H 1 (Ω )-stability of the L2-1σ FEM). The L2-1σ FEM solution unh of (25) satisfies

∥∇ unh ∥2 ≤ ∥∇ u0 ∥2 + Γ (1 − α )T α max ∥f j+σ ∥2 for n = 0, 1, . . . , N − 1. 0≤j≤N −1

n,σ

Proof. Fix n ∈ {0, 1, . . . , N − 1}. Multiplying (25) by −∆h uh

−(δ

n,σ tn+σ uh ∆h uh )

α

,

+∥

∆h unh,σ

2

∥ = −(Ph f

n+σ

,

∆h unh,σ )

and integrating over Ω gives

.

Using (12), Lemma 6.1 and Young’s inequality yields 1 2

δtαn+σ ∥∇ uh ∥2 + ∥∆h unh,σ ∥2 ≤

1 2

1

∥Ph f n+σ ∥2 + ∥∆h unh,σ ∥2 . 2

But ∥Ph f n+σ ∥ ≤ ∥f n+σ ∥, so it follows that δtαn+σ ∥∇ uh ∥2 ≤ ∥f n+σ ∥2 . Now Lemma 6.3 gives

∥∇ unh ∥2 ≤ ∥∇ u0h ∥2 + Γ (1 − α ) max

j=0,...,n

as ∥∇

u0h

{

tjα+σ δtαj+σ ∥∇ uh ∥2

}

≤ ∥∇ u0 ∥2 + Γ (1 − α )T α max ∥f j+σ ∥2 , 0≤j≤N −1

∥ = ∥∇ Rh u0 ∥ ≤ ∥∇ u0 ∥. □

We can now prove an error bound for the L2-1σ FEM solution. Similarly to (19), write un − unh = (Rh un − unh ) − (Rh un − un ) = ζ n − ρ n ,

(26)

where ζ n := Rh un − unh and ρ n := Rh un − un . From (1a), (13) and (25), for n = 0, 1, . . . , N − 1 one has

δtαn+σ ζ − ∆h ζ n,σ = (Rh δtαn+σ u − ∆h Rh un,σ ) − (δtαn+σ uh − ∆h unh,σ ) = (Rh − Ph )δtαn+σ u + Ph (δtαn+σ u − ∆un,σ ) − Ph f n+σ = Ph (Rh − I)δtαn+σ u + Ph (δtαn+σ u − ∆un,σ ) − Ph (Dαt un+σ − ∆un+σ ) = Ph δtαn+σ ρ + Ph (δtαn+σ u − Dαt un+σ ) + Ph (∆un+σ − ∆un,σ ) ( ) = Ph δtαn+σ ρ + ϕ n+σ + Rn+σ ,

(27)

where ϕ n+σ := δtαn+σ u − Dαt un+σ and Rn+σ := ∆un+σ − ∆un,σ . The next result establishes convergence of the L2-1σ FEM in L∞ (H 1 ). Theorem 6.6.

Suppose σ = 1 − α/2. Let un and unh be the solutions of (4) and (25), respectively. Assume that u ∈

L (0, T ; ∩ H k+1 (Ω )), Dαt u ∈ L∞ (0, T ; H01 (Ω ) ∩ H k+1 (Ω )), and ∥∂tl u∥3 ≲ 1 + t α−l for l = 0, 1, 2, 3. Then there exists a constant C such that ∞

H01 (Ω )

max ∥∇ un − ∇ unh ∥ ≤ C N − min{r α,2} + hk .

(

1≤n≤N

)

(28)

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

9

Proof. Fix n ∈ {0, 1, . . . , N }. Multiplying (27) by −∆h ζ n,σ then integrating over Ω , one has

) ( −(δtαn+σ ζ , ∆h ζ n,σ ) + ∥∆h ζ n,σ ∥2 = − Ph (δtαn+σ ρ + ϕ n+σ + Rn+σ ), ∆h ζ n,σ . Now (12), Lemma 6.1 and a Cauchy–Schwarz inequality yield 1 2

δtαn+σ ∥∇ζ ∥2 ≤ ∥∇ Ph (δtαn+σ ρ + ϕ n+σ + Rn+σ )∥ ∥∇ζ n,σ ∥.

By (10) and a triangle inequality, we get

) ( δtαn+σ ∥∇ζ ∥2 ≤ 2K ∥∇ (δtαn+σ ρ )∥ + ∥∇ϕ n+σ ∥ + ∥∇ Rn+σ ∥ max ∥∇ζ s ∥. 1≤s≤N

Hence Lemma 6.3 and

u0h

= Rh u0 give

} { ∥∇ζ n+1 ∥2 ≤ ∥∇ζ 0 ∥2 + Γ (1 − α ) max tjα+σ δtαj+σ ∥∇ζ ∥2 j=0,...,n )} { ( α ≤ C max tj+σ ∥∇ (δtαj+σ ρ )∥ + ∥∇ϕ j+σ ∥ + ∥∇ Rj+σ ∥ max ∥∇ζ s ∥. 1≤s≤N

0≤j≤N −1

Consequently max ∥∇ζ n ∥ ≤ C max 0≤j≤N −1

1≤n≤N

{

tjα+σ ∥∇δtαj+σ ρ∥ + ∥∇ϕ j+σ ∥ + ∥∇ Rj+σ ∥

(

)}

.

(29)

We bound separately the three terms in the right-hand side of (29). First, ∥∇ Rh ϕ j+σ ∥ ≤ ∥∇ϕ j+σ ∥ and (11) give

∥∇δtαj+σ ρ∥ = ∥∇δtαj+σ ρ − ∇ Dαt ρ j+σ + ∇ Dαt ρ j+σ ∥  ( ) ( )     = ∇ Dαt u(tj+σ ) − δtαj+σ u + ∇ Rh δtαj+σ u − Dαt u(tj+σ )  + ∇ Dαt ρ j+σ  ≤ ∥∇ϕ j+σ ∥ + ∥∇ Rh ϕ j+σ ∥ + Chk ∥Dαt uj+σ ∥k+1 ≤ 2∥∇ϕ j+σ ∥ + Chk ∥Dαt u∥L∞ (H k+1 (Ω )) .

(30)

Next, by Lemmas 6.4 and 6.2 we have max

{

0≤j≤N −1

tjα+σ ∥∇ϕ j+σ ∥

}

}} { { j+σ j,s ≤ CN − min{r α,3−α} . ≤ C max ∥ψu ∥1 + max ∥ψu ∥1

(31)

s=1,...,j

0≤j≤N −1

Finally, consider the term max0≤j≤N −1 tjα+σ ∥∇ Rj+σ ∥ in (29). When j = 0, the hypothesis that u ∈ L∞ (0, T ; H01 (Ω ) ∩

{

}

H k+1 (Ω )) yields tσα ∥∇ Rσ ∥ ≤ Ct1α ≤ CN −r α . From a Taylor expansion [13, Lemma 9] and the hypothesis that ∥∂t2 u(·, t)∥3 ≤

C (1 + t α−2 ) ≤ Ct α−2 , for each j ≥ 1 one has

tjα+σ ∥∇ Rj+σ ∥ ≤ Ctjα+σ τj2+1 tjα−2 ≤ C (j + 1)r α N −r α N −2r j2(r −1) jr(α−2) N −r(α−2) ≤ C (j/N)2r α−2 N −2 , where we used (6). Hence tjα+σ ∥∇ Rj+σ ∥

{ ≤

CN −2 C (1/N)2r α−2 N −2 = CN −2r α

for j = 1, 2, . . . , N − 1 if r ≥ 1/α, for j = 1, 2, . . . , N − 1 if 1 ≤ r < 1/α.

Combine this bound with the case j = 0 above to get max 0≤j≤N −1

tjα+σ ∥∇ Rj+σ ∥ ≤ CN − min{r α,2} .

{

}

(32)

Substituting (30), (31) and (32) into (29), one gets max1≤n≤N ∥ζ n ∥ ≤ C N − min{r α,2} + hk . This inequality, (11) and (26) yield (28). □

)

(

Corollary 6.7.

Suppose that r ≥ 2/α . Then the L2-1σ FEM solution unh satisfies

( ) ∥un − unh ∥H 1 (Ω ) ≤ C hk + N −2 for n = 0, 1, . . . , N . Proof. Theorem 6.6 and r ≥ 2/α imply that ∥∇ un − ∇ unh ∥ ≤ C hk + N −2 for each n. The result now follows from a Poincaré inequality. □

(

)

(

)

The O(N −2 ) temporal convergence rate of the L2-1σ FEM shown in Corollary 6.7 is superior to the O N −(2−α ) rate that, as we saw in Section 5, is attained by the L1 FEM. The O(hk ) convergence rate of both methods is optimal in H 1 (Ω ).

10

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435 Table 1 Example 7.1: (L∞ (H 1 ), N) errors and orders of convergence for L1 FEM with r = (2 − α )/α . M = N = 32

M = N = 64

M = N = 128

M = N = 256

M = N = 512

α = 0.4

2.0043E−3

7.3119E−4 1.4548

2.5834E−4 1.5009

8.9420E−5 1.5306

3.0529E−5 1.5504

α = 0.6

2.5834E−3

1.0224E−3 1.3373

3.9823E−4 1.3603

1.5356E−4 1.3747

5.8847E−5 1.3838

α = 0.8

3.9201E−3

1.7345E−3 1.1763

7.6222E−4 1.1862

3.3358E−4 1.1921

1.565E−4 1.1955

Remark 6.1. Our proof of Theorem 6.6 relies heavily on results from [13]. By inspection, one can verify easily that this proof remains valid if our graded mesh is replaced by the more general temporal mesh used in [13]. We now describe precisely this mesh and the analogue of Theorem 6.6 that one then obtains. ˆ be two positive integers with M ˆ ≤ M but M ˆ ≥ cM for some fixed constant c ∈ (0, 1). Divide the Let M and M ˆ ˆ ˆ interval [0, T ] into two subintervals [0, T ] and [T , T ], where T ∈ (0, T ] is arbitrary. Partition [0, Tˆ ] by the graded mesh ˆ r for j = 0, . . . , M, ˆ where the mesh grading parameter r ≥ 1 is chosen by the user. Partition [Tˆ , T ] by the tj = Tˆ (j/M) mesh Tˆ = tMˆ < tMˆ +1 < · · · < tM −1 < tM = T , where the mesh is quasiuniform on [tMˆ −1 , T ] (i.e., the ratio maxj τj /minj τj

ˆ − 1. (See [13, Remark 3] for the source of these is bounded by a fixed constant) and 3/4 ≤ τj+1 /τj ≤ 62 for j ≥ M inequalities.) Both these assumptions are satisfied if one takes a uniform mesh on [Tˆ , T ]. Then under the hypotheses of Theorem 6.6, when one uses the above temporal mesh of [13], there exists a constant C such that max ∥∇ un − ∇ unh ∥ ≤ C M − min{r α,2} + hk .

(

)

1≤n≤N

7. Numerical experiments We compute numerical solutions for two initial–boundary value problems (whose spatial domains have dimensions 1 and 2) of the form (1); near t = 0 the solutions of these problems behave as described in Lemma 2.1. Example 7.1.

Consider the following problem:

∂ u = f (x, t) for (x, t) ∈ (0, 1) × (0, 1], ∂ x2 u(0, t) = u(1, t) = 0 for t ∈ (0, 1], u(x, 0) = (ex − 1)(x − 1) for x ∈ [0, 1].

Dαt u −

2

(33)

The function f (x, t) in (33) is chosen such that the exact solution of the problem is u(x, t) = (Eα (−t α ) + t 3 )(ex − 1)(x − 1), ∑∞ j where Eα (z) = j=0 z /Γ (jα + 1) is the Mittag-Leffler function [14, Section 1.2] . This solution u displays typical layer behaviour near t = 0, i.e., its derivatives agree exactly with the bounds of Lemma 2.1. In this experiment, the solution unh is computed using the L1 FEM (9) and L2-1σ FEM (24). To estimate the L∞ (H 1 ) error, we compute

 ( ) ∥u − unh ∥L∞ (H 1 ),N := max ∇ u(x, tn ) − unh  , 0≤n≤N

where 3-point Gaussian quadrature is used on each spatial mesh interval to approximate ∥∇ u(x, tn ) − unh ∥ for each n. L1 FEM results: Theorem 5.2 predicts the rate of convergence O(hk + N − min{2−α,r α} ) for ∥u − unh ∥L∞ (H 1 ),N . Table 1 shows the ∥u − unh ∥L∞ (H 1 ),N errors for α = 0.4, 0.6, 0.8 with r = (2 − α )/α and k = 2. We take M = N so that the temporal error dominates the spatial error in the bound of Corollary 5.3. The orders of convergence displayed indicate that the rate of convergence is N −(2−α ) , as predicted by the corollary. Table 2 displays the ∥u − unh ∥H 1 errors and the associated orders of convergence at t = 1 when α = 0.4 and r = (2 − α )/α . Here we take N = 1000 so that the spatial error dominates the results and we observe O(hk ) convergence, as predicted by Corollary 5.3. L2-1σ FEM results: Theorem 6.6 predicts the rate of convergence O(hk + N − min{r α,2} ) for ∥u − unh ∥L∞ (H 1 ),N . Table 3 shows

(

)

the ∥u − unh ∥L∞ (H 1 ),N errors for α = 0.4, 0.6, 0.8 with r = 2/α and k = 2. We take M = N so that the temporal error dominates the spatial error. The orders of convergence displayed indicate that the rate of convergence is N −2 , as predicted by Corollary 6.7. Table 4 shows the ∥u − unh ∥H 1 errors and the associated orders of convergence at t = 1 when α = 0.4 and r = 2/α . Here N = 1000 is taken so that the spatial error dominates the results and O(hk ) convergence is observed, as predicted by Corollary 6.7.

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435

11

Table 2 Example 7.1: H 1 (Ω ) errors and orders of convergence at t = 1 for L1 FEM with α = 0.4. Polynomial P1

P2

M

| u − uh | 1

Order

5 10 20 40

5.0480E−1 2.3367E−1 1.0127E−1 4.3782E−2

1.1112 1.2058 1.2101

5 10 20 40

1.1506E−2 2.8381E−3 7.0689E−4 1.7658E−4

2.0194 2.0053 2.0011

Table 3 Example 7.1: (L∞ (H 1 ), N) errors and orders of convergence for L2-1σ FEM with r = 2/α . M = N = 32

M = N = 64

M = N = 128

M = N = 256

M = N = 512

α = 0.4

7.2649E−3

1.9215E−3 1.9186

4.9315E−4 1.9621

1.2481E−4 1.9822

3.1382E−5 1.9917

α = 0.6

4.5176E−3

1.1549E−3 1.9677

2.9153E−4 1.9860

7.3179E−5 1.9942

1.8321E−5 1.9979

α = 0.8

3.0216E−3

7.6079E−4 1.9897

1.9064E−4 1.9966

4.7674E−5 1.9996

1.1908E−5 2.0011

Table 4 Example 7.1: H 1 (Ω ) errors and orders of convergence at t = 1 for L2-1σ FEM with α = 0.4. Polynomial P1

P2

M

| u − uh | 1

Order

5 10 20 40

5.2477E−1 2.3899E−1 1.0238E−1 4.3937E−2

1.1346 1.2229 1.2193

5 10 20 40

1.1253E−2 2.8227E−3 7.0654E−4 1.7708E−4

1.9951 1.9982 1.9963

Table 5 Example 7.2: (L∞ (H 1 ), N) errors and orders of convergence for L1 FEM with r = (2 − α )/α .

Example 7.2.

N = 10

N = 20

N = 30

N = 40

α = 0.4

1.7446E−1

5.7675E−2 1.5969

3.0211E−2 1.5947

1.9068E−2 1.5997

α = 0.6

2.6843E−1

1.0416E−1 1.3657

5.9647E−2 1.3749

3.9879E−2 1.39948

α = 0.8

4.3637E−1

1.8862E−1 1.2100

1.1631E−1 1.1923

8.3082E−2 1.1696

We consider the following time-fractional diffusion problem on a two-dimensional spatial domain:

α

Dt u(x, y, t) − ∆u(x, y, t) = f (x, y, t)

for (x, t) ∈ Ω × (0, 1],

u(x, y, 0) = sin(π x) sin(π y) u|∂ Ω = 0

(34)

for (x, y) ∈ Ω ,

for t ∈ (0, 1],

where Ω = (0, 1) × (0, 1). The function f (x, y, t) in (34) is chosen such that the exact solution of the problem is u(x, t) = (1 + t α ) sin(π x) sin(π y), In this experiment, we partition Ω by a uniform square mesh with M +1 nodes in each spatial direction. Table 5 displays the errors ∥u − unh ∥L∞ (H 1 ),N of the L1 FEM and their associated orders of convergence, where M = ⌈N 2−α ⌉ is chosen so that the temporal error dominates the spatial error. The rates of convergence obtained agree with the theoretical rates of convergence stated in Corollary 5.3: they show that r = (2 − α )/α yields the optimal rate of convergence O(N −(2−α ) ). Table 6 shows the corresponding results for the L2-1σ FEM with r = 2/α , where M = N 2 is chosen so that the temporal error dominates the spatial error. The numerical results agree with Corollary 6.7: they show that r = 2/α yields the optimal rate of convergence O(N −2 ).

12

C. Huang and M. Stynes / Journal of Computational and Applied Mathematics 367 (2020) 112435 Table 6 Example 7.2: (L∞ (H 1 ), N) errors and orders of convergence for L2-1σ FEM with r = 2/α . N = 5

N = 10

N = 15

N = 20

α = 0.4

2.8265E−1

7.0504E−2 2.0032

3.1329E−2 2.0004

1.7623E−2 1.9999

α = 0.6

2.8102E−1

7.0089E−2 2.0034

3.1138E−2 2.0009

1.7513E−2 2.0004

α = 0.8

2.7982E−1

6.9864E−2 2.0019

3.1044E−2 2.0005

1.7461E−2 2.0002

Acknowledgement The research of Martin Stynes is supported in part by the National Natural Science Foundation of China under grant NSAF U1930402. References [1] Bangti Jin, Raytcho Lazarov, Zhi Zhou, Numerical methods for time-fractional evolution equations with nonsmooth data: a concise overview, Comput. Methods Appl. Mech. Engrg. 346 (2019) 332–358. [2] Hong-lin Liao, Dongfang Li, Jiwei Zhang, Sharp error estimate of the nonuniform L1 formula for linear reaction-subdiffusion equations, SIAM J. Numer. Anal. 56 (2) (2018) 1112–1133. [3] Martin Stynes, Eugene O’Riordan, José Luis Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal. 55 (2) (2017) 1057–1079. [4] Chaobao Huang, Martin Stynes, A direct discontinuous Galerkin method for a time-fractional diffusion equation with a Robin boundary condition, Appl. Numer. Math. 135 (2019) 15–29. [5] Natalia Kopteva, Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions, Math. Comp. 88 (319) (2019) 2135–2155. [6] Vidar Thomée, Galerkin Finite Element Methods for Parabolic Problems. 2nd Revised and Expanded ed., second revised and expanded ed., Springer, Berlin, 2006, p. xii + 370. [7] Jinchang Ren, Hong-lin Liao, Jiwei Zhang, Zhimin Zhang, Sharp H 1 -norm error estimates of two time-stepping schemes for reaction-subdiffusion problems. Preprint, arXiv:1811.08059, 2018. [8] Daniel Henry, Geometric Theory of Semilinear Parabolic Equations, in: Lecture Notes in Mathematics, vol. 840, Springer-Verlag, Berlin-New York, 1981, p. iv+348. [9] Kenichi Sakamoto, Masahiro Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, J. Math. Anal. Appl. 382 (1) (2011) 426–447. [10] James H. Bramble, Joseph E. Pasciak, Olaf Steinbach, On the stability of the L2 projection in H 1 (Ω ), Math. Comp. 71 (237) (2002) 147–156. [11] Chaobao Huang, Martin Stynes, Na An, Optimal L∞ (L2 ) error analysis of a direct discontinuous Galerkin method for a time-fractional reaction-diffusion problem, BIT 58 (3) (2018) 661–690. [12] Anatoly A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys. 280 (2015) 424–438. [13] Hu Chen, Martin Stynes, Error analysis of a second-order method on fitted meshes for a time-fractional diffusion problem, J. Sci. Comput. 79 (1) (2019) 624–647. [14] Igor Podlubny, Fractional Differential Equations, in: Mathematics in Science and Engineering, vol. 198, Academic Press, Inc., San Diego, CA, 1999, p. xxiv+340, An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications.