# Convergence in Positive Time for a Finite Difference Method Applied to a Fractional Convection-Diffusion Problem

José Luis Gracia , Eugene O’Riordan and Martin Stynes

## Abstract

A standard finite difference method on a uniform mesh is used to solve a time-fractional convection-diffusion initial-boundary value problem. Such problems typically exhibit a mild singularity at the initial time t=0. It is proved that the rate of convergence of the maximum nodal error on any subdomain that is bounded away from t=0 is higher than the rate obtained when the maximum nodal error is measured over the entire space-time domain. Numerical results are provided to illustrate the theoretical error bounds.

MSC 2010: 65M06; 65M12; 65M15

### 1 Introduction

In this paper, we examine the convergence rate of numerical approximations to a time-fractional convection-diffusion problem using a standard finite difference method on a uniform mesh. Initial-boundary value problems of this type, where the time derivative is fractional, have solutions that are mildly singular at the initial time t=0; that is, their temporal derivatives are unbounded on the closed space-time domain, but are bounded on any subdomain that is bounded away from t=0. It is shown in [9, 10] that the case of solutions with bounded temporal derivatives on the closed space-time domain is very special and that the weakly singular solutions examined here are much more typical of how solutions to this class of problems behave. As one would expect, the rate of convergence of the computed numerical approximations is affected adversely by the presence of large temporal derivatives at t=0.

This paper is a companion paper to [10], where it was shown that the convergence rate of the same finite difference scheme on a uniform mesh was O(N-α), where α(0,1) is the order of the fractional derivative and the mesh spacing in time is O(N-1). Results related to the main result in [10] are available in [7, 8], using a finite element framework. In contrast to [10], we shall prove here for the same scheme on a uniform mesh that the convergence rate of the numerical solution is O(N-1) on any subdomain that is bounded away from t=0.

Our analysis is carried out in the discrete L norm; an analogous convergence result in the L2 norm was derived in [4]. Using an alternative formulation of the continuous problem, the phenomenon of higher-order convergence at some fixed distance away from the initial singularity is examined in [6] for a homogeneous version of (2.1) in the case of non-smooth initial data.

#### Notation.

In this paper C denotes a generic constant that depends on the data of the boundary value problem (2.1) but is independent of T and of any mesh used to solve (2.1) numerically. Note that C can take different values in different places. For all x the ceiling function x is the smallest integer greater than or equal to x. For any continuous function z:Q with Q2 and any mesh function zmn with n=0,1,,N and m=0,1,,M, we set

z:=max(x,t)Q¯|z(x,t)|andzn:=max0mM|zmn|.

### 2 The Continuous Problem

Consider the initial-boundary value problem

(2.1)

(2.1a)Lu:=Dtαu-p(x)2ux2+q(x)ux+r(x)u=f(x,t)
for (x,t)Q:=(0,l)×(0,T], with initial and boundary conditions
(2.1b)u(0,t)=u(l,t)=0for t(0,T],
(2.1c)u(x,0)=ϕ(x)for x[0,l].
Here 0<α<1p(x)p0>0 on [0,l], the functions p(x),q(x) and r(x) are smooth on [0,l] and are assumed to satisfy the constraint
(2.1d)r(x)-q(x)+p′′(x)20for (x,t)Q.

The initial condition ϕ is also smooth on [0,l] and the function f is smooth on Q¯. Furthermore, in (2.1a) Dtα denotes the Caputo fractional derivative which is defined [1] by

Dtαg(x,t):=[J1-α(gt)](x,t)for 0xl, 0<tT,

where

(J1-αg)(x,t):=[1Γ(1-α)s=0t(t-s)-αg(x,s)𝑑s]for (x,t)Q¯

is the Riemann–Liouville fractional integral operator of order 1-α.

There is no loss of generality in assuming homogeneous boundary conditions in (2.1b), because inhomogeneous boundary conditions are easily made homogeneous by a simple change of variable.

Under the transformation

y(x,t):=u(x,t)p(0)p(x)e-120xq(s)p(s)𝑑s,

problem (2.1) becomes

(2.2)

