A posteriori error bounds for fully-discrete hp-discontinuous Galerkin timestepping methods for parabolic problems

We consider fully discrete time-space approximations of abstract linear parabolic partial differential equations (PDEs) consisting of an hp-version discontinuous Galerkin (DG) time stepping scheme in conjunction with standard (conforming) Galerkin discretizations in space. We derive abstract computable a posteriori error bounds resulting, for instance, in concrete bounds in - and -type norms when I is the temporal and the spatial domain for the PDE. We base our methodology for the analysis on a novel space-time reconstruction approach. Our approach is flexible as it works for any type of elliptic error estimator and leaves their choice to the user. It also exhibits mesh-change estimators in a clear and concise way. We also show how our approach allows the derivation of such bounds in the norm.


Introduction
Adaptive numerical methods have been shown to provide accurate and efficient numerical treatment of evolution PDEs thanks to their properties for localized mesh resolution especially in the context of moving fronts, interfaces, singularities, or layers (both boundary and interior).Such numerical methods predominantly admit spatial discretizations of variational type, e.g., finite element methods (FEMs), which allow for general, possibly unstructured, dynamic mesh modification.FEMs are also ideally suited for deriving mathematically rigorous a posteriori bounds, owing to their variational nature.
Some of the classical works on adaptive finite element methods for parabolic problems [EJ91,EJ95a,EJ95b,EJ95c,EJL98] are based on discontinuous Galerkin (DG) time stepping combined with FEM in space, and proving a posteriori bounds in various norms using duality techniques.The key motivation in using DG in time, which is also of variational type, is that it naturally allows for spatiallylocal-time stepping, i.e. different time step sizes in different parts of the spatial domain [Jam78,EJT85,EJ91,MB97].This classical, but as of yet undeveloped in full, concept of local adaptivity in both space and time has the potential of delivering substantial computational savings and even complexity reduction.
In addition to the ability of Galerkin time marching schemes to employ locally different time step sizes, their variational character also allows for arbitrary variations in the local approximation orders.They can therefore be cast naturally into the framework of hp-approximation schemes.In the context of parabolic PDEs, hp-version time marching methods can be used, for instance, to resolve an initial layer in the (otherwise smooth) solution at high algebraic or even exponential rates of convergence, see, e.g., the works [SS00,SS01,WGSS01] on linear parabolic PDEs, and also [MSW05,MSW06] which employ a combination of hp-version time stepping with suitable wavelet spatial discretizations to yield a log-linear complexity algorithm for nonlocal evolution processes involving pseudo-differential operators.Additionally, we note the numerical analysis of high-dimensional parabolic problems using sparse grids in space; see [vPS04].
More recent results on rigorous a posteriori bounds for parabolic problems have focused on extending the paradigm of the reliable and efficient a posteriori error analysis of elliptic problems to the parabolic case [Pic98,Ver98,Ver03].Such works typically involve basic low-order time stepping schemes combined with various types of FEM in space.A posteriori error bounds for DG time-stepping methods have also appeared in the last few years; we point to [MN06a, SW10,KMW18] which are based on the reconstruction technique, to [ESV17,ESV19] which employ an equilibrated flux approach, or to [GSKZ19] which presents a provably convergent adaptive algorithm for a residual-type a posteriori estimator.
In this paper, we present a posteriori error bounds for an hp-version DG-intime and conforming Galerkin discretization in space method for both L ∞ (I; L 2 )and L 2 (I; H 1 )-norm errors separately, allowing for what appears to be optimal order in each case.The key idea is the use of suitable reconstruction frameworks to derive a perturbed PDE for the reconstructed error of the numerical method; a posteriori error bounds are then deduced using PDE stability properties, cf.[MN03,MN06a,LM06a,AMN06].Our approach is based on new space and spacetime reconstructions which are built on the combination of respective ideas for DGtime stepping methods [MN06a, SW10] and elliptic reconstruction [MN03,LM06a] to the fully-discrete setting.To that end, the key challenge of constructing a globally time-continuous reconstruction in the presence of mesh modification between time-steps is addressed by first reconstructing onto the solution space with respect to the spatial variables via a novel elliptic reconstruction definition, given in (3.14), which is a modification of the one in [LM06a].In particular, the new proposed elliptic reconstruction takes into account the effect of mesh-change.
Our results are closely related, however, with important departures, to those of [ESV17,ESV19,GSKZ19].In particular, the new reconstructions defined below allow for the derivation of a posteriori upper error bounds for each of the following norms L 2 (I; X ) (Theorem 4.4) and H 1 (I; X ′ ) ( §6) separately; the Hilbert space X is the domain of a self-adjoint uniformly elliptic operator A (see §2 for details).
A key attribute of our approach is the flexibility in incorporating any a posteriori elliptic error estimators available, as the reconstruction-type approach, allows to separate the challenges in the a posteriori error estimation of elliptic and timeevolution errors.To facilitate a wide range of applications, we will present the theory within a Gelfand-triple-type abstract setting allowing, for instance, both secondand fourth-order spatial operators.This generality comes at the possible expense of different, yet quantitatively analogous, computable constants in the resulting a posteriori error estimators compared to the bounds in [ESV17,ESV19,GSKZ19].For instance, when the equilibrated flux elliptic error estimators from [BPS09] are used (with X = H 1 0 (Ω) and H = L 2 (Ω)), we recover similar estimators to the upper bounds derived in [ESV17,ESV19].Importantly, however, the work [ESV17] shows that these are also lower bounds for the "joint-norm" of H 1 (I; X ′ ) ∩ L 2 (I; X ), and the article [ESV19] does the same for the L 2 (I; X ) under the condition h 2 < cτ , relating the mesh-size h with the time-step τ for some constant c > 0. Also, in the present work, we are not concerned with the interesting question of convergence of adaptive algorithms as in [GSKZ19].Crucially, however, the novel space-time reconstruction, allows for the proof of an a posteriori error bound for the L ∞ (I; H )norm, which appears to be of optimal order; this result, to the best of our reading, is not captured in [ESV17,ESV19,GSKZ19].
Outline.The remainder of this work is structured as follows.In §2 we set up the abstract framework for the paper by introducing the model parabolic PDE problem and its DG-in-time and conforming Galerkin spatial discretization.Furthermore, in §3, we provide the necessary technical tools for the ensuing analysis, and state their essential properties.In §4, we derive a posteriori error bounds in the L 2 (I; X )norm using a time reconstruction and a novel elliptic reconstruction which includes mesh-change; this technical novelty is revealed to be crucial in our setting.In §5, upon defining a new space-time reconstruction, we derive a posteriori error bounds in the L ∞ (I; H )-norm. Finally, in §6 we briefly discuss an approach to arriving at H 1 (I; X ′ )-type a posteriori error estimates.

