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 . It is proved that the rate of convergence of the maximum nodal error on any subdomain that is bounded away from 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.
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 ; that is, their temporal derivatives are unbounded on the closed space-time domain, but are bounded on any subdomain that is bounded away from . 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 .
This paper is a companion paper to , where it was shown that the convergence rate of the same finite difference scheme on a uniform mesh was , where is the order of the fractional derivative and the mesh spacing in time is . Results related to the main result in  are available in [7, 8], using a finite element framework. In contrast to , we shall prove here for the same scheme on a uniform mesh that the convergence rate of the numerical solution is on any subdomain that is bounded away from .
Our analysis is carried out in the discrete norm; an analogous convergence result in the norm was derived in . 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  for a homogeneous version of (2.1) in the case of non-smooth initial data.
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 the ceiling function is the smallest integer greater than or equal to x. For any continuous function with and any mesh function with and , we set
2 The Continuous Problem
Consider the initial-boundary value problem
is the Riemann–Liouville fractional integral operator of order .
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
problem (2.1) becomes
Note that no first-order derivative in space appears in (2.2a), and (2.1d) implies that . Consequently, (2.2) belongs to the class of problems analysed in . In  it was assumed that was a positive constant, but the analysis of  can be extended to the case of a smooth variable positive coefficient . 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:
for all .
for all and some constant , where is a constant independent of t and is the norm associated with the vector space . This space is defined by
where is the inner product in the Hilbert space and are the eigenvalues and normalised eigenfunctions of the Sturm–Liouville two-point boundary value problem
3 The Discrete Problem
The solution of problem (2.1) is approximated by the solution of a finite difference scheme on a mesh , that is uniform in both space and time. Let M and N be positive integers. Set and for . Set for . The nodal approximation to the solution u computed at the mesh point is denoted by .
The first and second-order spatial derivatives are discretised using standard approximations:
The Caputo fractional derivative , which can be written as
is approximated by the classical L1 approximation
Thus, we approximate (2.1) by the discrete problem
This discretisation of (2.1) is standard.
To ensure the stability of the discrete operator by imposing the correct sign pattern in the associated matrix, we make the nonrestrictive assumption that N satisfies
for some constant C. In particular, the method has the low order of convergence in time when α is small. In the present paper we shall consider the rate of convergence in a subdomain , 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 . In this analysis the truncation error bound (4.2) indicates that the truncation error decreases as one moves further away from the initial time . 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
The truncation error in time is more tricky to estimate and this is done in the next lemma.
Assume that u satisfies (2.3). Then there exists a positive constant C such that for each mesh point one has
where for and we define the truncation error in the kth time cell to be
The following four bounds are established in [10, equations (5.9), (5.10), (5.11) and (5.14)]:
For , it is clear that and . Thus, we see that
by the Mean Value Theorem. Combine this bound with (4.4) to complete the proof. ∎
for where the positive weights are defined for by the recurrence relation
The coefficients satisfy for .
First, as for all , so the lemma is true when . The proof is completed by induction. Assume that for . We want to prove that . It is easy to check that for all k. Using this inequality and the inductive hypothesis, we require the inequality
which is established in [5, Lemma 3.2]. ∎
Let the parameter β satisfy . Then, for , one has
By Lemma 2, we have
But for and , one has
This completes the proof. ∎
We can now prove our main result.
for some constant C.
By (4.5) we then obtain
The bound in (4.8) implies that for any fixed one has
That is, on any subdomain that is bounded away from , we observe an improved rate of convergence in time compared with the rate of convergence (in time) of on 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 . In these numerical experiments we always take . Hence the bounds in (3.3) and (4.8) imply that the spatial error term will be dominated by the temporal error term or .
Consider the constant coefficient homogeneous problem
with initial condition , , and boundary conditions , . The exact solution of this problem is
where is the Mittag-Leffler function which is defined  by
For Example 5.1 we computed the maximum errors
and the orders of convergence
where can be the entire domain or the subdomain . The numerical results in (see Table 1) show that scheme (3.2) is convergent there (which agrees with [10, Theorem 5.2]), while it is convergent in the subdomain (see Table 2), which indicates that the error bound (4.9) is sharp.
Considering the convergence in time, identified by the factor in the error bound (4.8), the error is compared with the error bound in Figure 2, for and . These plots indicate that the exponent 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 , as increases.
Consider the variable coefficient inhomogeneous problem
Figure 3 displays the computed solution for and and we observe that the solution again exhibits an initial layer at .Figure 3
The exact solution of Example 5.2 is unknown and we shall estimate the order of convergence using the two-mesh principle . Let be the computed solution with scheme (3.2) on the mesh for , . To estimate the order of convergence, we compute a new approximation using the same scheme defined on the finer mesh for , , where and . We then compute the two-mesh differences
and hence the estimated orders of convergence
Tables 3 and 4 give the maximum two-mesh differences and their corresponding orders of convergence for Example 5.2 in the domain and the subdomain . The numerical results in both cases are again in agreement with Theorem 4: the order of convergence improves from on to on .
In , numerical results were given for the particular case of a fractional reaction-diffusion equation (i.e., with 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 are given in .
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.
 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
 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
 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
 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
 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
 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
 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
 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
© 2017 Walter de Gruyter GmbH, Berlin/Boston