On the time growth of the error of the DG method for advective problems

On the time growth of the error of the DG method for advective problems Abstract In this paper we derive a priori $$L^{\infty }(L^{2})$$ and L2(L2) error estimates for a linear advection–reaction equation with inlet and outlet boundary conditions. The goal is to derive error estimates for the discontinuous Galerkin method that do not blow up exponentially with respect to time, unlike the usual case when Gronwall’s inequality is used. While this is possible in special cases, such as divergence-free advection fields, we take a more general approach using exponential scaling of the exact and discrete solutions. Here we use a special scaling function, which corresponds to time taken along individual pathlines of the flow. For advection fields, where the time that massless particles carried by the flow spend inside the spatial domain is uniformly bounded from above by some $$\widehat{T}$$, we derive $$\mathcal{O}$$(hp+1/2) error estimates where the constant factor depends only on $$\widehat{T}$$, but not on the final time T. This can be interpreted as applying Gronwall’s inequality in the error analysis along individual pathlines (Lagrangian setting), instead of physical time (Eulerian setting). 1. Introduction The discontinuous Galerkin (DG) finite element method introduced in the study by Reed & Hill (1973) is an increasingly popular method for the numerical solution of partial differential equations. The DG method was first formulated for a neutron transport equation and such problems remain the major focus of the DG community. It is for problems of an advective or convective nature that the DG method is suited best and shows its strengths compared to other numerical methods for such problems. In this paper, we shall consider a scalar time-dependent linear advection–reaction equation of the form   \begin{align} \frac{\partial u}{\partial t}+a\cdotp\nabla\!\, u+cu=0. \end{align} (1.1)We will discretize the problem in space using the DG method with the upwind numerical flux on unstructured simplicial meshes in $$\mathbb{R}^{d}$$ which may contain hanging nodes. The first analysis of the DG method for a stationary advection–reaction problem with constant coefficients was made in the study by Lesaint & Raviart (1974), later improved in the study by Johnson & Pitkäranta (1986). In this paper we derive a priori estimates of the error eh = u − uh, where uh is the DG solution. The goal is to derive estimates of eh in the $$L^{\infty }\left (L^{2}\right )$$- and $$L^{2}\left (L^{2}\right )$$-norms of the order Chp+1/2, where the constant C does not depend exponentially on the final time $$T\to +\infty $$. Such results already exist in the literature; however, they are derived under the ellipticity condition   \begin{align} c-\tfrac{1}{2}\operatorname{div}a\geqslant\gamma_{0}>0 \end{align} (1.2)for some constant $$\gamma_0 $$. In the paper by Feistauer & Švadlenka (2004), $$\gamma_0=0 $$ is also admissible, which corresponds to the interesting case of a divergence-free advection field a. Without assuming (1.2), when one proceeds straightforwardly, at some point Gronwall’s inequality must be used in the proofs. This, however, results in exponential growth of the constant factor C in the error estimate with respect to time. Here, we will circumvent condition (1.2) by considering an exponential scaling transformation   \begin{align} u(x,t)=e^{\mu(x,t)}\tilde{u}(x,t),\quad u_{h}(x,t)=e^{\mu(x,t)}\tilde{u}_{h}(x,t) \end{align} (1.3)of the exact and discrete solutions. Substituting (1.3) into (1.1) and its discretization results in a new equation for $$\tilde{u}$$ with a modified reaction term $$\tilde{c}=\frac{\partial \mu }{\partial t}+a\cdotp \nabla \!\mu +c$$. The function μ can then be chosen in various ways in order to satisfy the ellipticity condition (1.2) for the new equation. After performing the error analysis, the resulting error estimates for $$\tilde{e}_{h}=\tilde{u}-\tilde{u}_{h}$$ can be transformed back via (1.3) to obtain estimates for eh. The use of transformations similar to (1.3) is not new. The trick known as exponential scaling, corresponding to taking μ(x, t) = αt for some sufficiently large constant α, is very well known from the partial differential equation and numerical analysis communities. However, the application of such a scaling transformation inherently leads to exponential dependence of the error on t after transforming from $$\tilde{e}_{h}$$ back to eh. Other choices of μ are possible, e.g. μ(x, t) = μ0⋅x for some constant vector μ0. This choice was used in the stationary case in the studies by Nävert (1982) and Johnson & Pitkäranta (1986). Taking μ(x) independent of t corresponds to the analysis of the stationary case from the study by Ayuso & Marini (2009). In the nonstationary case the use of this stationary transform would lead to too restrictive assumptions on the advection field. In this paper we construct the function μ from (1.3) using characteristics of the advection field. Namely, μ will be proportional to time along individual pathlines of the flow. If pathlines exist only for a finite time bounded by some $$\widehat{T}$$ for each particle entering and leaving the spatial domain $$\varOmega$$, we obtain error estimates exponential in $$\widehat{T}$$ and not the final physical time T. This is a result we would expect if we applied Gronwall’s lemma in the Lagrangian framework along individual pathlines (which exist only for a finite time uniformly bounded by $$\widehat{T}$$) instead of the usual application of Gronwall with respect to physical time in the Eulerian framework. The analysis can be carried out under mild assumptions that a(⋅, ⋅) satisfies the assumptions of the Picard–Lindelöf theorem and that there are no characteristic boundary points on the inlet (i.e. a ⋅ n is uniformly bounded away from zero on the inlet). The paper is organized as follows. After introducing the continuous problem in Section 2, we discuss the exponential scaling transform, its variants and application in the weak formulation in Section 3. In Section 4 we construct the scaling function μ and prove results on its regularity and other properties needed in the analysis. In Section 5 we introduce the DG formulation and its basic properties. In Sections 6 and 7 we analyse the DG advection and reaction forms and estimate the error of the method. In order to focus on the main ideas of the analysis, we postpone some of the technical estimates of the DG forms to the Appendix. We use (⋅, ⋅) to denote the L2($$\varOmega$$) scalar product and ∥⋅∥ for the L2($$\varOmega$$)-norm. To simplify the notation, we shall drop the argument $$\varOmega$$ in Sobolev norms, e.g. $$\|\cdot \|_{H^{p+1}}$$ denotes the Hp+1($$\varOmega$$)-norm. We will also denote the Bochner norms over the whole considered interval (0, T) in concise form, e.g. $$\|u\|_{L^{2}\left (H^{p+1}\right )}$$ denotes the $$L^{2}(0,T;H^{p+1}(\varOmega ))$$-norm of u. 2. Continuous problem Let $$\varOmega \subset \mathbb{R}^{d}$$, $$d\in \mathbb{N}$$ be a bounded polygonal (polyhedral) domain with Lipschitz boundary ∂$$\varOmega$$. Let $$0<T\leqslant +\infty $$ and let QT = $$\varOmega$$ × (0, T) be the space–time domain. We note that we admit the time interval to be infinite in the special case of $$T=+\infty $$. We seek $$u:Q_{T}\to \mathbb{R}$$ such that   \begin{alignat}{2} \frac{\partial u}{\partial t}+a\cdotp\nabla\!\, u+cu&=0 &&\textrm{in }Q_{T},\\ u&=u_{D} &&\textrm{on }\partial\varOmega^{-}\times(0,T),\nonumber\\ u(x,0)&=u^{0}(x),\quad && x\in\varOmega\nonumber. \end{alignat} (2.1)Here $$a:\overline{Q_{T}}\to \mathbb{R}^{d}$$ and $$c:\overline{Q_{T}}\to \mathbb{R}$$ are the given advective field and reaction coefficient, respectively. By ∂$$\varOmega$$− we denote the inflow part of the boundary, i.e.   \begin{align} \partial\varOmega^{-}=\left\{x\in\partial\varOmega; a(x,t)\cdot\mathbf{n}(x)<0\quad \forall\; t\in(0,T)\right\}\!, \nonumber \end{align}where n(x) is the unit outer normal to ∂$$\varOmega$$ at x. For simplicity we assume that the inflow boundary is independent of t, in other words we assume that a(x, t) ⋅ n(x) < 0 for all t and all x ∈ ∂$$\varOmega$$−, and a(x, t) ⋅ n(x) is non-negative on ∂$$\varOmega$$ ∖ ∂$$\varOmega$$− for all t. In general ∂$$\varOmega$$− can depend on t due to the nonstationary nature of a. However, we do not consider this case since it would only make the notation more complicated without changing the key points of the analysis itself. We assume that the reaction coefficient satisfies $$c\in C([0,T);L^{\infty }(\varOmega ))\cap L^{\infty }(Q_{T})$$. Furthermore, let $$a\in C([0,T);W^{1,\infty }(\varOmega ))$$ with a, ∇a uniformly bounded a.e. in QT. In other words, a (⋅, ⋅) satisfies the assumptions of the Picard–Lindelöf theorem: continuity with respect to x, t and uniform Lipschitz continuity with respect to x. The spatial Lipschitz constant of a will be denoted by La. The DG method for the stationary version of problem (2.1) was analysed in the studies by Nävert (1982), Johnson & Pitkäranta (1986) and Ayuso & Marini (2009). In this case the well-posedness of the problem is guaranteed by the assumption that the flow field possesses neither closed curves nor stationary points, cf. the study by Devinatz et al. (1974). For the nonstationary problem (2.1), the situation is simpler. Here, the problem can be solved using the method of characteristics, similarly to the method we use in Section 4, cf. Evans (2010). Along each characteristic curve (pathline), problem (2.1) reduces to a simple, ordinary differential equation, which can be solved in closed form. The result is exponential growth or decay of u along the pathline with the rate given by the reaction coefficient c. As the flow field a satisfies the assumptions of the Picard–Lindelöf theorem, we have existence and uniqueness of the pathlines. Since the Dirichlet boundary condition is prescribed only at the inflow boundary, this gives the well-posedness of problem (2.1). We note that also a ‘parabolic’ approach to (2.1) is possible to obtain ‘energy’ estimates of the time growth of the L2($$\varOmega$$)-norm of the solution. This is done by Gronwall’s lemma after testing the weak form of (2.1) with the solution itself and applying Green’s theorem in the advection terms; cf. e.g. Mizohata (1973) and Larsson & Thomée (2003). The result is an upper bound on the time growth of the solution that is exponential in the quantity $$c-\tfrac{1}{2}\textrm{div} a$$ if this quantity is non-negative; compare with the ellipticity condition (3.1). We note that from the theoretical and practical points of view, the case of systems of equations of the considered type is much more interesting; cf. e.g. Benzoni-Gavage & Serre (2007). 3. Exponential scaling The standard assumption on the coefficients a, c found throughout the numerical literature is   \begin{align} c-\tfrac{1}{2}\operatorname{div} a\geqslant\gamma_{0}>0 \quad\textrm{ on }Q_{T} \end{align} (3.1)for some constant $$\gamma_0 $$. This assumption comes from the requirement that the weak formulations of the advection and reaction terms give an elliptic bilinear form on the corresponding function space. In the study by Feistauer & Švadlenka (2004), this assumption is avoided by using exponential scaling in time, i.e. the transformation u(x, t) = eαtw(x, t) for $$\alpha \in \mathbb{R}$$. Substitution into (2.1) and division of the whole equation by the common positive factor eαt leads to a modified reaction coefficient $$\tilde c=c+\alpha $$ in the new equation for w. By choosing the constant α sufficiently large, (3.1) will be satisfied for the new equation. In the study by Feistauer & Švadlenka (2004), error estimates that grow linearly in time are then derived for the DG scheme. However, for α > 0, if one transforms the resulting estimates back to the original problem using the exponential scaling transformation, the result is an estimate that depends exponentially on T. Another possibility to avoid assumption (3.1) is a transformation similar to exponential scaling, but with respect to the spatial variable; cf. the studies by Nävert (1982), Johnson & Pitkäranta (1986) and Roos et al. (2008): let $$\mu _{0}\in \mathbb{R}^{d}$$, then we write   \begin{align} u(x,t)=e^{\mu_{0}\cdotp x}\tilde u(x,t). \end{align} (3.2)If we assume that u is sufficiently regular, as we shall in this paper, by substituting (3.2) into (2.1) and dividing the whole equation by the strictly positive function $$e^{\mu _{0}\cdotp x}$$, we get the new problem   \begin{align} \frac{\partial \tilde{u}}{\partial t}+a\cdotp\nabla\!\, \tilde{u}+(a\cdotp\mu_{0}+c)\tilde{u}=0. \end{align} (3.3)The condition corresponding to (3.1) now reads as follows: there exists $$\mu _{0}\in \mathbb{R}^{d}$$ such that   \begin{align} a\cdotp\mu_{0}+c-\tfrac{1}{2}\operatorname{div} a\geqslant\gamma_{0}>0 \quad\textrm{ on }Q_{T}. \end{align} (3.4) A possible generalization of the exponential scaling transformation (3.2) is taking a sufficiently smooth function $$\mu :\varOmega \to \mathbb{R}$$ and setting   \begin{align} u(x,t)=e^{\mu(x)}\tilde u(x,t). \end{align} (3.5)Again, substituting into (2.1) and dividing by eμ(x), we get the new problem   \begin{align} \frac{\partial \tilde{u}}{\partial t}+a\!\cdot\nabla \tilde{u}+(a\cdotp\nabla\!\mu+c)\tilde{u}=0. \nonumber \end{align}The condition corresponding to (3.1) and (3.4) now reads as follows: there exists $$\mu :\varOmega \to \mathbb{R}$$ such that   \begin{align} a\cdotp\nabla\!\mu+c-\tfrac{1}{2}\operatorname{div} a\geqslant\gamma_{0}>0 \quad\textrm{ on }Q_{T}. \end{align} (3.6)This is essentially the approach used in the study by Ayuso & Marini (2009) for a stationary advection–diffusion–reaction problem. As shown in the study by Devinatz et al. (1974), the existence of a function $$\mu :\varOmega \to \mathbb{R}$$ such that $$a\cdotp \nabla \!\mu \geqslant \gamma _{0}$$ is equivalent to the property that the advective field a possesses neither closed curves nor stationary points. The uniformly positive term $$a\cdotp \nabla \!\mu $$ can then be used to dominate the possibly negative term $$c-\frac{1}{2}\textrm{div} a$$ in order to satisfy condition (3.6). If we used the choice (3.5) in our analysis, we would need to assume the nonexistence of closed curves or stationary points of the flow field for all t. This assumption is too restrictive. Furthermore, in the study by Ayuso & Marini (2009), rather high regularity of μ is needed in the analysis, namely $$\mu \in W^{p+1,\infty }(\varOmega )$$, where p is the polynomial degree of the DG scheme. In this paper, we shall generalize the transformation (3.5) using a function $$\mu :Q_{T}\to \mathbb{R}$$ and set   \begin{align} u(x,t)=e^{\mu(x,t)}\tilde u(x,t). \end{align} (3.7)In our analysis we will need only minimum regularity of the scaling function μ, in contrast to the study by Ayuso & Marini (2009). Namely, we assume that μ is Lipschitz continuous both in space and time. This implies that the derivatives μt and ∇μ exist for almost all x and t. Hence, the following steps are valid for almost all x and t and particularly, they will hold in the integral (i.e. weak) sense. We note that the Lipschitz continuity in space and time will hold for the specific construction of μ from Section 4, as shown in Section 4.1. As in the previous cases, substituting (3.7) into (2.1) and dividing by eμ(x, t) gives the new problem   \begin{align} \frac{\partial \tilde{u}}{\partial t}+a\cdot\nabla\! \tilde{u}+\left(\frac{\partial \mu}{\partial t}+a\cdotp\nabla\!\mu+c\right)\tilde{u}=0. \end{align} (3.8)The condition corresponding to (3.1), (3.4) and (3.6) now reads as follows: there exists $$\mu :Q_{T}\to \mathbb{R}$$ such that   \begin{align} \frac{\partial \mu}{\partial t}+a\cdotp\nabla\!\mu+c-\tfrac{1}{2}\operatorname{div} a\ge\gamma_{0}>0 \quad\text{a.e. in }Q_{T}. \end{align} (3.9)This is the condition we will assume throughout the paper. Using the transformation (3.7) instead of (3.5) gives sub-exponential growth or even uniform boundedness of the constant factor in the error estimate with respect to time. If the original problem (2.1) does not satisfy (3.1), one is tempted to numerically solve the transformed equation (3.8) instead, using DG or any other method. This is done e.g. in Roos et al. (2008) for the case of linear μ, i.e. (3.2). However, then we are numerically solving a different problem and obtain different results, as the DG solutions of (2.1) and (3.8) are not related by the simple relation (3.7), unlike the exact solutions. However, as we will show in this paper, one can analyse the DG method for the original problem (2.1), while taking advantage of the weaker ellipticity condition (3.9) for the transformed problem (3.8). Since the DG scheme is based on a suitable weak formulation, the first step is to reformulate the transformation (3.7) within the weak, rather than strong formulation of (2.1). The key step in deriving (3.8) from (2.1) is dividing the whole equation by the common factor eμ(x, t). However, if we substitute (3.7) into the weak formulation   \begin{align} \int_{\varOmega}\frac{\partial u}{\partial t}v+a\cdotp\nabla\!\, uv+cuv\,\mathrm{d} x =0, \end{align} (3.10)it is not easy to divide the equation by eμ(x, t), since this is inside the integral. The solution is to test (3.10) by the new test function $$\hat{v}(x,t)=e^{-\mu (x,t)}v(x,t)$$. Due to the opposite signs in the exponents, the exponential factors cancel each other and we obtain the weak formulation of (3.8):   \begin{align} \int_{\varOmega}\frac{\partial \tilde{u}}{\partial t}v+a\cdotp\nabla\!\, \tilde{u} v+\left(\frac{\partial \mu}{\partial t} +a\cdotp\nabla\!\mu+c\right)\tilde{u} v\,\mathrm{d} x =0. \nonumber \end{align}We note that the transformations $$u\mapsto \tilde{u}$$ and $$v\mapsto \hat{v}$$ are bijections. Furthermore, since we assume that μ is Lipschitz continuous in space, then the factor eμ(⋅, t) lies in $$W^{1,\infty }(\varOmega )$$ and the transformation (3.7) is a bijection from the function space V into itself, where   \begin{align} V=\big\{u\in H^{1}(\varOmega); u|_{\partial\varOmega^{-}}=0\big\} \nonumber \end{align}is the space for weak solutions and test functions for (2.1) we consider. This is the case uD = 0, the general case can be treated by the standard Dirichlet lifting procedure. We note that in the general case, the appropriate function spaces for the solution of (2.1) would be e.g. L1($$\varOmega$$) or BV($$\varOmega$$) with respect to space. However, since we assume smooth solutions and the linear case, it is possible to consider the Hilbert setting of H1($$\varOmega$$) which makes the analysis easier in the finite element or DG framework. In the DG method, the above procedure using transform (3.7) cannot be directly applied, since if v ∈ Sh—the discrete space of discontinuous piecewise polynomial functions—then $$\hat{v}=e^{-\mu }v$$ will no longer lie in Sh and therefore cannot be used as a test function in the formulation of the method. The solution is to test with a suitable projection of $$\hat{v}$$ onto Sh and estimate the difference, as we shall do in Section 6. 4. Construction of the function μ In this section we show a construction of the function μ satisfying (3.9). Many different constructions of μ are possible depending on the assumptions on the vector field a. For example, we can always take μ(x, t) = αt for some $$\alpha \geqslant 0$$. This corresponds to standard exponential scaling and, as we have seen, this choice leads to exponential growth in time of the error estimate. Another possibility was the approach of Ayuso & Marini (2009) mentioned in Section 2, where a suitable function μ(x) exists if a, which is stationary, possesses no closed curves or stationary points. Here we show another possibility with an interesting interpretation. If $$c-\frac{1}{2}\textrm{div} a$$ is negative or changes sign frequently, we can use the expression μt + a ⋅∇μ to dominate this term everywhere. If we choose μ1 such that   \begin{align} \frac{\partial\mu_{1}}{\partial t}+a\cdot\nabla\!\mu_{1}=1\quad \textrm{on }Q_{T}, \end{align} (4.1)then by multiplying μ1 by a sufficiently large constant, we can satisfy the ellipticity condition (3.9) for a chosen $$\gamma_0 $$ > 0. Equation (4.1) can be explicitly solved using characteristics. We define pathlines of the flow, i.e. the family of curves $$S\left (t;x_{0},t_{0}\right )$$ by   \begin{align} S(t_{0};x_{0},t_{0})=x_{0}\in\overline{\varOmega},\quad \frac{\mathrm{d} S(t;x_{0},t_{0})}{\mathrm{d} t} =a\left(S(t;x_{0},t_{0}),t\right). \end{align} (4.2)This means that S(⋅; t0, x0) is the trajectory of a massless particle in the nonstationary flow field a passing through point x0 at time t0. It is convenient to choose the parameter t0 minimal for each pathline—then the pair (x0, t0) is the ‘origin’ of the pathline. In other words, for each $$x_{0}\in \overline{\varOmega }$$, there is a pathline S(⋅; x0, 0) corresponding to trajectories of particles present in $$\varOmega$$ at the initial time t0 = 0. Then there are trajectories of particles entering through the inlet part of ∂$$\varOmega$$: for each x0 ∈ ∂$$\varOmega$$− and all t0 there exists a pathline S(⋅; x0, t0) originating at (x0, t0) ∈ ∂$$\varOmega$$−× [0, T). Equation (4.1) can then be rewritten along the pathlines as   \begin{align} \frac{\mathrm{d}\, \mu_{1}(S(t;x_{0},t_{0}),t)}{\mathrm{d} t} = \left(\frac{\partial\mu_{1}}{\partial t}+a\cdot\nabla\!\mu_{1}\right) (S(t;x_{0},t_{0}),t)=1, \nonumber \end{align}therefore   \begin{align} \mu_{1}(S(t;x_{0},t_{0}),t) =t-t_{0}. \end{align} (4.3)Here any constant can be chosen instead of t0; however, this choice is the most convenient. With this choice, we have at the origin of a pathline   \begin{align} \mu_{1}(S(t_{0};x_{0},t_{0}),t_{0}) =0 \nonumber \end{align}and the value of μ1 along this pathline is simply the time elapsed since t0. In other words, μ1(x, t) is the time a particle carried by the flow passing through x ∈$$\varOmega$$ at time t has spent in $$\varOmega$$ up to the time t. In this paper we assume this quantity to be uniformly bounded in order for μ1 and μ to be uniformly bounded; cf. assumption (4.7). It is important to note that t0 depends on the specific trajectory considered and that t0 will in general be different for each (x, t). Therefore, even though t and t0 may be arbitrarily large if we consider an infinite time interval, the difference t − t0 = μ1, which measures the time elapsed since the entry of the corresponding particle into $$\varOmega$$, will remain bounded under this assumption. Now that we have constructed a function satisfying (4.1), we can choose $$\gamma_0 $$ and define e.g.   \begin{align} \mu(x,t)=\mu_{1}(x,t)\left(\left|\inf_{Q_{T}}\left(c-\tfrac{1}{2}\operatorname{div} a\right)^{-}\right|+\gamma_{0}\right), \end{align} (4.4)where $$f^{-}=\min \{0,f\}$$ is the negative part of f. Then   \begin{align} \frac{\partial\mu}{\partial t} +a\cdot\nabla\!\mu+c-\tfrac{1}{2}\operatorname{div} a = \left|\inf_{Q_{T}}\left(c-\tfrac{1}{2}\operatorname{div} a\right)^{-}\right|+\gamma_{0} +c-\tfrac{1}{2}\operatorname{div} a \ge\gamma_{0}. \nonumber \end{align} In the main result, of this paper, Theorem 7.1 and Corollary 7.2, this choice of μ leads to estimates of the following form for the DG error eh (cf. (7.1) and (7.10)):   \begin{align} \|e_{h}\|_{L^{\infty}\left(L^{2}\right)} +\sqrt{\gamma_{0}}\|e_{h}\|_{L^{2}\left(L^{2}\right)}\leqslant Ce^{\widehat{T}(|\inf_{Q_{T}}(c-\frac{1}{2}\operatorname{div} a)^{-}|+\gamma_{0})}h^{p+1/2}, \end{align} (4.5)where $$\widehat{T}=\sup _{Q_{T}}\mu _{1}$$, i.e. $$\widehat{T}$$ is the maximal particle ‘lifetime’, i.e. the maximal time any particle carried by the flow field a spends in $$\varOmega$$, since by (4.3), μ1(x, t) is the time the particle at (x, t) carried by the flow has spent in $$\varOmega$$. If we compare this to estimates obtained by a straightforward analysis using Gronwall’s lemma without any ellipticity assumption, we would expect   \begin{align} \|e_{h}\|_{L^{\infty}(L^{2})} \leqslant Ce^{T(\sup_{Q_{T}}|c-\frac{1}{2}\operatorname{div} a|)}h^{p+1/2}. \end{align} (4.6)Comparing (4.6) to (4.5), we see that the estimates are essentially of similar form, only the exponential dependence on global physical time T has been replaced by dependence on time $$\widehat{T}$$ along pathlines, which is bounded. Effectively, our analysis replaces the application of Gronwall’s lemma in the Eulerian framework with its application in the Lagrangian framework—along individual pathlines which have bounded length. Throughout the paper, we assume the particle lifetime $$\widehat{T}$$ to be bounded. In general, we could consider dependencies of the form   \begin{align} \widehat{T}(t)=\sup_{(x,\vartheta)\in\varOmega\times{(0,t)}} \mu_{1}(x,\vartheta), \nonumber \end{align}i.e. $$\widehat{T}(t)$$ is the maximal time any particle carried by the flow a spends in $$\varOmega$$ up to time t. From (4.5), we can expect exponential dependence of the $$L^\infty(0,t;L^2)$$ and $$L^2(0,t;L^2)$$ norms on $$\widehat{T}(t)$$. The result of the standard analysis (4.6) corresponds to the ‘worst case’ $$\widehat{T}(t)=t$$, i.e. that there is a particle that stays inside $$\varOmega$$ for all t ∈ [0, T). However, considering more general dependencies on time is possible, e.g. $$\widehat{T}=\sqrt{t}$$, leading to growth of the error that is exponential in $$\sqrt{t}$$. 4.1 Regularity of the function μ In our analysis of the DG scheme we will assume that μ satisfies   \begin{equation}\label{4.7}\begin{split} 0\leqslant\mu(x,t)\leqslant&\,\mu_{\max},\\ |\mu(x,t)-\mu(y,t)|\leqslant&\, L_{\mu}|x-y|,\end{split}\end{equation} (4.7)for all x, y ∈ $$\varOmega$$ and t ∈ (0, T). In other words, μ is non-negative, uniformly bounded and Lipschitz continuous in space, where the Lipschitz constant is uniformly bounded for all t. Since $$\varOmega$$ is a Lipschitz domain, hence quasiconvex, this means that $$\mu (t)\in W^{1,\infty }(\varOmega )$$ for all t, with its $$W^{1,\infty }(\varOmega )$$ seminorm uniformly bounded by Lμ for all t. Now we show when μ1 defined by (4.3) satisfies conditions (4.7), especially Lipschitz continuity in space which is necessary in the analysis. We note that since μ is obtained from μ1 by multiplication by a constant, properties derived for μ1 imply those for μ. An obvious example of when μ1 defined by (4.3) will not be Lipschitz continuous is when a vortex ‘touches’ ∂$$\varOmega$$− as in Fig. 1. Then if (x, t) lies on a pathline originating at (x0, t0) such that a(x0, t0) is tangent to ∂$$\varOmega$$, we can find $${\tilde{x}}$$ arbitrarily close to x such that the corresponding pathline is much longer, perhaps winding several times around the vortex. Then $$\mu _{1}(x,t)-\mu _{1}\left ({\tilde{x}},t\right )={\tilde{t}}_{0}-t_{0}$$ can be very large, while $$\|x-{\tilde{x}}\|$$ is arbitrarily small, hence μ1 is not Lipschitz continuous. In fact μ1 is discontinuous at (x, t). In the following lemma, we show that inlet points where a is tangent to ∂$$\varOmega$$ are the only troublemakers. Fig. 1. View largeDownload slide Left: proof of Lemma 4.1. Right: vortex touching ∂$$\varOmega$$− with $$t_{0}\gg{\tilde{t}}_{0}=0$$, i.e. $$t_{0}-{\tilde{t}}_{0}$$ is large, hence μ1 is not Lipschitz continuous at x. Fig. 1. View largeDownload slide Left: proof of Lemma 4.1. Right: vortex touching ∂$$\varOmega$$− with $$t_{0}\gg{\tilde{t}}_{0}=0$$, i.e. $$t_{0}-{\tilde{t}}_{0}$$ is large, hence μ1 is not Lipschitz continuous at x. We note that since $$\varOmega$$ is a bounded Lipschitz domain (hence quasiconvex), when proving Lipschitz continuity of μ1 in space it is sufficient to prove local Lipschitz continuity of μ1 in some neighborhood of each $$x\in \overline{\varOmega }$$ with a Lipschitz constant independent of x. Lemma 4.1 Let $$a\in L^{\infty }(Q_{T})$$ be continuous with respect to time and Lipschitz continuous with respect to space. Let there exist a constant $$a_{\min }>0$$ such that   \begin{align} -a(x,t)\cdot\mathbf{n}\geqslant a_{\min} \end{align} (4.8)for all x ∈ ∂$$\varOmega$$−, t ∈ [0, T). Let μ1 be defined by (4.3) on $$\overline{\varOmega }\times [0,T)$$. Let the time any particle carried by the flow field a(⋅, ⋅) spends in $$\varOmega$$ be uniformly bounded by $$\widehat{T}$$. Then μ1 satisfies assumption (4.7). Proof. By definition, we have $$\mu _{1}(x,t)\geqslant 0$$ for all x, t. By the above considerations, μ1 is bounded by the maximal time particles spend in $$\varOmega$$, which is uniformly bounded. This implies $$\mu _{1}(x,t)\leqslant \widehat{T}=\mu _{\max }$$ for some $$\mu _{\max }$$. Now we will prove Lipschitz continuity. Let t ∈ (0, T) be fixed and let $$x,{\tilde{x}}\in \overline{\varOmega }$$ such that $$|x-{\tilde{x}}|\leqslant \varepsilon $$, where ε will be chosen sufficiently small in the following. Due to the assumptions on a, the pathlines passing through x and $${\tilde{x}}$$ are uniquely determined and originate at some (x0, t0) and $$({\tilde{x}}_{0},{\tilde{t}}_{0})$$, respectively. In other words,   \begin{align} x=S(t;x_{0},t_{0}),\quad{\tilde{x}}=S\left(t;{\tilde{x}}_{0},{\tilde{t}}_{0}\right)\!. \nonumber \end{align}Without loss of generality, let $$t_{0}\ge{\tilde{t}}_{0}>0$$, hence $$x_{0},{\tilde{x}}_{0}\in \partial \varOmega ^{-}$$. The case when $${\tilde{t}}_{0}=0$$ can be treated similarly. Furthermore, we assume that x0 is not a vertex of ∂$$\varOmega$$−—this case will be treated at the end of the proof. Since a is uniformly bounded and Lipschitz continuous in space, there exists δ, such that if t ∈ [0, T) and $$\mathrm{dist}\{x,\partial \varOmega ^{-}\}\leqslant \delta $$ then   \begin{align} -a(x,t)\cdot\mathbf{n}\geqslant a_{\min}/2, \end{align} (4.9)due to (4.8). If ε is sufficiently small, then for the distance of the two considered pathlines at time t0 we have $$ |x_{0}-S (t_{0};{\tilde{x}}_{0},{\tilde{t}}_{0})| =|S(t_{0};x_{0},t_{0} )-S(t_{0};{\tilde{x}}_{0},{\tilde{t}}_{0})|\leqslant \delta $$. Since x0 ∈ ∂$$\varOmega$$− this means that $$S(t_{0};{\tilde{x}}_{0},{\tilde{t}}_{0})$$ is in the δ-neighborhood of ∂$$\varOmega$$− and by (4.9), $$\mathrm{dist}\{S(\vartheta ;\tilde{x}_{0},\tilde{t}_{0}),\partial \varOmega ^{-}\}$$ decreases as ϑ goes from t0 to $${\tilde{t}}_{0}$$ with a rate of at least $$a_{\min }/2$$ due to the uniformity of the bound (4.9). Therefore, $$S(\vartheta ;{\tilde{x}}_{0},{\tilde{t}}_{0})$$ stays in the δ-neighborhood of ∂$$\varOmega$$− for all $$\vartheta \in \left [{\tilde{t}}_{0},t_{0}\right ]$$ and   \begin{align} -a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\cdot\mathbf{n}\geqslant a_{\min}/2 \nonumber \end{align}for all $$\vartheta \in \left [{\tilde{t}}_{0},t_{0}\right ]$$. Moreover, since x0 lies in the interior of an edge of ∂$$\varOmega$$−, by choosing ε sufficiently small, we can ensure that $${\tilde{x}}_{0}$$ also lies on this edge (face). Now we estimate $$\left |\mu _{1}(x,t)-\mu _{1}\left ({\tilde{x}},t\right )\right |=\left |t_{0}-{\tilde{t}}_{0}\right |$$. We have   \begin{align*} x-x_{0} =&\,\int_{t_{0}}^{t}\frac{\mathrm{d} S}{\mathrm{d} t}(\vartheta;x_{0},t_{0})\mathrm{\ d}\vartheta =\int_{t_{0}}^{t} a\left(S(\vartheta;x_{0},t_{0}),\vartheta\right)\mathrm{\ d}\vartheta,\\{\tilde{x}}-{\tilde{x}}_{0} =&\,\int_{{\tilde{t}}_{0}}^{t}\frac{\mathrm{d} S}{\mathrm{d} t}\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right)\mathrm{\ d}\vartheta =\int_{{\tilde{t}}_{0}}^{t} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\mathrm{d}\vartheta. \end{align*}Subtracting these two identities gives us   \begin{align} \begin{split} x_{0}-{\tilde{x}}_{0} =x-{\tilde{x}} +\int_{{\tilde{t}}_{0}}^{t_{0}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\mathrm{\ d}\vartheta +\int_{t_{0}}^{t}a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right) -a\left(S(\vartheta;x_{0},t_{0}),\vartheta\right)\mathrm{d}\vartheta. \end{split} \end{align} (4.10)If we consider n, the normal to ∂$$\varOmega$$− at x0, we see that $$(x_{0}-{\tilde{x}}_{0})\cdot \mathbf{n}=0$$, as both x0 and $${\tilde{x}}_{0}$$ lie on the same edge on ∂$$\varOmega$$−. Therefore, if we multiply (4.10) by −n, we get   \begin{align} 0 =&\,(x-{\tilde{x}})\cdotp\!(-\mathbf{n}) +\int_{{\tilde{t}}_{0}}^{t_{0}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\cdotp\!(-\mathbf{n})\mathrm{\ d}\vartheta\nonumber\\ &+\int_{t_{0}}^{t}\left(a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right) -a\left(S\left(\vartheta;x_{0},t_{0}\right),\vartheta\right)\right)\cdotp\!(-\mathbf{n})\mathrm{\ d}\vartheta. \end{align} (4.11)Due to (4.9), we can estimate the first integral as   \begin{align} \int_{{\tilde{t}}_{0}}^{t_{0}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\cdotp\!(-\mathbf{n})\,\,\mathrm{d}\vartheta \geqslant \frac{a_{\min}}{2}\left(t_{0}-{\tilde{t}}_{0}\right). \nonumber \end{align}As for the second integral in (4.11), we have   \begin{align}&\left|\int_{t_{0}}^{t}\big(a\big( S \big(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\big),\vartheta \big) -a\left(S\left(\vartheta;x_{0},t_{0}\right),\vartheta\right)\big)\cdot \!(-\mathbf{n}) \mathrm{\ d}\vartheta \right|\nonumber\\ &\quad\leqslant \widehat{T}L_{a} \sup_{\vartheta\in(t_{0},t)}\left|S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right) -S\left(\vartheta;x_{0},t_{0}\right)\right|\leqslant \widehat{T}L_{a} \left(e^{\widehat{T}L_{a}}-1\right)\|x-{\tilde{x}}\|, \end{align} (4.12)where La is the Lipschitz constant of a with respect to x. The last inequality in (4.12) follows from standard results on ordinary differential equations, namely continuous dependence of the solution on the initial condition—here we consider the ordinary differential equations defining S(⋅;⋅, ⋅) backward in time on the interval t0, t with ‘initial’ conditions x and $${\tilde{x}}$$ at time t. Finally, $$|\left (x-{\tilde{x}}\right )\cdotp \!(-\mathbf{n})|\leqslant \|x-{\tilde{x}}\|$$. Therefore, we get from (4.11),   \begin{align} \frac{a_{\min}}{2}\left|\mu_{1}(x,t)-\mu_{1}\left({\tilde{x}},t\right)\right| =\frac{a_{\min}}{2}\left|t_{0}-{\tilde{t}}_{0}\right| =\frac{a_{\min}}{2}\left(t_{0}-{\tilde{t}}_{0}\right) \leqslant \left(1+\widehat{T}L_{a} \left(e^{\widehat{T}L_{a}}-1\right)\right)\|x-{\tilde{x}}\|. \end{align} (4.13)Dividing by $$a_{\min }/2>0$$ gives Lipschitz continuity of μ1(⋅, t). Now we return to the case when x0 is a vertex of ∂$$\varOmega$$−. Reasoning as above, by choosing ε sufficiently small, we can ensure that $${\tilde{x}}_{0}\in \partial \varOmega ^{-}$$ is sufficiently close to x0, i.e $${\tilde{x}}_{0}$$ lies on one of the edges adjoining x0. Then we can again multiply (4.11) by −n, the normal to ∂$$\varOmega$$− at $${\tilde{x}}_{0}$$. Hence, $$\left ({\tilde{x}}-{\tilde{x}}_{0}\right )\cdotp \!(\mathbf{n})=0$$ will also be satisfied and we can proceed as in the previous case. Remark 4.2 We note that the Lipschitz constant $$L_{\mu _{1}}$$ of μ1 can be estimated from (4.13) by   \begin{align} L_{\mu_{1}}\leqslant \frac{2}{a_{\min}}\left(1+\widehat{T}L_{a} \left(e^{\widehat{T}L_{a}}-1\right)\right)\!. \nonumber \end{align}The Lipschitz constant Lμ of μ can then be obtained by multiplying by the constant factor from (4.4). Having established uniform Lipschitz continuity of μ1 in space, we can prove its Lipschitz continuity with respect to time. This implies the differentiability of μ1 with respect to x and t a.e. in QT; therefore the left-hand side of (4.1) is well defined a.e. in QT and this expression is equal to the derivative of μ1 along pathlines. This is important, as all the following considerations are thus justified. Lemma 4.3 Let a satisfy the assumptions of Lemma 4.1. Then μ1 is uniformly Lipschitz continuous with respect to time t: there exists Lt ⩾ 0 such that   \begin{align} \left|\mu(x,t)-\mu\left(x,{\tilde{t}}\right)\right|\leqslant L_{t}\left|t-{\tilde{t}}\right| \nonumber \end{align}for all x ∈ $$\varOmega$$ and $$t,{\tilde{t}}\in [0,T)$$. Proof. Let e.g. $${\tilde{t}}>t$$ and let (x0, t0) and $$\left ({\tilde{x}}_{0},{\tilde{t}}_{0}\right )$$ denote the origin of the pathlines passing through (x, t) and $$\left (x,{\tilde{t}}\right )$$, respectively. Therefore $$x=S(t;x_{0},t_{0}) =S\left ({\tilde{t}};{\tilde{x}}_{0},{\tilde{t}}_{0}\right )$$. Denote $${\tilde{x}}=S\left (t;{\tilde{x}}_{0},{\tilde{t}}_{0}\right )$$. Then   \begin{align} |x-{\tilde{x}}| =\left|\int_{t}^{{\tilde{t}}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right) \mathrm{\ d}\vartheta\right| \leqslant \|a\|_{L^{\infty}(Q_{T})}|t-{\tilde{t}}|. \end{align} (4.14)Therefore, by (4.3),   \begin{align*}\left|\mu(x,t)-\mu\left(x,{\tilde{t}}\right)\right| \leqslant&\, |t-{\tilde{t}}|+|{\tilde{t}}_{0}-t_{0}| =\left|t-{\tilde{t}}\right|+\left|\mu(x,t)-\mu({\tilde{x}},t)\right| \leqslant \left|t-{\tilde{t}}\right|+L_{\mu}\left|x-{\tilde{x}}\right|\\ \leqslant&\, \left(1+L_{\mu} \|a\|_{L^{\infty}(Q_{T})}\right)\left|t-{\tilde{t}}\right| \end{align*}due to (4.14) and Lemma 4.1. This completes the proof. Example 4.4 In simple cases, the function μ1 can be explicitly written down. As a trivial example, we take the one-dimensional stationary flow field a(x, t) = x + 1 on $$\varOmega$$ = (0, 1) and the time interval $$(0,+\infty )$$. Then (4.2) can be easily solved to obtain   \begin{align} S(t;x_{0},t_{0})=(x_{0}+1)e^{t-t_{0}}-1. \nonumber \end{align}The (x,t)--plane is then separated into two regions separated by the pathline S(t;0, 0), which is the curve x = et − 1, i.e. $$t=\ln (x+1)$$. For points (x, t) beneath this curve, i.e. $$t\leqslant \ln (x+1)$$, we have t0 = 0, hence μ1(x, t) = t − t0 = t. For points (x, t) above the separation curve, we have x0 = 0, hence   \begin{align} x=S(t;0,t_{0})=e^{t-t_{0}}-1 \quad\Longrightarrow \quad \mu_{1}(x,t)=t-t_{0}=\ln(x+1). \nonumber \end{align}Altogether, we have   \begin{align} \mu_{1}(x,t)=\begin{cases} t & \text{if {$t\leqslant\ln(x+1)$},} \\ \ln(x+1)& \text{otherwise.} \end{cases} \nonumber \end{align}This function is globally bounded and globally Lipschitz continuous. We note that the standard exponential scaling trick which gives exponential growth of the error corresponds to taking μ1(x, t) = t for all x, t, which is an unbounded function. 5. DG method Having introduced the necessary tools for the analysis using the general exponential scaling argument, we can finally proceed to the DG method itself. Let $${\mathscr{T}}_{h}$$ be a triangulation of $$\varOmega$$, i.e. a partition of $$\overline{\varOmega }$$ into a finite number of closed simplices with mutually disjoint interiors. As usual with DG, $${\mathscr{T}}_{h}$$ need not be conforming, i.e. hanging nodes are allowed. For $$K\in{\mathscr{T}}_{h}$$ we set $$h_{K}=\textrm{diam}(K),\,\, h=\max _{K\in{\mathscr{T}}_{h}}h_{K}$$. For each $$K\in{\mathscr{T}}_{h}$$ we define its inflow and outflow boundaries by   \begin{align*} \partial K^{-}(t)&=\left\{x\in\partial K; a(x,t)\cdotp\mathbf{n}(x)<0\right\}\!,\\ \partial K^{+}(t)&=\{x\in\partial K; a(x,t)\cdotp\mathbf{n}(x)\ge 0\}, \nonumber \end{align*}where n(x) is the unit outer normal to ∂K. To simplify the notation, we shall usually omit the argument t in the following and write simply ∂K±. Here we remark that $${\mathscr{T}}_{h}$$ need not conform to ∂$$\varOmega$$−. In other words, the intersection of ∂$$\varOmega$$− and ∂K need not be an entire face (edge) of K. The triangulation $${\mathscr{T}}_{h}$$ defines the broken Sobolev space  \begin{align} H^{1}({\mathscr{T}}_{h})=\left\{v\in L^{2}(\varOmega);\, v|_{K} \in H^{1}(K)\quad \forall\, K\in{\mathscr{T}}_{h} \right\}\!. \nonumber \end{align}The approximate solution will be sought in the space of discontinuous piecewise polynomial functions   \begin{align} S_{h} = \left\{v_{h};\, v_{h}|_{K} \in P^{p}(K)\quad \forall\, K\in{\mathscr{T}}_{h}\right\}\!, \nonumber \end{align}where Pp(K) is the space of all polynomials on K of degree at most $$p\in \mathbb{N}$$. Given an element $$K\in{\mathscr{T}}_{h}$$, for vh ∈ Sh we define $$v_{h}^{-}$$ as the trace of vh on ∂K from the side of the element adjacent to K. Furthermore, on ∂K ∖ ∂$$\varOmega$$ we define the jump of vh as $$[v_{h}]=v_{h}-v_{h}^{-}$$, where vh is the trace from inside K. The DG formulation of (2.1) then reads as follows: we seek $$u_{h}\in C^{1}\left ([0,T);S_{h}\right )$$ such that $$u_{h}(0)={u_{h}^{0}}$$, an Sh-approximation of the initial condition u0, and for all t ∈ (0, T),   \begin{align} \Bigg(\frac{\partial u_{h}}{\partial t},v_{h}\Bigg)+b_{h}(u_{h},v_{h}) +c_{h}(u_{h},v_{h})=l_{h}(v_{h}) \quad\forall\, v_{h}\in S_{h}. \end{align} (5.1)Here bh, ch and lh are bilinear and linear forms, respectively, defined for $$u,v\in H^{1}({\mathscr{T}}_{h})$$ as follows. The bilinear advection form bh is   \begin{align} b_{h}(u,v)= \sum_{K\in{\mathscr{T}}_{h}}\int_{K} (a\cdotp \nabla\! u) v\,\mathrm{d} x -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[u]v\,\mathrm{d} S -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})uv\,\mathrm{d} S, \nonumber \end{align}the reaction form is defined by   \begin{align} c_{h}(u,v)= \int_{\varOmega} cu v\,\mathrm{d} x \nonumber \end{align}and lh is the right-hand side form  \begin{align} l_{h}(v)= -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})u_{D}v\,\mathrm{d} x . \nonumber \end{align}The definition of bh corresponds to the concept of upwinding which has already been applied in the context of DG in the original paper by Reed & Hill (1973). The right-hand side form lh is based on the weak enforcement of the Dirichlet boundary condition on ∂$$\varOmega$$− by penalization. For a detailed derivation of the specific forms of bh and lh as used in the paper, we refer e.g. to the study by Feistauer & Švadlenka (2004). 6. Analysis of the advection and reaction terms In this section, we prove the estimates of the bilinear forms bh and ch necessary for the error analysis. First, we need some standard results from approximation theory. 6.1 Auxiliary results In the following analysis we assume that the exact solution u is sufficiently regular, namely   \begin{align} u,u_{t}:=\frac{\partial u}{\partial t}\in L^{2}\left(H^{p+1}\right)\!. \nonumber \end{align}We consider a system $$\{{\mathscr{T}}_{h}\}_{h\in (0,\,h_{0})}$$, h0 > 0 of nonconforming triangulations of $$\varOmega$$ that are shape regular and satisfy the inverse assumption, cf. Ciarlet (1980). Under these assumptions we have the following standard results: Lemma 6.1 (Inverse inequality). There exists a constant CI > 0 independent of h, K such that for all $$K\in{\mathscr{T}}_{h}$$, and all v ∈ Pp(K),   \begin{align} |v|_{H^{1}(K)}\leqslant C_{I} h_{K}^{-1} \|v\|_{L^{2}(K)}. \nonumber \end{align} For v ∈ L2($$\varOmega$$) we denote by Πhv ∈ Sh the L2($$\varOmega$$)-projection of v onto Sh:   \begin{align} \left(\Pi_{h} v-v,\,\varphi_{h}\right)=0 \quad \forall\,\varphi_{h}\in S_{h}. \end{align} (6.1)Let ηh(t) = u(t) −Πhu(t) and ξh(t) =Πhu(t) − uh(t) ∈ Sh. Then we can write the error of the method as eh(t) := u(t) − uh(t) = ηh(t) + ξh(t). For simplicity, we shall usually drop the subscript h. We have the following standard approximation result; cf. Ciarlet (1980). Lemma 6.2 There exists a constant CΠ > 0 independent of h, K such that for all h ∈ (0, h0),   \begin{align*}\|\eta(t)\|_{L^{2}(K)}\leqslant&\, C_{\Pi} h^{p+1}|u(t)|_{{H^{p+1}}(K)},\\ |\eta(t)|_{H^{1}(K)}\leqslant&\, C_{\Pi} h^{p}|u(t)|_{{H^{p+1}}(K)},\\ \|\eta(t)\|_{L^{2}(\partial K)}\leqslant&\, C_{\Pi} h^{p+1/2}|u(t)|_{{H^{p+1}}(K)},\\ \Bigg\|\frac{\partial\eta(t)}{\partial t}\Bigg\|_{L^{2}(K)}\leqslant&\, C_{\Pi} h^{p+1}|u_{t}(t)|_{{H^{p+1}}(K)}. \end{align*} Throughout the paper, C will be a generic constant independent of h, t, T. In order to track the dependence of the estimates on the function μ, we will also assume that the generic constant C is independent of μ, particularly $$\mu _{\max }$$ and Lμ from (4.7), and we will explicitly track these dependencies. The dependence of the resulting estimates of Theorem 7.1 on μ will then be clarified in Section 7 for the specific construction of μ from Section 4. 6.2 Estimates of the advection and reaction terms In our analysis we will assume there exists a constant γ0 and a function $$\mu :Q_{T}\to \mathbb{R}$$ such that the ellipticity assumption (3.9) holds. Furthermore, we assume that μ satisfies (4.7), i.e. μ is non-negative, bounded and Lipschitz continuous in space. Such a function μ was constructed in Section 4. Similarly to Section 2, we wish to write $$\xi (x,t)=e^{\mu (x,t)}\tilde{\xi }(x,t)$$ and test the error equation with $$\phi (x,t)=e^{-\mu (x,t)}\tilde{\xi }(x,t)=e^{-2\mu (x,t)}\xi (x,t)$$ to obtain estimates for $$\tilde{\xi }$$. However, since ϕ(t) ∉ Sh this is not possible. One possibility is to test by Πhϕ(t) ∈ Sh and estimate the resulting difference Πhϕ(t) − ϕ(t). This is done in the stationary case in the study by Ayuso & Marini (2009) and the analysis is carried out under the assumption $$\mu \in W^{p+1,\infty }(\varOmega )$$. Such high regularity can be achieved by mollification of μ. However, this would be somewhat technical in the evolutionary case, as space–time smoothing would be required in which case the dependence of all constants on T must be carefully considered. Also QT is potentially an unbounded domain (for $$T=+\infty $$) which leads to technical difficulties. Here we carry out the analysis under the weaker assumption (4.7), i.e. $$\mu (t)\in W^{1,\infty }(\varOmega )$$. Lemma 6.3 Let μ satisfy assumptions (4.7). Let $$\phi (x,t)=e^{-\mu (x,t)}\tilde{\xi }(x,t)=e^{-2\mu (x,t)}\xi (x,t)$$, where ξ(t) ∈ Sh. Then there exists C independent of $$h,t,\xi ,\tilde{\xi }$$ and μ such that   \begin{align} \left\|\Pi_{h}\phi(t)-\phi(t)\right\|_{L^{2}(K)}\leqslant&\, CL_{\mu} e^{h_{K}L_{\mu}} h_{K} \max_{x\in K}e^{-\mu(x,t)} \|\tilde{\xi}(t)\|_{L^{2}(K)},\nonumber\\ \left\|\Pi_{h}\phi(t)-\phi(t)\right\|_{L^{2}(\partial K)}\leqslant&\, CL_{\mu} e^{h_{K}L_{\mu}} h^{1/2}_{K} \max_{x\in K}e^{-\mu(x,t)} \|\tilde{\xi}(t)\|_{L^{2}(K)}. \end{align} (6.2) Proof. Let xK be the centroid of K. On element K we introduce the constant μK(t) = μ(xK, t); then the function $$e^{-2\mu _{K}(t)}\xi (\cdot ,t)$$ lies in Pp(K), hence is fixed by the projection Πh. Therefore,   \begin{align*}\Pi_{h}\phi(t)-\phi(t)=&\, \Pi_{h}\left(e^{-2\mu(t)}\xi(t) -e^{-2\mu_{K}(t)}\xi(t)\right) -\left(e^{-2\mu(t)}\xi(t) -e^{-2\mu_{K}(t)}\xi(t)\right)\\ =&\,\Pi_{h} w(t)-w(t), \end{align*}where $$w(x,t)=(e^{-2\mu (x,t)} -e^{-2\mu (x_{K},t)})\xi (x,t)$$. Standard estimates of the interpolation error of Πh give   \begin{align} \left\|\Pi_{h}\phi(t)-\phi(t)\right\|_{L^{2}(K)}^{2} =\left\|\Pi_{h} w(t)-w(t)\right\|_{L^{2}(K)}^{2}\leqslant C{h_{K}^{2}}\left|w(t)\right|_{H^{1}(K)}^{2}. \end{align} (6.3)For the right-hand side seminorm we have   \begin{align*} |w(t)|_{H^{1}(K)}^{2}=&\,\int_{K}\left|\nabla\!\left(\left(e^{-2\mu(x,t)} -e^{-2\mu(x_{K},t)}\right)\xi(x,t)\right)\right|{}^{2}\,\mathrm{d} x \\ \leqslant&\, 2\int_{K}\left|\nabla\! e^{-2\mu(x,t)}\xi(x,t)\right| \mathrm{d} x +2\int_{K}\left|\left(e^{-2\mu(x,t)} -e^{-2\mu(x_{K},t)}\right)\nabla\!\xi(x,t)\right|{}^{2}\,\mathrm{d} x.\end{align*}If $$\mu (t)\in C^{1}(\overline{\varOmega })$$, by the mean value theorem   \begin{align} \left|e^{-2\mu(x_{K},t)}-e^{-2\mu(x,t)}\right|\leqslant h_{K} \left|\nabla\! e^{-2\mu(\zeta,t)}\right| = h_{K} e^{-2\mu(\zeta,t)}2\left|\nabla\mu(\zeta,t)\right| \nonumber \end{align}for some point ζ on the line between x and xK. Therefore, by the inverse inequality,   \begin{align*}|w(t)|_{H^{1}(K)}^{2}\leqslant&\, 2\int_{K} 4e^{-4\mu(x,t)} \left|\nabla\! \mu(x,t)\right|{}^{2}\left|\xi(x,t)\right|{}^{2}\,\mathrm{d} x +2\int_{K} {h_{K}^{2}}e^{-4\mu(\zeta,t)}4|\nabla\mu(\zeta,t)|^{2} |\nabla\!\xi(x,t)|^{2}\,\mathrm{d} x \\ \leqslant&\, 8 |\mu(t)|^{2}_{W^{1,\infty}}\int_{K} e^{-4\mu(x,t)}e^{2\mu(x,t)} |\tilde{\xi}(x,t)|^{2}\,\mathrm{d} x +8{C_{I}^{2}} \max_{x\in K}e^{-4\mu(x,t)}\left|\mu(t)\right|^{2}_{W^{1,\infty}} \|\xi(t)\|^{2}\\ \leqslant&\, 8 L_{\mu}^{2} \max_{x\in K} e^{-2\mu(x,t)} e^{2h_{K}L_{\mu}} \|\tilde{\xi}(t)\|^{2} +8{C_{I}^{2}} \max_{x\in K} e^{-2\mu(x,t)} e^{2h_{K}L_{\mu}} L_{\mu}^{2} \|\tilde{\xi}(t)\|^{2}.\end{align*}Substituting into (6.3) gives us the first inequality in (6.2) for $$\mu (t)\in C^{1}(\overline{\varOmega })$$. The case $$\mu (t)\in W^{1,\infty }(\varOmega )$$ follows by a density argument. The second inequality in (6.2) can be obtained similarly, only intermediately applying the trace inequality $$\|\xi (t)\|_{L^{2}(\partial K)}\leqslant Ch_{K}^{-1/2}\|\xi (t)\|_{L^{2}(K)}$$. Remark 6.4 The estimates of Lemma 6.3 are effectively $$\mathcal{O}(h_{K})\|\tilde{\xi }\|$$ and $$\mathcal{O}(h_{K}^{1/2})\|\tilde{\xi }\|$$ estimates, respectively. We note that the dependence of the estimates on Lμ is rather mild, effectively linear, due to the small factor hK in the exponent. Now we shall estimate individual terms in the DG formulation. Due to the consistency of the DG scheme, the exact solution u also satisfies (5.1). We subtract the formulations for u and uh to obtain the error equation   \begin{align} \begin{split} \left(\frac{\partial \xi}{\partial t},v_{h}\right) +\left(\frac{\partial \eta}{\partial t},v_{h}\right) + b_{h}(\xi,v_{h}) +b_{h}(\eta,v_{h}) + c_{h}(\xi,v_{h}) +c_{h}(\eta,v_{h}) =0 \end{split} \end{align} (6.4)for all vh ∈ Sh. As stated earlier we want to test (6.4) by $$\phi (x,t)=e^{-\mu (x,t)}\tilde{\xi }(x,t)$$; however, ϕ(t) ∉ Sh. We therefore set vh =Πhϕ(t) and estimate the difference using Lemma 6.3. We write (6.4) as   \begin{align}&\left(\frac{\partial \xi}{\partial t},\Pi_{h}\phi\right) +b_{h}(\xi,\phi) +b_{h}(\xi,\Pi_{h}\phi-\phi) +b_{h}(\eta,\phi) +b_{h}(\eta,\Pi_{h}\phi-\phi)\nonumber\\ &\qquad+\ c_{h}(\xi,\phi) +c_{h}(\xi,\Pi_{h}\phi-\phi) +c_{h}(\eta,\Pi_{h}\phi) +\left(\frac{\partial \eta}{\partial t},\Pi_{h}\phi\right)=0. \end{align} (6.5)We will estimate the individual terms of (6.5) in a series of lemmas. For this purpose, we introduce the following norm on a subset ω of ∂K or ∂$$\varOmega$$:   \begin{align} \|\,f\|_{a,\omega}=\left\|\sqrt{|a\cdotp\mathbf{n}|}\,f\right\|_{L^{2}(\omega)}, \nonumber \end{align}where n is the outer normal to ∂K or ∂$$\varOmega$$. We will usually omit the argument t to simplify the notation. In the following lemma we show that selected terms from (6.5) indeed give us the desired ellipticity similarly to (3.8) and (3.9) plus additional elliptic terms due to the use of the upwind flux. Lemma 6.5 Let $$\xi =e^{\mu }\tilde{\xi }, \phi =e^{-\mu }\tilde{\xi }$$ as above and let μ satisfy assumptions (3.9) and (4.7). Then   \begin{equation} \Bigg(\frac{\partial \xi}{\partial t},\Pi_{h}\phi\Bigg) +b_{h}(\xi,\phi)+c_{h}(\xi,\phi)\ge\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}\|^{2}+\gamma_{0}\|\tilde{\xi}\|^{2} +\frac{1}{2}\sum_{K\in{\mathscr{T}}_{h}}\bigg(\big\|[\tilde{\xi}]\big\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\big\|\tilde{\xi}\big \|_{a,\partial K\cap\partial\varOmega}^{2}\bigg). \end{equation} (6.6) Proof. Since ∂ξ/∂t ∈ Sh for each t, by definition (6.1) of Πh we have   \begin{align} \Bigg(\frac{\partial \xi}{\partial t},\Pi_{h}\phi\Bigg)= \Bigg(\frac{\partial \xi}{\partial t},\phi\Bigg) =\Bigg(e^{\mu}\frac{\partial \tilde{\xi}}{\partial t} +e^{\mu}\frac{\partial \mu}{\partial t}\tilde{\xi},e^{-\mu} \tilde{\xi}\Bigg) =\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}\|^{2} +\Bigg(\frac{\partial \mu}{\partial t}\tilde{\xi},\tilde{\xi}\Bigg). \end{align} (6.7)The reactive term satisfies   \begin{align} c_{h}(\xi,\phi)=\int_{\varOmega} c\xi\phi\,\mathrm{d} x =\int_{\varOmega} ce^{\mu}\tilde{\xi} e^{-\mu}\tilde{\xi}\,\mathrm{d} x =\int_{\varOmega} c\tilde{\xi}^{2}\,\mathrm{d} x . \end{align} (6.8)From the definition of bh, we get   \begin{align} b_{h}(\xi,\phi)=&\,\!\! \sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp(\nabla\!\mu\,\tilde{\xi}+\nabla\!\tilde{\xi})e^{\mu} e^{-\mu}\tilde{\xi}\,\mathrm{d} x -\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[e^{\mu}\tilde{\xi}]e^{-\mu}\tilde{\xi}\,\mathrm{d} S \nonumber\\ &-\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})e^{\mu}\tilde{\xi} e^{-\mu}\tilde{\xi}\,\mathrm{d} S \nonumber\\ =&\,\!\! \sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp(\nabla\!\mu\,\tilde{\xi}+\nabla\!\tilde{\xi})\tilde{\xi}\,\mathrm{d} x -\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\tilde{\xi}]\tilde{\xi}\,\mathrm{d} S -\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S .\end{align} (6.9)By Green’s theorem,   \begin{align} \int_{K} a\cdotp\!\nabla\!\tilde{\xi}\,\tilde{\xi}\,\mathrm{d} x =-\frac{1}{2}\int_{K} \operatorname{div} a\,\tilde{\xi}^{2}\,\mathrm{d} x +\frac{1}{2}\int_{\partial K}(a\cdotp\mathbf{n})\,\tilde{\xi}^{2}\,\mathrm{d} S . \end{align} (6.10)Splitting the last integral over the separate parts of ∂K, i.e. ∂K− ∖ ∂$$\varOmega$$, ∂K− ∩ $$ \partial \varOmega $$, ∂K+ ∖ ∂$$\varOmega$$ and ∂K+ ∩ $$ \partial \varOmega $$, and by substituting (6.10) into (6.9), we get   \begin{align}b_{h}(\xi,\phi)=&\,\!\! \sum_{K\in{\mathscr{T}}_{h}}\int_{K} \left(a\cdotp\nabla\!\mu-\tfrac{1}{2}\operatorname{div} a\right)\tilde{\xi}^{2}\,\mathrm{d} x \nonumber\\ &+\!\!\sum_{K\in{\mathscr{T}}_{h}}\bigg(-\frac{1}{2}\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})\left(\tilde{\xi}^{2}-2\tilde{\xi}\tilde{\xi}^{-}\right)\,\mathrm{d} S -\frac{1}{2}\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S \nonumber\\ &\qquad\qquad+\frac{1}{2}\int_{\partial K^{+}\setminus\partial\varOmega}(a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S +\frac{1}{2}\int_{\partial K^{+}\cap\partial\varOmega}(a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S \bigg).\end{align} (6.11)We note that   \begin{align} \sum_{K\in{\mathscr{T}}_{h}}\int_{\partial K^{+}\setminus\partial\varOmega}(a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S =-\sum_{K\in{\mathscr{T}}_{h}}\int_{\partial K^{-}\setminus\partial\varOmega}(a\cdotp \mathbf{n})(\tilde{\xi}^{-})^{2}\,\mathrm{d} S . \end{align} (6.12)Using the facts that $$\tilde{\xi }^{2}-2\tilde{\xi }\tilde{\xi }^{-}+(\tilde{\xi }^{-})^{2}= [\tilde{\xi }]^{2}$$ and − a ⋅ n = |a⋅n| on ∂K− and a ⋅ n = |a⋅n| on ∂K+, we get (6.6) by substituting (6.12) into (6.11) and applying assumption (3.9) in the resulting interior terms along with (6.7), (6.8). The following series of lemmas deals with estimating the remaining terms in (6.5). These estimates more or less follow standard procedures with several important differences pertaining to the specific choice of test functions ϕ and Πhϕ − ϕ. These include the inability to apply the inverse inequality to $$\tilde{\xi }\notin S_{h}$$ but only to ξ ∈ Sh and Lemma 6.3. In order to keep the main body of the text less cluttered, we leave the proofs of the following lemmas to the Appendix. Lemma 6.6 The advection terms satisfy   \begin{align} \hskip-30pt|b_{h}(\xi,\Pi_{h}\phi-\phi)|\leqslant C(1+L_{\mu})^{2} e^{4hL_{\mu}} h\|\tilde{\xi}\|^{2}+\frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\bigg(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\bigg), \end{align} (6.13)  \begin{align} {\hskip30pt}|b_{h}(\eta,\phi)|\leqslant C(1+L_{\mu})^{2}h\|\tilde{\xi}\|^{2} + Ch^{2p+1}|u(t)|_{H^{p+1}}^{2} +\frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\!\!\!\left(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\right), \end{align} (6.14)  \begin{align} |b_{h}(\eta,\Pi_{h}\phi-\phi)|\leqslant CL_{\mu}^{2} e^{2hL_{\mu}} h\|\tilde{\xi}\| +Ch^{2p+1}|u(t)|_{H^{p+1}}^{2}. \hskip150pt\end{align} (6.15) Lemma 6.7 The reaction terms satisfy   \begin{align*}|c_{h}(\xi,\Pi_{h}\phi-\phi)|\leqslant&\, CL_{\mu} e^{2hL_{\mu}}h\|\tilde{\xi}\|^{2},\\ |c_{h}(\eta,\Pi_{h}\phi)|\leqslant&\, C\big(1+L_{\mu} e^{hL_{\mu}}\big)^{2} h\|\tilde{\xi}\|^{2} +Ch^{2p+1}|u(t)|_{H^{p+1}}^{2}.\end{align*} Lemma 6.8 The evolutionary term satisfies   \begin{align} \Bigg|\Bigg(\frac{\partial \eta}{\partial t},\Pi_{h}\phi\Bigg)\Bigg| \leq C\big(1+L_{\mu} e^{hL_{\mu}}\big)^{2} h\|\tilde{\xi}\|^{2} +Ch^{2p+1}|u_{t}(t)|_{H^{p+1}}^{2}. \nonumber \end{align} 7. Error analysis Finally, we come to the error analysis. The starting point is the error identity (6.5) to which we apply the derived estimates of its individual terms. Theorem 7.1 Let there exist a function $$\mu :Q_{T}\to [0,\mu _{\max }]$$ for some constant $$\mu _{\max }$$, such that $$\mu (t)\in W^{1,\infty }(\varOmega )$$ and let there exist a constant γ0 > 0 such that the coefficients of (2.1) satisfy $$\mu _{t} +a\cdotp \nabla \!\mu +c-\tfrac{1}{2}\textrm{div} a\ge \gamma _{0}>0$$ a.e. in QT. Let the initial condition $${u_{h}^{0}}$$ satisfy $$\left \|\Pi _{h} u^{0}-{u_{h}^{0}}\right \|\leqslant Ch^{p+1/2}\left |u^{0}\right |_{H^{p+1}}$$. Then there exists a constant C independent of μ, h and T such that for h sufficiently small the error of the DG scheme (5.1) satisfies   \begin{align}&\max_{t\in[0,T]}\|e_{h}(t)\| +\sqrt{\gamma_{0}}\|e_{h}\|_{L^{2}(Q_{T})} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|[e_{h}(\vartheta)]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|e_{h}(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ce^{\mu_{\max}} h^{p+1/2}\left(|u^{0}|{}_{H^{p+1}}+|u|_{L^{2}\left({H^{p+1}}\right)} +\left|u_{t}\right|{}_{L^{2}\left({H^{p+1}}\right)} +h^{1/2}|u|_{L^{\infty}\left(H^{p+1}\right)}\right)\!.\end{align} (7.1) Proof. Applying Lemmas 6.5–6.8 to (6.5), multiplying by 2 for convenience and collecting similar terms, we get for all t,  \begin{align}&\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}(t)\|^{2} +2\gamma_{0}\|\tilde{\xi}(t)\|^{2} +\frac{1}{2}\sum_{K\in{\mathcal{T}}_{h}}\left(\|[\tilde{\xi}(t)]\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\tilde{\xi}(t)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\nonumber\\ &\quad\leqslant C_{1}(1+L_{\mu})^{2}e^{4hL_{\mu}}h\|\tilde{\xi}(t)\|^{2} +C_{2}h^{2p+1}\left(\left|u(t)\right|{}_{H^{p+1}}^{2}+|u_{t}(t)|_{H^{p+1}}^{2}\right)\!, \end{align} (7.2)where we have used $$C_{1} (1+L_{\mu })^{2}e^{4hL_{\mu }}$$ as a common upper bound for the constants from Lemmas 6.6–6.8 at the $$\mathcal{O}(h)\|\tilde{\xi }(t)\|^{2}$$ terms. Here C1 and C2 are independent of h, t, T, μ and $$\tilde{\xi }$$. Now if h is sufficiently small so that   \begin{align} C_{1}(1+L_{\mu})^{2}e^{4hL_{\mu}}h \leqslant\gamma_{0}, \end{align} (7.3)then the first right-hand side term can be absorbed by the left-hand side term $$2\gamma _{0}\|\tilde{\xi }(t)\|^{2}$$:   \begin{align*} &\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}(t)\|^{2} +\gamma_{0}\|\tilde{\xi}(t)\|^{2} +\frac{1}{2}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}(t)]\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\tilde{\xi}(t)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\\ &\quad\leqslant C_{2}h^{2p+1}\left(\left|u(t)\right|{}_{H^{p+1}}^{2}+|u_{t}(t)|_{H^{p+1}}^{2}\right)\!. \end{align*}Substituting t = ϑ and integrating over (0, t) give us   \begin{align}&\|\tilde{\xi}(t)\|^{2} +\gamma_{0}{\int_{0}^{t}}\|\tilde{\xi}(\vartheta)\|^{2}\mathrm{d} \vartheta +\frac{1}{2}{\int_{0}^{t}}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}(\vartheta)]\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\tilde{\xi}(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\nonumber\\ &\quad\leqslant \|\tilde{\xi}(0)\|^{2} +C_{2}h^{2p+1}\left(|u|_{L^{2}(0,t;{H^{p+1}})}^{2} +|u_{t}|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2}\right)\nonumber\\ &\quad\leqslant Ch^{2p+1}\left(|u^{0}|{}_{H^{p+1}}^{2} +|u|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2} +|u_{t}|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2}\right)\!, \end{align} (7.4)since the assumptions give   \begin{align} \|\tilde{\xi}(0)\|^{2}\leqslant \left\|\xi(0)\right\|^{2} =\|\Pi_{h} u^{0}-{u_{h}^{0}}\|^{2} \leqslant Ch^{2p+1}|u^{0}|{}_{H^{p+1}}^{2}. \nonumber \end{align}Now we reformulate estimate (7.4) as an estimate of ξ instead of $$\tilde{\xi }$$. Because $$\tilde{\xi }=e^{-\mu }\xi $$, we can estimate   \begin{align} \|\tilde{\xi}(t)\|^{2}\geqslant \min_{Q_{T}}e^{-2\mu(x,\vartheta)} \left\|\xi(t)\right\|^{2} =e^{-2\max_{Q_{T}}\mu(x,\vartheta)}\left\|\xi(t)\right\|^{2} =e^{-2\mu_{\max}} \|\xi(t)\|^{2} \nonumber \end{align}and similarly for the remaining left-hand side norms in (7.4). Substituting into (7.4) and multiplying by $$e^{2\mu _{\max }}$$ gives us   \begin{align*} &\|\xi(t)\|^{2} +\gamma_{0}\|\xi\|_{L^{2}\left(0,t;L^{2}\right)}^{2} +\frac{1}{2}{\int_{0}^{t}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|\left[\xi(\vartheta)\right]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\xi(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\nonumber\\ &\quad\leqslant Ch^{2p+1} e^{2\mu_{\max}}\left(|u^{0}|{}_{H^{p+1}}^{2} +|u|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2} +|u_{t}|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2}\right). \end{align*}By taking the square root of the last inequality and taking the maximum over all t ∈ [0, T], we get   \begin{align}&\max_{t\in[0,T]}\|\xi(t)\| +\sqrt{\gamma_{0}}\|\xi\|_{L^{2}(L^{2})} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\big\|[\xi(\vartheta)]\big\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\xi(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ch^{p+1/2} e^{\mu_{\max}}\left(\big|u^{0}\big|_{H^{p+1}} +|u|_{L^{2}({H^{p+1}})} +|u_{t}|_{L^{2}\left({H^{p+1}}\right)}\right) \end{align} (7.5)which is estimate (7.1) for the discrete part ξ of the error. Lemma 6.2 gives a similar inequality for η:   \begin{align} &\max_{t\in[0,T]}\|\eta(t)\| +\sqrt{\gamma_{0}}\|\eta\|_{L^{2}\left(L^{2}\right)} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|\left[\eta(\vartheta)\right]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\eta(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{\ d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ch^{p+1/2} \left(h^{1/2}|u|_{L^{\infty}\left(H^{p+1}\right)} +(h^{1/2}+1)|u|_{L^{2}\left({H^{p+1}}\right)}\right)\!. \end{align} (7.6)Since eh = ξ + η, the triangle inequality along with (7.5) and (7.6) gives us the statement of the theorem. We now interpret the results of Theorem 7.1 given the specific construction of the scaling function μ from Section 4. If the assumptions of Lemma 4.1 are satisfied, the function μ1 exists and has suitable regularity. If we choose γ0 (this can be chosen e.g. as γ0 = 1 for simplicity), then μ is constructed by the scaling (4.4). As μ1 is bounded by $$\widehat{T}$$, the exponential factor in (7.1) becomes   \begin{align} e^{\mu_{\max}}=e^{\widehat{T}(|\inf_{Q_{T}}(c-\frac{1}{2}\operatorname{div} a)^{-}|+\gamma_{0})}\!. \nonumber \end{align} Now we turn to the requirement ‘h is sufficiently small’ in Theorem 7.1. Specifically, (7.3) must be satisfied. If we denote $$A:=C_{1}\left (1+L_{\mu }\right )^{2}/\gamma _{0}$$ and B := 4Lμ, then condition (7.3) reads   \begin{align} Ae^{Bh}h \leqslant 1. \end{align} (7.7)It can be easily seen that if $$h\leqslant \tfrac{1}{2}\min \left \{A^{-1},B^{-1}\right \}$$, then   \begin{align} Ae^{Bh}h\leqslant \tfrac{1}{2}e^{1/2}\leqslant 1, \nonumber \end{align}i.e. (7.7) holds. Taking the definitions of A and B, we get the restriction on h,  \begin{align} h\leqslant \frac{1}{2}\min\left\{\frac{\gamma_{0}}{C_{1}\left(1+L_{\mu}\right)^{2}}, \frac{1}{4L_{\mu}}\right\}\!. \end{align} (7.8)If we take e.g. γ0 = 1 as mentioned and assume $$C_{1}\geqslant 1$$ (if not, C1 can be increased to 1 in (7.2)), condition (7.8) reduces to   \begin{align} h\leqslant \frac{1}{2C_{1}\left(1+L_{\mu}\right)^{2}}, \end{align} (7.9)where Lμ is the Lipschitz constant of μ from Lemma 4.1, which can be estimated as in Remark 4.2 using data from the continuous problem (2.1) such as La, the spatial Lipschitz constant of a. Therefore the condition (7.8) on the mesh size and the final error estimate can be formulated without μ, only in terms of $$\widehat{T}$$ and the data of the problem. Thus we can state the following corollary derived from Theorem 7.1 using the construction of μ from Section 4. Corollary 7.2 Let $$a\in L^{\infty }(Q_{T})$$ be continuous with respect to time and Lipschitz continuous with respect to space. Let there exist a constant $$a_{\min }>0$$ such that $$-a(x,t)\cdot \mathbf{n}\geqslant a_{\min }$$ for all x ∈ ∂$$\varOmega$$−, t ∈ [0, T). Let the time any particle carried by the flow field a(⋅, ⋅) spends in $$\varOmega$$ be uniformly bounded by $$\widehat{T}$$. Let the initial condition $${u_{h}^{0}}$$ satisfy $$\left \|\Pi _{h} u^{0}-{u_{h}^{0}}\right \|\leqslant Ch^{p+1/2}\left |u^{0}\right |_{H^{p+1}}$$. Then there exist constants c > 0 and C independent of h and T such that for $$h\leqslant c$$, the error of the DG scheme (5.1) satisfies   \begin{align}&\max_{t\in[0,T]}\left\|e_{h}(t)\right\| +\|e_{h}\|_{L^{2}(Q_{T})} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|[e_{h}(\vartheta)]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|e_{h}(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{\ d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ce^{\widehat{T}(|\inf_{Q_{T}}(c-\frac{1}{2}\operatorname{div} a)^{-}|+1)} h^{p+1/2}\left(|u^{0}|{}_{H^{p+1}}+|u|_{L^{2}\left({H^{p+1}}\right)} +|u_{t}|_{L^{2}\left({H^{p+1}}\right)} +h^{1/2}|u|_{L^{\infty}\left(H^{p+1}\right)}\right)\!. \end{align} (7.10) Remark 7.3 As mentioned in Section 4, estimate (7.10) is exponential in the maximal ‘lifetime’ of particles $$\widehat{T}$$, instead of the final time T, which would be expected when applying Gronwall’s lemma directly. Since we assume the maximal particle lifetime to be finite, the exponential factor in (7.10) remains uniformly bounded independent of time T, which can be infinite. 8. Conclusion and future work In this paper we derived a priori error estimates for a linear advection–reaction problem with inlet and outlet boundary conditions in the $$L^{\infty }\left (L^{2}\right )$$- and $$L^{2}\left (L^{2}\right )$$-norms. Unlike previous works, the analysis was performed without the usual ellipticity assumption $$c-\frac{1}{2}\textrm{div} a\geqslant 0$$. This is achieved by applying a general exponential scaling transformation in space and time to the exact and discrete solutions of the problem. We considered the case when the time spent by particles carried by the flow field inside the spatial domain $$\varOmega$$ is uniformly bounded by some $$\widehat{T}$$. The resulting error estimates are of the order Chp+1/2, where C depends exponentially on $$\widehat{T}$$ (which is a constant) and not on the final time $$T\leqslant +\infty $$, as would be expected from the use of Gronwall’s inequality. Effectively, due to the exponential scaling, we apply Gronwall’s lemma in the Lagrangian setting along pathlines, which exist only for time at most $$\widehat{T}$$, and not in the usual Eulerian sense. As for future work, we plan to extend the analysis to fully discrete DG schemes with discretization in time. Furthermore, we wish to extend the analysis to nonlinear convective problems, following the arguments of Zhang & Shu (2004) and Kučera (2014) to obtain error estimates without the exponential dependence on time of the error estimates. Funding The work of V. Kučera was supported by the J. William Fulbright Commission in the Czech Republic and research project No. 17-01747S of the Czech Science Foundation. The work of C.-W. Shu was supported by DOE grant DE-FG02-08ER25863 and NSF grants DMS-1418750 and DMS-1719410. References Ayuso, B. & Marini, L. D. ( 2009) Discontinuous Galerkin methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. , 47, 1391– 1420. Google Scholar CrossRef Search ADS   Benzoni-Gavage, S. & Serre, D. ( 2007) Multi-Dimensional Hyperbolic Partial Differential Equations: First-Order Systems and Applications. Oxford Mathematical Monographs . Oxford: Oxford University Press. Ciarlet, P. G. ( 1980) The finite element method for elliptic problems . Studies in Mathematics and its Applications . Amsterdam, New York: North-Holland. Devinatz, A., Ellis, R. & Friedman, A. ( 1974) The asymptotic behavior of the first real eigenvalue of second order elliptic operators with a small parameter in the highest derivatives, II. Indiana Univ. Math. J. , 23, 991– 1011. Google Scholar CrossRef Search ADS   Evans, L. C. ( 2010) Partial Differential Equations . Providence, RI: American Mathematical Society. Google Scholar CrossRef Search ADS   Feistauer, M. & Švadlenka, K. ( 2004) Discontinuous Galerkin method of lines for solving nonstationary singularly perturbed linear problems. J. Numer. Math. , 12, 97– 117. Google Scholar CrossRef Search ADS   Johnson, C. & Pitkäranta, J. ( 1986) An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Math. Comp. , 46, 1– 26. Google Scholar CrossRef Search ADS   Kučera, V. ( 2014) On diffusion-uniform error estimates for the DG method applied to singularly perturbed problems. IMA J. Numer. Anal. , 34, 820– 861. Google Scholar CrossRef Search ADS   Larsson, S. & Thomée, V. ( 2003) Partial Differential Equations with Numerical Methods. Texts in Applied Mathematics . Berlin: Springer. Lesaint, P. & Raviart, P. A. ( 1974) On a finite element method for solving the neutron transport equation. Mathematical Aspects of Finite Elements in Partial Differential Equations  ( C. de Boor ed). New York: Academic Press, pp. 89– 145. Google Scholar CrossRef Search ADS   Mizohata, S. ( 1973) The Theory of Partial Differential Equations . Cambridge: Cambridge University Press. Nävert, U. ( 1982) A finite element method for convection-diffusion problems. Ph.D. Thesis , Chalmers University of Technology. Reed, W. H. & Hill, T. ( 1973) Triangular mesh methods for the neutron transport equation. Los Alamos Report LA-UR-73–479. Roos, H.-G., Stynes, M. & Tobiska, L. ( 2008) Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems. Springer Series in Computational Mathematics , vol. 24. Berlin: Springer. Zhang, Q. & Shu, C.-W. ( 2004) Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws. SIAM J. Numer. Anal. , 42, 641– 666. Google Scholar CrossRef Search ADS   Appendix A.1 Proof of Lemma 6.6 A.1.1 Proof of estimate (6.13) We estimate the terms of bh(ξ, Πhϕ − ϕ) over the interiors and boundaries of elements separately. Let $${\Pi _{h}^{1}}$$ be the L2($$\varOmega$$)-projection onto the space of discontinuous piecewise linear polynomials on $${\mathcal{T}}_{h}$$. Since on each $$K\in{\mathcal{T}}_{h}$$ it holds that $$\nabla \!\tilde{\xi }|_{K}\in P^{p-1}(K)$$, then $${\Pi _{h}^{1}} a\cdotp \nabla \!\tilde{\xi }\in S_{h}$$. Hence, due to (6.1), we have   \begin{align} \sum_{K}\int_{K}{\Pi_{h}^{1}} a\cdotp\nabla\!\tilde{\xi}(\Pi_{h}\phi-\phi)\,\mathrm{d} x =0. \nonumber \end{align}Due to standard approximation results, we have $$\left \|a-{\Pi _{h}^{1}} a\right \|_{L^{\infty }(K)}\leqslant Ch_{K}|a|_{W^{1,\infty }(K)}$$, thus we can estimate the interior terms of bh(ξ, Πhϕ − ϕ) as   \begin{align} \sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp\nabla\!\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} x =&\,\sum_{K\in{\mathscr{T}}_{h}}\int_{K} \left(a-{\Pi_{h}^{1}} a\right)\cdotp\nabla\!\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} x \nonumber\\ \leqslant&\,\sum_{K\in{\mathscr{T}}_{h}} Ch_{K}|a|_{W^{1,\infty}}C_{I}h_{K}^{-1}\|\xi\|_{L^{2}(K)} \left\|\Pi_{h}\phi-\phi\right\|_{L^{2}(K)}\nonumber\\ \leqslant&\, C\sum_{K\in{\mathscr{T}}_{h}} \max_{x\in K}e^{\mu(x,t)} \|\tilde{\xi}\|_{L^{2}(K)} CL_{\mu} e^{hL_{\mu}} h \max_{x\in K}e^{-\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)}\nonumber\\\leqslant&\, CL_{\mu} e^{2hL_{\mu}}h\|\tilde{\xi}\|^{2},\end{align} (A.1)due to the inverse inequality, Lemma 6.3 and the estimate $$\max _{x\in K}e^{\mu (x,t)}\max _{x\in K}e^{-\mu (x,t)}\leqslant e^{L_{\mu } h_{K}}$$. For the boundary terms of bh(ξ, Πhϕ − ϕ) we get   \begin{align*} & -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\xi](\Pi_{h}\phi-\phi)\,\mathrm{d} S -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} S \\ &\ \leqslant\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} \!|a\cdotp \mathbf{n}|\max_{x\in K}(e^{\mu(x,t)})\big|[\tilde{\xi}]\big||\Pi_{h}\phi-\phi| \,\mathrm{d} S \!+\!\!\sum_{K\in{\mathscr{T}}_{h}}\! \int_{\partial K^{-}\cap\partial\varOmega} \!|a\cdotp \mathbf{n}|\max_{x\in K}(e^{\mu(x,t)})|\tilde{\xi}||\Pi_{h}\phi\!-\phi|\,\mathrm{d} S \\ &\ \leqslant \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} \!|a\cdotp \mathbf{n}|[\tilde{\xi}]^{2}\,\mathrm{d} S \!+\!\frac{1}{8}\!\sum_{K\in{\mathscr{T}}_{h}} \!\int_{\partial K^{-}\cap\partial\varOmega} |a\cdotp \mathbf{n}|\tilde{\xi}^{2}\,\mathrm{d} S\!+\!2\!\sum_{K\in{\mathscr{T}}_{h}}\!\max_{x\in K}e^{2\mu(x,t)}\!\!\int_{\partial K}\!|a\cdotp\mathbf{n}||\Pi_{h}\phi\!-\!\phi|^{2}\mathrm{d} S \\ &\ \leqslant \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\right) + CL_{\mu}^{2} e^{2hL_{\mu}}\sum_{K}\max_{x\in K}e^{2\mu(x,t)} h_{K} \max_{x\in K}e^{-2\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)}^{2} \end{align*}by Young’s inequality and Lemma 6.3. Again we estimate $$\max _{x\in K}e^{2\mu (x,t)}\max _{x\in K}e^{-2\mu (x,t)}\leqslant e^{2L_{\mu } h_{K}}$$ in the last inequality, which completes the proof after combining with (A.1). □ 8.1.2 Proof of estimate (6.14) We have by Green’s theorem   \begin{align} b_{h}(\eta,\phi)=&\, \sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K} a\cdotp\mathbf{n}\eta\phi\,\mathrm{d} S -\int_{K} (\operatorname{div} a)\eta\phi\,\mathrm{d} x -\int_{K} a\cdotp\nabla\!\phi\eta\,\mathrm{d} x \right.\nonumber\\ &\qquad\quad\left.-\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\eta]\phi\,\mathrm{d} S -\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\eta\phi\,\mathrm{d} S \right).\end{align} (A.2)The first integral over K can be estimated as   \begin{align} -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} (\operatorname{div} a)\eta\phi\,\mathrm{d} x \leqslant Ch^{p+1}|u(t)|_{H^{p+1}}\max_{x\in \varOmega}e^{-\mu(x,t)}\|\tilde{\xi}\| \leqslant Ch^{p+1}|u(t)|_{H^{p+1}}\|\tilde{\xi}\|, \nonumber \end{align}because $$\mu \geqslant 0$$. Since ϕ = e−2μξ, we get for the second integral over K in (A.2),   \begin{align} \begin{split} -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp\nabla\!\phi\,\eta\,\mathrm{d} x =\sum_{K\in{\mathscr{T}}_{h}}\int_{K} 2a\cdotp\nabla\!\mu\,e^{-2\mu}\xi\eta\,\mathrm{d} x -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} e^{-2\mu}a\cdotp\nabla\!\xi\,\eta\,\mathrm{d} x . \end{split} \end{align} (A.3)The first right-hand side term in (A.3) can be estimated by   \begin{align*} \sum_{K\in{\mathscr{T}}_{h}}\int_{K} 2a\cdotp\nabla\!\mu\,e^{-2\mu}\xi\eta\,\mathrm{d} x\! =&\,\!\sum_{K\in{\mathscr{T}}_{h}}\!\int_{K} \!2a\cdotp\!\nabla\!\mu\,e^{-\mu}\tilde{\xi}\eta\,\mathrm{d} x\!\leqslant\!\!\sum_{K\in{\mathscr{T}}_{h}} CL_{\mu}\max_{x\in K}e^{-\mu(x,t)}h^{p+1}|u(t)|_{{H^{p+1}}(K)} \|\tilde{\xi}\|_{L^{2}(K)} \nonumber\\ \leqslant&\, CL_{\mu}^{2} h\|\tilde{\xi}\|^{2} +C h^{2p+1}|u(t)|_{{H^{p+1}}}^{2}. \end{align*}The second right-hand side term in (A.3) can be estimated similarly to (A.1), by the definition of η:   \begin{align*} -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} e^{-2\mu}a\cdotp\nabla\!\xi\,\eta\,\mathrm{d} x =&\,\sum_{K\in{\mathscr{T}}_{h}}\int_{K} \left({\Pi_{h}^{1}}(e^{-2\mu}a)-e^{-2\mu}a\right) \cdotp\nabla\!\xi\eta\,\mathrm{d} x \\ \leqslant&\,\sum_{K\in{\mathscr{T}}_{h}} Ch_{K}|e^{-2\mu}a|_{W^{1,\infty}}C_{I}h_{K}^{-1}\|\xi\|_{L^{2}(K)} h^{p+1}|u(t)|_{H^{p+1}(K)}\\\leqslant&\, C(1+L_{\mu})\max_{x\in \varOmega}e^{-\mu(x,t)}h^{p+1}|u(t)|_{{H^{p+1}}}\|\tilde{\xi}\|\\ \leqslant&\, C(1+L_{\mu})^{2}h\|\tilde{\xi}\|^{2} +Ch^{2p+1}|u(t)|_{{H^{p+1}}}^{2}, \end{align*}since $$\left |e^{-2\mu }a\right |_{W^{1,\infty }}= \|e^{-2\mu }\textrm{div} a -2\nabla \!\mu \, e^{-2\mu }a \|_{L^{\infty }}\leqslant C(1+L_{\mu })$$. As for the boundary terms in (A.2), we can split the integral over ∂K into integrals over the separate parts ∂K− ∖ ∂$$\varOmega$$, ∂K− ∩ ∂$$\varOmega$$, ∂K+ ∖ ∂$$\varOmega$$ and ∂K+ ∩ ∂$$\varOmega$$, similarly to the proof of Lemma 6.5. Thus, several terms are canceled out:   \begin{align}&\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S -\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\eta]\phi\,\mathrm{d} S -\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\eta\phi\,\mathrm{d} S \right)\nonumber\\ &\quad=\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K^{+}\setminus\partial\varOmega} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S +\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})\eta^{-}\phi\,\mathrm{d} S \right)\nonumber\\ &\quad=\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp\mathbf{n})\eta^{-}[\phi]\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S \right)\end{align} (A.4)using a similar identity to (6.12). Finally, we can use Young’s inequality to estimate (A.4) further as   \begin{align*}\ldots\leqslant&\,\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K^{-}\setminus\partial\varOmega} |a\cdotp\mathbf{n}||\eta^{-}|e^{-\mu}\big|[\tilde{\xi}]\big|\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} |a\cdotp\mathbf{n}||\eta|e^{-\mu}|\tilde{\xi}|\,\mathrm{d} S \right)\\ \leqslant&\, \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\bigg(\int_{\partial K^{-}\setminus\partial\varOmega} |a\cdotp\mathbf{n}|[\tilde{\xi}]^{2}\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} |a\cdotp\mathbf{n}|\tilde{\xi}^{2}\,\mathrm{d} S \bigg) +2\sum_{K}\int_{\partial K}|a\cdotp\mathbf{n}|\eta^{2}e^{-2\mu}\,\mathrm{d} S \\ \leqslant&\, \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\right) +Ch^{2p+1}|u(t)|_{H^{p+1}}^{2}. \nonumber \end{align*}The proof is completed by gathering all the above estimates of the individual terms of bh(η, ϕ). □ 8.1.3 Proof of estimate (6.15) We use Lemmas 6.2 and 6.3 to estimate   \begin{align*} &b_{h}(\eta,\Pi_{h}\phi-\phi)\\ &\quad=\sum_{K\in{\mathscr{T}}_{h}}\bigg(\int_{K} (a\cdotp \nabla\! \eta) (\Pi_{h}\phi-\phi)\,\mathrm{d} x -\!\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\eta](\Pi_{h}\phi\!-\phi)\,\mathrm{d} S\!-\!\!\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\eta (\Pi_{h}\phi\!-\!\phi)\,\mathrm{d} S \bigg)\\ &\quad\leqslant Ch^{p}|u(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}} h \max_{x\in\varOmega}e^{-\mu(x,t)}\|\tilde{\xi}\|+Ch^{p+1/2}|u(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}}h^{1/2} \max_{x\in\varOmega}e^{-\mu(x,t)}\|\tilde{\xi}\|. \end{align*}The application of Young’s inequality completes the proof. □ 8.2 Proof of Lemma 6.7 Lemma 6.3 gives us   \begin{align*} \sum_{K\in{\mathscr{T}}_{h}}\int_{K} c\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} x \leqslant&\, C\sum_{K\in{\mathscr{T}}_{h}}\max_{x\in K}e^{\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)} L_{\mu} e^{hL_{\mu}}h_{K} \max_{x\in K}e^{-\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)}\\ \leqslant&\, CL_{\mu} e^{2hL_{\mu}} h\|\tilde{\xi}\|^{2}. \nonumber \end{align*}As for the second estimate, we write ch(η, Πhϕ) = ch(η, ϕ) + ch(η, Πhϕ − ϕ) and estimate by Lemmas 6.2 and 6.3:   \begin{align} \begin{split} c_{h}(\eta,\phi) =&\,\int_{\varOmega} c\eta\phi\,\mathrm{d} x =\int_{\varOmega} c\eta e^{-\mu}\tilde{\xi}\,\mathrm{d} x \leqslant Ch^{p+1}|u(t)|_{H^{p+1}}\|\tilde{\xi}\|,\\ c_{h}(\eta,\Pi_{h}\phi-\phi) =&\,\int_{\varOmega} c\eta (\Pi_{h}\phi-\phi)\,\mathrm{d} x \leqslant Ch^{p+1}|u(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}}h \|\tilde{\xi}\|. \nonumber \end{split} \end{align}Combining these two estimates with Young’s inequality gives the desired result. □ 8.3 Proof of Lemma 6.8 We use Lemmas 6.2, 6.3 and the definition of ϕ:   \begin{align} \begin{split} \left(\frac{\partial \eta}{\partial t},\Pi_{h}\phi\right)= \left(\frac{\partial \eta}{\partial t},\phi\right) +\left(\frac{\partial \eta}{\partial t},\Pi_{h}\phi-\phi\right) \!\leqslant\! Ch^{\,p+1}|u_{t}(t)|_{H^{p+1}}\|\tilde{\xi}\| +Ch^{p+1}|u_{t}(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}}h \|\tilde{\xi}\|. \nonumber \end{split} \end{align}Again, Young’s inequality yields the final estimate. □ © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

