Stabilization of Galerkin Finite Element Approximations to Transient Convection-Diffusion Problems
Abstract
A postprocessing technique to improve Galerkin finite element approximations to linear evolutionary convection-reaction-diffusion equations is considered. A steady convection-reactiondiffusion problem with data based on the computed standard Galerkin approximation is solved at any fixed time. The postprocessing approximation is obtained using the SUPG method over the same Galerkin finite element space. Error bounds for the method are obtained in the convectiondominated regime. The numerical experiments we present show a substantial reduction of spurious oscillations achieved by means of this procedure.
Full text
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. SIAM J. NUMER. ANAL.c 2010 Society for Industrial and Applied Mathematics Vol. 48, No. 3, pp. 953–979 STABILIZATION OF GALERKIN FINITE ELEMENT APPROXIMATIONS TO TRANSIENT CONVECTION-DIFFUSION PROBLEMS∗ JAVIER DE FRUTOS†,BOSCOGARC ´ IA-ARCHILLA‡,AND JULIA NOVO§ Abstract. A postprocessing technique to improve Galerkin finite element approximations to linear evolutionary convection-reaction-diffusion equations is considered. A steady convection-reaction- diffusion problem with data based on the computed standard Galerkin approximation is solved at any fixed time. The postprocessing approximation is obtained using the SUPG method over the same Galerkin finite element space. Error bounds for the method are obtained in the convectiondominated regime. The numerical experiments we present show a substantial reduction of spurious oscillations achieved by means of this procedure. Key words. stabilized methods, finite element methods, transient convection-diffusion problems AMS subject classifications. 65M60, 65M20, 65M15, 35M25 DOI. 10.1137/090763378 1. Introduction. A new technique to improve the accuracy of the spatial discretization of evolutionary convection-reaction-diffusion problems is studied. In this procedure a steady convection-reaction-diffusion problem, with data based on the computed standard Galerkin approximation, is solved at any time level where the output is desired. More precisely, we consider the problem ut−Δu+b·∇u+cu =fin Ω, u=0 in ∂Ω,(1.1) u(0,x)=u0(x)inΩ, where Ω is a bounded open domain with smooth boundary in Rn,n=1,2,3, band c are given functions, ≥0 is a constant diffusion coefficient, and u0a given initial data. Typically, in some applications, the size of the diffusion is much smaller than the size of the convective term and solutions develop sharp layers. In this case, it is well known that standard finite element methods perform poorly and develop nonphysical spurious oscillations, especially when using low order piecewise polynomials. Stabilization techniques (see, e.g., [23], [24], [25] and the references therein) are widely used in steady problems to suppress oscillations so that physically reasonable approximations can be obtained. The method we study is as follows. Let us denote by uh(t), t∈(0,T], the semidiscrete Galerkin finite element approximation to (1.1). Assume that the output is wanted at the final time T. Then we postprocess uh(T)togetanewapproximation˜uh by solving numerically the steady problem of finding v∈H1 0(Ω) such that L(v)= ∗Received by the editors June 30, 2009; accepted for publication (in revised form) April 20, 2010; published electronically June 29, 2010. The research of the first and third authors was supported by the Spanish MEC under grant MTM2007-60528. http://www.siam.org/journals/sinum/48-3/76337.html †Departamento de Matem´atica Aplicada, Universidad de Valladolid, 47011 Valladolid, Spain ([email protected]). ‡Departamento de Matem´atica Aplicada II, Universidad de Sevilla, 41002 Sevilla, Spain (bosco@ esi.us.es). This author’s research was supported by the Spanish MEC under grant MTM2006-0847. §Departamento de Matem´aticas, Universidad Aut´onoma de Madrid, Instituto de Ciencias Matem´aticas CSIC-UAM-UC3M-UCM, 28049 Madrid, Spain (julia.nov[email protected]). 953 Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 954 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO f−uh,t,whereL(v)=−Δv+b·∇v+cv. This steady problem is solved using the stabilized streamline-upwind Petrov–Galerkin (SUPG) method [14], [4]. In this paper we show that the new approximation ˜uhis essentially free of oscillations. This is most remarkable since the oscillation-free postprocessed approximation ˜uhis computed at a fixed time Tbased on the highly oscillatory uh,t(T) (which have been obtained using the plain Galerkin method along the full time interval (0,T]). This technique was first introduced in [11] for nonlinear convection-diffusion problems of evolution. It was proved in [11] that the new technique possesses a rate of convergence one unit higher up to a logarithmic term than that of the Galerkin method when is away from zero. In [12] the fully discrete case was also addressed. In this paper we carry out the error analysis of the spatial postprocessed semidiscretizations of (1.1) in the convection-dominated regime, that is, when tends to zero. We carry out the analysis when the SUPG method is used in the steady problem of the postprocessing step, although error bounds for other stabilized methods as the DWG (Douglas–Wang/Galerkin) or the GALS (Galerkin/least-squares) can be obtained in a similar way. The new technique is an alternative for the discretization of time-dependent problems with stabilized methods. These have been studied in [24] (see also [2], [3], [5], [6] [8], [9], [15], [16], [17], [18], [19], [20], [21]). Thus, in the present paper, we examine the possibility of computing with a standard method (without stabilization techniques) until a selected time and then use a stabilized method in a single steady problem. This procedure simplifies the task of stabilizing in evolutionary convection-diffusion problems and has the advantage of being extensible (see [11]) to more involved nonlinear problems. In this paper we obtain error bounds in the norm associated to the SUPG method that do not deteriorate when tends to zero. The rate of convergence we prove is suboptimal in the sense that it differs from the rate of convergence of the SUPG approximation to a steady problem by a factor of order h1/2. However, we remark that, to our knowledge, error bounds for the SUPG method in the evolutionary case are scarce or simply not available. We refer to [21], where error bounds for the GALS method are obtained in the fully discrete case using the θ-method assuming that the stabilization parameter is O(Δt), Δtbeing the time step; see also [5], where the limit case = 0 is studied. In practice, both the Galerkin uhand postprocessed approximation ˜uhcannot be computed exactly and, instead, approximations Un h≈uh(tn)and Un h≈˜u(tn) on time levels t0<t 1<···<t N=Tare computed. In this paper we get error bounds for the error in the fully discrete postprocessed approximation Un h. We prove that the temporal error of the fully discrete postprocessed approximation can be bounded in terms of the temporal error of the Galerkin approximation. The results are valid for any convergent time-stepping procedure. We include a refined analysis when the time integrator chosen is the popular midpoint rule. We also obtain improved error bounds in the one-dimensional case using linear finite elements. For this case, we analyze the postprocessed method in both the coercive and noncoercive cases and obtain a quasi-optimal error bound in the H1(0,x N−1)- norm (i.e., in the H1-norm but excluding the last interval). The analysis is carried out for nonuniform meshes and for variable coefficients in the convective and reactive term (i.e., general functions band cin (1.1)). Some numerical experiments are provided to show the reduction of the spurious oscillations that is achieved when using this postprocessing technique. The extension of the error analysis when tends to zero to the nonlinear problems considered in [11] will be the subject of future research. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 955 The outline of the paper is as follows. In section 2 we introduce the notation and state some preliminaries. In section 3 we carry out the error analysis. First, we consider the general case and analyze both the semidiscrete in space and the fully discrete cases. Then the one-dimensional case for linear finite elements is considered. In order to simplify the paper the variable coefficient case is left for the appendix. Finally, in section 4 we present some numerical experiments. 2. Preliminaries and notation. We will assume that band care sufficiently smooth functions of xand define μ(x)=c−1 2divb(x),μ 0=inf x∈Ωμ(x). We assume that Ω is a domain with smooth Crboundary. Let Th=(Kh i,φ h i)i∈Jh, h>0, be a family of shape-regular and quasi-uniform partitions of suitable domains Ωh,wherehis the maximum diameter of the elements Kh i∈T h,andφh iare the mappings of the reference simplex K0onto Kh i.Wedenote Vh,r(Ωh)=vh∈(C0(Ωh))n|vh◦φh i∈(Pr−1(K0)) ∀Kh i∈T h,v h=0in∂Ωh, where Pr−1(K0) is the space of polynomials of (total) degree less than or equal to r−1 over K0. We notice that when only linear elements are used, we may also assume that Ω is a convex polygonal or polyhedral domain and Ωh=Ω. Let us denote by Ah:Vh,r →Vh,r the positive self-adjoint operator defined by (Ahvh,w h)=(∇vh,∇wh)∀vh,w h∈Vh,r, where (·,·) denotes the standard inner product in L2(Ω). The standard L2(Ω) orthogonal projection onto Vh,r will be denoted by Ph. We will denote by πhu∈Vh,r the elliptic projection defined by (∇πhu, ∇ϕh)=(∇u, ∇ϕh)∀ϕh∈Vh,r. Assuming that the meshes are quasi-uniform, the following inverse inequality holds for each vh∈Vh,r (see, e.g., [7, Theorem 3.2.6]): vhWm,q (K)≤Chl−m−d(1 q−1 q)vhWl,q(K),(2.1) where 0 ≤l≤m≤1, 1 ≤q≤q≤∞,andKis an element in the partition Th. The following bound holds for any u∈H1 0(Ω) ∩Hr(Ω) (2.2) u−πhu0+hu−πhu1≤Chrur. Here and in the rest of the paper ·sdenotes the norm of the Sobolev space Hs(Ω) or, if s= 0, the norm of L2(Ω). We remark that for (2.2) to hold, either Ω = Ωhor δ(h)=max{dist(x, ∂Ω) |x∈∂Ωh}=O(h2(r−1)). Following an idea by Wahlbin [26] (see also [1]), we may assume that Ωh⊂Ω, and functions in Vh,r(Ωh) are extended by0toΩ,sothatVh,r(Ωh)⊂H1 0(Ω)n. We denote by uh:(0,T]→Vh,r the spatial semidiscrete Galerkin approximation to (1.1) satisfying uh(0,·)=u0 h∈Vh,r and (2.3) (uh,t,ϕ h)+(∇uh,∇ϕh)+(b·∇uh,ϕ h)+(cuh,ϕ h)=(f,ϕh)∀ϕh∈Vh,r. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 956 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO It is well known that the Galerkin approximation develops spurious oscillations in the advection-dominated regime. In this paper, we consider a procedure that is able to stabilize the Galerkin approximation at any time T. Let us fix any positive time T; we define the postprocessed approximation ˜uh=˜uh(T)∈Vh,r as the solution of the following stationary convection-reaction-diffusion problem: (2.4) (∇˜uh,∇ϕh)+(b·∇˜uh,ϕ h)+(c˜uh,ϕ h)=(f−uh,t,ϕ h) +(f−uh,t +Δ˜uh−b·∇˜uh−c˜uh,b·∇ϕh)h∀ϕh∈Vh,r, where all time-dependent functions are evaluated at the fixed time t=T.Hereand in what follows (·,·)hdenotes the broken inner product (2.5) (f,g)h= K∈Th δK(f,g)K, δKbeing the stabilization parameter and (·,·)Kthe standard inner product in L2(K). We denote by ·hits associated norm. Let us observe that we solve a stationary convection-diffusion problem with data based on the already computed Galerkin approximation using the SUPG method introduced by Hughes and Brooks [14], [4]. The new approximation belongs to the same finite element space as that of the Galerkin approximation. We will denote by aS(·,·) the bilinear form associated to the SUPG method, defined by (2.6) aS(vh,w h)=(∇vh,∇wh)+(b·∇vh,w h)+(cvh,w h) +(−Δvh+b·∇vh+cvh,b·∇wh)h∀vh,w h∈Vh,r. 3. Error analysis. 3.1. General case. In this section we obtain error bounds for the method in which the constants do not deteriorate when tends to zero. In what follows, for simplicity we will assume that uh(0) = πhu0. We notice, however, that as long as uh(0) −u(0)0+huh(0) −u(0)1≤Chru(0)r, the results that follow below are still valid. Lemma 1. Let ube the solution of (1.1) and let πhube its elliptic projection. Let uhbe the Galerkin approximation (2.3).Then,forT>0, there exists a constant C>0, independent of , such that the following bound holds: (3.1) max 0≤t≤T(uh−πhu)t0≤Chr−1u(0)r+max 0≤t≤T(utr+uttr−1). Proof. Let us denote by eh=uh−πhu. Let us call g1=(πhut−ut)+c(πhu−u),g 2=(b·∇(πhu−u)). It is easy to see that for all ϕh∈Vh,r (eh,t,ϕ h)+(∇eh,∇ϕh)+(b·∇eh,ϕ h)+(ceh,ϕ h)=(−g1−g2,ϕ h), Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 957 which we can write as (3.2) eh,t =−Gheh−Ph(g1+g2), where Gh:Vh,r →Vh,r is the operator defined by (3.3) Gh=Ah+Ph(b·∇+cI), Ibeing the identity operator. Since the solution of yh,t =−Ghyhcan be written as yh(t)=exp(−tGh)yh(0), standard energy arguments show that exp(−tGh)0≤e−μ0t. Taking derivatives with respect to tin (3.2) we have eh,tt =−Gheh,t −Ph(g1,t +g2,t), and, thus, eh,t(t)=exp(−tGh)eh,t(0) + t 0 exp(−(t−s)Gh)Ph(g1,t +g2,t)ds. Since by applying (2.2) we have (g1,t +g2,t)0≤Chr−1(utr+uttr−1), the proof will be finished if we show that eh,t(0)0can be bounded by the right-hand side of (3.1). In view of (3.2) and applying (2.2) we have that eh,t(0)0≤Gheh(0)0+Chr−1(u(0)r+ut(0)r−1), and since we are assuming that eh(0) = 0, the result follows. The following lemma states the coerciveness of the bilinear form associated to the SUPG method. The result is standard and can be found, for example, in [25, Lemma 10.3]. Lemma 2. Let us assume that μ0>0and, further, δKc2 K,∞≤μ0 2,δ K≤h2 K 2C2,(3.4) where Cis the constant in the inverse inequality (2.1). Then the bilinear form aS(·,·) associated to the SUPG method satisfies aS(uh,u h)≥1 2∇uh2 0+ K∈Th δKb·∇uh2 0,K +μ1/2uh2 0.(3.5) For linear elements the second condition in (3.4) can be omitted. In what follows we will assume δK=δ0hK/b∞,K whenever PeK=b∞,KhK/ (2)>1andδK=δ1h2 K/ for PeK≤1, where δ0and δ1are user-chosen positive constants. No precise general formula for an optimal value of the stabilization parameter δKis known; see [25, Remark 10.4]. Let us denote by wh∈Vh,r the SUPG approximation to the steady problem (3.6) −Δv+b·∇v+cv =g, subject to homogeneous Dirichlet boundary conditions. Then whsatisfies aS(wh,ϕ h)=(g,ϕh)+(g,b ·∇ϕh)h∀ϕh∈Vh,r. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 958 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO Under the conditions of Lemma 2 it is easy to obtain the following (well-known) error bound for the SUPG approximation (see, for example, [25, Theorem 10.5], [22]): ∇(v−wh)2 0+ K∈Th δkb·∇(v−wh)2 0,K +μ1/2(v−wh)2 01/2 ≤C(1/2+h1/2)hr−1vr.(3.7) We note that the solution uof the evolutionary problem (1.1) is also the solution of (3.6) for the particular case in which g=f−ut. Theorem 1. Let us fix T>0,letube the solution of (1.1),andlet˜uh∈Vh,r be the postprocessed approximation (2.4). Assume that μ0>0and condition (3.4) is satisfied. Then there exists a constant C>0that does not depend on such that the following bound holds: ∇(u−˜uh)2 0+ K∈Th δkb·∇(u−˜uh)2 0,K +μ1/2(u−˜uh)2 01/2 (3.8) ≤Chr−1(ur+utr+uttr−1)+C(1/2+h1/2)hr−1ur. Proof. Let us denote by whthe SUPG approximation to the steady problem (3.6) with g=f−ut. It is easy to see that this approximation satisfies (3.9) (∇wh,∇ϕh)+(b·∇wh,ϕ h)+(cwh,ϕ h)=(f−ut,ϕ h) +(f−ut+Δwh−b·∇wh−cwh,b·∇ϕh)h∀ϕh∈Vh,r. Let us decompose u−˜uh=(u−wh)+(wh−˜uh). To bound the first term we apply (3.7). Let us now get a bound for the second term. Let us denote ˜eh=˜uh−wh. Subtracting (3.9) from (2.4) we get (∇˜eh,∇ϕh)+(b·∇˜eh,ϕ h)+(c˜eh,ϕ h)+(−Δ˜eh+b·∇˜eh+c˜eh,b·∇ϕh)h =(ut−uh,t,ϕ h)+ K∈Th δK(ut−uh,t,b·∇ϕh)K. Taking ϕh=˜ehand applying (3.5) we get ∇˜eh2 0+b·∇˜eh2 h+μ1/2˜eh2 0≤2(ut−uh,t,˜eh)+2(ut−uh,t,b·∇˜eh)h ≤2 μ0ut−uh,t2 0+1 2μ1/2˜eh2 0 +2ut−uh,t2 h+1 2b·∇˜eh2 h. Now notice that for both the convection-dominated or the diffusion-dominated regime we have ut−uh,t2 h≤Chut−uh,t2 0. Thus, by writing ut−uh,t =(ut−πhut)+(πhut−uh,t) and applying (2.2) and Lemma 1 the proof is concluded. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 959 Remark 1. The case μ0= 0, which includes the interesting case c=0and∇·b=0, can be treated in the following way. The change of variables v=ue−αt, with α>0 arbitrary, transforms the equation into one of the same type satisfying μ0=α>0. We can apply the numerical method to this equation and then transform back to the original variables. The result of Theorem 1 is still valid with an O(1) multiplicative constant in the error bound if, for example, we choose α=1/T. A similar comment applies if μ0<0. Alternatively, one can argue as in [25, Remark 7.3]. Remark 2. Let us observe that while in the SUPG method for the evolutionary problem (1.1) the stabilization terms are computed along all the time integration, in our method we carry out the time integration using the standard Galerkin method and compute the stabilization terms only at a fixed time. This procedure, besides being easier to implement, can be easily extended to more complicated nonlinear problems (see [11]). On the contrary, the extension of the SUPG method to nonlinear evolutionary problems is not trivial (see, for example, the discussion in [13]). If we compare the error bound (3.8) that we have obtained for the new method with the error bound (3.5) for the SUPG approximation to a steady problem, we observe a difference of half an order in the rate of convergence. However, to the best of our knowledge, optimal error bounds for the semidiscrete SUPG method in the evolutionary case are not available, even in the linear case. In view of the results obtained in section 3.3 for the one-dimensional case using linear elements where we prove that the error in the energy norm in the last interval is only O(h1/2) (i.e., the same rate of convergence that provides the bound (3.8) for linear elements (r=2)in the L2-norm of the streamline derivative) we think that the error bound (3.8) cannot be improved in general. 3.2. Fully discrete case. Throughout this section we will assume that μ0>0. The general case can be treated as explained in Remark 1. We consider the case in which approximations Un h≈uh(tn) on time levels 0 = t0<t 1<···<t N=Tare computed by means of any convergent time integrator. Given an approximation d∗ tUn h to uh,t(tn) the fully discrete postprocessed approximation Un h∈Vh,r is obtained as the solution of the following problem: (3.10) aS( Un h,ϕ h)=(f−d∗ tUn h,ϕ h)+(f−d∗ tUn h,b·∇ϕh)h∀ϕh∈Vh,r. For the approximation d∗ tUn hto the derivative uh,t(tn)wepropose (3.11) d∗ tUn h=−AhUn h−Ph(b·∇Un h+cUn h)+Phf(tn). This is the approximation used in the numerical experiments of section 4 and the same that has been considered in [10], [12]. To estimate the error u(tn)− Un hwe write u(tn)− Un h=(u(tn)−˜uh(tn)) + ˜en, where ˜en=˜uh(tn)− Un h. The first term on the right-hand side above is the error in the spatial semidiscrete postprocessed approximation that has been bounded in Theorem 1. Next, we analyze the time discretization error ˜en. We estimate the size of ˜enin terms of en=uh(tn)−Un h, the temporal error of the fully discrete Galerkin approximation. We carry out the error analysis for any convergent time-stepping procedure satisfying lim k→0max 0≤n≤Nen0=0,lim k→0max 0≤n≤Nen1=0, where k=max{kn−1=tn−tn−1|1≤n≤N}. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 960 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO Theorem 2. Let us fix tn>0,andlet˜en=˜u(tn)−˜ Un hbe the temporal error of the fully discrete postprocessed approximation. Then there exists a constant C>0 that does not depend on such that the following bound holds: (3.12) ∇˜en2 0+ K∈Th δkb·∇˜en2 0,K +μ1/2˜en2 01/2 ≤Ch−1en1+b∞en1+c∞en0, where en=uh(tn)−Un his the temporal error of the Galerkin approximation. Proof. The proof is similar to the proof of Theorem 1. Subtracting (3.10) from (2.4) we get aS(˜en,ϕ h)=(d∗ tUn h−uh,t,ϕ h)+(d∗ tUn h−uh,t,b·∇ϕh)h∀ϕh∈Vh,r. Taking ϕh=˜ehand applying (3.4) we get ∇˜en2 0+b·∇˜en2 h+μ1/2˜en2 0≤2(d∗ tUn h−uh,t(tn),˜en) +2(d∗ tUn h−uh,t(tn),b·∇˜en)h ≤2 μ0d∗ tUn h−uh,t(tn)2 0+1 2μ1/2˜en2 0 +2d∗ tUn h−uh,t(tn)2 h+1 2b·∇˜en2 h. Now, since we have ·h≤Ch·0it is sufficient to get a bound for d∗ tUn h−uh,t(tn)0 to conclude the proof. Taking into account that uh,t =−Ahuh(tn)−Ph(b·∇uh(tn)+cuh(tn)) + Phf(tn) and using (3.11) we get (3.13) d∗ tUn h−uh,t(tn)0≤Ahen+Ph(b·∇en+cen)0. Applying the inverse inequality (2.1) we finally arrive at (3.14) d∗ tUn h−uh,t(tn)0≤Ch−1en1+b∞en1+c∞en0. In view of bound (3.12) we deduce that the temporal error of the fully discrete postprocessed approximation can be bounded in terms of the temporal error of the Galerkin approximation. Analogous results were obtained in [10] for an earlier postprocessed technique applied to the nonlinear evolutionary Navier–Stokes equations, and in [12], where the same postprocessing technique was applied to nonlinear evolutionary convection-diffusion equations in the diffusive regime. In [10] and [12] we prove that the temporal errors of the Galerkin finite element approximations satisfy en0+Ahen0≤Ckl0,1≤n≤N, where kis the time step, l0= 1 for the backward Euler method, and l0=2for the two-step backward differentiation formula. Let us point out that, as was already observed in [10] and [12], better bounds are obtained in the diffusion-dominated regime if Ahen0is used in (3.12) instead of h−1en1. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 961 Let us now consider the case in which the fully discrete approximations Un hare obtained by means of the trapezoidal rule. Then Un+1 h−Un h kn +GhUn+1 h+Un h 2=f(tn+1)+f(tn) 2,n=0,...,N −1, where, as in (3.3), Gh=Ah+Ph(b·∇+cI). We have that (3.15) vh2 G=(Ghvh,v h)≥vh2 1+μ0vh2 0. In the following theorem we bound the temporal error of the fully discrete postprocessed approximation for this particular case. Let us remark, though, that the result can be similarly proved for other A-stable time integrators. Theorem 3. Let us fix tn>0, let us assume that we integrate in time with the trapezoidal rule, and let ˜en=˜u(tn)−˜ Un hbe the temporal error of the fully discrete postprocessed approximation. Then there exists a constant C>0that does not depend on such that the following bound holds: (3.16) ∇˜en2 0+ K∈Th δkb·∇˜en2 0,K +μ1/2˜en2 01/2 ≤Ck2tn 0 d3 dt3f(t) 2 0 + d4 dt4uh(t) 2 0dt1/2 . Proof. Reasoning as in the proof of Theorem 2 but using (3.13) instead of (3.14) we arrive at ∇˜en2 0+ K∈Th δkb·∇˜en2 0,K +μ1/2˜en2 01/2 ≤CGhen0. In the rest of this proof we get a bound for Ghen0. For simplicity we will assume that (3.17) e0=0. We first observe that for n=0,1,...,N −1 (3.18) 1 kn (en+1 −en,ϕ h)+Ghen+1/2,ϕ h=(τn,ϕ h)∀ϕh∈Vh,r, where here and in the rest of the proof en+1/2=en+1 +en 2,n=0,1,...,N −1, and τnis the truncation error, that is, τn=uh(tn+1)−uh(tn) kn +Ghuh(tn+1)+uh(tn) 2−f(tn+1)+f(tn) 2. A simple calculation shows that (3.19) τn=1 kntn+1 tn uh,t(t)dt−uh,t(tn+1)+uh,t(tn) 2 =−1 kntn+1 tn (tn+1 −t)(t−tn) 2 d3 dt3uh(t)dt. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 968 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO In the next theorem we obtain the same error bound of Theorem 4 for this approximation. Since we are in the coercive case, Lemma 2 can be applied. The first condition on (3.4) must hold, whereas the second one is not required since we are dealing with linear elements. A simple calculation shows that this happens if h≤b/c. Theorem 5. Let ube the solution of (3.53) and ˜uhthe postprocessed approximation (3.54),andsetδKas in (3.25). Then there exists a constant Cthat does not depend on such that the following bounds hold for h≤b/c: u−˜uh1,[0,xN−1]≤Ch+ h1/2(u2+K),(3.55) u−˜uh1,[xN−1,xN]≤Ch1/2(u2+K),(3.56) where Kis the constant in (3.44). Proof. Let us denote by whthe SUPG approximation to the steady convectionreaction-diffusion equation −uxx +bux+cu =g, u(0) = u(1) = 0, for g=f−ut. Then the function whsolves (3.57) (wh,x,ϕ h,x)+(bwh,x +cwh,ϕ h)+(bwh,x +cwh,bϕ h,x)h =(g,ϕh)+(g,bϕh,x)h. Applying (3.7) we obtain u−wh0≤Ch3/2u2,u−wh1≤Chu2. To bound the error in the postprocessed approximation we decompose as usual u− ˜uh=(u−wh)+(wh−˜uh). Let us obtain a bound for the second term. We denote by ˜eh=˜uh−wh. Taking into account the equivalence (3.26) and applying Theorem 1, we obtain (3.58) ˜eh0≤CKh, ˜eh,x0≤CKh1/2. Let us first observe that from (3.58) we get (3.56) and as a consequence |DeN|≤C. To prove (3.55) we will apply Lemmas 3 and 4. Notice that subtracting (3.57) from (3.54) we get (˜eh,x,ϕ h,x)+(b˜eh,x +c˜eh,ϕ h)+(b˜eh,x +c˜eh,bϕ h,x)h =(ut−uh,t,ϕ h)+(ut−uh,t,bϕ h,x)h. Thus, taking vh=˜ehin Lemma 3, we have that sj=˜rjand Sj=˜r1+···+˜rj,where (3.59) ˜rj=((ut−uh,t)+c˜eh,ϕ j)+((ut−uh,t)+c˜eh,bϕ j,x)h,j=1,...,N −1. Observethatwehaveanestimateof˜eh0in (3.58). Also, as a consequence of (2.2) and Lemma 1, we have ut−uh,t0≤CKh. Thus, applying Lemma 4 (estimate (3.35)) to the residuals in (3.59), we have that, as in Theorem 4, (3.45) holds. Then, applying (3.27), we reach (3.55). Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 969 4. Numerical experiments. Next, we show a simple numerical experiment that illustrates the behavior of the postprocessed approximation. We consider (3.21) with forcing term f= 1, convection coefficient b= 1, and initial condition u0(x)=0. We compute the Galerkin approximation based on linear finite elements over a uniform partition of [0,1] of size h=1/N at time T=0.6. For the time integration we use the implicit Euler method with fixed time step k=0.001. To get the postprocessed approximation ˜uhwe solve (3.23). On the left of Figure 4.1 we have represented the Galerkin (solid line) and postprocessed (dashed line) approximations for =1e−4 and N= 80. We can observe that the postprocessing step annihilates the spurious oscillations of the Galerkin approximation. On the right of Figure 4.1 we have represented the Galerkin time derivative uh,t. The bound (3.52) indicates that the difference between two values in the postprocessed approximation is bounded in terms of the local L∞-norm of the error in the Galerkin time derivative uh,t. Indeed, if we look at a zoom of Figure 4.1 (see Figure 4.2) we can observe that the postprocessing step does not annihilate completely the Galerkin oscillations, although it considerably 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 ε=1e−4, N=80 x uh(x,1) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.80.9 1 −0.5 0 0.5 1 1.5 2 ε=1e−4, N=80 x (uh(x,1))t Fig. 4.1.On the left: Galerkin and postprocessed approximations. On the right: Galerkin time derivative. 0.9 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 0.565 0.57 0.575 0.58 0.585 0.59 0.595 ε=1e−4, N=80 x uh(x,1) Fig. 4.2.Zoom of Figure 4.1(left). Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 970 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 x ε=1e−5, N=80 uh(x,1) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.80.9 1 −1 −0.5 0 0.5 1 1.5 2 x ε=1e−5, N=80 (uh(x,1))t Fig. 4.3.On the left: Galerkin and postprocessed approximations. On the right: Galerkin time derivative. 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ε=1e−5, N=300 x uh(x,1) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.80.9 1 −0.5 0 0.5 1 1.5 2 ε=1e−5, N=300 x (uh(x,1))t Fig. 4.4.On the left: Galerkin and postprocessed approximations. On the right: Galerkin time derivative. diminishes them. To check the good behavior of our method when tends to zero we have computed the Galerkin and postprocessed approximations with N=80for =1e−5. In Figure 4.3 we have plotted the approximations on the left and the Galerkin time derivative on the right. We can observe in the figure that although the postprocessed approximation annihilates again the Galerkin oscillations, it is not accurate enough since the Galerkin approximation used in the postprocessing step (3.23) is completely inaccurate; observe the picture of the Galerkin time derivative on the right of Figure 4.3. This lack of accuracy can be solved by computing the Galerkin approximation over a refined mesh. In Figure 4.4 we show the results obtained using a partition of [0,1] into 300 subintervals. Now, even though the Galerkin approximation is still completely contaminated, the postprocessing step is able to produce an accurate and oscillation-free approximation. The lack of accuracy of ˜uhwhen the mesh is not fine enough, which is observed in Figure 4.3, can be a drawback of the method in practice. However, we can still take advantage of the nonoscillatory character of the postprocessed approximation. The main idea is to use the difference ηh=˜uh−uhbetween the postprocessed oscillation-free approximation ˜uhand the polluted Galerkin approximation uhas an a posteriori indicator, at each element of the mesh, of the oscillations presented in the Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 971 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.80.9 1 −1 −0.5 0 0.5 1 0 0 .1 0 .2 0 . 3 0 .4 0 . 5 0 . 6 0 .7 0 . 80 . 9 1 0 0.2 0.4 0.6 0.8 Fig. 4.5.Approximation obtained using the adaptive algorithm for =1e−6. Galerkin approximation. Using this indicator at each time step, we can detect the oscillations developed in the Galerkin approximation and consequently locally refine the mesh before these oscillations become excessively large and globally pollute the approximation. The adaptive procedure we now present provides a wider application of the postprocessing technique since it produces accurate and oscillation-free Galerkin approximations with a very small number of degrees of freedom. •Choose an initial subdivision of the interval [0,1]. •Compute the Galerkin approximation at the first time step. •Compute the postprocessed approximation. •Compute the error indicator as the difference between the postprocessed and the Galerkin approximations: ηh=˜uh−uh. •For every element Ij=[xj−1,x j] if the difference |ηh(xj)−ηh(xj−1)|is greater than a given tolerance TOL1,halvetheintervalIj. If the difference is less than a given tolerance TOL2<TOL1, suppress the point xj−1whenever the new interval does not exceed a maximum prescribed size. Interpolate the approximation and use it as initial condition for the next time step. •Continue with the procedure until the final time T. We now show a numerical experiment to illustrate the behavior of our algorithm. We consider the same experiment as before with a smaller value of , more precisely =1e−6, for which, in view of Figure 4.3, we cannot expect an accurate postprocessed approximation. The initial mesh for this experiment has 100 nodes and the maximum his set to 0.04. The parameters TOL1and TOL2were set to TOL1=0.01 and TOL2=TOL 1/100. In Figure 4.5 we show the approximation obtained at the final time T=0.6 (bottom) and the final mesh (top). Our algorithm ends with only 45 nodes from which 19 lie on the interval [0.9,1] and, as we can observe in the figure, produces an excellent approximation in which the boundary layer is perfectly solved. The extension of this adaptive procedure to more than one spatial dimension as well as to nonlinear problems will be the subject of future research. Next, we show a numerical experiment in a two-dimensional problem. Let us consider the equation ut−Δu+ux+uy=f Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 972 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x uh(0.6,x,x) ε=0.005 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x ~uh(0.6,x,x) ε=0.005 Fig. 4.6.Sections with the plane y=x. Exact solution shown as solid line. On the left, Galerkin approximations. On the right, postprocessed approximations. 00.2 0.4 0.6 0.81 0 0.5 1 −0.5 0 0.5 1 1.5 2 2.5 00.2 0.4 0.6 0.81 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Fig. 4.7.On the left, Galerkin approximation for h=1/10. On the right, postprocessed approximation. =0.001. in the domain Ω = [−1,1]×[−1,1] subject to homogeneous Dirichlet boundary conditions. We take f= 1 and as initial condition u0(x, y) = 0. Let us consider regular triangulations of Ω induced by the set of nodes (i/N, j/N),0≤i, j ≤N,whereN=1/h is an integer. We use linear finite elements. The final time chosen is T=0.6. All the experiments are carried out using MATLAB. For the time integration we use the midpoint rule with fixed time step. To compare the methods a reference approximation wascomputedwiththeGalerkinmethodonaveryfinemeshandwithsufficiently small time steps. In Figure 4.6 we have represented the sections of the Galerkin and postprocessed approximations along the plane y=xfor =0.005. The exact solution is plotted using a solid line. On the left, we plot the Galerkin approximations for h=1/10,1/20, and 1/40 using dashed, dotted, and dash-dotted lines, respectively. The same lines are used on the right for the postprocessed approximations. It can be observed, in agreement with the experiments shown before in one spatial dimension, that the postprocessing step considerably reduces the spurious oscillations for all the values of hin the figure. For the last value h=1/40, the postprocessed approximation does not oscillate at all and matches very precisely the section of the exact solution. To observe the behavior of our method when decreases, we have represented in Figures 4.7 and 4.8 the Galerkin (left) and postprocessed (right) approximations Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 973 00.2 0.4 0.6 0.81 0 0.5 1 −0.5 0 0.5 1 1.5 2 2.5 00.2 0.4 0.6 0.81 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Fig. 4.8.On the left, Galerkin approximation for h=1/20. On the right, postprocessed approximation. =0.001. for =0.001. In Figure 4.7 we represent the approximations for h=1/10 and in Figure 4.8 for h=1/20. We can observe that although the Galerkin approximation is completely contaminated by spurious oscillations all over the whole domain Ω, the postprocessed method provides quite accurate approximations with a coarse mesh of only N=10orN= 20 nodes for each variable. The postprocessed approximation with N= 10 still has some small oscillations away from the boundary layers, while in the case N= 20 the small oscillations remain only in the neighborhood of the boundary layer. Of course, the oscillations can be completely annihilated by increasing the number of degrees of freedom. However, the aim of these figures is to show that from a “completely wrong” Galerkin approximation this postprocessing procedure can recover enough information to compute quite accurate approximations. 5. Appendix: Variable coefficients. We analyze the one-dimensional case allowing for bto depend on x.Wedenote bM=max 0≤x≤1b(x),b 0=min 0≤x≤1b(x), and we will assume that b0>0. We will follow the results on the constant coefficient case, commenting on the differences. We start by noticing that the postprocessed approximation satisfies (3.23), and the SUPG approximation whto the solution uof (3.41) satisfies (3.42). We now set (5.1) δK=δj=hj/(2bj) for all elements K=[xj−1,x j], j=1,...,N, where (5.2) bj=1 2xj xj−1 b2(x)dxxj xj−1 b(x)ϕj−1(x)dx,j=1,...,N. With this choice of the parameter δKthe following relation holds: (5.3) (bvh,x,ϕ j)+(bvh,x,bϕ j,x)h=bj(vh(xj)−vh(xj−1)),j=1,...,N −1, where bj=1 hjxj xj−1 b(x)dx. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 974 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO Remark 5. The choice of δKhere is not restrictive. With other choices, like, for example, δK=δj=hj/(2 bL∞(xj−1,xj)), the right-hand side of (5.3) should be replaced by ˜ bj(vh(xj)−vh(xj−1)) −εj+1(vh(xj+1)−vh(xj)), with |bj−˜ bj|=O(hq j)and|εj+1|=O(hq j+1)forsomeq≥1. This makes the analysis much more cumbersome and lengthy, but the results are essentially those we state below. Observe also that by the mean value theorem, bj=b(ξj)forsomeξj∈(xj−1,x j), so that bj≥b0,j=1,...,N. Also, H¨older’s inequality and the mean value theorem shows that b2 0/bM<b j< b2 M/b0.Consequently,forj=1,...,N, b0 2b2 M h0≤b0 2b2 M hj≤δj≤bM 2b2 0 hj≤bM 2b2 0 h, so that the following relations follow: b0 2b2 M h0vh2 0<vh2 h≤bM 2b2 0 hvh2 0∀vh∈Vh. As in previous sections, we restrict ourselves to the convection-dominated regime, so that (5.4) b0h0 2>1, but we will need to further assume (5.5) b0h0 2>2exp 1 b01 0|b(x)|dx. Indeed, the bound that will be needed is (5.6) b0h0 2>2exp 1 b0 N−1 j=1 |bj+1 −bj|, but since, as commented above, bj=b(ξj)forsomeξj∈(xj−1,x j), we have |bj+1 −bj|=ξj+1 ξj b(x)dx≤ξj+1 ξj|b(x)|dx, and then (5.5) implies (5.6). We now state and prove the version of Lemma 3 for the variable coefficient case. In its proof, the following discrete Gronwall lemma will be needed, which can be easily proved by an induction argument. Lemma 5. Let (yn)∞ n=1 and (γn)∞ n=1 be sequences of nonnegative numbers and σ0>0such that yn≤γ0+ n−1 j=1 γjyj,n=1,2,.... Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 975 Then yn≤γ0exp(γ1+···+γn−1),n=1,2,.... Lemma 6. Assume that (5.4) and (5.6) hold, and for vh∈Vhand j=1,...,N−1, let sj=(vh,x,ϕ j,x)+(bvh,x,ϕ j)+(bvh,x,bϕ j,x)h and Sj=s1+···+sj. Then the following bounds hold: vh,x0,(x0,xN−1)≤2 b0h1/2 0s2 1+···+s2 N−1+|DvN|,(5.7) b0|DvN|≤ 1 hNs2 1+···+s2 N−1+eσN−1 hN max 1≤j≤N−1|Sj|,(5.8) bn|vn|≤eσn−1max 1≤j≤n|Sj|+2eσn−1 h1/2 0(vh)x0,(x0,xn+1) (5.9) for n=1,...,N −1,where σn=1 b0 n−1 j=1 |bj+1 −bj|,n=1,...,N −1. Proof. In view of (5.3) we have (5.10) sj=(+bjhj)Dvj−Dvj+1,j=1,...,N −1. Then, recalling (3.31), from (5.10) it follows that h1/2+b0h1/2 0vh,xL2(x0,xN−1)≤ h1/2 0vh,xL2(x1,xN−1)+|DeN| +N−1 j=1 |sj|21/2 , and since according to (5.4), /h1/2 0≤b0h1/2 0/2, the bound (5.7) follows. We will now prove (5.8). Summation in (5.10) from j=1toj=ngives (5.11) n j=1 bj(vj−vj−1)=(Dvn+1 −Dv1)+Sn,n=1,...,N −1. Summation by parts allows us to write (recall that v0=0) n j=1 bj(vj−vj−1)=bnvn− n−1 j=1 (bj+1 −bj)vj. Thus, we can rewrite (5.11) as (5.12) bnvn= n−1 j=1 (bj+1 −bj)vj+(Dvn+1 −Dv1)+Sn,n=1,...,N −1. Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. 976 J.DEFRUTOS,B.GARC ´ IA-ARCHILLA, AND J. NOVO Applying Lemma 5 we have (5.13) |vn|≤ 1 bn exp1 b0 n−1 j=1 |bj+1 −bj|max 1≤j≤n(Dvj+1 −Dv1)+Sj. On the other hand, since vN=0wehaveDvN=−vN−1/hN,so that setting n= N−1 in (5.12) and multiplying by −h−1 Nwe have bN−1DvN=− hN (DvN−Dv1)−1 hN SN−1−1 hN N−2 j=1 (bj+1 −bj)vj, that is, (5.14) bN−1+ hNDvN= hN Dv1−1 hN SN−1−1 hN N−2 j=1 (bj+1 −bj)vj. In order to finish the proof of (5.8) we need to express the first and third terms on the right-hand side above in terms of the previous bounds. We start with the third term. Observe that |bj+1 −bj||vj|≤|bj+1 −bj| b0||bjvj|=(σj+1 −σj)||bjvj|. Then, in view of (5.13), we have N−2 j=1 (bj+1 −bj)vj≤N−2 j=1 (σj+1 −σj)eσjmax 1≤j≤N−2(Dvj+1 −Dv1)+Sj. Notice that the sum above is a lower Riemann sum of the exponential function and, thus, smaller than the corresponding integral. Hence, N−2 j=1 (bj+1 −bj)vj≤(eσN−1−1) max 1≤j≤N−2(Dvj+1 −Dv1)+Sj. Recalling (3.34) we can write (5.15) max 1≤j≤N−2(Dvj+1 −Dv1)+Sj≤2 h1/2 0vh,xL2(x0,xN−1)+max 1≤j≤N−2Sj and, thus, N−2 j=1 (bj+1 −bj)vj≤(eσN−1−1) 2 h1/2 0vh,xL2(x0,xN−1)+max 1≤j≤N−2Sj. Then, going back to (5.14), we get bN−1+ hN|DvN|≤2eσN−1 hNh1/2 0vh,xL2(x0,xN−1)+eσN−1 hN max 1≤j≤N−1|Sj|, Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. TRANSIENT CONVECTION-DIFFUSION PROBLEMS 977 and applying (5.7) we obtain bN−1+ hN|DvN|≤4eσN−1 hNb0h0s2 1+···+s2 N−1+|DvN|eσN−1 hN max 1≤j≤N−1|Sj|, so that, rearranging terms, we have bN−1+ hN1−4eσN−1 b0h0|DvN|≤4eσN−1 hNb0h0s2 1+···+s2 N−1 eσN−1 hN max 1≤j≤N−1|Sj|. Since we are assuming that (5.6) holds, the left-hand side above can be bounded below by bN−1|DvN|, and the coefficient multiplying the square root on the right-hand side can be bounded by 1/hN. Hence, (5.8) follows. Finally, (5.9) follows from (5.13) and (5.15). Lemma 7. For δKas specified in (5.1),(5.2), the bounds (3.35–3.37) hold for v, w∈L2(0,1) and n=1,...,N −1. Proof. We follow the proof of Lemma 4. Observe that (3.35) follows from (3.38) and (3.40). In the present case (3.38) obviously holds, so that we are left to show that (3.40) also holds in the present case. To do this we notice that instead of (3.39) we now have (w,b(ϕj)x)h=1 2bjxj xj−1 b(x)w(x)dx−1 2bj+1 xj+1 xj b(x)w(x)dx. However, in view of the expression of bjin (5.2), applying H¨older’s inequality we have 1 bj≤ϕj−1L2(xj−1,xj) bL2(xj−1,xj) =h1/2 j √3bL2(xj−1,xj) , so that applying H¨older’s inequality again it follows that (5.16) 1 2bjxj xj−1 b(x)w(x)dx≤√3 6h1/2 jwL2(xj−1,xj),j=1,...,N, and thus |(w,b(ϕj)x)h|≤(√6/6)h1/2wL2(xj−1,xj+1), which implies (3.40). Thus, (3.35) follows. Also, thanks to (5.16), the same arguments used to prove (3.36) and (3.37) in Lemma 4 also apply in the present case. Lemma 8. Let ube the solution of (3.41) and let (5.4) and (5.6) hold. There exists a positive constant Cindependent of such that the SUPG approximation wh satisfies the error estimate (3.49). Proof. The proof can be obtained by reasoning as in the proof of Theorem 4, using Lemmas 6 and 7; see [22]. For the postprocessed approximation, we state Theorem 6 below, which is the variable coefficient version of Theorem 4. It can be proved by following the arguments in the proof of Theorem 4, provided that references to Lemmas 3, 4 and (3.49) are replaced by references to Lemmas 6, 7, and 8, respectively. Theorem 6. Let ube the solution of (3.21) and let (5.4) and (5.6) hold. Then there exists a positive constant Cthat does not depend on such that the postprocessed approximation ˜uhsolution of (3.23) satisfies estimates (3.47) and (3.48). Finally we consider the equivalent of Theorem 5 with variable coefficients, that is, we consider the problem (5.17) ut−uxx +bux+cu =f, 0<x<1,t>0; u(0,t)=u(1,t)=0 u(x, 0) = u0, Downloaded 04/18/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php