Pointwise-in-time a posteriori error control for time-fractional parabolic equations

For time-fractional parabolic equations with a Caputo time derivative of order $\alpha\in(0,1)$, we give pointwise-in-time a posteriori error bounds in the spatial $L_2$ and $L_\infty$ norms. Hence, an adaptive mesh construction algorithm is applied for the L1 method, which yields optimal convergence rates $2-\alpha$ in the presence of solution singularities.

where Γ(·) is the Gamma function, and ∂ t denotes the partial derivative in t.
Although there is a substantial literature on the a posteriori error estimation for classical parabolic equations, the pointwise-in-time a posteriori error control appears an open question for equations of type (1.1) (the few papers for similar problems give error estimates in global fractional Sobolev space norms [13]).
In this paper, we shall address this question by deriving pointwise-in-time a posteriori error bounds in the L 2 (Ω) and L ∞ (Ω) norms. Furthermore, explicit upper barriers on the residual will be obtained that guarantee that the error remains within a prescribed tolerance and within certain desirable pointwise-in-time error profiles. These residual barriers naturally lead to an adaptive mesh construction algorithm, which will be applied for the popular L1 method. It will be demonstrated that the constructed adaptive meshes successfully detect solution singularities and yield optimal convergence rates 2 − α, with the error profiles in remarkable agreement with the target.
The advantages of the proposed approach include: + no need to store past values of the sampled residual (even though the latter affects the local increments in the error non-locally); + applicability to wide classes of time discretizations; + low regularity assumptions; + the approach works seamlessly for arbitrarily large times. Notation. We use the standard inner product ·, · and the norm · in the space L 2 (Ω), as well as the standard spaces L ∞ (Ω), H 1 0 (Ω), L ∞ (0, t; L 2 (Ω)), and W 1,∞ (t , t ; L 2 (Ω)). The notation v + := max{0, v} is used for the positive part of a generic function v.

