Christine Bernardi
email: bernardi@ann.jussieu.fr

Frédéric Hecht
email: hecht@ann.jussieu.fr

Rüdiger Verfürth
email: ruediger.verfuerth@ruhr-uni-bochum.de

A Posteriori Error Analysis of the Method of Characteristics

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In this note we consider general linear parabolic differential equations of second order, equation (2.1), and their discretization by finite elements based on the method of characteristics, equations (3.1) - (3.3). Our focus is on deriving reliable, efficient and robust a posteriori error estimates for these discretizations. Here, as usual, robust means that the quality of the error indicators, i.e. the ratio of the upper and lower bounds, does not depend on the relative size of the diffusion to the convection or reaction.

The method of characteristics was initiated in [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF] as an alternative to standard space-time finite elements with or without upwinding or SUPG-stabilization. Its main advantages are that it implicitly includes an upwinding and thus stabilizing effect and that it requires only the solution of symmetric coercive problems. On the other hand it needs an efficient way for back-tracking the characteristics and re-interpolating finite element functions. Depending on the data structure of the implementation, this can be costly or not. The initial method was of order one but was soon extended to higher order in [START_REF] Boukir | A high order characteristics method for the incompressible Navier-Stokes equations[END_REF]. We refer, e.g., to [START_REF] Bermùdez | Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. Part I. Time discretization[END_REF][START_REF]Numerical analysis of convection-diffusion-reaction problems with higher order characteristics/finite elements. Part II. Fully discretized schemes and quadrature formulas[END_REF][START_REF] Boukir | A high order characteristics method for the incompressible Navier-Stokes equations[END_REF][START_REF]Finite element characteristics methods requiring no quadrature[END_REF][START_REF] Rui | A second order characteristic finite element scheme for convection-diffusion problems[END_REF][START_REF] Tourigny | The finite element method with nodes moving along the characteristics for convection-diffusion equations[END_REF] for a more detailed description of the method of characteristics including a priori error estimates.

On the other hand, a large amount of work has been done concerning the a posteriori analysis of parabolic equations aimed at time and space adaptivity, see e.g. [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF][START_REF]A posteriori error analysis of the fully discretized time-dependent Stokes equations[END_REF][START_REF] Eriksson | Adaptive finite element methods for parabolic problems. I. A linear model problem[END_REF][START_REF]Adaptive finite element methods for parabolic problems. IV. Nonlinear problems[END_REF][START_REF] Verfürth | A posteriori error estimates for nonlinear problems: L r (0, T ; W 1,ρ (Ω))-error estimates for finite element discretizations of parabolic equations[END_REF][START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF][START_REF]A posteriori error estimates for linear parabolic equations[END_REF][START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF] and the references therein. However, we are not aware of any a posteriori error analysis for the method of characteristics although, as we will see, it can be performed by adapting of the approach of [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF][START_REF]A posteriori error estimates for linear parabolic equations[END_REF][START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF].

The subsequent analysis follows the general framework of a posteriori error analysis and is performed in several steps. First, in Lemma 4.1, we prove that the energy norm of the error and the norm of its material derivative are equivalent to the dual norm of the residual (see section 2 for a detailed definition of these norms). This stability result is based on standard parabolic energy estimates and is independent of the discretization.

Next, in section 5, we decompose the residual into three contributions: a term representing temporal data oscillations, a temporal residual and a spatial residual. The first term can easily be estimated separately. For the estimation of two other terms it is crucial that the norm of their sum is equivalent to the sum of their norms, cf. Lemma 5.1. This result essentially is a strengthened Cauchy-Schwarz inequality. However, it is proved in a different way which is better adapted to dual norms.

The temporal residual is estimated in section 6. Here, the crucial step is Lemma 6.1 which allows to split the temporal residual into two contributions reflecting the effects of the convection and of the combination of the diffusion and reaction. In case of a dominant diffusion or reaction, the contribution of the convection can be controlled by the diffusion-reaction term. In case of a dominant convection, however, this leads to a non robust estimate. To recover robustness, as in [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF], the convection term must be controlled with the help of an auxiliary discrete diffusion-reaction problem.

The spatial residual is the term which reflects the difference between the method of characteristics and standard space-time finite element discretizations. It is a linear functional which admits an L 2representation involving contributions by the elements and faces of the partition. Thus it falls into the abstract class of functionals which are amenable to a posteriori error estimation [START_REF] Veeser | Explicit upper bounds for dual norms of residuals[END_REF]27]. However, it does not satisfy the Galerkin orthogonality, i.e. its kernel does not contain the lowest order conforming finite element space. Instead, its consistency error describes the difference between the method of characteristics and the backward Euler scheme. Contrary to the consistency errors of SUPG-or discontinuous Galerkin schemes, this term is crucial and cannot be controlled by the standard element and face residuals. Instead, in Lemma 7.1, we prove that this term can be controlled with the help of an auxiliary discrete diffusion-reaction problem. Thus its treatment is similar to the one of the convection in the convection-dominated regime.

In section 8 we collect all estimates and thus obtain our main result, Theorem 8.1. When compared with [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF] the error indicator of Theorem 8.1 requires, in each time step, the solution of an additional discrete diffusion-reaction problem corresponding to the lowest order conforming finite element space. This extra cost is compensated by the less costly discretization. We stress that all estimates are robust with respect to the relative size of the diffusion compared to the reaction or convection. Moreover, they hold for every conforming spatial discretization based on Lagrangian elements and every ODE-solver for the characteristic equations as long as the endpoints of the approximate characteristics remain within the closure of the domain.

In section 9 we present some numerical experiments both for convection and diffusion-dominated problems in a square. The numerical results confirm the robustness of our estimates.

The differential equation