Model problem and space-time discretization
We introduce most of the notation and technical background for the paper.In §2.1 we provide the functional analytic set-up for the abstract heat equation, a related concrete Example 2.2, and we present the fully discrete numerical scheme in §2.3.
2.1.Abstract setting.Throughout this work, Bochner spaces will be used.To that end, given an interval J ⊂ R, and a real Hilbert space Z with inner product (•, •) Z and induced norm • Z , we define (2.1) We write L p (J; Z ) to signify the space of measurable functions u : J → Z such that the corresponding norm is bounded.Note that L 2 (J; Z ) is a Hilbert space with inner product and induced norm given by and u L2(J;Z ) := (u, u) 1 /2 L2(J;Z ) , respectively.We also let H 1 (J; Z ) be the Sobolev space of all functions in L 2 (J; Z ) whose (weak temporal) derivative is bounded in L 2 (J; Z ), with the norm Finally, the space C 0 (J; V ) consists of all functions that are continuous on J, the closure of J, with values in Z , endowed with the standard maximum norm We consider henceforth two (real) Hilbert spaces X and H forming a Gelfand triple (2.5) where X ′ denotes the dual of X .The duality pairing • | • of X ′ and X can be seen as a continuous extension of the inner product (•, •) H .In particular, identifying H ′ ≃ H , for u ∈ H and v ∈ X , there holds see, e.g., [Rou13, §7.2] for details.Moreover, let (2.7) A : X → X ′ be a self-adjoint linear elliptic operator continuous and coercive in the sense that there exist constants β ≥ α > 0 such that (2.8) Given an initial value u 0 ∈ H , a final time T > 0, denoting henceforth the time interval (2.9) and given a source function f ∈ L 2 (I; X ′ ), we are interested in a Galerkin-type numerical approximation of the function which solves uniquely the linear parabolic initial value problem Incidentally, due to the continuous embedding it follows that u belongs to C 0 (I; H ) [Rou13, e.g., Lemma 7.3] and the initial condition in (2.11) makes sense.