(2.2a)Dtαy-p(x)2yx2+r1(x)y=f1(x,t)for (x,t)Q:=(0,l)×(0,T],
(2.2b)y(0,t)=y(l,t)=0for t(0,T],
(2.2c)y(x,0)=ϕ1(x):=ϕ(x)p(0)p(x)e-120xq(s)p(s)𝑑sfor x[0,l],
(2.2d)wherer1(x):=r(x)-(q(x)+p′′(x))2+(q(x)+p(x))24p,
(2.2e)f1(x,t):=f(x,t)p(0)p(x)e-120xq(s)p(s)𝑑s.

Note that no first-order derivative in space appears in (2.2a), and (2.1d) implies that r10. Consequently, (2.2) belongs to the class of problems analysed in [10]. In [10] it was assumed that p(x) was a positive constant, but the analysis of [10] can be extended to the case of a smooth variable positive coefficient p(x). Thus after placing suitable regularity and compatibility conditions on the data of (2.2), one can invoke [10, Theorem 2.1] to conclude that (2.2) has a unique solution y whose derivatives satisfy certain bounds. Transforming back to the original problem (2.1), under certain conditions on its data one obtains the following bounds on the derivatives of u:

(2.3)

(2.3a)|kuxk(x,t)|Cfor k=0,1,2,3,4,
(2.3b)|ut(x,t)|C(1+tα-)for =1,2,

for all (x,t)[0,l]×(0,T].

In [10, Theorem 2.1] the estimates in (2.3) are proved assuming that ϕ1D(5/2), (f1)(,t)D(5/2), (f1)t(,t) and (f1)tt(,t) are in D(1/2) for each t(0,T] and

(f1)(,t)5/2+(f1)t(,t)1/2+tρ(f1)tt(,t)1/2C1

for all t(0,T] and some constant ρ<1, where C1 is a constant independent of t and γ is the norm associated with the vector space D(γ). This space is defined by

D(γ):={gL2(0,l):i=1λi2γ|(g,ψi)|2<},γ0,

where (,) is the inner product in the Hilbert space L2(0,l) and {(λi,ψi):i=1,2,} are the eigenvalues and normalised eigenfunctions of the Sturm–Liouville two-point boundary value problem

ψi:=-pψi′′+cψi=λiψion (0,l),ψi(0)=ψi(l)=0.

### 3 The Discrete Problem

The solution of problem (2.1) is approximated by the solution of a finite difference scheme on a mesh {(xm,tn):m=0,1,,M,n=0,1,,N}, that is uniform in both space and time. Let M and N be positive integers. Set h=lM and xm:=mh for m=0,1,,M. Set tn=nτ=nTN for n=0,1,,N. The nodal approximation to the solution u computed at the mesh point (xm,tn) is denoted by umn.

The first and second-order spatial derivatives are discretised using standard approximations:

ux(xm,tn)Dx0umn:=um+1n-um-1n2h,
2ux2(xm,tn)δx2unm:=um+1n-2umn+um-1nh2.

The Caputo fractional derivative Dtαu, which can be written as

Dtαu(xm,tn)=1Γ(1-α)k=0n-1s=tktk+1(tn-s)-αu(xm,s)s𝑑s,

is approximated by the classical L1 approximation

(3.1)

DNαumn:=1Γ(1-α)k=0n-1umk+1-umkτs=tktk+1(tn-s)-α𝑑s
(3.1a)=1Γ(2-α)k=0n-1(umk+1-umk)dn-k,
where
(3.1b)dk:=k1-α-(k-1)1-α,k1.

Thus, we approximate (2.1) by the discrete problem

(3.2)

(3.2a)LM,Numn:=DNαumn-p(xm)δx2umn+q(xm)Dx0umn+r(xm)umn=f(xm,tn)for 1mM-1, 1nN,
(3.2b)u0n=uMn=0for 0<nN,
(3.2c)um0=ϕ(xm)for 0mM.

This discretisation of (2.1) is standard.

To ensure the stability of the discrete operator LM,N by imposing the correct sign pattern in the associated matrix, we make the nonrestrictive assumption that N satisfies

lq2p0<N.

After some minor modifications in the proof of [10, Theorem 5.2] to handle the term q(xm)Dx0umn, it follows that the solution umn of scheme (3.2) satisfies the error bound

(3.3)max(xm,tn)Q¯|u(xm,tn)-umn|C(h2+N-α)

for some constant C. In particular, the method has the low order of convergence O(N-α) in time when α is small. In the present paper we shall consider the rate of convergence in a subdomain [0,l]×[κ,T], where κ is a fixed positive value.