In what follows we consider general linear parabolic equations of second order:

∂ t u -div(D∇u) + a • ∇u + bu = f in Ω × (0, T ] u = 0 on Γ × (0, T ] u| t=0 = u 0 in Ω. (2.1)
Here, Ω ⊂ R d is a bounded polygonal two or three dimensional domain with a Lipschitz-continuous boundary Γ. The final time T is arbitrary, but kept fixed in the sequel. We assume that the coefficients are smooth enough:

D ∈ C(0, T ; L ∞ (Ω) d×d ), a ∈ C(0, T ; W 1,∞ (Ω) d ), b ∈ C(0, T ; L ∞ (Ω)),
and satisfy the following conditions:

• The diffusion D is symmetric, uniformly positive definite and uniformly isotropic, i.e.

ε = inf 0<t≤T,x∈Ω min z∈R d \{0} z T D(x, t)z z T z > 0
and the quantity

κ = ε -1 sup 0<t≤T,x∈Ω max z∈R d \{0} z T D(x, t)z z T z is of moderate size.
• The convection a is solenoidal and vanishes on the boundary:

div a = 0 in Ω × (0, T ], a = 0 on Γ × (0, T ].
• The reaction b is non-negative and does not vary too much, i.e.

there are two constants β ≥ 0 and c b ≥ 0 such that

β ≤ b ≤ c b β in Ω × (0, T ].
Moreover, the constant c b should be of moderate size. These assumptions guarantee that, for any data f in L 2 (Ω × (0, T )) and u 0 in L 2 (Ω), problem (2.1) is a well-posed parabolic problem, see e.g. [3, Sections II.4 and V.2] and [10, Theorems XVIII.3.1 and XVIII. 3.2]. The parameter κ is introduced in order to stress that the ratio of the constants in the error estimates depends on the condition number of the diffusion matrix. If this condition number is exceedingly large, length scales such as element diameters must be measured in a diffusion-dependent metric in order to recover robust estimates as described in [START_REF] Bernardi | Adaptive finite element methods for elliptic equations with non-smooth coefficients[END_REF] for elliptic equations. The assumptions cover a wide range of different regimes:

• dominant diffusion: a L ∞ (0,T ;W 1,∞ (Ω)) ≤ c c ε and β ≤ c ′ b ε with constants of moderate size; • dominant reaction: a L ∞ (0,T ;W 1,∞ (Ω)) ≤ c c ε and β ≫ ε with a constant c c of moderate size; • dominant convection: a L ∞ (0,T ;W 1,∞ (Ω)) ≫ ε.
When compared with [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF] we here consider a more general problem in that D = εI is assumed in [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF]. On the other hand, the assumptions of [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF] are slightly more general concerning the convection and the boundary conditions in that there the convection may not be solenoidal, may not vanish on the boundary and that Neumann boundary conditions are allowed. The present assumptions are mainly made in order to avoid unnecessary technical difficulties with the method of characteristics.

We conclude this section with some notation and technical results that will be used in the sequel. For any bounded open subset ω of Ω with Lipschitz boundary γ, we denote by H k (ω), k ∈ N, L 2 (ω) = H 0 (ω), and L 2 (γ) the usual Sobolev and Lebesgue spaces equipped with the standard norms [START_REF] Adams | Sobolev Spaces[END_REF]). Similarly, (•, •) ω and (•, •) γ denote the scalar products of L 2 (ω) and L 2 (γ), respectively. If ω = Ω, we omit the subscript Ω.

• k;ω = • H k (ω) and • γ = • L 2 (γ) (cf.
As usual, H 1 0 (Ω) is the Sobolev space of L 2 -functions having their derivatives in L 2 (Ω) and vanishing on the boundary Γ. It is equipped with the norm

|v | = ε ∇v 2 + β v 2 1 2 .
The dual space of H 1 0 (Ω) is denoted by H -1 (Ω) and is equipped with the norm

|ϕ | * = sup v∈H 1 0 (Ω)\{0} ϕ, v |v |
where •, • stands for the standard duality pairing.

For any separable Banach space V with norm • V and any two numbers a < b, we denote by L 2 (a, b; V ) and L ∞ (a, b; V ) the spaces of measurable functions u defined on (a, b) with values in V such that the mapping t → u(•, t) V is square integrable, respectively, essentially bounded. These are Banach spaces equipped with the norms

u L 2 (a,b;V ) = b a u(•, t) 2 V dt 1 2 , u L ∞ (a,b;V ) = ess. sup t∈(a,b) u(•, t) V (cf.
[10, Chap. XVIII, §1]). For abbreviation we introduce the space

X(a, b) = u ∈ L 2 (a, b; H 1 0 (Ω)) ∩ L ∞ (a, b; L 2 (Ω)) : ∂ t u + a • ∇u ∈ L 2 (a, b; H -1 (Ω))
and equip it with its graph norm

u X(a,b) = ess. sup t∈(a,b) u(•, t) 2 + b a |u(•, t) | 2 dt + b a |(∂ t u + a • ∇u)(•, t) | 2 * dt 1 2 .
Here the derivative ∂ t u has to be understood in the distributional sense [10, loc.cit.]. Note that ∂ t u+a•∇u is the material derivative associated with the differential operator of (2.1). Furthermore, according to [10, loc.cit.], every function u ∈ X(a, b) has point-values u(•, a) and u(•, b) in H -1 (Ω). With this notation, the variational formulation of problem (2.1) takes the form:

Find u ∈ X(0, T ) such that u(•, 0) = u 0 in H -1 (Ω) and for almost every t ∈ (0, T ) and all v ∈ H 1 0 (Ω) (∂ t u, v) + (D∇u, ∇v) + (a • ∇u, v) + (bu, v) = (f, v).
(2.2)

These assumptions imply that, for any data f in L 2 (0, T ; H -1 (Ω)) and The key idea of the method of characteristics relies on another formulation of problem (2.1). Since the function a is Lipschitz continuous with respect to the spatial variable and vanishes on the boundary Γ, for every (x * , t * ) ∈ Ω × (0, T ], standard global existence results for the flows of ordinary differential equations (cf., e.g., [START_REF] Amann | Ordinary differential equations[END_REF]Theorem II.7.6]) imply that the characteristic equation