Example (concrete elliptic operators).
A commonly encountered situation which can be cast in the above framework is the classical linear diffusion equation, i.e.

Time discontinuous and space conforming Galerkin approximation.
Given a (real) linear space Z , the space of all Z -valued polynomials of degree at most r, with r ∈ N 0 , on R is defined by In order to introduce the discontinuous Galerkin time stepping scheme for (2.11), we consider a finite sequence of time nodes and time steps, (2.16) 0 = t 0 < t 1 < . . .< t N = T, and τ n := t n − t n−1 for n = 1, . . ., N, as well as the corresponding time intervals (2.17) Thus, we have a partition I := {I n : n = 1, . . ., N } of the time interval I given in (2.9).
Given a I -piecewise continuous function g : I ⊆ R → Z , we define its time jump across t n , n = 0, . . ., N − 1, for given g(t − 0 ), by (2.18) where we introduce the one-sided limits g(t ± n ) = lim ǫ→0 + g(t n ± ǫ).Moreover, we associate with the finite sequence of time instants t 0 , . . ., t N a finite sequence of finite-dimensional conforming subspaces For a generic X -conforming Galerkin space X, we signify by π X the H -orthogonal projection from X ′ onto X: Note that, due to (2.6), for v ∈ H , we have When X is X n , for some n = 0, . . ., N , we write π n to indicate π X .In order to introduce the time semidiscrete and space-time fully discrete spaces, let r n ∈ N 0 , n = 1, . . ., N , be a polynomial degree.Then, consider the time semidiscrete Galerkin space

Reconstructions
We will next introduce some technical essentials.The main tools are the time lifting ( §3.1), the time reconstruction ( §3.3), and a new variant of the elliptic reconstruction from [LM06a] for fully discrete schemes ( §3.6).In §3.7 we postulate the availability of a posteriori error estimators for elliptic residuals, and we give some pointers to the relevant literature.In addition, we discuss various error estimates that measure the time reconstruction error; in particular, we state two identities which follow directly, respectively, from [SW10, Theorem 2] In particular, writing 1 A for the indicator function on a generic set A, and assuming w n ∈ X n for each n = 1, . . ., N , we have Proof.This result is a straightforward consequence of the explicit representation of χ n as described in [SW10, Lemma 6].
Equivalently, we note the following characterization of W in weak form on each I n , for each V ∈ P rn (I n ; H ), with the initial condition for n = 1, . . ., N ; cf.[MN06a,SW10].Evidently, the above construction carries over to any linear subspace W ⊂ H in an obvious way.