### 4 Error Analysis

The structure of our error analysis is the standard finite difference technique of estimating the truncation error at each mesh point, then invoking a stability argument to derive an error bound for the computed solution umn. In this analysis the truncation error bound (4.2) indicates that the truncation error decreases as one moves further away from the initial time t=0. The stability bound (4.5) shows that the error at any discrete time level depends on a weighted sum of the truncation errors at all the previous time levels.

The estimate of the truncation error in space is standard: using (2.3a), one gets

(4.1)

ux(xm,tn)=Dx0u(xm,tn)+O(h2),
2ux2(xm,tn)=δx2u(xm,tn)+O(h2).

The truncation error in time is more tricky to estimate and this is done in the next lemma.

### Lemma 1.

Assume that u satisfies (2.3). Then there exists a positive constant C such that for each mesh point (xm,tn)Q one has

(4.2)|(DNα-Dtα)u(xm,tn)|Cn-min{2-α,α+1}.

### Proof.

We modify the argument of [10, Lemma 5.1]. By (3.1a) and the definition of Dtα, for each mesh point (xm,tn)Q one has

(DNα-Dtα)u(xm,tn)=k=0n-1Tnk,

where for n=1,2,,N and k=0,1,,n-1 we define the truncation error in the kth time cell [tk,tk+1] to be

(4.3)Tnk:=1Γ(1-α)s=tktk+1(tn-s)-α[u(xm,tk+1)-u(xm,tk)τ-us(xm,s)]𝑑s.

The following four bounds are established in [10, equations (5.9), (5.10), (5.11) and (5.14)]:

(4.4)

(4.4a)k=1n2-1|Tnk|Cn-(α+1)for 1k<n-1,
(4.4b)k=n2n-2|Tnk|Cn-(2-α)for 1k<n-1,
(4.4c)|T10|C,
(4.4d)|Tn,n-1|Cn-(2-α).

It remains to bound |Tn0| for n>1. We sharpen the bound [10, (5.12)] of this term. An integration by parts in (4.3) yields

Tn0=-αΓ(1-α)s=0t1(tn-s)-α-1(ϕ-ψ)(xm,s)𝑑s,

where

ϕ(xm,s):=s[u(xm,t1)-u(xm,0)τ]andψ(xm,s):=u(xm,s)-u(xm,0).

For 0sτ, it is clear that |ϕ(xm,s)||u(xm,τ)-u(xm,0)| and |ψ(xm,s)|0τ|ut(xm,t)|𝑑t. Thus, we see that

|ϕ(xm,s)|+|ψ(xm,s)|20τ|ut(xm,t)|𝑑tC0τ(1+tα-1)𝑑tCτα.

Hence

|Tn0|Cταs=0t1(tn-s)-α-1𝑑sCτα[(tn-t1)-α-tn-α]=C[(n-1)-α-n-α]Cn-(α+1),

by the Mean Value Theorem. Combine this bound with (4.4) to complete the proof. ∎

Observe that min{2-α,α+1}>1 in (4.2) for all values of α(0,1); thus this bound is sharper than the truncation error bound of O(n-α) proved in [10, Lemma 5.1]. This improvement is critical in establishing our main result later.

Next, we derive some new information about the stability constants that appear in [10, Section 4]. It follows from [10, Lemma 4.2] that the computed solution umn of (3.2) satisfies