u 0 in L 2 (Ω),
d dt x(t; x * , t * ) = a(x(t; x * , t * ), t), t ∈ (0, t * ), x(t * ; x * , t * ) = x * (2.3)
has a unique solution x(•; x * , t * ) which exists for all t ∈ [0, t * ] and stays within Ω ∪ Γ. Hence, we may set U (x * , t) = u(x(t; x * , t * ), t). The total derivative d t U satisfies

d t U = ∂ t u + a • ∇u.
Therefore, the first line of problem (2.1) can equivalently be written as

d t U -div(D∇u) + bu = f in Ω × (0, T ). (2.4)
The discretization introduced in the next section relies on a separate treatment of equations (2.3) and (2.4).

For later use we finally note that integration by parts and the previous assumptions imply for all v, w ∈ H 1 0 (Ω) and all t ∈ (0, T ]:

(D∇v, ∇v) + (a • ∇v, v) + (bv, v) ≥ |v | 2 , (2.5) (D∇v, ∇w) + (bv, w) ≤ max{κ, c b } |v | |w |.
(2.6)

The finite element discretization

For the discretization of problem (2.1) we consider partitions I = (t n-1 , t n ] : 1 ≤ n ≤ N I of the time-interval [0, T ] into subintervals satisfying 0 = t 0 < . . . < t N I = T . For every n with 1 ≤ n ≤ N I we denote by I n = (t n-1 , t n ] the n-th subinterval and by τ n = t nt n-1 its length. With every intermediate time t n , 0 ≤ n ≤ N I , we associate a partition T n of Ω and a corresponding finite element space V n ⊂ H 1 0 (Ω) which have to satisfy the following assumptions:

• The closure of Ω is the union of all elements in T n .

• Every element in T n is either a simplex or a parallelepiped, i.e., it is the image of the d-dimensional reference simplex {x ∈ R d :

x 1 ≥ 0, . . . , x d ≥ 0, x 1 + . . . + x d ≤ 1} or of the d-dimensional
reference cube [0, 1] d under an affine mapping (affine equivalence). • Any two elements in T n are either disjoint or share a complete lower dimensional face of their boundaries (admissibility). • For every element K, the ratio of its diameter h K to the diameter ρ K of the largest ball inscribed into K is bounded uniformly with respect to all partitions T n and I (shape regularity). • For every n with 1 ≤ n ≤ N I there is an affine equivalent, admissible, and shape-regular partition T n such that it is a refinement of both T n and T n-1 and such that sup

1≤n≤N I sup K∈ e Tn sup K ′ ∈Tn K⊂K ′ h K ′ h K < ∞
uniformly with respect to all partitions I which are obtained by adaptive or uniform refinement of any initial partition of [0, T ] (transition condition). • Each V n consists of continuous functions which are polynomials on each element of T n , the degrees being at least one and being bounded uniformly with respect to all partitions T n and I (degree condition). • For every n there is a set V n of nodes such that every function in V n is uniquely defined by its values in the nodes (Lagrange condition). When compared with [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF] the Lagrange condition is new. It simplifies the formulation of the method of characteristics. The transition condition is needed for the a posteriori error estimates of Lemmas 6.2 and 7.1.

Every set V n splits into the subsets V n,Ω and V n,Γ of all nodes inside Ω respective on Γ. For every intermediate time t n , 1 ≤ n ≤ N I , and every node z ∈ V n,Ω we compute an approximation x n-1 z to x(t n-1 ; z, t n ) by applying an arbitrary but fixed ODE-solver such as e.g. the explicit Euler scheme to the characteristic equation (2.3) with (x * , t * ) = (z, t n ). We assume that the time-step τ n and the ODE-solver are chosen such that x n-1 z lies within Ω ∪ Γ for every n ∈ {1, . . . , N I } and every z ∈ V n,Ω . The assumptions on the convection a imply that this condition is satisfied for a single explicit Euler step if

τ n < 1/ a(•, t n ) W 1,∞ (Ω) [18, Proposition 1]. Note that x n-1 z = z -τ n a(z, t n )
in this case. We stress that our results hold for every ODE-solver which satisfies the previous condition on the x n-1 z .

As standard, all functions in V n will be denoted with an index h. For every n ∈ {0, . . . , N I } we denote by π n the L 2 -projection operator onto V n and set

D n = D(•, t n ), a n = a(•, t n ), b n = b(•, t n ), f n = f (•, t n ).
With this notation the discretization of problem (2.1) is given by:

Set u 0 h = π 0 u 0 . (3.1) For n = 1, . . . , N I successively compute u n-1 h ∈ V n such that u n-1 h (z) = u n-1 h (x n-1 z ) if z ∈ V n,Ω , 0 if z ∈ V n,Γ , (3.2)
and find

u n h ∈ V n such that u n h -u n-1 h τ n , v h + (D n ∇u n h , ∇v h ) + (b n u n h , v h ) = (f n , v h ) (3.3)
holds for all v h ∈ V n . The Lagrange condition and the property x n-1 z ∈ Ω ∪ Γ for all n ∈ {1, . . . , N I } and all z ∈ V n,Ω imply that the u n-1 h are well-defined. Similarly, inequality (2.5) implies that, for every n, problem (3.3) admits a unique solution u n h . Hence the above discretization yields a unique sequence (u n h ) 0≤n≤N I . With this we associate the function u I which is continuous and linear with respect to time on each interval I n and which equals u n h at time t n , 0 ≤ n ≤ N I . Note that