Proposition (time-reconstruction error identities). Consider any (real)
Hilbert space W , and W | In ∈ P rn (I n ; W ), n = 1, . . ., N , with W defined from W through (3.6).Then, for given W (t − 0 ) ∈ W , the following approximation identities hold where , and Proof.The identity (3.9) was first proven in [MN06a, Lemma 2.2], and extended to this exact form in [SW10, Theorem 2] accounting for the dependence on the polynomial degree explicitly.The second equality (3.11) follows directly by combining the explicit representation formula derived in [SW10, Eq. ( 33)] with [HW18, Lemma 1].
3.5.Remark (continuity of the time-reconstruction). Owing to [MN06b, Lemma 2.1], the semidiscrete (spatially exact) time-reconstruction (3.6) originally defined in [MN06b] and [SW10] is a continuous function in time.In particular, the time-reconstruction U of the fully discrete solution U, defined in (2.24) and (2.25), is still continuous across the time nodes t 0 , . . ., t N −1 , despite having π n U = U, on I n , when the spatial mesh changes across t n−1 in a non-hierarchical fashion.
3.6.Elliptic reconstruction.Let X ⊂ X be a generic conforming Galerkin space.Then, given the elliptic operator A from (2.7), we define the discrete elliptic operator A X : X → X, for each w ∈ X, as A X w ∈ X such that From the ellipticity of A , it follows that A X : X → X is invertible.Note that the discrete elliptic operator's domain may be extended from X to all of X ; indeed, this may be convenient in some cases where we are ready to give up its invertibility.If X is one of X n s, for some n = 0, . . ., N , we denote A Xn by A n .
To optimize on the structure of the mesh-change indicator below, we use a nonstandard elliptic reconstruction on each time interval I n .To that end, for each t ∈ I n , we define the elliptic reconstruction U (t) ∈ X , by where U is the time-reconstruction of the fully discrete solution U from (2.24) and (2.25).It follows that U ∈ Y and may be written implicitly, for any n = 1, . . ., N , as the solution of the t-dependent elliptic problem The initial value of U is given by Upon restricting the test functions in (3.13) to X n , we obtain cf. (3.12).Evidently this identity implies (3.17) for each V ∈ P rn (I n ; X n ), n = 1, . . ., N .
3.7.Assumption (elliptic a posteriori error estimates).For given g ∈ H , consider the abstract elliptic problem of finding w ∈ X such that A w = g.Moreover, let X ⊆ X be a generic X -conforming Galerkin space, and let w ∈ X be w's Galerkin approximation in X, defined implicitly as the solution of A X w = π X g.Then some a posteriori error bound holds, viz., with a suitable a posteriori error estimator E Z ,X , which we assume to be available for Z representing any of the spaces X , H or X ′ .Recalling (3.14), assumption (3.18) allows, for instance, to get a posteriori error control of the elliptic reconstruction error in the Z -norm, i.e. (3.19) for the selection of spaces Z = X , H , or X ′ .Details on such a posteriori error estimates can be found, e.g., in [AO00, Bra07, BS07, DM99].It is worth mentioning that w does not need to belong to X for (3.18) to hold; instead, it is usually enough that (w − w) is A -orthogonal to X in order to derive elliptic a posteriori error estimates.
3.8.Pointwise form.For n = 1, . . ., N , we denote by the time-local fully discrete L 2 (I n ; H )-orthogonal projection defined by (3.21) and its time-global counterpart Let V ∈ Y, and note that using Definition 3.3 and identity (3.17), the fully discrete local DG formulation (2.25) transforms into (3.23) Thus, for each n = 1, . . ., N , we have noting that Π n U ′ = π n U ′ from Fubini's Theorem.Hence, from the elliptic reconstruction (3.13), we deduce the pointwise form 3.9.Remark.Upon noting the trivial identity for t ∈ I n , we observe that the above definition of the elliptic reconstruction constitutes a high order extension of the zero-th order version from [CGM14, Definition 6.1].

L 2 (I; X )-norm a posteriori error analysis
To highlight the potential generality of the reconstruction approach, we first embark on proving a posteriori error bounds in the L 2 (I; X )-norm.As we will see below, the resulting bounds are qualitatively closely related to the bounds from [ESV17, ESV19, GSKZ19] when X = H 1 (Ω).A key attribute of the approach presented below is that the resulting a posteriori bounds are flexible with respect to the choice of respective elliptic bounds, such as residual, recovery, or flux-reconstruction ones.Recalling that our ultimate goal is to find a posteriori estimates for e, and noting that such estimates for σ and σ are provided by Proposition 3.4 and Assumption 3.7, respectively, it will suffice to estimate e, or e, in terms of σ or σ.

