A direct discontinuous Galerkin method for a time-fractional diffusion equation with a Robin boundary condition

Applied Numerical Mathematics 135 (2019) 15–29

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A direct discontinuous Galerkin method for a time-fractional diffusion equation with a Robin boundary condition Chaobao Huang, Martin Stynes ∗ Applied and Computational Mathematics Division, Beijing Computational Science Research Center, Beijing 100193, PR China

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 3 May 2018 Received in revised form 10 August 2018 Accepted 10 August 2018 Available online 15 August 2018

A time-fractional reaction–diffusion initial-boundary value problem with Robin boundary condition is considered on the domain  × [0, T ], where  = (0, l) ⊂ R. The coefficient of the zero-order reaction term is not required to be non-negative, which complicates the analysis. In general the unknown solution will have a weak singularity at the initial time t = 0. Existence and uniqueness of the solution and pointwise bounds on some of its derivatives are derived. A fully discrete numerical method for computing an approximate solution is investigated; it uses the well-known L1 discretisation on a graded mesh in time and a direct discontinuous Galerkin (DDG) finite element method on a uniform mesh in space. Discrete stability of the computed solution is proved. Its error is bounded in the L 2 () and H 1 () norms at each discrete time level tn by means of a non-trivial projection of the unknown solution into the finite element space. The L 2 () bound is optimal for all tn ; the H 1 () bound is optimal for tn not close to t = 0. An optimal grading of the temporal mesh can be deduced from these bounds. Numerical results show that our analysis is sharp. © 2018 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: DDG method Fractional reaction–diffusion equation Robin boundary condition Graded mesh

1. Introduction Fractional derivatives are now used frequently in the modelling of physical processes. In this paper we shall examine the most popular class of fractional-derivative problems, which has received much attention in the research literature. Thus, consider the time-fractional initial-boundary value problem

∂ 2u + c (x)u = f (x, t ) ∀(x, t ) ∈  × (0, T ], ∂ x2 σ0 u (0, t ) − u x (0, t ) = 0 and σ1 u (l, t ) + u x (l, t ) = 0 for t ∈ (0, T ], D tα u −

u (x, 0) = φ(x) for x ∈ .

(1a) (1b) (1c)