u I = u n-1 h + t -t n-1 τ n (u n h -u n-1 h ) (3.4)
and

∂ t u I = u n h -u n-1 h τ n on every time-interval I n .
We conclude this section with some notation which will be useful for the formulation of the a posteriori error estimates. For every n ∈ {0, . . . , N I } we denote by N n and E n the collection of all vertices respective (d -1)-dimensional faces corresponding to the partition T n . A subscript Ω or Γ to any of these sets indicates that only those vertices or faces are considered which are contained in the respective set. Similarly, N K and E K denote the vertices and faces of a given element K and N E is the collection of the vertices of a given face E.

With every face E we associate a unit vector n E which is perpendicular to E and which coincides with the exterior normal of Ω if E is contained in Γ. For every piece-wise continuous function v and every face E which is not contained in Γ, we denote by J E (v) the jump of v in direction n E . It depends on the orientation of n E , but quantities of the form J E (n E D∇v) are independent thereof.

Finally, we denote by D τ , a τ , b τ and f τ functions which are piece-wise constant on the time-intervals and which, on every interval (t n-1 , t n ], equal D n , a n , b n and f n , respectively. Similarly, D τ,h , a τ,h , b τ,h and f τ,h are piece-wise constant functions on the time-intervals which, on every interval (t n-1 , t n ], equal the L 2 -projection of D n , a n , b n and f n respectively onto the space of piece-wise constant functions corresponding to T n .

The equivalence of error and residual