Lemma (time reconstruction error bounds).
For n = 1, . . ., N , we have the explicitly computable bounds Proof.The proof follows immediately from Proposition 3.4.

Lemma (space reconstruction error bounds).
Assume the availability of elliptic a posteriori error bounds as per Assumption 3.7.Then, for n = 1, . . ., N , we have the explicitly computable bound We can now prove the first main result of this work.
4.5.Remark.For the particular case A = −∆, X = H 1 0 (Ω) and H = L 2 (Ω), cf.Example 2.2, note that the constants in (2.8) can be chosen to be α = β = 1.Especially, selecting ϑ 2 = 1 /2, we have c α,β ϑ = 1 in (4.8).4.6.Remark (mesh change error via elliptic reconstruction).The elliptic reconstruction (3.13) features a mesh-change type term.This is an important departure from the (standard) elliptic reconstruction proposed in [LM06b]; cf. also [CGM14, Definition 6.1].For instance, if E X ,Xn from (3.18) is the standard residual energy norm a posteriori error estimator, the element residual includes the expression π n U ′ − U ′ , which is effectively a mesh-change term.Indeed, from (3.6) we have on I n : from the commutativity of the spatial projection π n and operations on the time variable.Let now κ n : I n → R be the polynomial of degree r n representing the time lifting χ n from (3.2), i.e. such that for all t ∈ I n and W ∈ Y ; we refer to [SW10, Lemmas 6 & 7] or [HW18, Remark 1] for an explicit formula for κ n .Then, we can conclude i.e. the arbitrary order DG-analogue to the classical mesh-change indicator.We note that the representation (4.17) is particularly relevant in implementation as it can be used to efficiently realize the space reconstruction error bounds in Lemma 4.3.Finally, observing that κ n takes its maximum at t n−1 , which follows from [SW10, Lemma 6], and recalling (3.2), it is possible to compute the maximum value of κ n on I n : 4.7.Remark (reconstruction vs direct approach).As noted by one referee, it is possible to spare the use of elliptic reconstruction in the proof of the L 2 (I; X )norm a posteriori bound, for particular cases of operators A and of Gelfand triples (2.5), as discussed in [ESV17,ESV19,GSKZ19] for the specific case A = −∆, X = H 1 0 (Ω) and H = L 2 (Ω); cf.Example 2.2.As a result, different a posteriori error estimators arise with, possibly, slightly different constants multiplying common terms in the estimator.Elliptic reconstruction, nonetheless, offers the ability to use various elliptic estimators from the literature in the bound: this feature may become important for multiscale operators A , e.g., singular perturbations.In general, elliptic reconstruction allows for, crucial in some cases, flexibility in the handling for more complicated spatial operators A , e.g., nonlinear or singularly perturbed operators [CGM14, GM14, CGS20a].