¯ ), f ∈ C ( Q¯ ) where ¯ = [0, l] and σ0 and σ1 are positive constants. We assume that c ∈ C ( Here 0 < α < 1 with  = (0, l),  ¯ ). (Stronger regularity hypotheses will be placed on these functions later.) In (1a), D tα u Q¯ := [0, l] × [0, T ], and φ ∈ C ( denotes the Caputo fractional derivative, which is defined by

*

Corresponding author. E-mail address: [email protected] (M. Stynes).

https://doi.org/10.1016/j.apnum.2018.08.006 0168-9274/© 2018 IMACS. Published by Elsevier B.V. All rights reserved.

16

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

D tα u (x, t ) =

1

(1 − α )

t

(t − s)−α

∂ u (x, s) ds. ∂s

0

Problems similar to (1) have been considered in many numerical analysis papers; see [10,4,13] and their references. More general boundary conditions σ0 u (0, t ) − u x (0, t ) = g 1 (t ) and σ1 u (0, t ) + u x (0, t ) = g 2 (t ), where g 1 (t ) and g 2 (t ) are smooth, are easily reduced to the homogeneous case of (1b) by a standard linear change of variable. ¯ ) implies that − M 1 ≤ c (·) ≤ M 2 for some positive constants M 1 and M 2 . The hypothesis that c ∈ C ( Remark 1.1. For classical parabolic initial-boundary value problems (i.e., α = 1 in (1a)), it is well known that one can make a simple change of variable that yields an equivalent problem where c ≥ 0. This facilitates the use of both maximum principle and energy-norm analyses. But for (1a) no such change of variable is available, so our analysis must deal with the awkward possibility that c may take negative values. Remark 1.2. If a convection term q(x)u x were present in (1a), it could be removed by a simple change of variable; see [10]. Our aim in this paper is to analyse a numerical method that uses the standard L1 discretisation in time — on a graded mesh — to approximate the fractional time derivative, and a direct discontinuous Galerkin (DDG) finite element method on a uniform mesh in space to approximate the spatial derivative. The DDG permits high-order approximation of the solution in the spatial direction, where the solutions are typically smooth; it has a practical advantage over other discontinuous Galerkin methods since no new (artificial) variables are introduced in the weak formulation of the problem, but its analysis poses difficulties that are outlined in [4], where an initial-boundary value problem with Dirichlet boundary conditions was considered — in fact [4] gives the first analysis of the DDG method for a problem with non-periodic boundary conditions. In the current paper we continue our development of the analysis of the DDG for fractional-in-time problems by examining the Robin boundary conditions (1b). The structure of our paper is as follows. In Section 2 bounds are derived on derivatives of the solution that are needed for the subsequent numerical analysis. In Section 3 we present the details of our numerical method. The analysis of the method is given in Section 4, in the realistic setting where u has a layer at t = 0: first, we derive a discrete stability result in L 2 (), then we invoke the carefully-constructed projection P of [4, equation (33)] to handle the Robin boundary conditions, finally obtaining our convergence results in L 2 () and H 1 (). Numerical results are presented in Section 5 to demonstrate the sharpness of our theoretical convergence results. Notation We use C to denote a generic constant that is independent of the mesh; it can take different values in different places. We write  ·  p ,ω for the norm in L p (ω), where 1 ≤ p ≤ ∞ and ω is any measurable subset of [0, l]. If ω = , we drop it from this notation. When p = 2 and ω = , we write  ·  for the norm in L 2 (). Let N = {1, 2, 3, . . .} and N0 = {0, 1, 2, . . .}. The notation H k (), for k ∈ N, is used for the standard Sobolev spaces. 2. A priori bounds on the solution of (1) Imitating [12,13], we use separation of variables to construct an infinite series that is a classical solution u (x, t ) of (1). ˜ i , ψi ) : i = 1, 2, . . . } be the eigenvalues and eigenfunctions for the two-point boundary value problem Let {(λ

i ψi on (0, l), σ0 ψi (0) − ψi (0) = σ1 ψi (l) + ψi (l) = 0, Lψi := −ψi

+ c ψi = λ

(2)

where the eigenfunctions are normalised by requiring that ψi 2 = 1 for all i. Now (2) can be rewritten so the eigenvalues and eigenfunctions (λi , ψi ) for the Sturm–Liouville two-point boundary value problem

−ψi

+ (c + M 1 )ψi = λi ψi on (0, l), σ0 ψi (0) − ψi (0) = σ1 ψi (l) + ψi (l) = 0,

(3)

i + M 1 , with ψi the same in both (2) and (3). Since c + M 1 ≥ 0, it is well known that ψi ∞ ≤ C [1, p. 335] and satisfy λi = λ λi ≈ (i π )2 /l2 [1, p. 415] for all i. Thus we have λi ≈ (i π )2 /l2 − M 1 , and it follows that only a finite number of eigenvalues i can be negative. As in [13, (2.2)], a standard separation of variables technique leads formally to λ u (x, t ) =

∞  

 i t α ) + J i (t ) ψi (x), (φ, ψi ) E α ,1 (−λ

i =1

where

t J i (t ) := 0

i sα ) f i (t − s) ds with f i (t ) := ( f (·, t ), ψi (·)), sα −1 E α ,α (−λ

(4)

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

17

and the generalised Mittag-Leffler function [11, Section 1.2] is defined by

E ρ ,β ( z) =

∞  k =0

zk

(ρ k + β)

.

Note that when β = 1, one usually writes E ρ ,1 as E ρ . For 0 < ρ < 2, β ∈ R, and r ≥ 0, by [11, Theorems 1.5 and 1.6] there exists a constant C such that

C

E ρ ,β (r ) ≤ C (1 + r )(1−β)/ρ exp(r 1/ρ ) +

1+r

and E ρ ,β (−r ) ≤

C 1+r

.

Define the space C W 1 [0, T ] := C [0, T ] ∩ W t1 (0, T ], where W t1 (0, T ] is the space of functions φ ∈ C 1 (0, T ] such that

φ ∈ L 1 (0, T ).

Next, we show that the problem (1) satisfies a maximum principle in the special case where c ≥ 0.

Lemma 2.1. Let u (·, t ) ∈ C 2 [0, l] for each t, and u (x, ·) ∈ C W 1 [0, T ] for each x. Assume that f ≤ 0 on Q and c ≥ 0 on . Assume that u satisfies the boundary conditions σ0 u (0, t ) − u x (0, t ) ≤ 0 and σ1 u (l, t ) + u x (l, t ) ≤ 0, t ∈ (0, T ]. Then





max u (x, t ) ≤ max max u (x, 0), 0 .

(5)

¯ x∈

(x,t )∈ Q¯

Proof. Suppose that (5) is false, i.e., the function u attains its global positive maximum at (x∗ , t ∗ ) ∈ [0, l] × (0, T ]. If x∗ = 0, then u x (0, t ∗ ) ≤ 0, as otherwise one does not have a maximum at (x∗ , t ∗ ). Consequently σ0 u (0, t ∗ ) − u x (0, t ∗ ) > 0 which contradicts the hypothesis σ0 u (0, t ∗ ) − u x (0, t ∗ ) ≤ 0. Similarly x∗ = l is impossible. Hence (x∗ , t ∗ ) ∈ (0, l) × (0, T ]. The argument is completed by imitating the proof of [8, Theorem 2]. 2 Continue for a moment to assume that c ≥ 0 on . If

f ≥ 0 on Q , with σ0 u (0, t ) − u x (0, t ) ≥ 0 and

¯ , (6) σ1 u (l, t ) + u x (l, t ) ≥ 0 for t ∈ (0, T ], and u (x, 0) ≥ 0 for x ∈ 

then u ≥ 0 on Q . This follows on changing u to −u in Lemma 2.1. When c ≥ 0 does not hold, one can nevertheless combine this result for c ≥ 0 with the iterative argument of [9, Theorem 1], modified to accommodate our Robin boundary conditions, to show that (1) still satisfies the same result: when the inequalities of (6) are satisfied, then u ≥ 0 on Q . This implies the uniqueness of any solution of (1). For each non-negative γ ∈ R, the fractional power Lγ of the operator L (see, e.g., [3]) has domain

 D (Lγ ) := g ∈ L 2 (0, l) :

∞ 

 2 2γ λi ( g , ψi ) < ∞ .

i =1

Lemma 2.2. Assume that φ ∈ D (L7/2 ), f (·, t ) ∈ D (L5/2 ), f t (·, t ) and f tt (·, t ) are in D (L1/2 ) for each t ∈ (0, T ] with

 f (·, t )L5/2 +  f t (·, t )L1/2 + t δ  f tt (·, t )L1/2 ≤ C 1 , for all t ∈ (0, T ] and some constant δ < 1, where C 1 is a constant independent of t. Assume also that c ∈ C 2 [0, 1]. Then (1) has a unique solution u that satisfies (1a), (1b) and (1c) pointwise, and there exists a constant C such that

p ∂ u ≤ C for p = 0, 1, 2, 3, 4, ( x , t ) ∂ xp q ∂ u α −q ) for q = 0, 1, 2, ∂ t q (x, t ) ≤ C (1 + t α p ∂ ∂ u ≤ C for p = 0, 1, 2, 3, 4, ( x , t ) ∂tα ∂ xp for all (x, t ) ∈ [0, l] × (0, T ]. Proof. Imitate the proof of [4, Lemma 2.2].

2

18

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

3. The numerical method First we describe a semi-discrete DDG method where we discretise only in space. Let M be a positive integer. Set h = l/ M and xm := mh for m = 0, 1, . . . , M. Define the mesh intervals I m = [xm−1 , xm ] and xm− 1 = (xm−1 + xm )/2 for m = 1, 2, . . . , M. The corresponding weak formulation of (1) is: Find u (·, t ) ∈ H 1 () for each 2

t ∈ (0, T ], such that for m = 1, 2, . . . , M one has

⎧ ⎪ ⎪ ⎨  ⎪ ⎪ ⎩



xm

D tα uv dx − u x v x

Im

Im

m −1

 u (x, 0) v dx =

+

Im



u x v x dx +

 Im

c (x)uv dx =

Im

f v dx ∀ v ∈ H 1 (), (7)

∀ v ∈ H 1 (),

φ(x) v dx Im

where u x (x0 ) = σ0 u (x0 ) and u x (x M ) = −σ1 u (x M ) are used for m = 1 and m = M, respectively. Fix k ≥ 1. Define the finite element space on the spatial mesh by







V h = v ∈ L 2 () : v I ∈ P k ( I m ), m = 1, 2, . . . , M , m

d

where P ( I m ) denotes the space of polynomials of degree at most d (≥ 0) defined on I m . We adopt the following notation, which is frequently used in our later analyses:



+ − u ± m := u (xm ± 0, t ), [u ]m := um − um , {u }m :=

+ + u− um m

2

for each m.

From (7), we construct the following DDG semi-discrete scheme: Find u h (·, t ) ∈ V h for each t ∈ (0, T ], such that for m = 1, . . . , M one has

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪  ⎪ ⎪ ⎪ ⎪ ⎩



xm D tα u h v h dx − ( u h )x v h x

m −1

Im 

=

Im

Im



1

+ Im

−

(uh )x ( v h )x dx + ([uh ]( v h )x )m + 2

1 2

+ ([uh ]( v h )x )m −1



+

c (x)u h v h dx Im

f v h dx ∀ v h ∈ V h ,



u (x, 0) v h dx =

Im

φ(x) v h dx ∀ v h ∈ V h , (8)

where we define [u h ]m to be 0 at m = 0 and m = M. Furthermore, for each v h ∈ V h the numerical flux ( v h )x at each point xm is defined by

⎧ v h )x x := β0 h−1 [ v h ]m + {( v h )x }m + β1 h[∂x2 v h ]m for m = 1, 2, . . . , M − 1, ( ⎪ ⎪ m ⎪ ⎨ ( v h )x x := σ0 v h (x0 ), 0 ⎪ ⎪ ⎪ ⎩  ( v h )x x := −σ1 v h (xM ),

(9)

M

where β0 and β1 are user-chosen positive constants; we will discuss their possible values in Section 4. Now we discretise (8) in the temporal variable. 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.) Set τn = tn − tn−1 for n = 0, 1, . . . , N. The Caputo fractional derivative D tα u (xm , tn ) for n ≥ 1 can be written as

D tα u (xm , tn ) =

n −1 

1

(1 − α )

ti+1 ∂ u (xm , s) (tn − s)−α ds. ∂s

i =0 t i

This is approximated by the standard L1 formula

D αN unm :=

=

n −1 

1

(1 − α ) 1

(2 − α )

i =0

n −1  i =0

i +1 i um − um

τ i +1 i +1 um

i  − um

τ i +1

t i+1 (tn − s)−α ds ti

(tn − t i )1−α − (tn − t i +1 )1−α



C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

dn,1

=

(2 − α )

unm −

dn,n

(2 − α )

1

0 um +

(2 − α )

n −1 

19

unm−i (dn,i +1 − dn,i ),

(10)

i =1

where dn,i := [(tn − tn−i )1−α − (tn − tn−i +1 )1−α ]/τn−i +1 for i ≥ 1; in particular, dn,1 = τn−α . Invoking the mean value theorem, one sees easily that

dn,i +1 ≤ dn,i for i ≥ 1.

ϕ n (x) is defined by

The temporal truncation error

ϕ

n

(x) := D tα u (x, tn ) − D αN u (x, tn ) =

n −1  i =0

  t i+1 ∂u u (x, t i +1 ) − u (x, tk ) (tn − s)−α (x, s) − ds. (1 − α ) ∂s τ i +1 1

ti

Lemma 3.1. [13, Lemma 5.1] There exists a constant C such that for all (xm , tn ) ∈ Q one has

n ϕ (xm ) ≤ Cn− min{2−α ,r α } . Now we apply the L1 scheme (10) to (8) for the temporal discretisation:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪  ⎪ ⎪ ⎪ ⎪ ⎩



xm

D αN unh v h dx − ( unh )x v h x

m −1

Im 

+

Im

1

1

2

2

+ − (unh )x ( v h )x dx + ([unh ]( v h )x )m + ([unh ]( v h )x )m −1 +

 Im

c (x)unh v h dx

f n v h dx ∀ v h ∈ V h for m = 1, . . . , M and n = 1, 2, . . . , N ,

= Im

Im



u h0 v h dx

=

Im

φ(x) v h dx ∀ v h ∈ V h , (11)

where ∈ V h is the numerical solution at time t = tn and f (·) := f (·, tn ) for n = 1, 2, . . . , N. Summing (11) for m = 1, . . . , M yields the fully discrete L1-DDG method: unh

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪  ⎪ ⎪ ⎪ ⎩

n



D αN unh v h dx + n 1 (u h v h )|x M

+σ

M  

unh )x ( v h )x dx +