With the function u I defined from the solution of problems (3.1), (3.2) and (3.3) we associate the residual R(u

I ) ∈ L 2 (0, T ; H -1 (Ω)) by setting for every v ∈ H 1 0 (Ω) R(u I ), v = (f, v) -(∂ t u I , v) -(D∇u I , ∇v) -(a • ∇u I , v) -(bu I , v). ( 4
R(u I ) L 2 (t n-1 ,tn;H -1 (Ω)) ≤ √ 2 max{1, κ, c b } u -u I X(t n-1 ,tn)
and from above by

u -u I X(0,tn) ≤ 2(1 + max{κ, c b } 2 ) u 0 -π 0 u 0 2 + 2(2 + max{κ, c b } 2 ) R(u I ) 2 L 2 (0,tn;H -1 (Ω)) 1 2
where κ and c b are introduced in section 2.

Proof. The proof of Lemma 4.1 is based on standard parabolic energy estimates and is a straightforward modification of the proof of Lemma 4.1 in [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF]. We therefore only sketch the crucial steps. Equations (2.2) and (4.1) yield for every

v ∈ H 1 0 (Ω) R(u I ), v = (∂ t (u -u I ), v) + (D∇(u -u I ), ∇v) + (a • ∇(u -u I ), v) + (b(u -u I ), v). ( 4
v ∈ H 1 0 (Ω) R(u I ), v ≤ |(∂ t (u -u I ) + a • ∇(u -u I ))(•, t) | * |v | + max{κ, c b } |(u -u I )(•, t) | |v |.
Taking into account the definition of X(t n-1 , t n ) and of its norm, this estimate proves the first part of the lemma. To prove the second part, we insert v = (uu I )(•, t) in (4.2). Taking into account inequality (2.5) this gives 1 2

d dt (u -u I )(•, t) 2 + |(u -u I )(•, t) | 2 ≤ |R(u I )(•, t) | * |(u -u I )(•, t) |
and thus

d dt (u -u I )(•, t) 2 + |(u -u I )(•, t) | 2 ≤ |R(u I )(•, t) | 2 * .
Integrating this estimate with respect to time yields for every n

u -u I 2 L ∞ (0,tn;L 2 (Ω)) ≤ u 0 -π 0 u 0 2 + R(u I ) 2 L 2 (0,tn;H -1 (Ω))
and

u -u I 2 L 2 (0,tn;H 1 0 (Ω)) ≤ u 0 -π 0 u 0 2 + R(u I ) 2 L 2 (0,tn;H -1 (Ω))
. Equation (4.2) and inequality (2.6), on the other hand, imply

|∂ t (u -u I ) + a • ∇(u -u I ) | * ≤ |R(u I ) | * + max{κ, c b } |u -u I |.
Taking the square of this inequality, integrating the result from 0 to t n and combining this with the previous two estimates completes the proof of the lemma.

Decomposition of the residual

We may write the residual R(u I ) as the sum of three contributions: a term R osc (u I ) representing oscillations of the data and the coefficients with respect to time, a temporal residual R τ (u I ) and a spatial residual R h (u I ). For every v ∈ H 1 0 (Ω) and every interval (t n-1 , t n ] these are defined by

R osc (u I ), v = (f -f n , v) + ((D n -D)∇u I , ∇v) + ((a n -a) • ∇u I , v) + ((b n -b)u I , v), R τ (u I ), v = (D n ∇(u n h -u I ), ∇v) + (a n • ∇(u n h -u I ), v) + (b n (u n h -u I ), v), R h (u I ), v = (f n , v) -(∂ t u I , v) -(D n ∇u n h , ∇v) (5.1) -(a n • ∇u n h , v) -(b n u n h , v
). Since R osc (u I ) only involves the known discrete solution u I and temporal oscillations of the data and coefficients of (2.1), the task of deriving upper and lower bounds for the norms R(u I ) L 2 (t n-1 ,tn;H -1 (Ω)) reduces to the estimation of the norms R τ (u I ) + R h (u I ) L 2 (t n-1 ,tn;H -1 (Ω)) . The following lemma shows that this can be achieved by estimating the contributions of R τ (u I ) and R h (u I ) separately. It is implicitly established in the proof of [25, Lemma 7.1]. Lemma 5.1. For every n ∈ {1, . . . , N I } we have

1 5 R τ (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) + R h (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) 1 2 ≤ R τ (u I ) + R h (u I ) L 2 (t n-1 ,tn;H -1 (Ω)) ≤ R τ (u I ) L 2 (t n-1 ,tn;H -1 (Ω)) + R h (u I ) L 2 (t n-1 ,tn;H -1 (Ω)) .
Proof. The second inequality is the triangle inequality. Thus, we only have to prove the first inequality. Since R h (u

I ) is constant on (t n-1 , t n ], we have R h (u I ) L 2 (t n-1 ,tn;H -1 (Ω)) = √ τ n |R h (u I ) | * . (5.2) Define r n ∈ H -1 (Ω) by setting for every v ∈ H 1 0 (Ω) r n , v = (D n ∇(u n h -u n-1 h ), ∇v) + (a n • ∇(u n h -u n-1 h ), v) + (b n (u n h -u n-1 h ), v).
(5.

3)

The definitions of R τ (u I ) and u I then yield

R τ (u I ), v = t n -t τ n r n , v .
Hence, we have

R τ (u I ) L 2 (t n-1 ,tn;H -1 (Ω)) = τ n 3 |r n | * . (5.4)
Next, owing to the Riesz isomorphism between the Hilbert spaces H 1 0 (Ω) and H -1 (Ω), there exist functions v and w in

H 1 0 (Ω) such that |v | = |R h (u I ) | * , R h (u I ), v = |R h (u I ) | 2 * , |w | = |r n | * , r n , w = |r n | 2 * .
For every α > 0 we obtain

(α + 1) t -t n-1 τ n α v + t n -t τ n w L 2 (t n-1 ,tn;H 1 (Ω)) ≤ α + 1 √ 2α + 1 √ τ n |v | + τ n 3 |w | = α + 1 √ 2α + 1 √ τ n |R h (u I ) | * + τ n 3 |r n | * ≤ α 2 + 4α + 2 2α + 1 R h (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) + R τ (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) 1 2 .
We also have

tn t n-1 R h (u I ) + R τ (u I ), (α + 1) t -t n-1 τ n α v + t n -t τ n w dt = τ n R h (u I ), v + τ n 2 R h (u I ), w + r n , v tn t n-1 t n -t τ n (α + 1) t -t n-1 τ n α dt + τ n 3 r n , w .
Computing the last integral by integration by parts yields

tn t n-1 R h (u I ) + R τ (u I ), (α + 1) t -t n-1 τ n α v + t n -t τ n w dt = τ n |R h (u I ) | 2 * + τ n 2 R h (u I ), w + τ n α + 2 r n , v + τ n 3 |r n | 2 * .
Next we use the inequality ab ≥ -

√ 3 2 a 2 -1 2 √ 3 b 2 twice, which gives tn t n-1 R h (u I ) + R τ (u I ), (α + 1) t -t n-1 τ n α v + t n -t τ n w dt ≥ 1 - √ 3 4 - √ 3 2α + 4 R h (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) + R τ (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) . Combining these estimates we arrive at R τ (u I ) + R h (u I ) L 2 (t n-1 ,tn;H -1 (Ω)) ≥ 2α + 1 α 2 + 4α + 2 1 - √ 3 4 - √ 3 2α + 4 R h (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) + R τ (u I ) 2 L 2 (t n-1 ,tn;H -1 (Ω)) 1 2 .
Since α > 0 is arbitrary, we may in particular insert α = 2 in this inequality. Since

5 14 1 -3 √ 3 8 ≥ 1
5 , this proves the first inequality of the lemma.

Estimation of the temporal residual

Equation (5.4) shows that the estimation of the temporal residual requires a reliable, efficient and robust error indicator for the functionals r n of (5.3). The following lemma is a first step in this direction. It is a slight modification of [START_REF]Robust a posteriori error estimates for stationary convection-diffusion equations[END_REF]Lemma 3.1] in that we consider a more general diffusion here. Lemma 6.1. For every n ∈ {1, . . . , N I } the dual norm of the functional r n defined in equation (5.3) can be bounded from above and from below by

1 2 + max{κ, c b } |u n h -u n-1 h | + |a n • ∇(u n h -u n-1 h ) | * ≤ |r n | * ≤ max{κ, c b } |u n h -u n-1 h | + |a n • ∇(u n h -u n-1 h ) | *
, where κ and c b are introduced in section 2.

Proof. Since the proof of Lemma 6.1 is very similar to the one of [START_REF]Robust a posteriori error estimates for stationary convection-diffusion equations[END_REF]Lemma 3.1], we only sketch the crucial steps. The second inequality immediately follows from equation (5.3) and inequality (2.6). To prove the first inequality, we observe that equation (5.3) and inequality (2.5) yield

|r n | * ≥ r n , u n h -u n-1 h |u n h -u n-1 h | ≥ |u n h -u n-1 h |.
On the other hand, owing to the definition of the dual norm, using the Riesz isomorphism implies that there exists a function v

∈ H 1 0 (Ω) with |v | = 1 and (a n • ∇(u n h -u n-1 h ), v) = |a n • ∇(u n h -u n-1 h ) | * . We insert v in (5.
3) and obtain with the help of (2.6)

|r n | * ≥ r n , v ≥ |a n • ∇(u n h -u n-1 h ) | * -max{κ, c b } |u n h -u n-1 h |.
The first inequality of the lemma now follows by considering all convex combinations of both estimates for |r n | * and taking the supremum with respect to the convex weight.