A posteriori error estimates in the L 2 (Ω) norm
Given a solution approximation u h such that u h = u for t = 0 and on ∂Ω, we shall use its residual R h (·, t) : for t > 0, as well as the operator [10, (2.11)] that (2.1) gives a solution of the equation The main results of the paper are as follows.
. Suppose a unique solution u of (1.1) and its approximation u h are in L ∞ (0, t; L 2 (Ω)) ∩ W 1,∞ ( , t; L 2 (Ω)) for any 0 < < t ≤ T , and also in H 1 0 (Ω) for any t > 0, while u h (·, 0) = u 0 and R h (·, t) = (D α t + L)u h (·, t) − f (·, t). Then The above corollary may seem to imply that one can get any desirable pointwise-intime error profile E(t) on demand. The tricky part is to ensure that (D α t + λ)E(t) > 0 for t > 0, which is not true for a general positive E. Two possible error profiles will be described by the following result.
Remark 2.6 (R 0 v R 1 ). If uniform-in-time accuracy is targeted, then the first bound in (2.3a), with the residual barrier R 0 , is to be employed. The second bound, with R 1 , is less intuitive. It may be viewed as an a posteriori analogue of pointwise-in-time a priori error bounds of type [8, (3.2)] and [6, (4.2)] on graded meshes {T (j/M ) r } M j=0 . Let q denote the order of the method (with q = 2 for the L1 method). The latter bounds show (for three discretizations) that the error behaves as t α−1 M −r for 1 ≤ r ≤ q − α (with a logarithmic factor for r = q − α), while the optimal convergence rate q − α in positive time is attained if r ≈ q − α. Hence, it is reasonable to expect that an adaptive algorithm using residual barriers R 0 and R 1 will respectively yield optimal convergence rates q − α globally or in positive time. This agrees, and remarkably well, with the numerical results in §4 for the L1 method, and in [3] for a number of higher-order methods.
Remark 2.7 (Lu 0 ∈ L 2 (Ω)). If u 0 is not sufficiently smooth (see, e.g., test problem C in §4.2), then (depending on the interpolation of u h in time) the residual R h (·, t) on the first time interval (0, t 1 ) may fail to be in L 2 (Ω). One way to rectify this is to With this modification, all above results become applicable. Importantly, all changes in u h need to be reflected when computing its residual R h ; in particular, as u h has been made discontinuous at t = 0, Corollary 2.5 is to be employed.
The remainder of this section is devoted to the proofs of the above results. The key role will be played by the following auxiliary lemma, a discrete version of which has been useful in the a priori error analysis; see, e.g., in [5, (3.4)].
Remark 2.9. One may consider (2.4) an alternative definition of D α t (with an obvious modification for the case v(·, 0) = 0; see also [1]), which can be applied to less smooth functions, including functions discontinuous at t = 0. Consider E 0 (t) := 1 for t > 0 with 3 E 0 (0) := 0. Then, a calculation using The same result may be obtained using the original definition (1.2) combined with ∂ t E 0 (t) = δ(t), the Dirac delta-function, or representing E 0 as the limit of a sequence of continuous piecewise-linear functions (similarly to [7,Remark 2.4]).
For the first bound in (2.3a), recall from Remark 2.9 that for the function E 0 (t) := 1 for t > 0 with E 0 (0) := 0 one has (D α t +λ)E 0 (t) = R 0 (t). So an application of Corollary 2.3 with E(t) := E 0 (t) yields the first desired bound e(·, t) ≤ E 0 (t) = 1 for t > 0. To evaluate D α t E 1 (t), setτ := τ /t, and note that E Proof of Corollary 2.5. Theorem 2.2 and its two corollaries immediately apply to u h once it is reset to u 0 at t = 0 (after which, it is worth noting, u h becomes right-discontinuous at t = 0). However, this modification of u h needs to be reflected in the computation of the residual R h as follows. Given u h , continuous in t for t ≥ 0, setū h := u h , and then reset u h (·, 0) := u 0 (soū h is continuous for t ≥ 0, while u h is continuous and equal toū h for t > 0). Now, the residual R h of u h for t > 0 is computed using Remark 2.9. In other words, to compensate for u h (·, 0 + ) = u 0 , one needs to add [u h (·, 0 3. Generalization for the L ∞ (Ω) norm , with sufficiently smooth coefficients {a k }, {b k } and c in C(Ω), for which we assume that a k > 0 inΩ, and also c ≥ λ ≥ 0 (while Lv, v ≥ λ v 2 is not required in this section).
Proof. This result is given in [9, Theorem 2] under a stronger condition that v(x, ·) ∈ C 1 (0, T ]∩W 1,1 (0, T ). An inspection of the proof shows that this condition is only required to apply [9, Theorem 1] (the maximum principle for D α t ). The proof of the latter relies on the representation of type (2.4) and remains valid under our weaker assumptions. (A similar, but not identical, result is also given in [1, Theorem 4.1].) Theorem 3.2. Under the above assumptions on L, let a unique solution u of (1.1) and its approximation u h be in C(Ω × [0, t]) ∩ W 1,∞ ( , t; L ∞ (Ω)) for any 0 < < t ≤ T , and also in C 2 (Ω) for any t > 0. Then the error bounds of Theorem 2.2 and Corollaries 2.3 and 2.4 remain true with · = · L2(Ω) replaced by · L∞(Ω) .

4.1.
Numerical results for a test without spatial derivatives Test problem A. We start our numerical experiments with a version of (1.1) without spatial derivatives, with L := 3 and the exact solution u = u(t) = t α − t 2 (which exhibits a typical singularity at t = 0) for t ∈ (0, 1]. For this problem, a straightforward adaptive algorithm (see §4.3) was employed, motivated by (2.3), and so constructing a temporal mesh such that R h (·, t) ≤ TOL · R p (t), p = 0, 1, with τ := t 1 in R 1 . The errors and rates of convergence obtained using residual barriers R 0 (t) and R 1 (t) are presented in Fig. 1 & 2. For R 0 , the errors on the adaptive meshes were compared with the errors on the optimal graded meshes {t j = T (j/M ) r } M j=0 with r = (2 − α)/α [5,8,12] for the same values of M . We observe that in both cases the optimal global rates of convergence 2−α are attained. Furthermore, not only the adaptive meshes successfully detect the solution singularity, but they slightly outperform the optimal graded meshes. For R 1 , we observe the optimal rates of convergence 2 − α at terminal time t = 1, which is consistent with the error bound [8, (3.2)] for a mildly graded mesh (see Remark 2.6).

Numerical results for fractional parabolic test problems
Test problem B. Next, we consider (1.1) for (x, t) ∈ (0, π) × (0, 1] with L = −∂ 2 x and the exact solution u := (t α − t 2 ) sin(x 2 /π), so we set λ := 1. The same adaptive algorithm was employed with R 0 (t) from (2.3) to generate temporal meshes, while in space the problem was discretized on the uniform mesh with 10 4 intervals using standard finite differences (equivalent to lumped-mass linear finite elements). The numerical results are given on Fig. 3 (left, centre) are similar to those on Fig. 1 for test problem A.

The residual becomes
for t > t 1 , where σ(t) := t −α −(1−α) −1 [t 1−α −(t−t 1 ) 1−α ]/t 1 . A fixed mesh with 10 5 subintervals was used in space. The reference solution was computed on a finer mesh. The numerical results, given on Fig. 3 (right) indicate that our adaptive algorithm provides adequate error control for piecewise-linear initial data, as well as for more typical solution singularities at initial time. For a further numerical study of this approach, we refer the reader to [3].