On the time growth of the error of the DG method for advective problems

Loading next page...
 
/lp/ou_press/on-the-time-growth-of-the-error-of-the-dg-method-for-advective-1msSvtMPIV
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/dry013
Publisher site
See Article on Publisher Site

Abstract

Abstract In this paper we derive a priori $$L^{\infty }(L^{2})$$ and L2(L2) error estimates for a linear advection–reaction equation with inlet and outlet boundary conditions. The goal is to derive error estimates for the discontinuous Galerkin method that do not blow up exponentially with respect to time, unlike the usual case when Gronwall’s inequality is used. While this is possible in special cases, such as divergence-free advection fields, we take a more general approach using exponential scaling of the exact and discrete solutions. Here we use a special scaling function, which corresponds to time taken along individual pathlines of the flow. For advection fields, where the time that massless particles carried by the flow spend inside the spatial domain is uniformly bounded from above by some $$\widehat{T}$$, we derive $$\mathcal{O}$$(hp+1/2) error estimates where the constant factor depends only on $$\widehat{T}$$, but not on the final time T. This can be interpreted as applying Gronwall’s inequality in the error analysis along individual pathlines (Lagrangian setting), instead of physical time (Eulerian setting). 1. Introduction The discontinuous Galerkin (DG) finite element method introduced in the study by Reed & Hill (1973) is an increasingly popular method for the numerical solution of partial differential equations. The DG method was first formulated for a neutron transport equation and such problems remain the major focus of the DG community. It is for problems of an advective or convective nature that the DG method is suited best and shows its strengths compared to other numerical methods for such problems. In this paper, we shall consider a scalar time-dependent linear advection–reaction equation of the form   \begin{align} \frac{\partial u}{\partial t}+a\cdotp\nabla\!\, u+cu=0. \end{align} (1.1)We will discretize the problem in space using the DG method with the upwind numerical flux on unstructured simplicial meshes in $$\mathbb{R}^{d}$$ which may contain hanging nodes. The first analysis of the DG method for a stationary advection–reaction problem with constant coefficients was made in the study by Lesaint & Raviart (1974), later improved in the study by Johnson & Pitkäranta (1986). In this paper we derive a priori estimates of the error eh = u − uh, where uh is the DG solution. The goal is to derive estimates of eh in the $$L^{\infty }\left (L^{2}\right )$$- and $$L^{2}\left (L^{2}\right )$$-norms of the order Chp+1/2, where the constant C does not depend exponentially on the final time $$T\to +\infty $$. Such results already exist in the literature; however, they are derived under the ellipticity condition   \begin{align} c-\tfrac{1}{2}\operatorname{div}a\geqslant\gamma_{0}>0 \end{align} (1.2)for some constant $$\gamma_0 $$. In the paper by Feistauer & Švadlenka (2004), $$\gamma_0=0 $$ is also admissible, which corresponds to the interesting case of a divergence-free advection field a. Without assuming (1.2), when one proceeds straightforwardly, at some point Gronwall’s inequality must be used in the proofs. This, however, results in exponential growth of the constant factor C in the error estimate with respect to time. Here, we will circumvent condition (1.2) by considering an exponential scaling transformation   \begin{align} u(x,t)=e^{\mu(x,t)}\tilde{u}(x,t),\quad u_{h}(x,t)=e^{\mu(x,t)}\tilde{u}_{h}(x,t) \end{align} (1.3)of the exact and discrete solutions. Substituting (1.3) into (1.1) and its discretization results in a new equation for $$\tilde{u}$$ with a modified reaction term $$\tilde{c}=\frac{\partial \mu }{\partial t}+a\cdotp \nabla \!\mu +c$$. The function μ can then be chosen in various ways in order to satisfy the ellipticity condition (1.2) for the new equation. After performing the error analysis, the resulting error estimates for $$\tilde{e}_{h}=\tilde{u}-\tilde{u}_{h}$$ can be transformed back via (1.3) to obtain estimates for eh. The use of transformations similar to (1.3) is not new. The trick known as exponential scaling, corresponding to taking μ(x, t) = αt for some sufficiently large constant α, is very well known from the partial differential equation and numerical analysis communities. However, the application of such a scaling transformation inherently leads to exponential dependence of the error on t after transforming from $$\tilde{e}_{h}$$ back to eh. Other choices of μ are possible, e.g. μ(x, t) = μ0⋅x for some constant vector μ0. This choice was used in the stationary case in the studies by Nävert (1982) and Johnson & Pitkäranta (1986). Taking μ(x) independent of t corresponds to the analysis of the stationary case from the study by Ayuso & Marini (2009). In the nonstationary case the use of this stationary transform would lead to too restrictive assumptions on the advection field. In this paper we construct the function μ from (1.3) using characteristics of the advection field. Namely, μ will be proportional to time along individual pathlines of the flow. If pathlines exist only for a finite time bounded by some $$\widehat{T}$$ for each particle entering and leaving the spatial domain $$\varOmega$$, we obtain error estimates exponential in $$\widehat{T}$$ and not the final physical time T. This is a result we would expect if we applied Gronwall’s lemma in the Lagrangian framework along individual pathlines (which exist only for a finite time uniformly bounded by $$\widehat{T}$$) instead of the usual application of Gronwall with respect to physical time in the Eulerian framework. The analysis can be carried out under mild assumptions that a(⋅, ⋅) satisfies the assumptions of the Picard–Lindelöf theorem and that there are no characteristic boundary points on the inlet (i.e. a ⋅ n is uniformly bounded away from zero on the inlet). The paper is organized as follows. After introducing the continuous problem in Section 2, we discuss the exponential scaling transform, its variants and application in the weak formulation in Section 3. In Section 4 we construct the scaling function μ and prove results on its regularity and other properties needed in the analysis. In Section 5 we introduce the DG formulation and its basic properties. In Sections 6 and 7 we analyse the DG advection and reaction forms and estimate the error of the method. In order to focus on the main ideas of the analysis, we postpone some of the technical estimates of the DG forms to the Appendix. We use (⋅, ⋅) to denote the L2($$\varOmega$$) scalar product and ∥⋅∥ for the L2($$\varOmega$$)-norm. To simplify the notation, we shall drop the argument $$\varOmega$$ in Sobolev norms, e.g. $$\|\cdot \|_{H^{p+1}}$$ denotes the Hp+1($$\varOmega$$)-norm. We will also denote the Bochner norms over the whole considered interval (0, T) in concise form, e.g. $$\|u\|_{L^{2}\left (H^{p+1}\right )}$$ denotes the $$L^{2}(0,T;H^{p+1}(\varOmega ))$$-norm of u. 2. Continuous problem Let $$\varOmega \subset \mathbb{R}^{d}$$, $$d\in \mathbb{N}$$ be a bounded polygonal (polyhedral) domain with Lipschitz boundary ∂$$\varOmega$$. Let $$0<T\leqslant +\infty $$ and let QT = $$\varOmega$$ × (0, T) be the space–time domain. We note that we admit the time interval to be infinite in the special case of $$T=+\infty $$. We seek $$u:Q_{T}\to \mathbb{R}$$ such that   \begin{alignat}{2} \frac{\partial u}{\partial t}+a\cdotp\nabla\!\, u+cu&=0 &&\textrm{in }Q_{T},\\ u&=u_{D} &&\textrm{on }\partial\varOmega^{-}\times(0,T),\nonumber\\ u(x,0)&=u^{0}(x),\quad && x\in\varOmega\nonumber. \end{alignat} (2.1)Here $$a:\overline{Q_{T}}\to \mathbb{R}^{d}$$ and $$c:\overline{Q_{T}}\to \mathbb{R}$$ are the given advective field and reaction coefficient, respectively. By ∂$$\varOmega$$− we denote the inflow part of the boundary, i.e.   \begin{align} \partial\varOmega^{-}=\left\{x\in\partial\varOmega; a(x,t)\cdot\mathbf{n}(x)<0\quad \forall\; t\in(0,T)\right\}\!, \nonumber \end{align}where n(x) is the unit outer normal to ∂$$\varOmega$$ at x. For simplicity we assume that the inflow boundary is independent of t, in other words we assume that a(x, t) ⋅ n(x) < 0 for all t and all x ∈ ∂$$\varOmega$$−, and a(x, t) ⋅ n(x) is non-negative on ∂$$\varOmega$$ ∖ ∂$$\varOmega$$− for all t. In general ∂$$\varOmega$$− can depend on t due to the nonstationary nature of a. However, we do not consider this case since it would only make the notation more complicated without changing the key points of the analysis itself. We assume that the reaction coefficient satisfies $$c\in C([0,T);L^{\infty }(\varOmega ))\cap L^{\infty }(Q_{T})$$. Furthermore, let $$a\in C([0,T);W^{1,\infty }(\varOmega ))$$ with a, ∇a uniformly bounded a.e. in QT. In other words, a (⋅, ⋅) satisfies the assumptions of the Picard–Lindelöf theorem: continuity with respect to x, t and uniform Lipschitz continuity with respect to x. The spatial Lipschitz constant of a will be denoted by La. The DG method for the stationary version of problem (2.1) was analysed in the studies by Nävert (1982), Johnson & Pitkäranta (1986) and Ayuso & Marini (2009). In this case the well-posedness of the problem is guaranteed by the assumption that the flow field possesses neither closed curves nor stationary points, cf. the study by Devinatz et al. (1974). For the nonstationary problem (2.1), the situation is simpler. Here, the problem can be solved using the method of characteristics, similarly to the method we use in Section 4, cf. Evans (2010). Along each characteristic curve (pathline), problem (2.1) reduces to a simple, ordinary differential equation, which can be solved in closed form. The result is exponential growth or decay of u along the pathline with the rate given by the reaction coefficient c. As the flow field a satisfies the assumptions of the Picard–Lindelöf theorem, we have existence and uniqueness of the pathlines. Since the Dirichlet boundary condition is prescribed only at the inflow boundary, this gives the well-posedness of problem (2.1). We note that also a ‘parabolic’ approach to (2.1) is possible to obtain ‘energy’ estimates of the time growth of the L2($$\varOmega$$)-norm of the solution. This is done by Gronwall’s lemma after testing the weak form of (2.1) with the solution itself and applying Green’s theorem in the advection terms; cf. e.g. Mizohata (1973) and Larsson & Thomée (2003). The result is an upper bound on the time growth of the solution that is exponential in the quantity $$c-\tfrac{1}{2}\textrm{div} a$$ if this quantity is non-negative; compare with the ellipticity condition (3.1). We note that from the theoretical and practical points of view, the case of systems of equations of the considered type is much more interesting; cf. e.g. Benzoni-Gavage & Serre (2007). 3. Exponential scaling The standard assumption on the coefficients a, c found throughout the numerical literature is   \begin{align} c-\tfrac{1}{2}\operatorname{div} a\geqslant\gamma_{0}>0 \quad\textrm{ on }Q_{T} \end{align} (3.1)for some constant $$\gamma_0 $$. This assumption comes from the requirement that the weak formulations of the advection and reaction terms give an elliptic bilinear form on the corresponding function space. In the study by Feistauer & Švadlenka (2004), this assumption is avoided by using exponential scaling in time, i.e. the transformation u(x, t) = eαtw(x, t) for $$\alpha \in \mathbb{R}$$. Substitution into (2.1) and division of the whole equation by the common positive factor eαt leads to a modified reaction coefficient $$\tilde c=c+\alpha $$ in the new equation for w. By choosing the constant α sufficiently large, (3.1) will be satisfied for the new equation. In the study by Feistauer & Švadlenka (2004), error estimates that grow linearly in time are then derived for the DG scheme. However, for α > 0, if one transforms the resulting estimates back to the original problem using the exponential scaling transformation, the result is an estimate that depends exponentially on T. Another possibility to avoid assumption (3.1) is a transformation similar to exponential scaling, but with respect to the spatial variable; cf. the studies by Nävert (1982), Johnson & Pitkäranta (1986) and Roos et al. (2008): let $$\mu _{0}\in \mathbb{R}^{d}$$, then we write   \begin{align} u(x,t)=e^{\mu_{0}\cdotp x}\tilde u(x,t). \end{align} (3.2)If we assume that u is sufficiently regular, as we shall in this paper, by substituting (3.2) into (2.1) and dividing the whole equation by the strictly positive function $$e^{\mu _{0}\cdotp x}$$, we get the new problem   \begin{align} \frac{\partial \tilde{u}}{\partial t}+a\cdotp\nabla\!\, \tilde{u}+(a\cdotp\mu_{0}+c)\tilde{u}=0. \end{align} (3.3)The condition corresponding to (3.1) now reads as follows: there exists $$\mu _{0}\in \mathbb{R}^{d}$$ such that   \begin{align} a\cdotp\mu_{0}+c-\tfrac{1}{2}\operatorname{div} a\geqslant\gamma_{0}>0 \quad\textrm{ on }Q_{T}. \end{align} (3.4) A possible generalization of the exponential scaling transformation (3.2) is taking a sufficiently smooth function $$\mu :\varOmega \to \mathbb{R}$$ and setting   \begin{align} u(x,t)=e^{\mu(x)}\tilde u(x,t). \end{align} (3.5)Again, substituting into (2.1) and dividing by eμ(x), we get the new problem   \begin{align} \frac{\partial \tilde{u}}{\partial t}+a\!\cdot\nabla \tilde{u}+(a\cdotp\nabla\!\mu+c)\tilde{u}=0. \nonumber \end{align}The condition corresponding to (3.1) and (3.4) now reads as follows: there exists $$\mu :\varOmega \to \mathbb{R}$$ such that   \begin{align} a\cdotp\nabla\!\mu+c-\tfrac{1}{2}\operatorname{div} a\geqslant\gamma_{0}>0 \quad\textrm{ on }Q_{T}. \end{align} (3.6)This is essentially the approach used in the study by Ayuso & Marini (2009) for a stationary advection–diffusion–reaction problem. As shown in the study by Devinatz et al. (1974), the existence of a function $$\mu :\varOmega \to \mathbb{R}$$ such that $$a\cdotp \nabla \!\mu \geqslant \gamma _{0}$$ is equivalent to the property that the advective field a possesses neither closed curves nor stationary points. The uniformly positive term $$a\cdotp \nabla \!\mu $$ can then be used to dominate the possibly negative term $$c-\frac{1}{2}\textrm{div} a$$ in order to satisfy condition (3.6). If we used the choice (3.5) in our analysis, we would need to assume the nonexistence of closed curves or stationary points of the flow field for all t. This assumption is too restrictive. Furthermore, in the study by Ayuso & Marini (2009), rather high regularity of μ is needed in the analysis, namely $$\mu \in W^{p+1,\infty }(\varOmega )$$, where p is the polynomial degree of the DG scheme. In this paper, we shall generalize the transformation (3.5) using a function $$\mu :Q_{T}\to \mathbb{R}$$ and set   \begin{align} u(x,t)=e^{\mu(x,t)}\tilde u(x,t). \end{align} (3.7)In our analysis we will need only minimum regularity of the scaling function μ, in contrast to the study by Ayuso & Marini (2009). Namely, we assume that μ is Lipschitz continuous both in space and time. This implies that the derivatives μt and ∇μ exist for almost all x and t. Hence, the following steps are valid for almost all x and t and particularly, they will hold in the integral (i.e. weak) sense. We note that the Lipschitz continuity in space and time will hold for the specific construction of μ from Section 4, as shown in Section 4.1. As in the previous cases, substituting (3.7) into (2.1) and dividing by eμ(x, t) gives the new problem   \begin{align} \frac{\partial \tilde{u}}{\partial t}+a\cdot\nabla\! \tilde{u}+\left(\frac{\partial \mu}{\partial t}+a\cdotp\nabla\!\mu+c\right)\tilde{u}=0. \end{align} (3.8)The condition corresponding to (3.1), (3.4) and (3.6) now reads as follows: there exists $$\mu :Q_{T}\to \mathbb{R}$$ such that   \begin{align} \frac{\partial \mu}{\partial t}+a\cdotp\nabla\!\mu+c-\tfrac{1}{2}\operatorname{div} a\ge\gamma_{0}>0 \quad\text{a.e. in }Q_{T}. \end{align} (3.9)This is the condition we will assume throughout the paper. Using the transformation (3.7) instead of (3.5) gives sub-exponential growth or even uniform boundedness of the constant factor in the error estimate with respect to time. If the original problem (2.1) does not satisfy (3.1), one is tempted to numerically solve the transformed equation (3.8) instead, using DG or any other method. This is done e.g. in Roos et al. (2008) for the case of linear μ, i.e. (3.2). However, then we are numerically solving a different problem and obtain different results, as the DG solutions of (2.1) and (3.8) are not related by the simple relation (3.7), unlike the exact solutions. However, as we will show in this paper, one can analyse the DG method for the original problem (2.1), while taking advantage of the weaker ellipticity condition (3.9) for the transformed problem (3.8). Since the DG scheme is based on a suitable weak formulation, the first step is to reformulate the transformation (3.7) within the weak, rather than strong formulation of (2.1). The key step in deriving (3.8) from (2.1) is dividing the whole equation by the common factor eμ(x, t). However, if we substitute (3.7) into the weak formulation   \begin{align} \int_{\varOmega}\frac{\partial u}{\partial t}v+a\cdotp\nabla\!\, uv+cuv\,\mathrm{d} x =0, \end{align} (3.10)it is not easy to divide the equation by eμ(x, t), since this is inside the integral. The solution is to test (3.10) by the new test function $$\hat{v}(x,t)=e^{-\mu (x,t)}v(x,t)$$. Due to the opposite signs in the exponents, the exponential factors cancel each other and we obtain the weak formulation of (3.8):   \begin{align} \int_{\varOmega}\frac{\partial \tilde{u}}{\partial t}v+a\cdotp\nabla\!\, \tilde{u} v+\left(\frac{\partial \mu}{\partial t} +a\cdotp\nabla\!\mu+c\right)\tilde{u} v\,\mathrm{d} x =0. \nonumber \end{align}We note that the transformations $$u\mapsto \tilde{u}$$ and $$v\mapsto \hat{v}$$ are bijections. Furthermore, since we assume that μ is Lipschitz continuous in space, then the factor eμ(⋅, t) lies in $$W^{1,\infty }(\varOmega )$$ and the transformation (3.7) is a bijection from the function space V into itself, where   \begin{align} V=\big\{u\in H^{1}(\varOmega); u|_{\partial\varOmega^{-}}=0\big\} \nonumber \end{align}is the space for weak solutions and test functions for (2.1) we consider. This is the case uD = 0, the general case can be treated by the standard Dirichlet lifting procedure. We note that in the general case, the appropriate function spaces for the solution of (2.1) would be e.g. L1($$\varOmega$$) or BV($$\varOmega$$) with respect to space. However, since we assume smooth solutions and the linear case, it is possible to consider the Hilbert setting of H1($$\varOmega$$) which makes the analysis easier in the finite element or DG framework. In the DG method, the above procedure using transform (3.7) cannot be directly applied, since if v ∈ Sh—the discrete space of discontinuous piecewise polynomial functions—then $$\hat{v}=e^{-\mu }v$$ will no longer lie in Sh and therefore cannot be used as a test function in the formulation of the method. The solution is to test with a suitable projection of $$\hat{v}$$ onto Sh and estimate the difference, as we shall do in Section 6. 4. Construction of the function μ In this section we show a construction of the function μ satisfying (3.9). Many different constructions of μ are possible depending on the assumptions on the vector field a. For example, we can always take μ(x, t) = αt for some $$\alpha \geqslant 0$$. This corresponds to standard exponential scaling and, as we have seen, this choice leads to exponential growth in time of the error estimate. Another possibility was the approach of Ayuso & Marini (2009) mentioned in Section 2, where a suitable function μ(x) exists if a, which is stationary, possesses no closed curves or stationary points. Here we show another possibility with an interesting interpretation. If $$c-\frac{1}{2}\textrm{div} a$$ is negative or changes sign frequently, we can use the expression μt + a ⋅∇μ to dominate this term everywhere. If we choose μ1 such that   \begin{align} \frac{\partial\mu_{1}}{\partial t}+a\cdot\nabla\!\mu_{1}=1\quad \textrm{on }Q_{T}, \end{align} (4.1)then by multiplying μ1 by a sufficiently large constant, we can satisfy the ellipticity condition (3.9) for a chosen $$\gamma_0 $$ > 0. Equation (4.1) can be explicitly solved using characteristics. We define pathlines of the flow, i.e. the family of curves $$S\left (t;x_{0},t_{0}\right )$$ by   \begin{align} S(t_{0};x_{0},t_{0})=x_{0}\in\overline{\varOmega},\quad \frac{\mathrm{d} S(t;x_{0},t_{0})}{\mathrm{d} t} =a\left(S(t;x_{0},t_{0}),t\right). \end{align} (4.2)This means that S(⋅; t0, x0) is the trajectory of a massless particle in the nonstationary flow field a passing through point x0 at time t0. It is convenient to choose the parameter t0 minimal for each pathline—then the pair (x0, t0) is the ‘origin’ of the pathline. In other words, for each $$x_{0}\in \overline{\varOmega }$$, there is a pathline S(⋅; x0, 0) corresponding to trajectories of particles present in $$\varOmega$$ at the initial time t0 = 0. Then there are trajectories of particles entering through the inlet part of ∂$$\varOmega$$: for each x0 ∈ ∂$$\varOmega$$− and all t0 there exists a pathline S(⋅; x0, t0) originating at (x0, t0) ∈ ∂$$\varOmega$$−× [0, T). Equation (4.1) can then be rewritten along the pathlines as   \begin{align} \frac{\mathrm{d}\, \mu_{1}(S(t;x_{0},t_{0}),t)}{\mathrm{d} t} = \left(\frac{\partial\mu_{1}}{\partial t}+a\cdot\nabla\!\mu_{1}\right) (S(t;x_{0},t_{0}),t)=1, \nonumber \end{align}therefore   \begin{align} \mu_{1}(S(t;x_{0},t_{0}),t) =t-t_{0}. \end{align} (4.3)Here any constant can be chosen instead of t0; however, this choice is the most convenient. With this choice, we have at the origin of a pathline   \begin{align} \mu_{1}(S(t_{0};x_{0},t_{0}),t_{0}) =0 \nonumber \end{align}and the value of μ1 along this pathline is simply the time elapsed since t0. In other words, μ1(x, t) is the time a particle carried by the flow passing through x ∈$$\varOmega$$ at time t has spent in $$\varOmega$$ up to the time t. In this paper we assume this quantity to be uniformly bounded in order for μ1 and μ to be uniformly bounded; cf. assumption (4.7). It is important to note that t0 depends on the specific trajectory considered and that t0 will in general be different for each (x, t). Therefore, even though t and t0 may be arbitrarily large if we consider an infinite time interval, the difference t − t0 = μ1, which measures the time elapsed since the entry of the corresponding particle into $$\varOmega$$, will remain bounded under this assumption. Now that we have constructed a function satisfying (4.1), we can choose $$\gamma_0 $$ and define e.g.   \begin{align} \mu(x,t)=\mu_{1}(x,t)\left(\left|\inf_{Q_{T}}\left(c-\tfrac{1}{2}\operatorname{div} a\right)^{-}\right|+\gamma_{0}\right), \end{align} (4.4)where $$f^{-}=\min \{0,f\}$$ is the negative part of f. Then   \begin{align} \frac{\partial\mu}{\partial t} +a\cdot\nabla\!\mu+c-\tfrac{1}{2}\operatorname{div} a = \left|\inf_{Q_{T}}\left(c-\tfrac{1}{2}\operatorname{div} a\right)^{-}\right|+\gamma_{0} +c-\tfrac{1}{2}\operatorname{div} a \ge\gamma_{0}. \nonumber \end{align} In the main result, of this paper, Theorem 7.1 and Corollary 7.2, this choice of μ leads to estimates of the following form for the DG error eh (cf. (7.1) and (7.10)):   \begin{align} \|e_{h}\|_{L^{\infty}\left(L^{2}\right)} +\sqrt{\gamma_{0}}\|e_{h}\|_{L^{2}\left(L^{2}\right)}\leqslant Ce^{\widehat{T}(|\inf_{Q_{T}}(c-\frac{1}{2}\operatorname{div} a)^{-}|+\gamma_{0})}h^{p+1/2}, \end{align} (4.5)where $$\widehat{T}=\sup _{Q_{T}}\mu _{1}$$, i.e. $$\widehat{T}$$ is the maximal particle ‘lifetime’, i.e. the maximal time any particle carried by the flow field a spends in $$\varOmega$$, since by (4.3), μ1(x, t) is the time the particle at (x, t) carried by the flow has spent in $$\varOmega$$. If we compare this to estimates obtained by a straightforward analysis using Gronwall’s lemma without any ellipticity assumption, we would expect   \begin{align} \|e_{h}\|_{L^{\infty}(L^{2})} \leqslant Ce^{T(\sup_{Q_{T}}|c-\frac{1}{2}\operatorname{div} a|)}h^{p+1/2}. \end{align} (4.6)Comparing (4.6) to (4.5), we see that the estimates are essentially of similar form, only the exponential dependence on global physical time T has been replaced by dependence on time $$\widehat{T}$$ along pathlines, which is bounded. Effectively, our analysis replaces the application of Gronwall’s lemma in the Eulerian framework with its application in the Lagrangian framework—along individual pathlines which have bounded length. Throughout the paper, we assume the particle lifetime $$\widehat{T}$$ to be bounded. In general, we could consider dependencies of the form   \begin{align} \widehat{T}(t)=\sup_{(x,\vartheta)\in\varOmega\times{(0,t)}} \mu_{1}(x,\vartheta), \nonumber \end{align}i.e. $$\widehat{T}(t)$$ is the maximal time any particle carried by the flow a spends in $$\varOmega$$ up to time t. From (4.5), we can expect exponential dependence of the $$L^\infty(0,t;L^2)$$ and $$L^2(0,t;L^2)$$ norms on $$\widehat{T}(t)$$. The result of the standard analysis (4.6) corresponds to the ‘worst case’ $$\widehat{T}(t)=t$$, i.e. that there is a particle that stays inside $$\varOmega$$ for all t ∈ [0, T). However, considering more general dependencies on time is possible, e.g. $$\widehat{T}=\sqrt{t}$$, leading to growth of the error that is exponential in $$\sqrt{t}$$. 4.1 Regularity of the function μ In our analysis of the DG scheme we will assume that μ satisfies   \begin{equation}\label{4.7}\begin{split} 0\leqslant\mu(x,t)\leqslant&\,\mu_{\max},\\ |\mu(x,t)-\mu(y,t)|\leqslant&\, L_{\mu}|x-y|,\end{split}\end{equation} (4.7)for all x, y ∈ $$\varOmega$$ and t ∈ (0, T). In other words, μ is non-negative, uniformly bounded and Lipschitz continuous in space, where the Lipschitz constant is uniformly bounded for all t. Since $$\varOmega$$ is a Lipschitz domain, hence quasiconvex, this means that $$\mu (t)\in W^{1,\infty }(\varOmega )$$ for all t, with its $$W^{1,\infty }(\varOmega )$$ seminorm uniformly bounded by Lμ for all t. Now we show when μ1 defined by (4.3) satisfies conditions (4.7), especially Lipschitz continuity in space which is necessary in the analysis. We note that since μ is obtained from μ1 by multiplication by a constant, properties derived for μ1 imply those for μ. An obvious example of when μ1 defined by (4.3) will not be Lipschitz continuous is when a vortex ‘touches’ ∂$$\varOmega$$− as in Fig. 1. Then if (x, t) lies on a pathline originating at (x0, t0) such that a(x0, t0) is tangent to ∂$$\varOmega$$, we can find $${\tilde{x}}$$ arbitrarily close to x such that the corresponding pathline is much longer, perhaps winding several times around the vortex. Then $$\mu _{1}(x,t)-\mu _{1}\left ({\tilde{x}},t\right )={\tilde{t}}_{0}-t_{0}$$ can be very large, while $$\|x-{\tilde{x}}\|$$ is arbitrarily small, hence μ1 is not Lipschitz continuous. In fact μ1 is discontinuous at (x, t). In the following lemma, we show that inlet points where a is tangent to ∂$$\varOmega$$ are the only troublemakers. Fig. 1. View largeDownload slide Left: proof of Lemma 4.1. Right: vortex touching ∂$$\varOmega$$− with $$t_{0}\gg{\tilde{t}}_{0}=0$$, i.e. $$t_{0}-{\tilde{t}}_{0}$$ is large, hence μ1 is not Lipschitz continuous at x. Fig. 1. View largeDownload slide Left: proof of Lemma 4.1. Right: vortex touching ∂$$\varOmega$$− with $$t_{0}\gg{\tilde{t}}_{0}=0$$, i.e. $$t_{0}-{\tilde{t}}_{0}$$ is large, hence μ1 is not Lipschitz continuous at x. We note that since $$\varOmega$$ is a bounded Lipschitz domain (hence quasiconvex), when proving Lipschitz continuity of μ1 in space it is sufficient to prove local Lipschitz continuity of μ1 in some neighborhood of each $$x\in \overline{\varOmega }$$ with a Lipschitz constant independent of x. Lemma 4.1 Let $$a\in L^{\infty }(Q_{T})$$ be continuous with respect to time and Lipschitz continuous with respect to space. Let there exist a constant $$a_{\min }>0$$ such that   \begin{align} -a(x,t)\cdot\mathbf{n}\geqslant a_{\min} \end{align} (4.8)for all x ∈ ∂$$\varOmega$$−, t ∈ [0, T). Let μ1 be defined by (4.3) on $$\overline{\varOmega }\times [0,T)$$. Let the time any particle carried by the flow field a(⋅, ⋅) spends in $$\varOmega$$ be uniformly bounded by $$\widehat{T}$$. Then μ1 satisfies assumption (4.7). Proof. By definition, we have $$\mu _{1}(x,t)\geqslant 0$$ for all x, t. By the above considerations, μ1 is bounded by the maximal time particles spend in $$\varOmega$$, which is uniformly bounded. This implies $$\mu _{1}(x,t)\leqslant \widehat{T}=\mu _{\max }$$ for some $$\mu _{\max }$$. Now we will prove Lipschitz continuity. Let t ∈ (0, T) be fixed and let $$x,{\tilde{x}}\in \overline{\varOmega }$$ such that $$|x-{\tilde{x}}|\leqslant \varepsilon $$, where ε will be chosen sufficiently small in the following. Due to the assumptions on a, the pathlines passing through x and $${\tilde{x}}$$ are uniquely determined and originate at some (x0, t0) and $$({\tilde{x}}_{0},{\tilde{t}}_{0})$$, respectively. In other words,   \begin{align} x=S(t;x_{0},t_{0}),\quad{\tilde{x}}=S\left(t;{\tilde{x}}_{0},{\tilde{t}}_{0}\right)\!. \nonumber \end{align}Without loss of generality, let $$t_{0}\ge{\tilde{t}}_{0}>0$$, hence $$x_{0},{\tilde{x}}_{0}\in \partial \varOmega ^{-}$$. The case when $${\tilde{t}}_{0}=0$$ can be treated similarly. Furthermore, we assume that x0 is not a vertex of ∂$$\varOmega$$−—this case will be treated at the end of the proof. Since a is uniformly bounded and Lipschitz continuous in space, there exists δ, such that if t ∈ [0, T) and $$\mathrm{dist}\{x,\partial \varOmega ^{-}\}\leqslant \delta $$ then   \begin{align} -a(x,t)\cdot\mathbf{n}\geqslant a_{\min}/2, \end{align} (4.9)due to (4.8). If ε is sufficiently small, then for the distance of the two considered pathlines at time t0 we have $$ |x_{0}-S (t_{0};{\tilde{x}}_{0},{\tilde{t}}_{0})| =|S(t_{0};x_{0},t_{0} )-S(t_{0};{\tilde{x}}_{0},{\tilde{t}}_{0})|\leqslant \delta $$. Since x0 ∈ ∂$$\varOmega$$− this means that $$S(t_{0};{\tilde{x}}_{0},{\tilde{t}}_{0})$$ is in the δ-neighborhood of ∂$$\varOmega$$− and by (4.9), $$\mathrm{dist}\{S(\vartheta ;\tilde{x}_{0},\tilde{t}_{0}),\partial \varOmega ^{-}\}$$ decreases as ϑ goes from t0 to $${\tilde{t}}_{0}$$ with a rate of at least $$a_{\min }/2$$ due to the uniformity of the bound (4.9). Therefore, $$S(\vartheta ;{\tilde{x}}_{0},{\tilde{t}}_{0})$$ stays in the δ-neighborhood of ∂$$\varOmega$$− for all $$\vartheta \in \left [{\tilde{t}}_{0},t_{0}\right ]$$ and   \begin{align} -a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\cdot\mathbf{n}\geqslant a_{\min}/2 \nonumber \end{align}for all $$\vartheta \in \left [{\tilde{t}}_{0},t_{0}\right ]$$. Moreover, since x0 lies in the interior of an edge of ∂$$\varOmega$$−, by choosing ε sufficiently small, we can ensure that $${\tilde{x}}_{0}$$ also lies on this edge (face). Now we estimate $$\left |\mu _{1}(x,t)-\mu _{1}\left ({\tilde{x}},t\right )\right |=\left |t_{0}-{\tilde{t}}_{0}\right |$$. We have   \begin{align*} x-x_{0} =&\,\int_{t_{0}}^{t}\frac{\mathrm{d} S}{\mathrm{d} t}(\vartheta;x_{0},t_{0})\mathrm{\ d}\vartheta =\int_{t_{0}}^{t} a\left(S(\vartheta;x_{0},t_{0}),\vartheta\right)\mathrm{\ d}\vartheta,\\{\tilde{x}}-{\tilde{x}}_{0} =&\,\int_{{\tilde{t}}_{0}}^{t}\frac{\mathrm{d} S}{\mathrm{d} t}\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right)\mathrm{\ d}\vartheta =\int_{{\tilde{t}}_{0}}^{t} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\mathrm{d}\vartheta. \end{align*}Subtracting these two identities gives us   \begin{align} \begin{split} x_{0}-{\tilde{x}}_{0} =x-{\tilde{x}} +\int_{{\tilde{t}}_{0}}^{t_{0}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\mathrm{\ d}\vartheta +\int_{t_{0}}^{t}a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right) -a\left(S(\vartheta;x_{0},t_{0}),\vartheta\right)\mathrm{d}\vartheta. \end{split} \end{align} (4.10)If we consider n, the normal to ∂$$\varOmega$$− at x0, we see that $$(x_{0}-{\tilde{x}}_{0})\cdot \mathbf{n}=0$$, as both x0 and $${\tilde{x}}_{0}$$ lie on the same edge on ∂$$\varOmega$$−. Therefore, if we multiply (4.10) by −n, we get   \begin{align} 0 =&\,(x-{\tilde{x}})\cdotp\!(-\mathbf{n}) +\int_{{\tilde{t}}_{0}}^{t_{0}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\cdotp\!(-\mathbf{n})\mathrm{\ d}\vartheta\nonumber\\ &+\int_{t_{0}}^{t}\left(a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right) -a\left(S\left(\vartheta;x_{0},t_{0}\right),\vartheta\right)\right)\cdotp\!(-\mathbf{n})\mathrm{\ d}\vartheta. \end{align} (4.11)Due to (4.9), we can estimate the first integral as   \begin{align} \int_{{\tilde{t}}_{0}}^{t_{0}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right)\cdotp\!(-\mathbf{n})\,\,\mathrm{d}\vartheta \geqslant \frac{a_{\min}}{2}\left(t_{0}-{\tilde{t}}_{0}\right). \nonumber \end{align}As for the second integral in (4.11), we have   \begin{align}&\left|\int_{t_{0}}^{t}\big(a\big( S \big(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\big),\vartheta \big) -a\left(S\left(\vartheta;x_{0},t_{0}\right),\vartheta\right)\big)\cdot \!(-\mathbf{n}) \mathrm{\ d}\vartheta \right|\nonumber\\ &\quad\leqslant \widehat{T}L_{a} \sup_{\vartheta\in(t_{0},t)}\left|S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right) -S\left(\vartheta;x_{0},t_{0}\right)\right|\leqslant \widehat{T}L_{a} \left(e^{\widehat{T}L_{a}}-1\right)\|x-{\tilde{x}}\|, \end{align} (4.12)where La is the Lipschitz constant of a with respect to x. The last inequality in (4.12) follows from standard results on ordinary differential equations, namely continuous dependence of the solution on the initial condition—here we consider the ordinary differential equations defining S(⋅;⋅, ⋅) backward in time on the interval t0, t with ‘initial’ conditions x and $${\tilde{x}}$$ at time t. Finally, $$|\left (x-{\tilde{x}}\right )\cdotp \!(-\mathbf{n})|\leqslant \|x-{\tilde{x}}\|$$. Therefore, we get from (4.11),   \begin{align} \frac{a_{\min}}{2}\left|\mu_{1}(x,t)-\mu_{1}\left({\tilde{x}},t\right)\right| =\frac{a_{\min}}{2}\left|t_{0}-{\tilde{t}}_{0}\right| =\frac{a_{\min}}{2}\left(t_{0}-{\tilde{t}}_{0}\right) \leqslant \left(1+\widehat{T}L_{a} \left(e^{\widehat{T}L_{a}}-1\right)\right)\|x-{\tilde{x}}\|. \end{align} (4.13)Dividing by $$a_{\min }/2>0$$ gives Lipschitz continuity of μ1(⋅, t). Now we return to the case when x0 is a vertex of ∂$$\varOmega$$−. Reasoning as above, by choosing ε sufficiently small, we can ensure that $${\tilde{x}}_{0}\in \partial \varOmega ^{-}$$ is sufficiently close to x0, i.e $${\tilde{x}}_{0}$$ lies on one of the edges adjoining x0. Then we can again multiply (4.11) by −n, the normal to ∂$$\varOmega$$− at $${\tilde{x}}_{0}$$. Hence, $$\left ({\tilde{x}}-{\tilde{x}}_{0}\right )\cdotp \!(\mathbf{n})=0$$ will also be satisfied and we can proceed as in the previous case. Remark 4.2 We note that the Lipschitz constant $$L_{\mu _{1}}$$ of μ1 can be estimated from (4.13) by   \begin{align} L_{\mu_{1}}\leqslant \frac{2}{a_{\min}}\left(1+\widehat{T}L_{a} \left(e^{\widehat{T}L_{a}}-1\right)\right)\!. \nonumber \end{align}The Lipschitz constant Lμ of μ can then be obtained by multiplying by the constant factor from (4.4). Having established uniform Lipschitz continuity of μ1 in space, we can prove its Lipschitz continuity with respect to time. This implies the differentiability of μ1 with respect to x and t a.e. in QT; therefore the left-hand side of (4.1) is well defined a.e. in QT and this expression is equal to the derivative of μ1 along pathlines. This is important, as all the following considerations are thus justified. Lemma 4.3 Let a satisfy the assumptions of Lemma 4.1. Then μ1 is uniformly Lipschitz continuous with respect to time t: there exists Lt ⩾ 0 such that   \begin{align} \left|\mu(x,t)-\mu\left(x,{\tilde{t}}\right)\right|\leqslant L_{t}\left|t-{\tilde{t}}\right| \nonumber \end{align}for all x ∈ $$\varOmega$$ and $$t,{\tilde{t}}\in [0,T)$$. Proof. Let e.g. $${\tilde{t}}>t$$ and let (x0, t0) and $$\left ({\tilde{x}}_{0},{\tilde{t}}_{0}\right )$$ denote the origin of the pathlines passing through (x, t) and $$\left (x,{\tilde{t}}\right )$$, respectively. Therefore $$x=S(t;x_{0},t_{0}) =S\left ({\tilde{t}};{\tilde{x}}_{0},{\tilde{t}}_{0}\right )$$. Denote $${\tilde{x}}=S\left (t;{\tilde{x}}_{0},{\tilde{t}}_{0}\right )$$. Then   \begin{align} |x-{\tilde{x}}| =\left|\int_{t}^{{\tilde{t}}} a\left(S\left(\vartheta;{\tilde{x}}_{0},{\tilde{t}}_{0}\right),\vartheta\right) \mathrm{\ d}\vartheta\right| \leqslant \|a\|_{L^{\infty}(Q_{T})}|t-{\tilde{t}}|. \end{align} (4.14)Therefore, by (4.3),   \begin{align*}\left|\mu(x,t)-\mu\left(x,{\tilde{t}}\right)\right| \leqslant&\, |t-{\tilde{t}}|+|{\tilde{t}}_{0}-t_{0}| =\left|t-{\tilde{t}}\right|+\left|\mu(x,t)-\mu({\tilde{x}},t)\right| \leqslant \left|t-{\tilde{t}}\right|+L_{\mu}\left|x-{\tilde{x}}\right|\\ \leqslant&\, \left(1+L_{\mu} \|a\|_{L^{\infty}(Q_{T})}\right)\left|t-{\tilde{t}}\right| \end{align*}due to (4.14) and Lemma 4.1. This completes the proof. Example 4.4 In simple cases, the function μ1 can be explicitly written down. As a trivial example, we take the one-dimensional stationary flow field a(x, t) = x + 1 on $$\varOmega$$ = (0, 1) and the time interval $$(0,+\infty )$$. Then (4.2) can be easily solved to obtain   \begin{align} S(t;x_{0},t_{0})=(x_{0}+1)e^{t-t_{0}}-1. \nonumber \end{align}The (x,t)--plane is then separated into two regions separated by the pathline S(t;0, 0), which is the curve x = et − 1, i.e. $$t=\ln (x+1)$$. For points (x, t) beneath this curve, i.e. $$t\leqslant \ln (x+1)$$, we have t0 = 0, hence μ1(x, t) = t − t0 = t. For points (x, t) above the separation curve, we have x0 = 0, hence   \begin{align} x=S(t;0,t_{0})=e^{t-t_{0}}-1 \quad\Longrightarrow \quad \mu_{1}(x,t)=t-t_{0}=\ln(x+1). \nonumber \end{align}Altogether, we have   \begin{align} \mu_{1}(x,t)=\begin{cases} t & \text{if {$t\leqslant\ln(x+1)$},} \\ \ln(x+1)& \text{otherwise.} \end{cases} \nonumber \end{align}This function is globally bounded and globally Lipschitz continuous. We note that the standard exponential scaling trick which gives exponential growth of the error corresponds to taking μ1(x, t) = t for all x, t, which is an unbounded function. 5. DG method Having introduced the necessary tools for the analysis using the general exponential scaling argument, we can finally proceed to the DG method itself. Let $${\mathscr{T}}_{h}$$ be a triangulation of $$\varOmega$$, i.e. a partition of $$\overline{\varOmega }$$ into a finite number of closed simplices with mutually disjoint interiors. As usual with DG, $${\mathscr{T}}_{h}$$ need not be conforming, i.e. hanging nodes are allowed. For $$K\in{\mathscr{T}}_{h}$$ we set $$h_{K}=\textrm{diam}(K),\,\, h=\max _{K\in{\mathscr{T}}_{h}}h_{K}$$. For each $$K\in{\mathscr{T}}_{h}$$ we define its inflow and outflow boundaries by   \begin{align*} \partial K^{-}(t)&=\left\{x\in\partial K; a(x,t)\cdotp\mathbf{n}(x)<0\right\}\!,\\ \partial K^{+}(t)&=\{x\in\partial K; a(x,t)\cdotp\mathbf{n}(x)\ge 0\}, \nonumber \end{align*}where n(x) is the unit outer normal to ∂K. To simplify the notation, we shall usually omit the argument t in the following and write simply ∂K±. Here we remark that $${\mathscr{T}}_{h}$$ need not conform to ∂$$\varOmega$$−. In other words, the intersection of ∂$$\varOmega$$− and ∂K need not be an entire face (edge) of K. The triangulation $${\mathscr{T}}_{h}$$ defines the broken Sobolev space  \begin{align} H^{1}({\mathscr{T}}_{h})=\left\{v\in L^{2}(\varOmega);\, v|_{K} \in H^{1}(K)\quad \forall\, K\in{\mathscr{T}}_{h} \right\}\!. \nonumber \end{align}The approximate solution will be sought in the space of discontinuous piecewise polynomial functions   \begin{align} S_{h} = \left\{v_{h};\, v_{h}|_{K} \in P^{p}(K)\quad \forall\, K\in{\mathscr{T}}_{h}\right\}\!, \nonumber \end{align}where Pp(K) is the space of all polynomials on K of degree at most $$p\in \mathbb{N}$$. Given an element $$K\in{\mathscr{T}}_{h}$$, for vh ∈ Sh we define $$v_{h}^{-}$$ as the trace of vh on ∂K from the side of the element adjacent to K. Furthermore, on ∂K ∖ ∂$$\varOmega$$ we define the jump of vh as $$[v_{h}]=v_{h}-v_{h}^{-}$$, where vh is the trace from inside K. The DG formulation of (2.1) then reads as follows: we seek $$u_{h}\in C^{1}\left ([0,T);S_{h}\right )$$ such that $$u_{h}(0)={u_{h}^{0}}$$, an Sh-approximation of the initial condition u0, and for all t ∈ (0, T),   \begin{align} \Bigg(\frac{\partial u_{h}}{\partial t},v_{h}\Bigg)+b_{h}(u_{h},v_{h}) +c_{h}(u_{h},v_{h})=l_{h}(v_{h}) \quad\forall\, v_{h}\in S_{h}. \end{align} (5.1)Here bh, ch and lh are bilinear and linear forms, respectively, defined for $$u,v\in H^{1}({\mathscr{T}}_{h})$$ as follows. The bilinear advection form bh is   \begin{align} b_{h}(u,v)= \sum_{K\in{\mathscr{T}}_{h}}\int_{K} (a\cdotp \nabla\! u) v\,\mathrm{d} x -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[u]v\,\mathrm{d} S -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})uv\,\mathrm{d} S, \nonumber \end{align}the reaction form is defined by   \begin{align} c_{h}(u,v)= \int_{\varOmega} cu v\,\mathrm{d} x \nonumber \end{align}and lh is the right-hand side form  \begin{align} l_{h}(v)= -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})u_{D}v\,\mathrm{d} x . \nonumber \end{align}The definition of bh corresponds to the concept of upwinding which has already been applied in the context of DG in the original paper by Reed & Hill (1973). The right-hand side form lh is based on the weak enforcement of the Dirichlet boundary condition on ∂$$\varOmega$$− by penalization. For a detailed derivation of the specific forms of bh and lh as used in the paper, we refer e.g. to the study by Feistauer & Švadlenka (2004). 6. Analysis of the advection and reaction terms In this section, we prove the estimates of the bilinear forms bh and ch necessary for the error analysis. First, we need some standard results from approximation theory. 6.1 Auxiliary results In the following analysis we assume that the exact solution u is sufficiently regular, namely   \begin{align} u,u_{t}:=\frac{\partial u}{\partial t}\in L^{2}\left(H^{p+1}\right)\!. \nonumber \end{align}We consider a system $$\{{\mathscr{T}}_{h}\}_{h\in (0,\,h_{0})}$$, h0 > 0 of nonconforming triangulations of $$\varOmega$$ that are shape regular and satisfy the inverse assumption, cf. Ciarlet (1980). Under these assumptions we have the following standard results: Lemma 6.1 (Inverse inequality). There exists a constant CI > 0 independent of h, K such that for all $$K\in{\mathscr{T}}_{h}$$, and all v ∈ Pp(K),   \begin{align} |v|_{H^{1}(K)}\leqslant C_{I} h_{K}^{-1} \|v\|_{L^{2}(K)}. \nonumber \end{align} For v ∈ L2($$\varOmega$$) we denote by Πhv ∈ Sh the L2($$\varOmega$$)-projection of v onto Sh:   \begin{align} \left(\Pi_{h} v-v,\,\varphi_{h}\right)=0 \quad \forall\,\varphi_{h}\in S_{h}. \end{align} (6.1)Let ηh(t) = u(t) −Πhu(t) and ξh(t) =Πhu(t) − uh(t) ∈ Sh. Then we can write the error of the method as eh(t) := u(t) − uh(t) = ηh(t) + ξh(t). For simplicity, we shall usually drop the subscript h. We have the following standard approximation result; cf. Ciarlet (1980). Lemma 6.2 There exists a constant CΠ > 0 independent of h, K such that for all h ∈ (0, h0),   \begin{align*}\|\eta(t)\|_{L^{2}(K)}\leqslant&\, C_{\Pi} h^{p+1}|u(t)|_{{H^{p+1}}(K)},\\ |\eta(t)|_{H^{1}(K)}\leqslant&\, C_{\Pi} h^{p}|u(t)|_{{H^{p+1}}(K)},\\ \|\eta(t)\|_{L^{2}(\partial K)}\leqslant&\, C_{\Pi} h^{p+1/2}|u(t)|_{{H^{p+1}}(K)},\\ \Bigg\|\frac{\partial\eta(t)}{\partial t}\Bigg\|_{L^{2}(K)}\leqslant&\, C_{\Pi} h^{p+1}|u_{t}(t)|_{{H^{p+1}}(K)}. \end{align*} Throughout the paper, C will be a generic constant independent of h, t, T. In order to track the dependence of the estimates on the function μ, we will also assume that the generic constant C is independent of μ, particularly $$\mu _{\max }$$ and Lμ from (4.7), and we will explicitly track these dependencies. The dependence of the resulting estimates of Theorem 7.1 on μ will then be clarified in Section 7 for the specific construction of μ from Section 4. 6.2 Estimates of the advection and reaction terms In our analysis we will assume there exists a constant γ0 and a function $$\mu :Q_{T}\to \mathbb{R}$$ such that the ellipticity assumption (3.9) holds. Furthermore, we assume that μ satisfies (4.7), i.e. μ is non-negative, bounded and Lipschitz continuous in space. Such a function μ was constructed in Section 4. Similarly to Section 2, we wish to write $$\xi (x,t)=e^{\mu (x,t)}\tilde{\xi }(x,t)$$ and test the error equation with $$\phi (x,t)=e^{-\mu (x,t)}\tilde{\xi }(x,t)=e^{-2\mu (x,t)}\xi (x,t)$$ to obtain estimates for $$\tilde{\xi }$$. However, since ϕ(t) ∉ Sh this is not possible. One possibility is to test by Πhϕ(t) ∈ Sh and estimate the resulting difference Πhϕ(t) − ϕ(t). This is done in the stationary case in the study by Ayuso & Marini (2009) and the analysis is carried out under the assumption $$\mu \in W^{p+1,\infty }(\varOmega )$$. Such high regularity can be achieved by mollification of μ. However, this would be somewhat technical in the evolutionary case, as space–time smoothing would be required in which case the dependence of all constants on T must be carefully considered. Also QT is potentially an unbounded domain (for $$T=+\infty $$) which leads to technical difficulties. Here we carry out the analysis under the weaker assumption (4.7), i.e. $$\mu (t)\in W^{1,\infty }(\varOmega )$$. Lemma 6.3 Let μ satisfy assumptions (4.7). Let $$\phi (x,t)=e^{-\mu (x,t)}\tilde{\xi }(x,t)=e^{-2\mu (x,t)}\xi (x,t)$$, where ξ(t) ∈ Sh. Then there exists C independent of $$h,t,\xi ,\tilde{\xi }$$ and μ such that   \begin{align} \left\|\Pi_{h}\phi(t)-\phi(t)\right\|_{L^{2}(K)}\leqslant&\, CL_{\mu} e^{h_{K}L_{\mu}} h_{K} \max_{x\in K}e^{-\mu(x,t)} \|\tilde{\xi}(t)\|_{L^{2}(K)},\nonumber\\ \left\|\Pi_{h}\phi(t)-\phi(t)\right\|_{L^{2}(\partial K)}\leqslant&\, CL_{\mu} e^{h_{K}L_{\mu}} h^{1/2}_{K} \max_{x\in K}e^{-\mu(x,t)} \|\tilde{\xi}(t)\|_{L^{2}(K)}. \end{align} (6.2) Proof. Let xK be the centroid of K. On element K we introduce the constant μK(t) = μ(xK, t); then the function $$e^{-2\mu _{K}(t)}\xi (\cdot ,t)$$ lies in Pp(K), hence is fixed by the projection Πh. Therefore,   \begin{align*}\Pi_{h}\phi(t)-\phi(t)=&\, \Pi_{h}\left(e^{-2\mu(t)}\xi(t) -e^{-2\mu_{K}(t)}\xi(t)\right) -\left(e^{-2\mu(t)}\xi(t) -e^{-2\mu_{K}(t)}\xi(t)\right)\\ =&\,\Pi_{h} w(t)-w(t), \end{align*}where $$w(x,t)=(e^{-2\mu (x,t)} -e^{-2\mu (x_{K},t)})\xi (x,t)$$. Standard estimates of the interpolation error of Πh give   \begin{align} \left\|\Pi_{h}\phi(t)-\phi(t)\right\|_{L^{2}(K)}^{2} =\left\|\Pi_{h} w(t)-w(t)\right\|_{L^{2}(K)}^{2}\leqslant C{h_{K}^{2}}\left|w(t)\right|_{H^{1}(K)}^{2}. \end{align} (6.3)For the right-hand side seminorm we have   \begin{align*} |w(t)|_{H^{1}(K)}^{2}=&\,\int_{K}\left|\nabla\!\left(\left(e^{-2\mu(x,t)} -e^{-2\mu(x_{K},t)}\right)\xi(x,t)\right)\right|{}^{2}\,\mathrm{d} x \\ \leqslant&\, 2\int_{K}\left|\nabla\! e^{-2\mu(x,t)}\xi(x,t)\right| \mathrm{d} x +2\int_{K}\left|\left(e^{-2\mu(x,t)} -e^{-2\mu(x_{K},t)}\right)\nabla\!\xi(x,t)\right|{}^{2}\,\mathrm{d} x.\end{align*}If $$\mu (t)\in C^{1}(\overline{\varOmega })$$, by the mean value theorem   \begin{align} \left|e^{-2\mu(x_{K},t)}-e^{-2\mu(x,t)}\right|\leqslant h_{K} \left|\nabla\! e^{-2\mu(\zeta,t)}\right| = h_{K} e^{-2\mu(\zeta,t)}2\left|\nabla\mu(\zeta,t)\right| \nonumber \end{align}for some point ζ on the line between x and xK. Therefore, by the inverse inequality,   \begin{align*}|w(t)|_{H^{1}(K)}^{2}\leqslant&\, 2\int_{K} 4e^{-4\mu(x,t)} \left|\nabla\! \mu(x,t)\right|{}^{2}\left|\xi(x,t)\right|{}^{2}\,\mathrm{d} x +2\int_{K} {h_{K}^{2}}e^{-4\mu(\zeta,t)}4|\nabla\mu(\zeta,t)|^{2} |\nabla\!\xi(x,t)|^{2}\,\mathrm{d} x \\ \leqslant&\, 8 |\mu(t)|^{2}_{W^{1,\infty}}\int_{K} e^{-4\mu(x,t)}e^{2\mu(x,t)} |\tilde{\xi}(x,t)|^{2}\,\mathrm{d} x +8{C_{I}^{2}} \max_{x\in K}e^{-4\mu(x,t)}\left|\mu(t)\right|^{2}_{W^{1,\infty}} \|\xi(t)\|^{2}\\ \leqslant&\, 8 L_{\mu}^{2} \max_{x\in K} e^{-2\mu(x,t)} e^{2h_{K}L_{\mu}} \|\tilde{\xi}(t)\|^{2} +8{C_{I}^{2}} \max_{x\in K} e^{-2\mu(x,t)} e^{2h_{K}L_{\mu}} L_{\mu}^{2} \|\tilde{\xi}(t)\|^{2}.\end{align*}Substituting into (6.3) gives us the first inequality in (6.2) for $$\mu (t)\in C^{1}(\overline{\varOmega })$$. The case $$\mu (t)\in W^{1,\infty }(\varOmega )$$ follows by a density argument. The second inequality in (6.2) can be obtained similarly, only intermediately applying the trace inequality $$\|\xi (t)\|_{L^{2}(\partial K)}\leqslant Ch_{K}^{-1/2}\|\xi (t)\|_{L^{2}(K)}$$. Remark 6.4 The estimates of Lemma 6.3 are effectively $$\mathcal{O}(h_{K})\|\tilde{\xi }\|$$ and $$\mathcal{O}(h_{K}^{1/2})\|\tilde{\xi }\|$$ estimates, respectively. We note that the dependence of the estimates on Lμ is rather mild, effectively linear, due to the small factor hK in the exponent. Now we shall estimate individual terms in the DG formulation. Due to the consistency of the DG scheme, the exact solution u also satisfies (5.1). We subtract the formulations for u and uh to obtain the error equation   \begin{align} \begin{split} \left(\frac{\partial \xi}{\partial t},v_{h}\right) +\left(\frac{\partial \eta}{\partial t},v_{h}\right) + b_{h}(\xi,v_{h}) +b_{h}(\eta,v_{h}) + c_{h}(\xi,v_{h}) +c_{h}(\eta,v_{h}) =0 \end{split} \end{align} (6.4)for all vh ∈ Sh. As stated earlier we want to test (6.4) by $$\phi (x,t)=e^{-\mu (x,t)}\tilde{\xi }(x,t)$$; however, ϕ(t) ∉ Sh. We therefore set vh =Πhϕ(t) and estimate the difference using Lemma 6.3. We write (6.4) as   \begin{align}&\left(\frac{\partial \xi}{\partial t},\Pi_{h}\phi\right) +b_{h}(\xi,\phi) +b_{h}(\xi,\Pi_{h}\phi-\phi) +b_{h}(\eta,\phi) +b_{h}(\eta,\Pi_{h}\phi-\phi)\nonumber\\ &\qquad+\ c_{h}(\xi,\phi) +c_{h}(\xi,\Pi_{h}\phi-\phi) +c_{h}(\eta,\Pi_{h}\phi) +\left(\frac{\partial \eta}{\partial t},\Pi_{h}\phi\right)=0. \end{align} (6.5)We will estimate the individual terms of (6.5) in a series of lemmas. For this purpose, we introduce the following norm on a subset ω of ∂K or ∂$$\varOmega$$:   \begin{align} \|\,f\|_{a,\omega}=\left\|\sqrt{|a\cdotp\mathbf{n}|}\,f\right\|_{L^{2}(\omega)}, \nonumber \end{align}where n is the outer normal to ∂K or ∂$$\varOmega$$. We will usually omit the argument t to simplify the notation. In the following lemma we show that selected terms from (6.5) indeed give us the desired ellipticity similarly to (3.8) and (3.9) plus additional elliptic terms due to the use of the upwind flux. Lemma 6.5 Let $$\xi =e^{\mu }\tilde{\xi }, \phi =e^{-\mu }\tilde{\xi }$$ as above and let μ satisfy assumptions (3.9) and (4.7). Then   \begin{equation} \Bigg(\frac{\partial \xi}{\partial t},\Pi_{h}\phi\Bigg) +b_{h}(\xi,\phi)+c_{h}(\xi,\phi)\ge\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}\|^{2}+\gamma_{0}\|\tilde{\xi}\|^{2} +\frac{1}{2}\sum_{K\in{\mathscr{T}}_{h}}\bigg(\big\|[\tilde{\xi}]\big\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\big\|\tilde{\xi}\big \|_{a,\partial K\cap\partial\varOmega}^{2}\bigg). \end{equation} (6.6) Proof. Since ∂ξ/∂t ∈ Sh for each t, by definition (6.1) of Πh we have   \begin{align} \Bigg(\frac{\partial \xi}{\partial t},\Pi_{h}\phi\Bigg)= \Bigg(\frac{\partial \xi}{\partial t},\phi\Bigg) =\Bigg(e^{\mu}\frac{\partial \tilde{\xi}}{\partial t} +e^{\mu}\frac{\partial \mu}{\partial t}\tilde{\xi},e^{-\mu} \tilde{\xi}\Bigg) =\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}\|^{2} +\Bigg(\frac{\partial \mu}{\partial t}\tilde{\xi},\tilde{\xi}\Bigg). \end{align} (6.7)The reactive term satisfies   \begin{align} c_{h}(\xi,\phi)=\int_{\varOmega} c\xi\phi\,\mathrm{d} x =\int_{\varOmega} ce^{\mu}\tilde{\xi} e^{-\mu}\tilde{\xi}\,\mathrm{d} x =\int_{\varOmega} c\tilde{\xi}^{2}\,\mathrm{d} x . \end{align} (6.8)From the definition of bh, we get   \begin{align} b_{h}(\xi,\phi)=&\,\!\! \sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp(\nabla\!\mu\,\tilde{\xi}+\nabla\!\tilde{\xi})e^{\mu} e^{-\mu}\tilde{\xi}\,\mathrm{d} x -\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[e^{\mu}\tilde{\xi}]e^{-\mu}\tilde{\xi}\,\mathrm{d} S \nonumber\\ &-\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})e^{\mu}\tilde{\xi} e^{-\mu}\tilde{\xi}\,\mathrm{d} S \nonumber\\ =&\,\!\! \sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp(\nabla\!\mu\,\tilde{\xi}+\nabla\!\tilde{\xi})\tilde{\xi}\,\mathrm{d} x -\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\tilde{\xi}]\tilde{\xi}\,\mathrm{d} S -\!\!\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S .\end{align} (6.9)By Green’s theorem,   \begin{align} \int_{K} a\cdotp\!\nabla\!\tilde{\xi}\,\tilde{\xi}\,\mathrm{d} x =-\frac{1}{2}\int_{K} \operatorname{div} a\,\tilde{\xi}^{2}\,\mathrm{d} x +\frac{1}{2}\int_{\partial K}(a\cdotp\mathbf{n})\,\tilde{\xi}^{2}\,\mathrm{d} S . \end{align} (6.10)Splitting the last integral over the separate parts of ∂K, i.e. ∂K− ∖ ∂$$\varOmega$$, ∂K− ∩ $$ \partial \varOmega $$, ∂K+ ∖ ∂$$\varOmega$$ and ∂K+ ∩ $$ \partial \varOmega $$, and by substituting (6.10) into (6.9), we get   \begin{align}b_{h}(\xi,\phi)=&\,\!\! \sum_{K\in{\mathscr{T}}_{h}}\int_{K} \left(a\cdotp\nabla\!\mu-\tfrac{1}{2}\operatorname{div} a\right)\tilde{\xi}^{2}\,\mathrm{d} x \nonumber\\ &+\!\!\sum_{K\in{\mathscr{T}}_{h}}\bigg(-\frac{1}{2}\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})\left(\tilde{\xi}^{2}-2\tilde{\xi}\tilde{\xi}^{-}\right)\,\mathrm{d} S -\frac{1}{2}\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S \nonumber\\ &\qquad\qquad+\frac{1}{2}\int_{\partial K^{+}\setminus\partial\varOmega}(a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S +\frac{1}{2}\int_{\partial K^{+}\cap\partial\varOmega}(a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S \bigg).\end{align} (6.11)We note that   \begin{align} \sum_{K\in{\mathscr{T}}_{h}}\int_{\partial K^{+}\setminus\partial\varOmega}(a\cdotp \mathbf{n})\tilde{\xi}^{2}\,\mathrm{d} S =-\sum_{K\in{\mathscr{T}}_{h}}\int_{\partial K^{-}\setminus\partial\varOmega}(a\cdotp \mathbf{n})(\tilde{\xi}^{-})^{2}\,\mathrm{d} S . \end{align} (6.12)Using the facts that $$\tilde{\xi }^{2}-2\tilde{\xi }\tilde{\xi }^{-}+(\tilde{\xi }^{-})^{2}= [\tilde{\xi }]^{2}$$ and − a ⋅ n = |a⋅n| on ∂K− and a ⋅ n = |a⋅n| on ∂K+, we get (6.6) by substituting (6.12) into (6.11) and applying assumption (3.9) in the resulting interior terms along with (6.7), (6.8). The following series of lemmas deals with estimating the remaining terms in (6.5). These estimates more or less follow standard procedures with several important differences pertaining to the specific choice of test functions ϕ and Πhϕ − ϕ. These include the inability to apply the inverse inequality to $$\tilde{\xi }\notin S_{h}$$ but only to ξ ∈ Sh and Lemma 6.3. In order to keep the main body of the text less cluttered, we leave the proofs of the following lemmas to the Appendix. Lemma 6.6 The advection terms satisfy   \begin{align} \hskip-30pt|b_{h}(\xi,\Pi_{h}\phi-\phi)|\leqslant C(1+L_{\mu})^{2} e^{4hL_{\mu}} h\|\tilde{\xi}\|^{2}+\frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\bigg(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\bigg), \end{align} (6.13)  \begin{align} {\hskip30pt}|b_{h}(\eta,\phi)|\leqslant C(1+L_{\mu})^{2}h\|\tilde{\xi}\|^{2} + Ch^{2p+1}|u(t)|_{H^{p+1}}^{2} +\frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\!\!\!\left(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\right), \end{align} (6.14)  \begin{align} |b_{h}(\eta,\Pi_{h}\phi-\phi)|\leqslant CL_{\mu}^{2} e^{2hL_{\mu}} h\|\tilde{\xi}\| +Ch^{2p+1}|u(t)|_{H^{p+1}}^{2}. \hskip150pt\end{align} (6.15) Lemma 6.7 The reaction terms satisfy   \begin{align*}|c_{h}(\xi,\Pi_{h}\phi-\phi)|\leqslant&\, CL_{\mu} e^{2hL_{\mu}}h\|\tilde{\xi}\|^{2},\\ |c_{h}(\eta,\Pi_{h}\phi)|\leqslant&\, C\big(1+L_{\mu} e^{hL_{\mu}}\big)^{2} h\|\tilde{\xi}\|^{2} +Ch^{2p+1}|u(t)|_{H^{p+1}}^{2}.\end{align*} Lemma 6.8 The evolutionary term satisfies   \begin{align} \Bigg|\Bigg(\frac{\partial \eta}{\partial t},\Pi_{h}\phi\Bigg)\Bigg| \leq C\big(1+L_{\mu} e^{hL_{\mu}}\big)^{2} h\|\tilde{\xi}\|^{2} +Ch^{2p+1}|u_{t}(t)|_{H^{p+1}}^{2}. \nonumber \end{align} 7. Error analysis Finally, we come to the error analysis. The starting point is the error identity (6.5) to which we apply the derived estimates of its individual terms. Theorem 7.1 Let there exist a function $$\mu :Q_{T}\to [0,\mu _{\max }]$$ for some constant $$\mu _{\max }$$, such that $$\mu (t)\in W^{1,\infty }(\varOmega )$$ and let there exist a constant γ0 > 0 such that the coefficients of (2.1) satisfy $$\mu _{t} +a\cdotp \nabla \!\mu +c-\tfrac{1}{2}\textrm{div} a\ge \gamma _{0}>0$$ a.e. in QT. Let the initial condition $${u_{h}^{0}}$$ satisfy $$\left \|\Pi _{h} u^{0}-{u_{h}^{0}}\right \|\leqslant Ch^{p+1/2}\left |u^{0}\right |_{H^{p+1}}$$. Then there exists a constant C independent of μ, h and T such that for h sufficiently small the error of the DG scheme (5.1) satisfies   \begin{align}&\max_{t\in[0,T]}\|e_{h}(t)\| +\sqrt{\gamma_{0}}\|e_{h}\|_{L^{2}(Q_{T})} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|[e_{h}(\vartheta)]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|e_{h}(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ce^{\mu_{\max}} h^{p+1/2}\left(|u^{0}|{}_{H^{p+1}}+|u|_{L^{2}\left({H^{p+1}}\right)} +\left|u_{t}\right|{}_{L^{2}\left({H^{p+1}}\right)} +h^{1/2}|u|_{L^{\infty}\left(H^{p+1}\right)}\right)\!.\end{align} (7.1) Proof. Applying Lemmas 6.5–6.8 to (6.5), multiplying by 2 for convenience and collecting similar terms, we get for all t,  \begin{align}&\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}(t)\|^{2} +2\gamma_{0}\|\tilde{\xi}(t)\|^{2} +\frac{1}{2}\sum_{K\in{\mathcal{T}}_{h}}\left(\|[\tilde{\xi}(t)]\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\tilde{\xi}(t)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\nonumber\\ &\quad\leqslant C_{1}(1+L_{\mu})^{2}e^{4hL_{\mu}}h\|\tilde{\xi}(t)\|^{2} +C_{2}h^{2p+1}\left(\left|u(t)\right|{}_{H^{p+1}}^{2}+|u_{t}(t)|_{H^{p+1}}^{2}\right)\!, \end{align} (7.2)where we have used $$C_{1} (1+L_{\mu })^{2}e^{4hL_{\mu }}$$ as a common upper bound for the constants from Lemmas 6.6–6.8 at the $$\mathcal{O}(h)\|\tilde{\xi }(t)\|^{2}$$ terms. Here C1 and C2 are independent of h, t, T, μ and $$\tilde{\xi }$$. Now if h is sufficiently small so that   \begin{align} C_{1}(1+L_{\mu})^{2}e^{4hL_{\mu}}h \leqslant\gamma_{0}, \end{align} (7.3)then the first right-hand side term can be absorbed by the left-hand side term $$2\gamma _{0}\|\tilde{\xi }(t)\|^{2}$$:   \begin{align*} &\frac{\mathrm{d}}{\mathrm{d} t}\|\tilde{\xi}(t)\|^{2} +\gamma_{0}\|\tilde{\xi}(t)\|^{2} +\frac{1}{2}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}(t)]\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\tilde{\xi}(t)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\\ &\quad\leqslant C_{2}h^{2p+1}\left(\left|u(t)\right|{}_{H^{p+1}}^{2}+|u_{t}(t)|_{H^{p+1}}^{2}\right)\!. \end{align*}Substituting t = ϑ and integrating over (0, t) give us   \begin{align}&\|\tilde{\xi}(t)\|^{2} +\gamma_{0}{\int_{0}^{t}}\|\tilde{\xi}(\vartheta)\|^{2}\mathrm{d} \vartheta +\frac{1}{2}{\int_{0}^{t}}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}(\vartheta)]\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\tilde{\xi}(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\nonumber\\ &\quad\leqslant \|\tilde{\xi}(0)\|^{2} +C_{2}h^{2p+1}\left(|u|_{L^{2}(0,t;{H^{p+1}})}^{2} +|u_{t}|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2}\right)\nonumber\\ &\quad\leqslant Ch^{2p+1}\left(|u^{0}|{}_{H^{p+1}}^{2} +|u|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2} +|u_{t}|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2}\right)\!, \end{align} (7.4)since the assumptions give   \begin{align} \|\tilde{\xi}(0)\|^{2}\leqslant \left\|\xi(0)\right\|^{2} =\|\Pi_{h} u^{0}-{u_{h}^{0}}\|^{2} \leqslant Ch^{2p+1}|u^{0}|{}_{H^{p+1}}^{2}. \nonumber \end{align}Now we reformulate estimate (7.4) as an estimate of ξ instead of $$\tilde{\xi }$$. Because $$\tilde{\xi }=e^{-\mu }\xi $$, we can estimate   \begin{align} \|\tilde{\xi}(t)\|^{2}\geqslant \min_{Q_{T}}e^{-2\mu(x,\vartheta)} \left\|\xi(t)\right\|^{2} =e^{-2\max_{Q_{T}}\mu(x,\vartheta)}\left\|\xi(t)\right\|^{2} =e^{-2\mu_{\max}} \|\xi(t)\|^{2} \nonumber \end{align}and similarly for the remaining left-hand side norms in (7.4). Substituting into (7.4) and multiplying by $$e^{2\mu _{\max }}$$ gives us   \begin{align*} &\|\xi(t)\|^{2} +\gamma_{0}\|\xi\|_{L^{2}\left(0,t;L^{2}\right)}^{2} +\frac{1}{2}{\int_{0}^{t}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|\left[\xi(\vartheta)\right]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\xi(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\nonumber\\ &\quad\leqslant Ch^{2p+1} e^{2\mu_{\max}}\left(|u^{0}|{}_{H^{p+1}}^{2} +|u|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2} +|u_{t}|_{L^{2}\left(0,t;{H^{p+1}}\right)}^{2}\right). \end{align*}By taking the square root of the last inequality and taking the maximum over all t ∈ [0, T], we get   \begin{align}&\max_{t\in[0,T]}\|\xi(t)\| +\sqrt{\gamma_{0}}\|\xi\|_{L^{2}(L^{2})} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\big\|[\xi(\vartheta)]\big\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\xi(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ch^{p+1/2} e^{\mu_{\max}}\left(\big|u^{0}\big|_{H^{p+1}} +|u|_{L^{2}({H^{p+1}})} +|u_{t}|_{L^{2}\left({H^{p+1}}\right)}\right) \end{align} (7.5)which is estimate (7.1) for the discrete part ξ of the error. Lemma 6.2 gives a similar inequality for η:   \begin{align} &\max_{t\in[0,T]}\|\eta(t)\| +\sqrt{\gamma_{0}}\|\eta\|_{L^{2}\left(L^{2}\right)} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|\left[\eta(\vartheta)\right]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|\eta(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{\ d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ch^{p+1/2} \left(h^{1/2}|u|_{L^{\infty}\left(H^{p+1}\right)} +(h^{1/2}+1)|u|_{L^{2}\left({H^{p+1}}\right)}\right)\!. \end{align} (7.6)Since eh = ξ + η, the triangle inequality along with (7.5) and (7.6) gives us the statement of the theorem. We now interpret the results of Theorem 7.1 given the specific construction of the scaling function μ from Section 4. If the assumptions of Lemma 4.1 are satisfied, the function μ1 exists and has suitable regularity. If we choose γ0 (this can be chosen e.g. as γ0 = 1 for simplicity), then μ is constructed by the scaling (4.4). As μ1 is bounded by $$\widehat{T}$$, the exponential factor in (7.1) becomes   \begin{align} e^{\mu_{\max}}=e^{\widehat{T}(|\inf_{Q_{T}}(c-\frac{1}{2}\operatorname{div} a)^{-}|+\gamma_{0})}\!. \nonumber \end{align} Now we turn to the requirement ‘h is sufficiently small’ in Theorem 7.1. Specifically, (7.3) must be satisfied. If we denote $$A:=C_{1}\left (1+L_{\mu }\right )^{2}/\gamma _{0}$$ and B := 4Lμ, then condition (7.3) reads   \begin{align} Ae^{Bh}h \leqslant 1. \end{align} (7.7)It can be easily seen that if $$h\leqslant \tfrac{1}{2}\min \left \{A^{-1},B^{-1}\right \}$$, then   \begin{align} Ae^{Bh}h\leqslant \tfrac{1}{2}e^{1/2}\leqslant 1, \nonumber \end{align}i.e. (7.7) holds. Taking the definitions of A and B, we get the restriction on h,  \begin{align} h\leqslant \frac{1}{2}\min\left\{\frac{\gamma_{0}}{C_{1}\left(1+L_{\mu}\right)^{2}}, \frac{1}{4L_{\mu}}\right\}\!. \end{align} (7.8)If we take e.g. γ0 = 1 as mentioned and assume $$C_{1}\geqslant 1$$ (if not, C1 can be increased to 1 in (7.2)), condition (7.8) reduces to   \begin{align} h\leqslant \frac{1}{2C_{1}\left(1+L_{\mu}\right)^{2}}, \end{align} (7.9)where Lμ is the Lipschitz constant of μ from Lemma 4.1, which can be estimated as in Remark 4.2 using data from the continuous problem (2.1) such as La, the spatial Lipschitz constant of a. Therefore the condition (7.8) on the mesh size and the final error estimate can be formulated without μ, only in terms of $$\widehat{T}$$ and the data of the problem. Thus we can state the following corollary derived from Theorem 7.1 using the construction of μ from Section 4. Corollary 7.2 Let $$a\in L^{\infty }(Q_{T})$$ be continuous with respect to time and Lipschitz continuous with respect to space. Let there exist a constant $$a_{\min }>0$$ such that $$-a(x,t)\cdot \mathbf{n}\geqslant a_{\min }$$ for all x ∈ ∂$$\varOmega$$−, t ∈ [0, T). Let the time any particle carried by the flow field a(⋅, ⋅) spends in $$\varOmega$$ be uniformly bounded by $$\widehat{T}$$. Let the initial condition $${u_{h}^{0}}$$ satisfy $$\left \|\Pi _{h} u^{0}-{u_{h}^{0}}\right \|\leqslant Ch^{p+1/2}\left |u^{0}\right |_{H^{p+1}}$$. Then there exist constants c > 0 and C independent of h and T such that for $$h\leqslant c$$, the error of the DG scheme (5.1) satisfies   \begin{align}&\max_{t\in[0,T]}\left\|e_{h}(t)\right\| +\|e_{h}\|_{L^{2}(Q_{T})} +\left(\frac{1}{2}{\int_{0}^{T}}\sum_{K\in{\mathscr{T}}_{h}}\left(\left\|[e_{h}(\vartheta)]\right\|_{a,\partial K^{-}\setminus\partial\varOmega}^{2} +\|e_{h}(\vartheta)\|_{a,\partial K\cap\partial\varOmega}^{2}\right)\mathrm{\ d}\vartheta\right)^{1/2}\nonumber\\ &\quad\leqslant Ce^{\widehat{T}(|\inf_{Q_{T}}(c-\frac{1}{2}\operatorname{div} a)^{-}|+1)} h^{p+1/2}\left(|u^{0}|{}_{H^{p+1}}+|u|_{L^{2}\left({H^{p+1}}\right)} +|u_{t}|_{L^{2}\left({H^{p+1}}\right)} +h^{1/2}|u|_{L^{\infty}\left(H^{p+1}\right)}\right)\!. \end{align} (7.10) Remark 7.3 As mentioned in Section 4, estimate (7.10) is exponential in the maximal ‘lifetime’ of particles $$\widehat{T}$$, instead of the final time T, which would be expected when applying Gronwall’s lemma directly. Since we assume the maximal particle lifetime to be finite, the exponential factor in (7.10) remains uniformly bounded independent of time T, which can be infinite. 8. Conclusion and future work In this paper we derived a priori error estimates for a linear advection–reaction problem with inlet and outlet boundary conditions in the $$L^{\infty }\left (L^{2}\right )$$- and $$L^{2}\left (L^{2}\right )$$-norms. Unlike previous works, the analysis was performed without the usual ellipticity assumption $$c-\frac{1}{2}\textrm{div} a\geqslant 0$$. This is achieved by applying a general exponential scaling transformation in space and time to the exact and discrete solutions of the problem. We considered the case when the time spent by particles carried by the flow field inside the spatial domain $$\varOmega$$ is uniformly bounded by some $$\widehat{T}$$. The resulting error estimates are of the order Chp+1/2, where C depends exponentially on $$\widehat{T}$$ (which is a constant) and not on the final time $$T\leqslant +\infty $$, as would be expected from the use of Gronwall’s inequality. Effectively, due to the exponential scaling, we apply Gronwall’s lemma in the Lagrangian setting along pathlines, which exist only for time at most $$\widehat{T}$$, and not in the usual Eulerian sense. As for future work, we plan to extend the analysis to fully discrete DG schemes with discretization in time. Furthermore, we wish to extend the analysis to nonlinear convective problems, following the arguments of Zhang & Shu (2004) and Kučera (2014) to obtain error estimates without the exponential dependence on time of the error estimates. Funding The work of V. Kučera was supported by the J. William Fulbright Commission in the Czech Republic and research project No. 17-01747S of the Czech Science Foundation. The work of C.-W. Shu was supported by DOE grant DE-FG02-08ER25863 and NSF grants DMS-1418750 and DMS-1719410. References Ayuso, B. & Marini, L. D. ( 2009) Discontinuous Galerkin methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. , 47, 1391– 1420. Google Scholar CrossRef Search ADS   Benzoni-Gavage, S. & Serre, D. ( 2007) Multi-Dimensional Hyperbolic Partial Differential Equations: First-Order Systems and Applications. Oxford Mathematical Monographs . Oxford: Oxford University Press. Ciarlet, P. G. ( 1980) The finite element method for elliptic problems . Studies in Mathematics and its Applications . Amsterdam, New York: North-Holland. Devinatz, A., Ellis, R. & Friedman, A. ( 1974) The asymptotic behavior of the first real eigenvalue of second order elliptic operators with a small parameter in the highest derivatives, II. Indiana Univ. Math. J. , 23, 991– 1011. Google Scholar CrossRef Search ADS   Evans, L. C. ( 2010) Partial Differential Equations . Providence, RI: American Mathematical Society. Google Scholar CrossRef Search ADS   Feistauer, M. & Švadlenka, K. ( 2004) Discontinuous Galerkin method of lines for solving nonstationary singularly perturbed linear problems. J. Numer. Math. , 12, 97– 117. Google Scholar CrossRef Search ADS   Johnson, C. & Pitkäranta, J. ( 1986) An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. Math. Comp. , 46, 1– 26. Google Scholar CrossRef Search ADS   Kučera, V. ( 2014) On diffusion-uniform error estimates for the DG method applied to singularly perturbed problems. IMA J. Numer. Anal. , 34, 820– 861. Google Scholar CrossRef Search ADS   Larsson, S. & Thomée, V. ( 2003) Partial Differential Equations with Numerical Methods. Texts in Applied Mathematics . Berlin: Springer. Lesaint, P. & Raviart, P. A. ( 1974) On a finite element method for solving the neutron transport equation. Mathematical Aspects of Finite Elements in Partial Differential Equations  ( C. de Boor ed). New York: Academic Press, pp. 89– 145. Google Scholar CrossRef Search ADS   Mizohata, S. ( 1973) The Theory of Partial Differential Equations . Cambridge: Cambridge University Press. Nävert, U. ( 1982) A finite element method for convection-diffusion problems. Ph.D. Thesis , Chalmers University of Technology. Reed, W. H. & Hill, T. ( 1973) Triangular mesh methods for the neutron transport equation. Los Alamos Report LA-UR-73–479. Roos, H.-G., Stynes, M. & Tobiska, L. ( 2008) Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems. Springer Series in Computational Mathematics , vol. 24. Berlin: Springer. Zhang, Q. & Shu, C.-W. ( 2004) Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws. SIAM J. Numer. Anal. , 42, 641– 666. Google Scholar CrossRef Search ADS   Appendix A.1 Proof of Lemma 6.6 A.1.1 Proof of estimate (6.13) We estimate the terms of bh(ξ, Πhϕ − ϕ) over the interiors and boundaries of elements separately. Let $${\Pi _{h}^{1}}$$ be the L2($$\varOmega$$)-projection onto the space of discontinuous piecewise linear polynomials on $${\mathcal{T}}_{h}$$. Since on each $$K\in{\mathcal{T}}_{h}$$ it holds that $$\nabla \!\tilde{\xi }|_{K}\in P^{p-1}(K)$$, then $${\Pi _{h}^{1}} a\cdotp \nabla \!\tilde{\xi }\in S_{h}$$. Hence, due to (6.1), we have   \begin{align} \sum_{K}\int_{K}{\Pi_{h}^{1}} a\cdotp\nabla\!\tilde{\xi}(\Pi_{h}\phi-\phi)\,\mathrm{d} x =0. \nonumber \end{align}Due to standard approximation results, we have $$\left \|a-{\Pi _{h}^{1}} a\right \|_{L^{\infty }(K)}\leqslant Ch_{K}|a|_{W^{1,\infty }(K)}$$, thus we can estimate the interior terms of bh(ξ, Πhϕ − ϕ) as   \begin{align} \sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp\nabla\!\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} x =&\,\sum_{K\in{\mathscr{T}}_{h}}\int_{K} \left(a-{\Pi_{h}^{1}} a\right)\cdotp\nabla\!\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} x \nonumber\\ \leqslant&\,\sum_{K\in{\mathscr{T}}_{h}} Ch_{K}|a|_{W^{1,\infty}}C_{I}h_{K}^{-1}\|\xi\|_{L^{2}(K)} \left\|\Pi_{h}\phi-\phi\right\|_{L^{2}(K)}\nonumber\\ \leqslant&\, C\sum_{K\in{\mathscr{T}}_{h}} \max_{x\in K}e^{\mu(x,t)} \|\tilde{\xi}\|_{L^{2}(K)} CL_{\mu} e^{hL_{\mu}} h \max_{x\in K}e^{-\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)}\nonumber\\\leqslant&\, CL_{\mu} e^{2hL_{\mu}}h\|\tilde{\xi}\|^{2},\end{align} (A.1)due to the inverse inequality, Lemma 6.3 and the estimate $$\max _{x\in K}e^{\mu (x,t)}\max _{x\in K}e^{-\mu (x,t)}\leqslant e^{L_{\mu } h_{K}}$$. For the boundary terms of bh(ξ, Πhϕ − ϕ) we get   \begin{align*} & -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\xi](\Pi_{h}\phi-\phi)\,\mathrm{d} S -\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} S \\ &\ \leqslant\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} \!|a\cdotp \mathbf{n}|\max_{x\in K}(e^{\mu(x,t)})\big|[\tilde{\xi}]\big||\Pi_{h}\phi-\phi| \,\mathrm{d} S \!+\!\!\sum_{K\in{\mathscr{T}}_{h}}\! \int_{\partial K^{-}\cap\partial\varOmega} \!|a\cdotp \mathbf{n}|\max_{x\in K}(e^{\mu(x,t)})|\tilde{\xi}||\Pi_{h}\phi\!-\phi|\,\mathrm{d} S \\ &\ \leqslant \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}} \int_{\partial K^{-}\setminus\partial\varOmega} \!|a\cdotp \mathbf{n}|[\tilde{\xi}]^{2}\,\mathrm{d} S \!+\!\frac{1}{8}\!\sum_{K\in{\mathscr{T}}_{h}} \!\int_{\partial K^{-}\cap\partial\varOmega} |a\cdotp \mathbf{n}|\tilde{\xi}^{2}\,\mathrm{d} S\!+\!2\!\sum_{K\in{\mathscr{T}}_{h}}\!\max_{x\in K}e^{2\mu(x,t)}\!\!\int_{\partial K}\!|a\cdotp\mathbf{n}||\Pi_{h}\phi\!-\!\phi|^{2}\mathrm{d} S \\ &\ \leqslant \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\right) + CL_{\mu}^{2} e^{2hL_{\mu}}\sum_{K}\max_{x\in K}e^{2\mu(x,t)} h_{K} \max_{x\in K}e^{-2\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)}^{2} \end{align*}by Young’s inequality and Lemma 6.3. Again we estimate $$\max _{x\in K}e^{2\mu (x,t)}\max _{x\in K}e^{-2\mu (x,t)}\leqslant e^{2L_{\mu } h_{K}}$$ in the last inequality, which completes the proof after combining with (A.1). □ 8.1.2 Proof of estimate (6.14) We have by Green’s theorem   \begin{align} b_{h}(\eta,\phi)=&\, \sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K} a\cdotp\mathbf{n}\eta\phi\,\mathrm{d} S -\int_{K} (\operatorname{div} a)\eta\phi\,\mathrm{d} x -\int_{K} a\cdotp\nabla\!\phi\eta\,\mathrm{d} x \right.\nonumber\\ &\qquad\quad\left.-\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\eta]\phi\,\mathrm{d} S -\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\eta\phi\,\mathrm{d} S \right).\end{align} (A.2)The first integral over K can be estimated as   \begin{align} -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} (\operatorname{div} a)\eta\phi\,\mathrm{d} x \leqslant Ch^{p+1}|u(t)|_{H^{p+1}}\max_{x\in \varOmega}e^{-\mu(x,t)}\|\tilde{\xi}\| \leqslant Ch^{p+1}|u(t)|_{H^{p+1}}\|\tilde{\xi}\|, \nonumber \end{align}because $$\mu \geqslant 0$$. Since ϕ = e−2μξ, we get for the second integral over K in (A.2),   \begin{align} \begin{split} -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} a\cdotp\nabla\!\phi\,\eta\,\mathrm{d} x =\sum_{K\in{\mathscr{T}}_{h}}\int_{K} 2a\cdotp\nabla\!\mu\,e^{-2\mu}\xi\eta\,\mathrm{d} x -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} e^{-2\mu}a\cdotp\nabla\!\xi\,\eta\,\mathrm{d} x . \end{split} \end{align} (A.3)The first right-hand side term in (A.3) can be estimated by   \begin{align*} \sum_{K\in{\mathscr{T}}_{h}}\int_{K} 2a\cdotp\nabla\!\mu\,e^{-2\mu}\xi\eta\,\mathrm{d} x\! =&\,\!\sum_{K\in{\mathscr{T}}_{h}}\!\int_{K} \!2a\cdotp\!\nabla\!\mu\,e^{-\mu}\tilde{\xi}\eta\,\mathrm{d} x\!\leqslant\!\!\sum_{K\in{\mathscr{T}}_{h}} CL_{\mu}\max_{x\in K}e^{-\mu(x,t)}h^{p+1}|u(t)|_{{H^{p+1}}(K)} \|\tilde{\xi}\|_{L^{2}(K)} \nonumber\\ \leqslant&\, CL_{\mu}^{2} h\|\tilde{\xi}\|^{2} +C h^{2p+1}|u(t)|_{{H^{p+1}}}^{2}. \end{align*}The second right-hand side term in (A.3) can be estimated similarly to (A.1), by the definition of η:   \begin{align*} -\sum_{K\in{\mathscr{T}}_{h}}\int_{K} e^{-2\mu}a\cdotp\nabla\!\xi\,\eta\,\mathrm{d} x =&\,\sum_{K\in{\mathscr{T}}_{h}}\int_{K} \left({\Pi_{h}^{1}}(e^{-2\mu}a)-e^{-2\mu}a\right) \cdotp\nabla\!\xi\eta\,\mathrm{d} x \\ \leqslant&\,\sum_{K\in{\mathscr{T}}_{h}} Ch_{K}|e^{-2\mu}a|_{W^{1,\infty}}C_{I}h_{K}^{-1}\|\xi\|_{L^{2}(K)} h^{p+1}|u(t)|_{H^{p+1}(K)}\\\leqslant&\, C(1+L_{\mu})\max_{x\in \varOmega}e^{-\mu(x,t)}h^{p+1}|u(t)|_{{H^{p+1}}}\|\tilde{\xi}\|\\ \leqslant&\, C(1+L_{\mu})^{2}h\|\tilde{\xi}\|^{2} +Ch^{2p+1}|u(t)|_{{H^{p+1}}}^{2}, \end{align*}since $$\left |e^{-2\mu }a\right |_{W^{1,\infty }}= \|e^{-2\mu }\textrm{div} a -2\nabla \!\mu \, e^{-2\mu }a \|_{L^{\infty }}\leqslant C(1+L_{\mu })$$. As for the boundary terms in (A.2), we can split the integral over ∂K into integrals over the separate parts ∂K− ∖ ∂$$\varOmega$$, ∂K− ∩ ∂$$\varOmega$$, ∂K+ ∖ ∂$$\varOmega$$ and ∂K+ ∩ ∂$$\varOmega$$, similarly to the proof of Lemma 6.5. Thus, several terms are canceled out:   \begin{align}&\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S -\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\eta]\phi\,\mathrm{d} S -\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\eta\phi\,\mathrm{d} S \right)\nonumber\\ &\quad=\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K^{+}\setminus\partial\varOmega} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S +\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})\eta^{-}\phi\,\mathrm{d} S \right)\nonumber\\ &\quad=\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp\mathbf{n})\eta^{-}[\phi]\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} (a\cdotp\mathbf{n})\eta\phi\,\mathrm{d} S \right)\end{align} (A.4)using a similar identity to (6.12). Finally, we can use Young’s inequality to estimate (A.4) further as   \begin{align*}\ldots\leqslant&\,\sum_{K\in{\mathscr{T}}_{h}}\left(\int_{\partial K^{-}\setminus\partial\varOmega} |a\cdotp\mathbf{n}||\eta^{-}|e^{-\mu}\big|[\tilde{\xi}]\big|\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} |a\cdotp\mathbf{n}||\eta|e^{-\mu}|\tilde{\xi}|\,\mathrm{d} S \right)\\ \leqslant&\, \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\bigg(\int_{\partial K^{-}\setminus\partial\varOmega} |a\cdotp\mathbf{n}|[\tilde{\xi}]^{2}\,\mathrm{d} S +\int_{\partial K^{+}\cap\partial\varOmega} |a\cdotp\mathbf{n}|\tilde{\xi}^{2}\,\mathrm{d} S \bigg) +2\sum_{K}\int_{\partial K}|a\cdotp\mathbf{n}|\eta^{2}e^{-2\mu}\,\mathrm{d} S \\ \leqslant&\, \frac{1}{8}\sum_{K\in{\mathscr{T}}_{h}}\left(\|[\tilde{\xi}]\|_{a,\partial K^{-}\!\setminus\partial\varOmega}^{2} +\|\tilde{\xi}\|_{a,\partial K\cap\partial\varOmega}^{2}\right) +Ch^{2p+1}|u(t)|_{H^{p+1}}^{2}. \nonumber \end{align*}The proof is completed by gathering all the above estimates of the individual terms of bh(η, ϕ). □ 8.1.3 Proof of estimate (6.15) We use Lemmas 6.2 and 6.3 to estimate   \begin{align*} &b_{h}(\eta,\Pi_{h}\phi-\phi)\\ &\quad=\sum_{K\in{\mathscr{T}}_{h}}\bigg(\int_{K} (a\cdotp \nabla\! \eta) (\Pi_{h}\phi-\phi)\,\mathrm{d} x -\!\int_{\partial K^{-}\setminus\partial\varOmega} (a\cdotp \mathbf{n})[\eta](\Pi_{h}\phi\!-\phi)\,\mathrm{d} S\!-\!\!\int_{\partial K^{-}\cap\partial\varOmega} (a\cdotp \mathbf{n})\eta (\Pi_{h}\phi\!-\!\phi)\,\mathrm{d} S \bigg)\\ &\quad\leqslant Ch^{p}|u(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}} h \max_{x\in\varOmega}e^{-\mu(x,t)}\|\tilde{\xi}\|+Ch^{p+1/2}|u(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}}h^{1/2} \max_{x\in\varOmega}e^{-\mu(x,t)}\|\tilde{\xi}\|. \end{align*}The application of Young’s inequality completes the proof. □ 8.2 Proof of Lemma 6.7 Lemma 6.3 gives us   \begin{align*} \sum_{K\in{\mathscr{T}}_{h}}\int_{K} c\xi(\Pi_{h}\phi-\phi)\,\mathrm{d} x \leqslant&\, C\sum_{K\in{\mathscr{T}}_{h}}\max_{x\in K}e^{\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)} L_{\mu} e^{hL_{\mu}}h_{K} \max_{x\in K}e^{-\mu(x,t)}\|\tilde{\xi}\|_{L^{2}(K)}\\ \leqslant&\, CL_{\mu} e^{2hL_{\mu}} h\|\tilde{\xi}\|^{2}. \nonumber \end{align*}As for the second estimate, we write ch(η, Πhϕ) = ch(η, ϕ) + ch(η, Πhϕ − ϕ) and estimate by Lemmas 6.2 and 6.3:   \begin{align} \begin{split} c_{h}(\eta,\phi) =&\,\int_{\varOmega} c\eta\phi\,\mathrm{d} x =\int_{\varOmega} c\eta e^{-\mu}\tilde{\xi}\,\mathrm{d} x \leqslant Ch^{p+1}|u(t)|_{H^{p+1}}\|\tilde{\xi}\|,\\ c_{h}(\eta,\Pi_{h}\phi-\phi) =&\,\int_{\varOmega} c\eta (\Pi_{h}\phi-\phi)\,\mathrm{d} x \leqslant Ch^{p+1}|u(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}}h \|\tilde{\xi}\|. \nonumber \end{split} \end{align}Combining these two estimates with Young’s inequality gives the desired result. □ 8.3 Proof of Lemma 6.8 We use Lemmas 6.2, 6.3 and the definition of ϕ:   \begin{align} \begin{split} \left(\frac{\partial \eta}{\partial t},\Pi_{h}\phi\right)= \left(\frac{\partial \eta}{\partial t},\phi\right) +\left(\frac{\partial \eta}{\partial t},\Pi_{h}\phi-\phi\right) \!\leqslant\! Ch^{\,p+1}|u_{t}(t)|_{H^{p+1}}\|\tilde{\xi}\| +Ch^{p+1}|u_{t}(t)|_{H^{p+1}} CL_{\mu} e^{hL_{\mu}}h \|\tilde{\xi}\|. \nonumber \end{split} \end{align}Again, Young’s inequality yields the final estimate. □ © The Author(s) 2018. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) For permissions, please e-mail: journals. permissions@oup.com

Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Apr 6, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off