In contrast to |u

n h -u n-1 h |, the quantity |a n •∇(u n h -u n-1 h
) | * cannot be used as an error indicator since it involves a dual norm. In order to bound this quantity by a computable one denote by

c Ω = sup v∈H 1 0 (Ω)\{0}
v ∇v the optimal constant in the Poincaré-Friedrichs inequality for Ω. Then, since a n is divergence-free, we have

|a n • ∇(u n h -u n-1 h ) | * ≤ ε -1 2 min{β -1 2 , ε -1 2 c Ω } a n ∞ |u n h -u n-1 h |.
In the case of a dominant diffusion or dominant reaction, using once more (2.5) and (2.6), we therefore obtain

1 2 + max{κ, c b } |u n h -u n-1 h | ≤ |r n | * ≤ max{κ, c b }{1 + c c c Ω } |u n h -u n-1 h |,
where the constants κ, c b and c c are introduced in section 2. If the convection is dominant, however, none of these estimates is robust.

In this case we have to proceed as in [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF]Lemma 8.1] (note that the present situation corresponds to the case θ = 1, T h,n = T n in this reference). For every n ∈ {1, . . . , N I } denote by Y n the space of continuous piecewise linear finite element functions vanishing on Γ associated with T n and by v n h ∈ Y n the unique solution of the discrete diffusion-reaction problem

ε(∇v n h , ∇ϕ h ) + β(v n h , ϕ h ) = (a n • ∇(u n h -u n-1 h ), ϕ h ) (6.1) for all ϕ h ∈ Y n . For every S ∈ T n ∪ E n set α S = min{ε -1 2 diam(S), β -1 2 }
and define the error indicator η n τ by

η n τ = K∈Tn α 2 K a n • ∇(u n h -u n-1 h ) + ε∆v n h -βv n h 2 K + E∈En ε -1 2 α E J E (n E • ∇v n h ) 2 E 1 2 . (6.2)
The transition condition ensures that u n h and u n-1 h , which are defined on different meshes, are piece-wise polynomials on the common refinement which has a local mesh-size comparable to the one of T n . With the same arguments as in the proof of [START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF]Lemma 8.1] we therefore obtain: Lemma 6.2. There are two constants c * and c * which depend only on the ratios h K /ρ K such that the following estimates are valid

c * |v n h | + η n τ ≤ |a n • ∇(u n h -u n-1 h ) | * ≤ c * |v n h | + η n τ .

Estimation of the spatial residual

Equation (5.2) shows that the estimation of the spatial residual requires, on each time-interval I n separately, a reliable, efficient and robust error indicator for the restriction of R h (u I ) to I n . Using integration by parts element-wise we easily conclude from the definition (5.1) of R h (u I ) that it admits an L 2 -representation consisting of element residuals and inter-element jumps. Thus it falls into the class of functionals amenable to a posteriori error estimates (cf., e.g., [START_REF] Veeser | Explicit upper bounds for dual norms of residuals[END_REF]27]). A comparison of equations (3.3) and (5.1), however, shows that R h (u I ) does not satisfy the Galerkin orthogonality, i.e. its kernel does not contain the space Y n introduced in Section 6. Instead, for every n ∈ {1, . . . , N I } and every v h ∈ V n , we have

R h (u I ), v h = u n-1 h -u n-1 h τ n -a n • ∇u n h , v h . (7.1)
This term describes the difference between the method of characteristics and the backward Euler scheme and must be taken into account by the error indicator. The following lemma shows a possible way to achieve this. The underlying idea is the same as for Lemma 6.2. Making use of the space Y n and of the piece-wise constant functions f τ,h , D τ,h , a τ,h and b τ,h defined at the end of section 3, we introduce the unique solution of the discrete diffusion-reaction problem

ε(∇w n h , ∇ϕ h ) + β(w n h , ϕ h ) = u n-1 h -u n-1 h τ n -a n • ∇u n h , ϕ h (7.2)
for all ϕ h ∈ Y n and define the error indicator η n h by

η n h = K∈Tn α 2 K f τ,h - u n h -u n-1 h τ n + div(D τ,h ∇u n h ) -a τ,h • ∇u n h -b τ,h u n h + ε∆w n h -βw n h 2 K + E∈En ε -1 2 α E J E (ε -1 n E • (D τ,h ∇u n h ) -n E • ∇w n h ) 2 E 1 2 (7.3)
Lemma 7.1. There are two constants c † and c † which depend only on the ratios h K /ρ K such that the following estimates are valid for every n ∈ {1, . . . , N I }

|w n h | + c † η n h -Θ n h ≤ |R h (u I ) | * ≤ |w n h | + c † η n h + Θ n h ,
where the term Θ n h , defined by

Θ n h = K∈Tn α 2 K f n -f τ,h + div((D n -D τ,h )∇u n h ) -(a n -a τ,h ) • ∇u n h -(b n -b τ,h )u n h 2 K + E∈En ε -1 2 α E J E (ε -1 n E • ((D n -D τ,h )∇u n h )) 2 E 1 2 , (7.4)
only depends on the oscillations of the data and coefficients.

Proof. Since Y n is contained in V n we conclude from (7.1) that ε(∇w n h , ∇ϕ h ) + β(w n h , ϕ h ) = R h (u I ), ϕ h holds for all ϕ h ∈ Y n . This in particular implies

|w n h | + |w n -w n h | ≤ |R h (u I ) | * ≤ |w n h | + |w n -w n h | .
Integration by parts element-wise yields the following L 2 -representation for the residual of