(4.5)u(xm,tn)-umnταΓ(2-α)j=1nσn-j(LM,N(u(xm,tj))j-uj

for n=1,2,,N, where the positive weights σi are defined for i=0,1,2,,n-1 by the recurrence relation

(4.6)σ0:=1,σi:=k=1i(dk-dk+1)σi-kfor i=1,2,.

Note that when the mesh is uniform, the weights θn,j in [10, Lemma 4.2] are the same as the weights σn-j defined in (4.6).

### Lemma 2.

The coefficients σi satisfy σi<(i+1)α-1 for i=1,2,.

### Proof.

First, σ1=(d1-d2)σ0=2-21-α<2-1+α as 0.5w+2w-1-2>0 for all w(1,2), so the lemma is true when i=1. The proof is completed by induction. Assume that σj<(j+1)α-1 for j=1,2,,i-1. We want to prove that σi<(i+1)α-1. It is easy to check that dk-dk+1>0 for all k. Using this inequality and the inductive hypothesis, we require the inequality

k=1i(dk-dk+1)(i+1-k)α-1<(i+1)α-1,

which is established in [5, Lemma 3.2]. ∎

The next result, which is a variant of [10, Lemma 4.3], bounds a weighted sum of the σn-j that will be used in the proof of Theorem 4.

### Lemma 3.

Let the parameter β satisfy β>1. Then, for n=1,2,,N, one has

ταj=1nj-βσn-jCN-1Ttnα-1.

### Proof.

By Lemma 2, we have

(4.7)ταj=1nj-βσn-jταj=1nj-β(n+1-j)α-1τα[(n2)α-1j=1n2j-β+(n2)-(β-α)j=n2+1nj-α(n+1-j)α-1].

But for sj-1 and jn, one has

(n+1-s)α-1(n+2-j)α-12α-1(n+1-j)α-1.

Hence,

j=n2+1nj-α(n+1-j)α-1j=n2+1n(n+1-j)α-1s=j-1js-α𝑑s
j=n2+1n21-αs=j-1j(n+1-s)α-1s-α𝑑s
21-αs=0n+1(n+1-s)α-1s-α𝑑s
=21-αΓ(α)Γ(1-α),

by [1, Theorem D.6]. Substituting this inequality into (4.7) and using tn=nτ and β>1, we get

ταj=1nj-βσn-jCταnα-1j=1n2j-β+Cταn-(β-α)Cταnα-1+Cταn-(β-α)=Ctnα-1(τ+τn-(β-1))Ctnα-1N-1T.

This completes the proof. ∎

We can now prove our main result.

### Theorem 4.

Assume that u satisfies (2.3). Then, for n=1,2,3,,N, the solution umn of scheme (3.2) satisfies

(4.8)max0mM|u(xm,tn)-umn|C(Tαh2+TN-1tnα-1)

for some constant C.

### Proof.

Fix (xm,tn)Q. By (4.1) and Lemma 1, the truncation error at (xm,tn) satisfies

(LM,N(u(xm,tn))n-un)C(h2+n-min{2-α,1+α}).

By (4.5) we then obtain

max0mM|u(xm,tn)-umn|Cταj=1nh2σn-j+Cταj=1nj-min{2-α,1+α}σn-j.

Invoking Lemma 3 (with β=min{2-α,1+α}) for the j-min{2-α,1+α} term and [10, Lemma 4.3] (with β=0) for the term involving h2, we obtain (4.8). ∎

The bound in (4.8) implies that for any fixed κ>0 one has

(4.9)max(xm,tn)Q¯{tnκ>0}|(xm,tn)-umn|CTα(h2+N-1).

That is, on any subdomain that is bounded away from t=0, we observe an improved rate of convergence in time compared with the rate of convergence (in time) of N-α on Q¯ that is given by (3.3).

### 5 Numerical Results

In this section we give numerical results for the numerical method (3.2) applied to two particular examples from the problem class (2.1). In the first example the exact solution of the problem is known; in the second example it is unknown, so we estimate the order of convergence using the double-mesh principle [2]. In these numerical experiments we always take N=M. Hence the bounds in (3.3) and (4.8) imply that the spatial error term Ch2 will be dominated by the temporal error term CN-α or CN-1.

### Example 5.1.

Consider the constant coefficient homogeneous problem

Dtαu-2ux2+ux=0for (x,t)Q=(0,π)×(0,1],

with initial condition u(x,0)=ex2sinx, 0<x<π, and boundary conditions u(0,t)=u(π,t)=0, 0t1. The exact solution of this problem is

u(x,t)=Eα(-1.25tα)ex2sinx,

where Eα is the Mittag-Leffler function which is defined [1] by

Eα(z):=k=0zkΓ(αk+1).

In Figure 1 we display the computed solutions with scheme (3.2) for α=0.4,0.8 and N=M=32 and we observe that the solution has an initial layer at t=0, which becomes sharper as the parameter α decreases.

Figure 1

Example 5.1: Computed solutions with scheme (3.2) for N=M=32.

α=0.4

### (b)

α=0.8

For Example 5.1 we computed the maximum errors

eM,N:=max(xm,tn)Q|u(xm,tn)-umn|

and the orders of convergence

pM,N:=log2(eM,Ne2M,2N),

where Q can be the entire domain Q¯ or the subdomain Q¯*:=[0,π]×[0.1,1]. The numerical results in Q¯ (see Table 1) show that scheme (3.2) is O(N-α) convergent there (which agrees with [10, Theorem 5.2]), while it is O(N-1) convergent in the subdomain Q¯* (see Table 2), which indicates that the error bound (4.9) is sharp.

Considering the convergence in time, identified by the factor tnα-1 in the error bound (4.8), the error maxm|u(xm,tn)-umn| is compared with the error bound N-1tnα-1 in Figure 2, for α=0.2 and α=0.8. These plots indicate that the exponent α-1 in the error bound is sharp for small values of α. However, in the case of larger values of α close to one, the maximum error decreases at a faster rate than α-1, as tn increases.

### Table 1

Example 5.1: Maximum errors and orders of convergence for scheme (3.2) in the domain Q¯.

 α N=M=128 N=M=256 N=M=512 N=M=1024 N=M=2048 0.4 8.438E-2 6.714E-2 5.282E-2 4.120E-2 3.191E-2 0.330 0.346 0.359 0.368 0.6 3.759E-2 2.512E-2 1.672E-2 1.109E-2 7.342E-3 0.581 0.588 0.592 0.595 0.8 1.121E-2 6.401E-3 3.666E-3 2.102E-3 1.206E-3 0.809 0.804 0.803 0.802

### Table 2

Example 5.1: Maximum errors and orders of convergence for scheme (3.2) in the subdomain Q¯*.

 α N=M=128 N=M=256 N=M=512 N=M=1024 N=M=2048 0.4 1.024E-2 4.966E-3 2.436E-3 1.214E-3 6.050E-4 1.044 1.027 1.005 1.005 0.6 1.300E-2 6.432E-3 3.190E-3 1.595E-3 7.965E-4 1.015 1.012 1.000 1.002 0.8 9.844E-3 5.123E-3 2.644E-3 1.361E-3 6.963E-4 0.942 0.954 0.959 0.966
Figure 2

Example 5.1: Log-log plot of the error bound N-1tnα-1 () and the maximum errors max0mM|u(xm,tn)-umn| () at each time level t=tn,n=1,2,,N, generated by scheme (3.2) for N=M=100.

α=0.2

α=0.8

### Example 5.2.

Consider the variable coefficient inhomogeneous problem

(5.1)

(5.1a)Dtαu-2ux2+(1+x2)ux+(1+x)u=4π2x(π-x)for (x,t)Q,
with Q=(0,π)×(0,1] and
(5.1b)u(x,0)=0for 0<x<π,u(0,t)=u(π,t)=0for 0t1.

Figure 3 displays the computed solution for α=0.4,0.8 and N=M=64 and we observe that the solution again exhibits an initial layer at t=0.

Figure 3

Example 5.2: Computed solutions with scheme (3.2) for N=M=64.

α=0.4

### (b)

α=0.8

The exact solution of Example 5.2 is unknown and we shall estimate the order of convergence using the two-mesh principle [2]. Let umn be the computed solution with scheme (3.2) on the mesh {(xm,tn)} for m=0,1,,M, n=0,1,,N. To estimate the order of convergence, we compute a new approximation zm/2n2 using the same scheme defined on the finer mesh {(xm/2,tn/2)} for m=0,1,,2M, n=0,1,,2N, where xm+1/2:=12(xm+1+xm) and tn+1/2:=12(tn+1+tn)/2. We then compute the two-mesh differences

dM,N:=max(xm,tn)Q|umn-zmn|

and hence the estimated orders of convergence

qM,N:=log2(dM,Nd2M,2N).

Tables 3 and 4 give the maximum two-mesh differences and their corresponding orders of convergence for Example 5.2 in the domain Q¯ and the subdomain Q¯*. The numerical results in both cases are again in agreement with Theorem 4: the order of convergence improves from O(N-α) on Q¯ to O(N-1) on Q¯*.

### Table 3

Example 5.2: Maximum two-mesh differences and orders of convergence for scheme (3.2) in the domain Q¯.

 α N=M=128 N=M=256 N=M=512 N=M=1024 N=M=2048 0.4 1.031E-2 8.673E-3 7.123E-3 5.740E-3 4.558E-3 0.250 0.284 0.311 0.333 0.6 4.935E-3 3.338E-3 2.234E-3 1.486E-3 9.857E-4 0.564 0.579 0.588 0.593 0.8 1.661E-3 9.441E-4 5.368E-4 3.060E-4 1.748E-4 0.815 0.815 0.811 0.808

### Table 4

Example 5.2: Maximum two-mesh differences and orders of convergence for scheme (3.2) in the subdomain Q¯*.

 α N=M=128 N=M=256 N=M=512 N=M=1024 N=M=2048 0.4 5.849E-4 2.783E-4 1.351E-4 6.711E-5 3.337E-5 1.072 1.042 1.010 1.008 0.6 1.148E-3 5.457E-4 2.628E-4 1.291E-4 6.356E-5 1.073 1.054 1.025 1.023 0.8 1.335E-3 6.752E-4 3.387E-4 1.703E-4 8.531E-5 0.984 0.995 0.992 0.997

In [10], numerical results were given for the particular case of a fractional reaction-diffusion equation (i.e., with q0 in (2.1)) showing that the scheme also converges with order α when the whole domain is considered. Additional numerical results that illustrate the improved rate of convergence away from t=0 are given in [3].

Funding source: National Natural Science Foundation of China

Award Identifier / Grant number: 91430216

Award Identifier / Grant number: NSAF U1530401

Funding source: Ministerio de Economía y Competitividad

Award Identifier / Grant number: MTM2016-75139-R

Funding statement: The research of José Luis Gracia was partly supported by the Institute of Mathematics and Applications (IUMA), the project MTM2016-75139-R funded by the Spanish Ministry of Economy and Competitiveness (MINECO) and the Diputación General de Aragón. The research of Martin Stynes was supported in part by the National Natural Science Foundation of China under grants 91430216 and NSAF U1530401.

#### References

[1] K. Diethelm, The Analysis of Fractional Differential Equations, Lecture Notes in Math. 2004, Springer, Berlin, 2010. 10.1007/978-3-642-14574-2Search in Google Scholar

[2] P. A. Farrell, A. F. Hegarty, J. J. H. Miller, E. O’Riordan and G. I. Shishkin, Robust Computational Techniques for Boundary Layers, Appl. Math. Math. Comput. 16, Chapman & Hall/CRC, Boca Raton, 2000. 10.1201/9781482285727Search in Google Scholar

[3] J. L. Gracia, E. O’Riordan and M. Stynes, Convergence outside the initial layer for a numerical method for the time-fractional heat equation, Numerical Analysis and Its Applications, Lecture Notes in Comput. Sci. 10187, Springer, Cham (2017), 82–94. 10.1007/978-3-319-57099-0_8Search in Google Scholar

[4] B. Jin, R. Lazarov and Z. Zhou, An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data, IMA J. Numer. Anal. 36 (2016), no. 1, 197–221. 10.1093/imanum/dru063Search in Google Scholar

[5] B. Jin and Z. Zhou, An analysis of Galerkin proper orthogonal decomposition for subdiffusion, ESAIM Math. Model. Numer. Anal. 51 (2017), no. 1, 89–113. 10.1051/m2an/2016017Search in Google Scholar

[6] W. McLean and K. Mustapha, Time-stepping error bounds for fractional diffusion problems withnon-smooth initial data, J. Comput. Phys. 293 (2015), 201–217. 10.1016/j.jcp.2014.08.050Search in Google Scholar

[7] K. Mustapha, An implicit finite-difference time-stepping method for a sub-diffusion equation, with spatial discretization by finite elements, IMA J. Numer. Anal. 31 (2011), no. 2, 719–739. 10.1093/imanum/drp057Search in Google Scholar

[8] K. Mustapha and W. McLean, Piecewise-linear, discontinuous Galerkin method for a fractional diffusion equation, Numer. Algorithms 56 (2011), no. 2, 159–184. 10.1007/s11075-010-9379-8Search in Google Scholar

[9] M. Stynes, Too much regularity may force too much uniqueness, Fract. Calc. Appl. Anal. 19 (2016), no. 6, 1554–1562. 10.1515/fca-2016-0080Search in Google Scholar

[10] M. Stynes, E. O’Riordan and J. L. Gracia, Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation, SIAM J. Numer. Anal. 55 (2017), no. 2, 1057–1079. 10.1137/16M1082329Search in Google Scholar