m= 1

Im

+

c (x)unh v h dx





u h (x, 0) v h dx = 





M −1  

  n  [unh ]{( v h )x } + ( u h )x [ v h ] + σ0 (unh v h )|x0

m =1

=

m

f v h dx ∀ v h ∈ V h ,

(12)



φ(x) v h dx ∀ v h ∈ V h . 

Define the operator A and the L 2 () inner product (·, ·) by

A (ζ, η) =

M   m =1 I



m

ζx ηx dx +

M −1  m =1



 [ζ ]{ηx } + ζx [η] m + σ0 (ζ η)|x0 + σ1 (ζ η)|xM ∀ ζ, η ∈ V h ,

v w dx, ∀ v , w ∈ L 2 ().

(v , w ) = 

Then the fully discrete L1-DDG scheme (12) can be written compactly as

⎧ α n ⎨ ( D N uh , v h ) + A (unh , v h ) + (cunh , v h ) = ( f n , v h ) ∀ v h ∈ V h , ⎩

(uh0 , v h ) = (φ, v h ) ∀ v h ∈ V h .

This formulation is convenient when analysing our numerical method. 4. Stability and convergence of the L1-DDG method Imitating [7, Definition 2.1], we adopt the following criterion, which will ensure the stability of our method.

(13)

20

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