L ∞ (I; H )-norm a posteriori error analysis
We continue by proving an a posteriori error bound in the L ∞ (I; H )-norm, which appears to be of higher order than the L 2 (I; X )-norm one presented above.
A key difference with respect to the above proof of the L 2 (I; X )-norm bound is the use of a combined space-time reconstruction defined below.5.1.Error-residual relation.Set w := U , i.e. the time reconstruction of the elliptic reconstruction, noting that Definition 3.3 is still valid for functions on X in space.Now, introducing the time error which is continuous on I = [0, T ], and subtracting the PDE (2.11) from (3.25), we have: our aim being to deduce an evolution equation for ρ.To this end, we have from the linearity of the time reconstruction.We define T n ∧ T n−1 to be the finest common coarsening of the two meshes, with the associated largest common finite element subspace given by X 5.2.Lemma (time reconstruction error bound).Let t ∈ I n , n = 0, 1 . . ., N .Then, we have the abstract estimate (5.4) A (w − U ) L2(In;X ′ ) ≤ C τn,rn η time n,1 , with Assume further that, for every v ∈ X , there exists a v ∈ X ⊖ n , such that for some s > 0, and for generic constants C ap , C stab > 0, independent of v and of h ⊖ n .Then, we also have the bound (5.7) Therefore, setting η time n := min{η time n,1 , η time n,2 }, we have the combined estimate (5.9)A (w − U ) L2(In;X ′ ) ≤ C τn,rn η time n .
Proof.Recalling (3.6) and using the fact that the elliptic operator A is time independent, we immediately observe that (5.10) Therefore, by Proposition 3.4 and of Remark 3.9, for n = 1, . . ., N − 1, we conclude that , which is (5.4).Furthermore, exploit the orthogonality identity (5.12) which follows from the definition of U and the properties of Π. From the latter, together with Remark 3.9, and the temporal independence of A , we have (5.13) From (5.6) and the continuity of A , therefore, we deduce (5.14) ; this is particularly pertinent when the dimension of X ⊖ n is comparable to those of X n and X n−1 .In cases, however, in which (5.6) may not be effective, or even known, we may be forced to revert to (5.4); see [CGS20b] for two such settings, where one involves virtual element methods, and the other incorporates moving mesh methods.The dual norm in η time n,1 typically requires a global inversion to be evaluated, in the absence of a natural Galerkin orthogonality property like (5.12) above.
Although solving an elliptic problem to determine/estimate the time estimator η time n,1 may appear to be computationally demanding, such global solves may be typically performed in a lower dimensional space than X n itself, through an a posteriori error controlled fashion, as we now demonstrate.Let ϕ ∈ X be the solution to the problem A ϕ = w (with w := Πf − U ′ n−1 for instance).Let also an approximation Φ ∈ Xn to ϕ, for some Galerkin space Xn , be given by (5.15) Now, the continuity of A implies w X ′ ≤ β ϕ X .In addition, from coercivity, continuity and the Galerkin orthogonality A (ϕ−Φ) | Φ = 0, we have, respectively, (5.16) X , using also the self-adjointness of A .The above, together with the assumed availability of an a posteriori error indicator in the energy norm give the computable bound (5.17) or, for the specific case w = Πf − U ′ n−1 , (5.18) 5.4.Lemma (space reconstruction error bound).For t ∈ I n , n = 1, . . ., N , we have the bound (5.20) with κ n : I n → R from (4.16) in Remark 4.6, and the estimator is explicitly computable, via Assumption 3.7.
It remains to bound the last term on the right-hand side of (5.33).To that end, applying (3.9), and involving Lemmata 5.2 and 5.4, results in (5.34) The triangle inequality now gives u − U L∞(0,tn;H ) ≤ ρ L∞(0,tn;H ) + w − U L∞(0,tn;H ) + U − U L∞(0,tn;H ) .(5.35) To estimate the second and third term on the right-hand side of (5.35), we apply (3.11), and the triangle inequality, to obtain Also, with the aid of (3.19), (5.37) Combining the last two estimates, we conclude (5.41)Now, using Hölder's inequality, recalling (5.10), and applying Proposition 3.4 and Remark 3.9, we have, respectively, (5.42) (5.45) this, in conjunction with (5.38) now yields an alternative a posteriori error bound which may be superior to the one given in Theorem 5.5 for small final time t n .
The second term on the right-hand side of the last estimate is a computable quantity, while the first can be immediately bounded using Lemma 5.4, through the linearity of the time reconstruction.

Conclusions
In this article we have presented a posteriori error estimates for fully discrete hp-DG-in-time and general Galerkin-in-space discretizations of abstract linear parabolic PDE.Our approach is based on a novel combination of temporal and elliptic reconstructions, the latter allowing for arbitrary elliptic error spatial estimators.Our main results include computable L ∞ (H )-and L 2 (X )-a posteriori error estimates; some remarks concerning H 1 (X ′ )-norm error estimation are given as well.Finally, we note that a series of numerical experiments showcasing the optimality of the a posteriori error estimators derived above are given in [CGS20a].
and the result already follows.5.3.Remark (Practical upper bound for η timen,1 ).In standard settings, e.g., in the canonical example of the heat equation, cf.Example 2.2, for which (5.6) holds, we can use η time n,2 for η time n