w n h R h (u I ), ϕ -ε(∇w n h , ∇ϕ) -β(w n h , ϕ) = K∈Tn K f n - u n h -u n-1 h τ n + div(D n ∇u n h ) -a n • ∇u n h -b n u n h + ε∆w n h -βw n h ϕ + E∈En E J E (n E • (D n ∇u n h ) -εn E • ∇w n h )ϕ.
The transition condition once more ensures that u n h and u n-1 h , which are defined on different meshes, are piece-wise polynomials on the common refinement which has a local mesh-size comparable to the one of T n . Standard techniques [START_REF]Robust a posteriori error estimators for a singularly perturbed reaction-diffusion equation[END_REF] therefore yield

c † η n h -Θ n h ≤ |w n -w n h | ≤ c † η n h + Θ n h
and thus prove the lemma.

The a posteriori error estimates

We may now collect the results of the preceding sections and formulate our main result.

u -u I X(0,T ) ≤ c ♯ u 0 -u 0 h 2 + f -f τ + div((D -D τ )∇u I -(a -a τ ) • ∇u I -(b -b τ )u I 2 L 2 (0,T ;H -1 (Ω)) + N I n=1 τ n |u n h -u n-1 h | 2 + |v n h | 2 + (η n τ ) 2 + N I n=1 τ n |w n h | 2 + (η n h ) 2 + (Θ n h ) 2 1 2
and, for every n ∈ {1, . . . , N I }, from below by

τ 1 2 n |u n h -u n-1 h | + |v n h | + η n τ + |w n h | + η n h ≤ c ♯ u -u I 2 X(t n-1 ,tn) + f -f τ + div((D -D τ )∇u I -(a -a τ ) • ∇u I -(b -b τ )u I 2 L 2 (t n-1 ,tn;H -1 (Ω)) + τ n (Θ n h ) 2 1 2 .
The functions v n h and w n h and the quantities η n τ , η n h and Θ n h are defined in equations (6.1), (7.2), (6.2), (7.3) and (7.4) respectively. The constants c ♯ and ♯ depend on the ratios h K /ρ K but are independent of any relative size of the diffusion with respect to the convection or reaction. If the diffusion or reaction are dominant, the quantities |v n h | and η n τ can be omitted.

Remark 8.2. According to Section 6, the terms τ

1 2 n |u n h -u n-1 h | + |v n h | + η n τ control the restriction to (t n-1 , t n ] of the temporal resid- ual R τ (u I ).
Hence, they can be used for adapting the time-step τ n . Evidently, they split into element-contributions:

τ 1 2 n |u n h -u n-1 h | + |v n h | + η n τ = τ n K∈Tn |u n h -u n-1 h | 2 K + |v n h | 2 K + η n τ,K 2 1 2 with η n τ,K = α 2 K a n • ∇(u n h -u n-1 h ) + ε∆v n h -βv n h 2 K + 1 2 E∈E K,Ω ε -1 2 α E J E (n E • ∇v n h ) 2 E 1 2 .
Thus, the terms τ

1 2 n |u n h -u n-1 h | K + |v n h | K + η n τ,K are
natural candidates for an eventual local (in space) time-step control. We once more stress that the terms |v n h | and η n τ and their local analogues can be omitted in the case of dominant reaction or diffusion. If, however, the convection is dominant these terms are mandatory for conserving the full robustness of the error estimates. n Θ n h . Evidently, these terms split into elementcontributions:

τ 1 2 n |w n h | + η n h = τ n K∈Tn |w n h | 2 K + η n h,K 2 1 2 , with η n h,K = α 2 K f τ,h - u n h -u n-1 h τ n + div(D τ,h ∇u n h ) -a τ,h • ∇u n h -b τ,h u n h + ε∆w n h -βw n h 2 K + 1 2 E∈E K,Ω ε -1 2 α E J E (ε -1 n E • (D τ,h ∇u n h ) -n E • ∇w n h ) 2 E 1 2 .
Taking into account that the spatial residual is piece-wise constant on the time intervals, combining the arguments used in the proof of Lemma 7.1 with the techniques of [START_REF]Robust a posteriori error estimators for a singularly perturbed reaction-diffusion equation[END_REF] for establishing lower error bounds and defining Θ n h,K in an obvious way as the local analogue of Θ n h , we can easily prove that

|w n h | K + η n h,K ≤ c sup v∈H 1 0 (e ω K )\{0} R h (u I ), v |v | K + Θ n h,K
holds for every interval (t n-1 , t n ] and every element K ∈ T n . Here, ω K denotes the union of all elements in T n that share at least a vertex with K. The constant c only depends on the shape parameter of T n but is independent of any relative size of the diffusion with respect to the convection or reaction. This shows that the terms

|w n h | K + η n h,K
can efficiently be used for adapting the spatial meshes.

Remark 8.4. We finally want to mention that the approach of [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF][START_REF]A posteriori error analysis of the fully discretized time-dependent Stokes equations[END_REF] could also be adapted to the present problem. Instead of splitting the residual as in Lemma 5.1, this would mean to split the error in the form uu τ X(0,T ) + u τu I X(0,T ) by introducing a semi-discrete approximation u τ as follows: Set u 0 = π 0 u 0 . For n = 1, . . . , N I define u n-1 by u n-1 (y) = u n-1 (x(t n-1 ; y, t n )) for all y ∈ Ω and find

u n in H 1 0 (Ω) such that u n -u n-1 τ n , v + (D n ∇u n , ∇v) + (b n u n , v) = (f n , v)
holds for all v in H 1 0 (Ω). Then, u τ is the continuous, piece-wise linear (with respect to time) function which coincides with u n at time t n , 0 ≤ n ≤ N I . Standard techniques can then be used to obtain upper and lower bounds for the terms u-u τ X(0,T ) and u τ -u I X(0,T ) separately. At first sight this approach is appealing since it yields an obvious uncoupling of spatial and temporal error contributions. But it has the severe drawback that it leads to error estimates which are not fully robust with respect to the relative sizes of diffusion, convection and reaction. Essentially, this is due to the fact that uu I X(0,T ) can be bounded from above by uu τ X(0,T ) + u τu I X(0,T ) but not from below. The present approach does not suffer from this drawback since Lemma 5.1 provides two-sided bounds.

As a conclusion, when replacing D n , a n , b n and f n by their finite element approximations D τ,h , a τ,h , b τ,h and f τ,h , the quantities v n h , w n h , η n τ and η n h (or their local analogues η n τ,K and η n h,K ) are very easy to compute. Hence, they are an efficient tool for mesh and time adaptivity, for all types of diffusion-convection-reaction equations.

Some numerical experiments

The numerical simulation that we now present has been performed with the code FreeFem++ due to F. Hecht and O. Pironneau, see [START_REF] Hecht | Freefem++[END_REF]. We only consider a model domain, the square Ω = (-1, 1) 2 . The space V n is the space of continuous piece-wise linear finite element functions vanishing on Γ associated with T n . Thus, it coincides with the space Y n introduced in section 6.

The diffusion coefficient D is equal to ε times the identity, and we intend to work with different values of ε. The coefficient b is taken equal to zero and the convection coefficient a is constant in time, given by

a(x) = curl(1 -x 2 1 ) 2 (1 -x 2 2 ) 2 = 4(1 -x 2 1 )(1 -x 2 2 ) -(1-x 2 1 )x 2 x 1 (1-x 2 2 )
.

The data f and u 0 are chosen such that the exact solution is

u(x, t) = v( 1 ε |x -z(t)| 2 ),
with

z(t) = 1 2 cos t 1 2 sin t , v(s) = (1 + 2s)(1 -s) 2 if 0 ≤ s ≤ 1 0 if s > 1.
The final time T is equal to 4. For simplicity, we do not perform time adaptivity and take all the τ n equal to 10 -1 . Mesh adaptivity is performed according to the following simple strategy where ℓ ∈ {0, 1, 2, 3}:

• Given the triangulation T 10ℓ , for n ∈ {10ℓ, . . . , 10ℓ+9}, we compute the solution u n h , the indicators η n τ,K and η n h,K , the solutions v n h and w n h together with their local norms |v n h | K and |w n h | K and set

H n K = |u n h -u n-1 h | K + η n τ,K + |v n h | K + η n h,K + |w n h | K .
• Next, we define H ℓ K as the maximum of the H n K , for n ∈ {10ℓ, . . . , 10ℓ + 9}.

• Then we construct a new triangulation T 10(ℓ+1) such that the diameter of any element in T 10(ℓ+1) which contains or is contained in an element K of T 10ℓ is proportional to the diameter of K times the ratio H * /H ℓ K , where H * is a fixed threshold.

• Finally, we perform a new computation of the solution u n h , n ∈ {10ℓ, . . . , 10ℓ + 9}, on the new mesh. We refer to [START_REF] Hecht | Mesh adaptation for time dependent simulation and optimization[END_REF] for more details on this algorithm and to [START_REF] Frey | Maillages, applications aux éléments finis[END_REF]Chap. 21] for the way of constructing such a mesh. Note also that this adaptation step can be iterated several times when necessary.

. 1 )Lemma 4 . 1 .

 141 The following lemma shows that R(u I ) and the error u-u I are equivalent. It is a slight modification of[START_REF]Robust a posteriori error estimates for nonstationary convectiondiffusion equations[END_REF] Lemma 4.1], since a more general diffusion term is considered in the present work. Recall that H 1 0 (Ω) and its dual space H -1 (Ω) are equipped with the energy-norm |• | and the dual norm |• | * respectively. For every n ∈ {1, . . . , N I } the error can be bounded from below by