Definition 4.1. A numerical flux ( v h )x of the form (9) is said to be admissible if there exist constants (0, 1] such that

γ1

M   m =1 I

( v h )2x dx +

M −1  2  [ v h ]m v h )x [ v h ] ≥ γ2 for all v h ∈ V h . [ v h ]{( v h )x } + (

M −1  

m

m =1

m

γ1 ∈ (0, 1) and γ2 ∈

h

m =1

Recall that V h contains piecewise polynomials of degree k. Sufficient conditions for admissibility of the numerical flux

v h )x are given in [7, Lemma 2.1]. We do not discuss these conditions here since this is already done in [4]; instead, we ( merely point out that the choices β0 = 4 and β1 = 1/12 that are used with k = 1, 2 in the numerical results of Section 5 below yield an admissible numerical flux. In the analysis that follows, we assume that our numerical flux is admissible. Define the L 2 (), H 1 () and discrete energy norms by

 v 2 =

M   m =1 I

v 2 dx,

| v |21 =

m =1 I

m

 v 2E = | v |21 +

M  

M −1 

β0 h

m =1

v 2x dx ∀ v ∈ V h ,

m

2 [ v ]m ∀ v ∈ V h.

From the admissibility criterion Definition 4.1, we have

A(v , v ) =

M   m =1 I

v 2x dx +

m =1

m

M  

≥ (1 − γ1 )

m =1 I

m =1 I



 [ v ]{ v x } + vx [ v ] m + σ0 v 2 |x0 + σ1 v 2 |xM

v 2x dx + γ2

M −1 

2 [ v ]m

m =1

m

M  

= (1 − γ1 ) Hence there exists

M −1 

v 2x dx +

m

h

M −1 γ 2  β0

β0

m =1

h

+ σ0 v 2 |x0 + σ1 v 2 |xM

2 [ v ]m + σ0 v 2 |x0 + σ1 v 2 |xM .

γ ∈ (0, 1) such that

A ( v , v ) ≥ γ  v 2E , ∀ v ∈ V h .

(14)

To analyse the numerical stability of the L1 scheme (10), we are unable to invoke the results of [13] (as was done in [4]) because there one had c ≥ 0 but this property does not hold in the present paper. Thus instead we shall use some results from [5] that do not require c ≥ 0. (n) Define the coefficients an−i := dn,n−i +1 /(2 − α ) for i = 1, . . . , n and

⎧ n ⎪ ⎨(2 − α )τ α  (a( j ) − a( j ) ) Q (n) (n) i j − i −1 j −i n− j Q n−i := j = i +1 ⎪ ⎩(2 − α )τ α n

Then the L1 scheme (10) can be written as D αN un =

when i = 1, 2, . . . , n − 1,

n

(n)

(15)

when i = n. (n) k k=1 an−k ∇t u ,

where ∇t uk := uk − uk−1 . The next result bounds a

weighted sum of the Q n−i that will be essential in our analysis. Set ωβ (t ) = t β /(1 + β) for t

> 0 and β > −1.

Lemma 4.1. [5, Lemma 2.1] Let the nonnegative integer m satisfy 0 ≤ m ≤ 1/α  (the largest integer that is not greater than 1/α ). Then n  i =1

(n)

Q n−i ωmα −α (t i ) ≤ ωmα (tn ) for n = 1, . . . , N .

Next, we state a discrete fractional Gronwall inequality.

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

21

n−1

Lemma 4.2. [5, Lemma 2.2] Let the κs be nonnegative constants with 0 < s=0 κs ≤ κ for n ≥ 1, where κ is some positive constant independent of n. Suppose that the nonnegative sequences {ξn , ηn | n ≥ 1} are bounded and the grid function { v n | n ≥ 0} satisfies

D αN v n ≤

n 

κn−s v s + ξn + ηn for n ≥ 1.

s =1



If the nonuniform grid satisfies the maximum time-step criterion τ N ≤ 2(2 − α )κ

 v k ≤ 2E α (2λtkα ) v 0 + max

j 

1≤ j ≤k

s =1

−1/α

, then



( j)

Q j −s ξ s + ωα (tk ) max

1≤s≤k

η s for k = 1, . . . , N ,

( j)

where Q j −s is defined in (15) and E α (·) is the Mittag-Leffler function of Section 2. Our stability result will be presented in a general framework. Suppose that for each n, the function



μn ∈ V h satisfies

( D αN μn , v h ) + A (μn , v h ) + (c μn , v h ) = ( g n , v h ) ∀ v h ∈ V h ,

(16)

(μ0 , v h ) = (φ, v h ) ∀ v h ∈ V h , with

σ0 μ(0, t n ) − μx (0, t n ) = σ1 μ(l, t n ) + μx (l, t n ) = 0. A stability bound for (16) is presented in the next lemma.

Lemma 4.3. If the graded mesh satisfies the maximum time-step condition discrete problem (16) satisfies



τN ≤ 2(2 − α ) M 1

−1/α

, then the solution

  j  ( j) μn  ≤ 2E α (2M 1 tnα ) μ0  + max Q j −s  g s  1≤ j ≤n

μn of the

(17a)

s =1

and

γ |μn |21 + γ

M  β0 m =1

h

2 [μn ]m ≤

dn,1

max μi 2 +

(2 − α )

(2 − α ) 0≤i ≤n−1

dn,1

 g n 2

(17b)

for n = 1, . . . , N, where C is a constant independent of n and N. Proof. Fix n ∈ {1, 2, . . . , N }. Letting v h = μn in equation (16) gives

  ( D αN μn , μn ) + A (μn , μn ) + c μn , μn = ( g n , μn ).

(18)

Hence, invoking (14), we obtain

  ( D αN μn , μn ) + c μn , μn ≤ ( g n , μn ).

(19)

Recall the definition (10) of D αN μn ; here Cauchy–Schwarz inequalities imply that

( D αN μn , μn ) =

≥

dn,1

(2 − α ) dn,1

(2 − α )

μn 2 −

μn 2 −

dn,n

(2 − α ) dn,n

(2 − α )

(μ0 , μn ) +

n −1 

1

(2 − α )

μ0  μn  +

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

i =1

1

(2 − α )

n −1 

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

i =1

  = D αN μn  μn . Insert this bound in (19), then by − M 1 ≤ c (x) ≤ M 2 and a Cauchy–Schwarz inequality we get

D αN μn  μn  ≤ M 1 μn 2 +  g n  μn , which is equivalent to

D αN μn  ≤ M 1 μn  +  g n . Now apply Lemma 4.2 with v n = μn , ξn =  g n ,

ηn = 0, κ0 = M 1 and κs = 0 for s ≥ 1. This gives (17a).

22

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

Next, we prove (17b). As the mesh width τn is an increasing function of n, the hypothesis on τ N is also satisfied by τn for each n = 1, . . . , N, and this is equivalent to M 1 ≤ (τn )−α /[2(2 − α )] = dn,1 /[2(2 − α )]. Using this inequality with − M 1 ≤ c (x) and (14), and once again appealing to (18), by (14) we get

( D αN μn , μn ) + γ |μn |21 + γ

M  β0

h

m =1

2 [μn ]m ≤ ( g n , μn ) +

dn,1 2(2 − α )

μn 2 .

Consequently (10) and the inequality ab ≤ a2 + (b2 /4) for a, b ≥ 0 give

dn,1

(2 − α ) ≤

≤

≤

(2 − α )





dn,n μ0 2 +

(2 − α ) dn,1 dn,1

 g n 2 +

μ  +

(2 − α ) dn,1

 g n 2 +

1 4

2(2 − α )

μn 2



−1  n   1 μn 2 + (dn,i − dn,i +1 ) μn−i 2 + μn 2

4

i =1

dn,1 4(2 − α )

μn 2 +

dn,1 2(2 − α )

μn 2



1

n 2

4(2 − α )

dn,1

n

i =1

(2 − α )

+

2 [μn ]m

 n −1  n −i n dn,n (μ , μ ) + (dn,i − dn,i +1 )(μ , μ ) + ( g n , μn ) + 0

1

h

m =1



1

+

M  β0

μn 2 + γ |μn |21 + γ

 n −1  n −i 2 dn,n μ  + (dn,i − dn,i +1 )μ  0 2

(2 − α )

i =1

3dn,1 4(2 − α )

μn 2 .

Thus,

γ |μn |21 + γ

M  β0 m =1

h

2 [μn ]m ≤

≤



1

(2 − α )

dn,n μ0 2 +

n −1 

 (2 − α ) n 2 (dn,i − dn,i +1 )μn−i 2 + g  dn,1

i =1

dn,1

max μi 2 +

(2 − α ) 0≤i ≤n−1

(2 − α ) dn,1

 g n 2 ,

2

as desired.

Theorem 4.4 (Stability of the L1-DDG method). Let μn be the solution of (16). If the maximum time-step of the graded mesh satisfies



τN ≤ 2(2 − α ) M 1

−1/α

, then for n = 1, 2, . . . , N, one has

 n

μ



 ≤ 2E α (2M 1 tnα ) μ0  + ωα (tn ) max  g j  ,

(20a)

1≤ j ≤n

|μn |21 +

M  β0 m =1

h

  2 [μn ]m ≤ C τn−α μ0 2 + max  g j 2 .

Proof. By Lemma 4.1 with m = 1, we have prove (20b), by (17b) and (20a) we have n 2 |1

γ |μ

+γ

M  β0 m =1

(20b)

1≤ j ≤n

h

j

s=1



n 2 ]m

[μ

≤

( j)

Q j −s ≤ ωα (t j ), and the inequality (17a) reduces immediately to (20a). To

2

2 2E α (2M 1 tnα−1 ) dn,1

(2 − α )



 0 2

2

μ  + ωα (tn−1 )

j 2

max  g 

1≤ j ≤n−1

  ≤ C τn−α μ0 2 + max  g j 2 , 1≤ j ≤n

for some constant C , using E α (2M 1 tnα−1 ) ≤ E α (2M 1 T α ) and w α (tn−1 ) ≤ T α /(1 + α ).

2

+

(2 − α ) dn,1

 g n 2

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

23

We now define a special projection P into V h which differs from the projection used in [6, (29b)–(29(c)]. Given a smooth function w, define P w ∈ V h by



( P w − w ) v dx = 0 ∀ v ∈ P k−2 ( I m ) for m = 1, . . . , M ,

(21a)

Im

−1 [ P w ] + {∂x ( P w )} + β1 h[∂x2 ( P w )] x = ∂x w (xm ) for m = 1, . . . , M − 1, ∂ x ( P w ) := β0 h m

( P w )(x0 ) = w (x0 ), ( P w )(xM ) = w (xM ), { P w } x = w (xm ) for m = 1, . . . , M − 1. m

(21b) (21c)

(In the case k = 1, the condition (21a) is ignored.) For a piecewise smooth function w with w | I m ∈ H k+1 ( I m ), the above x and { w }, respectively. definition should be modified by replacing the right-hand sides of (21b) and (21c) by w Provided that



β0 > k

2

2

1 − β1 (k − 1) +

β12 3



2

2

(k − 1) ,

it is proved in [4, Lemma 5.5] that the projection P defined in (21) is well defined. Lemma 4.5 (Projection error estimates). If w ∈ H k+1 (), then we have: M  m =1 M  m =1 M  m =1

 P w − w 20, I m ≤ Ch2k+2 | w |k2+1, , 2 | P w − w |m ≤ Ch2k+1 | w |k2+1, ,

| P w − w |2∞, I m ≤ Ch2k | w |k2+1, ,

where C is a constant independent of I m and w. Proof. Imitate the proof of [6, Theorem 7.1].

2

Next, we state a technical result which will be used in our error analysis. Lemma 4.6. For any integer j ∈ {1, 2, . . . , N }, one has j  s =1







Q j −s hk+1 + s− min{2−α ,r α } ≤ C hk+1 + N − min{2−α ,r α } ( j)



for some constant C that is independent of j and N, Proof. Taking m = 1 in Lemma 4.1 gives j  s =1

as

j

s=1

( j)

Q j −s ≤ ωm (t j ), so

Q j −s hk+1 ≤ ωα (t j )hk+1 ≤ Chk+1 , ( j)

(22)

ωα (t j ) = t αj /(1 + α ) ≤ T α /(1 + α ). Taking m = 0 in Lemma 4.1 gives j  s =1

Q j −s s− min{2−α ,r α } ≤ ( j)

j  s =1

( j)

Q j −s ω−α (t s ) max

1≤i ≤n

i − min{2−α ,r α }

ω−α (t i )

≤ max (1 − α )t iα i − min{2−α ,r α } 1≤i ≤n

= (1 − α ) max

1≤i ≤n

 α rα T i

i − min{2−α ,r α } N rα



j

s=1

( j)

Q j −s ω−α (t s ) ≤ 1. Hence

24

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

= (1 − α ) T α N − min{2−α ,r α } max



1≤i ≤n

i

r α −min{2−α ,r α } 

N

≤ C N − min{2−α ,r α } . Combining (22) and (23), the Lemma is proved.

(23)

2

We can now give our main result. Theorem 4.7 (Error estimate for L1-DDG method). Assume the hypotheses of Lemma 2.2. Assume also that φ ∈ H k+1 (). Assume that



the graded mesh satisfies the maximum time-step condition τ N ≤ 2(2 − α ) M 1 respectively at time t = tn . Then for n = 1, . . . , N, one has

−1/α

. Let un and unh be the solutions of (7) and (13)

  un − unh  ≤ C E α (2M 1 tnα ) hk+1 + N − min{2−α ,r α } ,

(24a)

M  n   β0 n 2 u − un 2 + [u − unh ]m ≤ C τn−α h2k + N −2 min{2−α ,r α } , h 1

(24b)

m =1

h

where the constant C is independent of h and n. Proof. Recall the definition of P u in (21). Write

un − unh = ( P un − unh ) − ( P un − un ) = en − εn ,

(25)

where u and represent the solution of (7) and (13) respectively at t = tn , and e := P u − ε := P u − u . Here can be estimated immediately by Lemma 4.5, but the analysis of en is more difficult and we examine it now. Combining (8) and (13), we obtain the error equation n

unh

n

n

unh ,

n

( D αN en , v h ) + A (en , v h ) + (cen , v h ) = ( D αN εn , v h ) + A (εn , v h ) + (c εn , v h ) − (ϕ n , v h ) ∀ v h ∈ V h .

n

n

εn

(26)

Take v h = en in (26). Then

( D αN en , en ) + A (en , en ) + (cen , en ) = ( D αN εn , en ) + A (εn , en ) + (c εn , en ) − (ϕ n , en ). n

(27)

n

To deal with the right-hand side of (27), we first write out all the terms in A (ε , e ):

A (εn , en ) =

M   m =1 I

=−

εxn enx dx +

 [εn ]{enx } + εxn [en ] + σ0 (εn en )|x0 + σ1 (εn en )|xM m

m =1

m

M   m =1 I

M −1  

enxx εn dx +

M −1  

 − [enx εn ] + [εn ]{enx } + εxn [en ] + (εn enx )|x0 + (εn enx )|xM m

m =1

m

+ σ0 (εn en )|x0 + σ1 (εn en )|xM =−

M   m =1 I

enxx n dx +

ε

m

M −1  

 n n n n  εx [e ] − {ε }[e x ] + (εn enx )|x0 + (εn enx )|xM

m =1

m

+ σ0 (εn en )|x0 + σ1 (εn en )|xM . By the definition of P in (21), we have and at xm one has



Im

enxx εn dx = 0 for each m because the enxx are polynomials of degree at most k − 2,

n n n n n εxn = ∂ x P u − ∂ x u = 0, { ε } = { P u } − u = 0,

εn |x0 = ( P un − un )|x0 = 0, εn |xM = ( P un − un )|xM = 0. Consequently A (εn , en ) = 0. Thus (27) simplifies to

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

25

    ( D αN en , en ) + A (en , en ) + c (x)en , en = χ n , en ,

(28)

where χ n = D αN εn + c (x)εn − ϕ n . To complete the proof, we shall apply Lemma 4.3 to (28). Thus we must estimate the terms on the right-hand sides of (17a) and (17b) that correspond to (28).   First, by [4, equation (53)], one has χ n  ≤ C hk+1 + n− min{2−α ,r α } . Now take μn = en and g n = χ n in (16) and invoke   (17a) from Lemma 4.3. By χ s  ≤ C hk+1 + s− min{2−α ,r α } for s = 1, . . . , n and Lemma 4.6 we obtain

 e  ≤ E α (2M 1 tnα ) e 0  + max

j 

n

1≤ j ≤n

s =1



( j) Q j −s χ s 

  j   ( j )  k +1 α 0 − min{2−α ,r α } Q j −s h +s ≤ E α (2M 1 tn ) e  + C max 1≤ j ≤n

s =1

  α k +1 − min{2−α ,r α } , ≤ C E α (2M 1tn ) h +N

(29)

where we used e 0  ≤ ε 0  ≤ Chk+1 , which follows from (e 0 , v h ) = (ε 0 , v h ) ∀ v h ∈ V h (take v h = e 0 ) and Lemma 4.5. If instead we invoke (17b) from Lemma 4.3, then using (29) and e 0  ≤ Chk+1 we get

2

γ en 1 + γ

M  β0 m =1

h

dn,1 (2 − α ) n 2 max e i 2 + χ  (2 − α ) 0≤i ≤n−1 dn,1  2 ≤ C τn−α hk+1 + N − min{2−α ,r α } + C τnα (h2k+2 + n−2 min{2−α ,r α } )

2 [en ]m ≤

−α

≤ C τn

−α

≤ C τn

 h

 h

2k+2

2k+2

+N

−2 min{2−α ,r α }

+N

−2 min{2−α ,r α }

2α −2 min{2−α ,r α } n n



+τ

+T

2α

 n 2r α N

n

−2 min{2−α ,r α }



   n 2r α −2 min{2−α ,r α } N −2 min{2−α ,r α } ≤ C τn−α h2k+2 + N −2 min{2−α ,r α } + T 2α N

  ≤ C τn−α h2k+2 + N −2 min{2−α ,r α } ,

(30)

where τn ≤ tn = T (n/ N )r was used. Combining the bounds (29), (30) and Lemma 4.5 with (25) yields the estimates claimed in the theorem. 2 The L 2 () bound of Theorem 4.7 exhibits optimal orders of convergence for all values of n. From the definition of the 1− 1



τn ≈ C N −1 tn r = Cnr −1 N −r for all n. Consequently the mesh condition τN ≤ 2(2 − −1/α α)M1 of the theorem is mild unless α is close to 0. The result indicates that one obtains the highest orders of convergence when r ≥ (2 − α )/α . On the other hand, the H 1 () bound of the theorem is suboptimal because of the factor τn−α . We have been unable graded mesh, one sees easily that

to improve this bound. The derivation of optimal H 1 () bounds has not yet been achieved in the research literature for the L1 scheme, and many published papers (in all areas of numerical analysis) contain suboptimal error estimates; thus we included our H 1 () bound since it does give some information about the performance of our method.

Remark 4.1. The constants C in our error estimates will blow up as α → 1. This is undesirable; but at far as we know, the same weakness appears in every published analysis of the L1 scheme for time-fractional problems resembling our problem (1), because — as in the derivation of (23) — a term like (1 − α ) usually creeps into the analysis. Otherwise, we believe that our constants C are well-behaved as functions of the data of the problem and of the mesh grading r. Remark 4.2. If the Robin boundary conditions (1b) are replaced by Dirichlet boundary conditions, and one still has c < 0, then the stability and convergence results obtained in this section remain valid.

26

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

Table 1 L ∞ ( L 2 ) errors and orders of convergence on temporal direction for Example 5.1 with r = (2 − α )/α . M = N = 32

M = N = 64

M = N = 128

M = N = 256

M = N = 512

M = N = 1024

α = 0.4

5.7930E-3

α = 0.6

6.8536E-3

α = 0.8

9.2792E-3

2.3344E-3 1.3112 2.9032E-3 1.2391 4.2941E-3 1.1116

8.7485E-4 1.4159 1.1764E-3 1.3032 1.9331E-3 1.1514

3.1426E-4 1.4770 4.6468E-4 1.3401 8.5748E-4 1.1727

1.0995E-4 1.5150 1.8075E-4 1.3622 3.7724E-4 1.1845

3.7812E-5 1.5400 6.9644E-5 1.3759 1.6520E-4 1.1912

Table 2 L ∞ ( H 1 ) errors and orders of convergence on temporal direction for Example 5.1 with r = (2 − α )/α . M = N = 32

M = N = 64

M = N = 128

M = N = 256

M = N = 512

M = N = 1024

α = 0.4

3.9718E-3

α = 0.6

5.0238E-3

α = 0.8

7.6285E-3

1.4349E-3 1.4688 1.9721E-3 1.3490 3.3614E-3 1.1823

5.0435E-4 1.5084 7.6486E-4 1.3665 1.4747E-3 1.1885

1.7402E-4 1.5351 2.9423E-4 1.3782 6.4513E-4 1.1928

5.9297E-5 1.5532 1.1259E-4 1.3858 2.8169E-4 1.1954

2.0028E-5 1.5658 4.2935E-5 1.3908 1.2285E-4 1.1971

5. Numerical experiments In this section, we present two numerical examples to check the theoretical results of Theorem 4.7. In both examples we take β0 = 4 and β1 = 1/12 in the L1-DDG method (12). We take r = (2 − α )/α to get the best possible rates of convergence in Theorem 4.7. The L ∞ ( L 2 ) and L ∞ ( H 1 ) errors in the computed solutions are defined by

u (x, tn ) − unh L ∞ ( L 2 ) = max u (x, tn ) − unh ,

u (x, tn ) − unh L ∞ ( H 1 ) = max u x (x, tn ) − (uh )nx .

0≤n≤ N

0≤n≤ N

Example 5.1. Consider the problem

⎧ α for (x, t ) ∈ (0, 2) × (0, 1], ⎪ ⎨ D t u − u xx + 2(x − 1)u = f (x, t ) u (0, t ) − u x (0, t ) = 0, u (2, t ) + u x (2, t ) = 0 for t ∈ (0, 1], ⎪ ⎩ u (x, 0) = 13 x3 − x2 + 13 x + 13 for x ∈ [0, 2], where the function f (x, t ) is chosen such that the exact solution of the problem is



u (x, t ) = E α (−t α ) + t 3

1 3

3

2

x −x +

1 3

x+

1 3



,

where the Mittag-Leffler function E α (·) was defined in Section 1. The derivatives of this solution are precisely in accordance with the bounds stated in Lemma 2.2. In Tables 1–5 we give numerical experiments for this example, with α = 0.4. Tables 1 and 2 display the L ∞ ( L 2 ) and L ∞ ( H 1 ) errors and the associated orders of convergence; in these tables we have taken M = N so that the temporal error N − min{2−α ,r α } = N −(2−α ) dominates the results. In Tables 3–5 we present the L 2 () error, the H 1 () error and the error

 M

m=1 [u

− unh ]m /h

1/2

, and their associated orders of convergence, at the time level t N = 1. Here we took N = 10000 so that the spatial errors dominate the results. The numerical experiments agree precisely with the theoretical rate of convergence of Theorem 4.7 for the L 2 errors: Table 1 shows orders O ( N −(2−α ) ) and Table 3 shows O (hk+1 ). For the H 1 -type errors of (24b), it appears that our theoretical bound is suboptimal by the factor τn−α , because Table 2 shows O ( N −(2−α ) ) convergence while Tables 4 and 5 show O (hk ) convergence. n

Example 5.2. The second numerical example demonstrates the effect of negative c (·) in (1). We solve the problem

⎧ α for (x, t ) ∈ (0, 2) × (0, 1], ⎪ ⎪ D t u − u xx − b(2 − x)u = 0 ⎨ u (0, t ) − u x (0, t ) = 0, u (2, t ) + u x (2, t ) = 0 for t ∈ (0, 1], ⎪ ⎪ ⎩ u (x, 0) = sin( 12 π x) for x ∈ [0, 2], λi in where b is a positive constant. One has M 1 = 2b, where M 1 was defined in Section 1. A finite number of eigenvalues  (4) will be negative, so the solution will grow rapidly as t increases.

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

27

Table 3 L 2 errors and orders of convergence on space direction for Example 5.1. Polynomial

M

u − u h  L 2

Order

P1

5 10 20 40 5 10 20 40

2.8213E-2 7.5687E-3 1.9236E-3 4.7501E-4 1.4342E-3 1.8203E-4 2.2851E-5 2.9060E-6

– 1.8982 1.9762 2.0177 – 2.9779 2.9938 2.9752

P2

Table 4 H 1 errors and orders of convergence at t = 1 for Example 5.1 with

α = 0.4.

Polynomial

M

|u − uh |1

Order

P1

5 10 20 40 5 10 20 40

2.8978E-1 1.4363E-1 7.1457E-2 3.5669E-2 2.6846E-2 6.7027E-3 1.6752E-3 4.1882E-4

– 1.0125 1.0072 1.0023 – 2.0018 2.0003 1.9999

P2

Table 5  M

m=1 [u

n

− unh ]m /h

1/2

errors and orders of convergence at t = 1 for Example 5.1 with

Polynomial

M

P1

5 10 20 40 5 10 20 40

P2

 M

m=1 [u

n

α = 0.4.

1/2 − unh ]m /h

8.7583E-2 2.4681E-2 6.4915E-3 1.6616E-4 4.3261E-4 2.8230E-5 1.7810E-6 1.1067E-7

Fig. 1. The numerical solution uh for b = 1 (left) and b = 1.5 (right) of Example 5.2 with

Order – 1.8272 1.9267 1.9659 – 3.9377 3.9863 4.0084

α = 0.6.

Figs. 1–2 show the numerical solution for b = 1, 1.5, 2, 3 with α = 0.6, N = 200 and r = (2 − α )/α . The exact solution of the Example 5.2 is unknown, so the two-mesh principle [2] is used to estimate the orders of convergence in our computed solutions. We now describe this technique. Let unh with 0 ≤ n ≤ N and 1 ≤ m ≤ M be the solution computed by our scheme (11). Then consider a second mesh that is uniform in the spatial direction with 2M mesh intervals, and in the temporal direction is defined by

28

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

Fig. 2. The numerical solution u h for b = 2 (left) and b = 3 (right) of Example 5.2 with

α = 0.6.

Table 6 L 2 errors and orders of convergence on temporal direction for Example 5.2 with r = (2 − α )/α . M = N = 64

M = N = 128

M = N = 256

M = N = 512

α = 0.4

6.1699E-1

α = 0.6

5.0511E-2

α = 0.8

1.9274E-2

2.1466E-1 1.5231 1.9723E-3 1.3566 8.2407E-3 1.2258

7.4776E-2 1.5214 7.6628E-4 1.3639 3.5375E-3 1.2200

2.5802E-2 1.5350 2.9575E-4 1.3735 1.5215E-3 1.2171



r

tn = n/(2N )

for 0 ≤ n ≤ 2N .

The solution computed on this mesh is denoted by znh with 0 ≤ n ≤ 2N and 1 ≤ m ≤ 2M. Now the L 2 error at the final time T = 1 is taken to be

D M , N := u hN − zh2N . Then the rates of convergence are computed by

 log2

 . 2M ,2N

D M ,N D

These L () errors and their corresponding estimated orders of convergence for Example 5.2 with b = 2 are given in Table 6. The numerical results agree with Theorem 4.7: they show that r = (2 − α )/α yields the optimal rate of convergence O ( N −(2−α ) ). 2

Acknowledgements The research of Martin Stynes is supported in part by the National Natural Science Foundation of China under grants 91430216 and NSAF U1530401. References [1] R. Courant, D. Hilbert, Methods of Mathematical Physics. Vol. I, Interscience Publishers, Inc., New York, NY, 1953. [2] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Applied Mathematics (Boca Raton), vol. 16, Chapman & Hall/CRC, Boca Raton, FL, 2000. [3] Daniel Henry, Geometric Theory of Semilinear Parabolic Equations, Lecture Notes in Mathematics, vol. 840, Springer-Verlag, Berlin–New York, 1981. [4] Chaobao Huang, Martin Stynes, Na An, Optimal L ∞ ( L 2 ) error analysis of a direct discontinuous Galerkin method for a time-fractional reaction–diffusion equation, BIT Numer. Math. (2018), https://doi.org/10.1007/s10543-018-0707-z, in press. [5] 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.

C. Huang, M. Stynes / Applied Numerical Mathematics 135 (2019) 15–29

29

[6] Hailiang Liu, Optimal error estimates of the direct discontinuous Galerkin method for convection–diffusion equations, Math. Comput. 84 (295) (2015) 2263–2295. [7] Hailiang Liu, Jue Yan, The direct discontinuous Galerkin (DDG) method for diffusion with interface corrections, Commun. Comput. Phys. 8 (3) (2010) 541–564. [8] Yury Luchko, Maximum principle for the generalized time-fractional diffusion equation, J. Math. Anal. Appl. 351 (1) (2009) 218–223. [9] Yuri Luchko, Masahiro Yamamoto, On the maximum principle for a time-fractional diffusion equation, Fract. Calc. Appl. Anal. 20 (5) (2017) 1131–1145. [10] José Luis Gracia, Eugene O’Riordan, Martin Stynes, Convergence in positive time for a finite difference method applied to a fractional convection– diffusion problem, Comput. Methods Appl. Math. 18 (1) (2018) 33–42. [11] Igor Podlubny, Fractional Differential Equations. An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, Mathematics in Science and Engineering, vol. 198, Academic Press, Inc., San Diego, CA, 1999. [12] 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. [13] 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.