. 2 )

 2 This identity, inequality (2.6) and the definition of the norms |• | and |• | * yield for all 0 < t < T and all

  |w n h | ≤ |R h (u I ) | * . Denote by w n ∈ H 1 0 (Ω) the unique solution of the weak diffusion-reaction equation ε(∇w n , ∇ϕ) + β(w n , ϕ) = R h (u I ), ϕ for all ϕ ∈ H 1 0 (Ω). Then we have |w n | = |R h (u I ) | * and therefore 1 3

Theorem 8 . 1 .

 81 Denote by u the solution of the variational formulation (2.2) of the differential equation (2.1) and by u I its discrete approximation defined by problems (3.1), (3.2), (3.3) and equation(3.4). Then the error uu I can be bounded from above by

Remark 8 . 3 . 2 n

 832 According to Section 7, the terms τ1 |w n h | + η n h control the restriction to (t n-1 , t n ] of the spatial residual R h (u I ) up to the data-oscillation τ 1 2

Figures 9 .

 9 1 -9.3 show the adapted mesh and the corresponding solution for several values of ε. In all cases the singularity of the solution is well detected and resolved.

Figure 9 . 1 .

 91 Figure 9.1. The adapted mesh and the solution at time T for ε = 1

Figure 9 . 2 .

 92 Figure 9.2. The adapted mesh and the solution at time T for ε = 10 -1

Figure 9 . 3 .

 93 Figure 9.3. The adapted mesh and the solution at time T for ε = 10 -2

  problem (2.2) admits a unique solution, see e.g. [3, Sections II.4 and V.2] or [10, Theorems XVIII.3.1 and XVIII.3.2].

Acknowledgements: We are very grateful to Brigitte Métivet for a first discussion on the subject of this paper and to Martin Vohralík for his kind invitation of the third author to visit the Laboratoire Jacques-Louis Lions.