Postprocessing the Galerkin method: the finite-element case

POSTPROCESSING THE GALERKIN METHOD: THE FINITE-ELEMENT CASE∗ BOSCO GARC´ IA-ARCHILLA†AND EDRISS S. TITI‡ SIAM J. NUMER. ANAL.c °2000 Society for Industrial and Applied Mathematics Vol. 37, No. 2, pp. 470–499 Abstract. A postprocessing technique, developed earlier for spectral methods, is extended here to Galerkin finite-element methods for dissipative evolution partial differential equations. The postprocessing amounts to solving a linear elliptic problem on a finer grid (or higher-order space) once the time integration on the coarser mesh is completed. This technique increases the convergence rate of the finite-element method to which it is applied, and this is done at almost no additional computational cost. The numerical experiments presented here show that the resulting postprocessed method is computationally more efficient than the method to which it is applied (say, quadratic finite elements) as well as standard methods of similar order of convergence as the postprocessed one (say, cubic finite elements). The error analysis of the new method is performed in L2and in L∞norms. Key words. dissipative equations, postprocessing finite-element methods, multilevel methods, nonlinear Galerkin methods AMS subject classifications. 65P25, 65M70, 65M15, 65M20 PII. S0036142998335893 1. Introduction. Finite-element methods do not seem to have benefited as much as spectral methods from some of the recent advances in the dynamical systems approach to partial differential equations (PDEs) like those in [15], [14], [13], and [41]. Since these advances are mostly based on spectral decompositions, they are readily adapted to spectral methods, but do not seem to have a clear-cut translation to finite elements. In [21] and based on the results of [15] and [13], an inexpensive novel technique to increase the accuracy and computational efficiency of Fourier spectral methods was developed. In this paper, we present a general technique to improve the convergence rate of Galerkin methods which applies to finite-element methods (and in the particular case of Fourier–Galerkin methods it coincides with that in [21]; see also [37]). We do this with no direct address to [13]. Let Ω ⊂Rdbe a domain with a smooth boundary. We consider dissipative PDEs (see, for instance, [8], [24], [25], [41]) which can be written as evolution systems of the form du dt +νAu +F(u)=0,(1) in the Hilbert space H(in this work H=L2(Ω)), with solutions determined uniquely by the initial condition u(·,0) = u0.(2) ∗Received by the editors March 19, 1998; accepted for publication (in revised form) March 8, 1999; published electronically January 5, 2000. This work was partially financed by DGICYT project PB95-216 and the National Science Foundation and was completed while the second author was the Orson Anderson Visiting Scholar at the Institute of Geophysics and Planetary Physics (IGPP) at the Los Alamos National Laboratory. http://www.siam.org/journals/sinum/37-2/33589.html †Departamento de Matem´aticas, Universidad Aut´onoma de Madrid, 28049 Madrid, Spain (b[email protected]). Current address: Depto. Matematica Aplicada II, Universidad de Sevilla, Escuela Superior de Ingenieror, Camino de lor Descubrimientor S/N, 41092 Sevilla (Spain) ‡Department of Mathematics and Department of Mechanical and Aerospace Engineering, University of California, Irvine, CA 92697-3875 ([email protected]). 470 Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 471 Here, ν>0 is a scalar, A:D(A)⊂H→His a densely defined, unbounded, self-adjoint, and positive operator with compact inverse. For simplicity, we treat here only the case A=−∆, where ∆ stands for the Laplacian operator subject to homogeneous Dirichlet boundary conditions, but our results also apply to other elliptic operators in divergence form with sufficiently smooth coefficients and (with adequate modifications) to other boundary conditions. The (nonlinear) operator F that we consider is one of the following two types: 1. reaction-diffusion equations, F(u)=g(u),(3) for some smooth function g:R→R; 2. reaction-diffusion equations plus nonlinear convection, F(u)=g(u)+b(u)·∇u,(4) for some smooth mapping b:R→Rd, and gis as above. For simplicity, we concentrate here on the two types of nonlinearity above, but both gand bcan be made dependent on xand tas well. Also, although we concentrate on scalar equations, the results here can be extended to systems. We treat here the case where Ω ⊂Rdwith d>1; the case d= 1 is particular and much simpler to analyze. Let Th=(τ h i,φ h i) i∈I h,h>0, be a family of partitions of suitable domains Ωh, where the parameter his the maximum diameter of the elements τh iin Th, and φh i are the mappings of the reference simplex τ0onto τh i.Forr≥2 we consider the finite-element spaces Sh,r =nχh∈C(Ωh)|χh|τh i◦φh i∈Pr−1(τ0),χ h (x)=0 ∀x∈∂Ω h o, where Pr−1(τ0) denotes the space of polynomials of degree at most r−1onτ 0 . Let us denote by a(·,·) the positive definite, bilinear form induced by A, that is, a(u, v)=(A 1/2 u, A1/2u)=(∇u, ∇v) for u, v ∈D(A1/2)=H 1 0(Ω), where (·,·) denotes the standard inner product in L2(Ω) (or in L2(Ω)d). In a similar manner let ah(·,·) be the corresponding bilinear form on Sh,r (in our case ah(χh,ψ h)=(∇χ h ,∇ψ h ) h , χ h ,ψ h∈S h,r, where (·,·)hstands for the standard inner product in L2(Ωh)); and let us denote by Ahthe associated positive, self-adjoint operator in Sh,r, that is, ah(χh,ψ h)=(A h χ h ,ψ h) h=(χ h ,A hψ h) h∀χ h,ψ h∈S h,r. Let us also denote by Phthe standard L2orthogonal projection onto Sh,r, and by Rhthe elliptic projection onto Sh,r (also known as the H1projection) which, for u∈H1 0(Ω) is defined by ah(Rhu, χh)h=ah(u, χh)=(∇u, ∇χh)h∀χh∈Sh,r. The discrepancy between Ω and Ωhis usually solved by an adequate selection of quadrature rules, although other alternatives are also possible (see section 2.1). The Galerkin approximation to the solution uof (1)–(2) is described as follows: Find uh:[0,T]→S h,r such that d dtuh+νAhuh+PhF(uh)=0,(5) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 472 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI uh(0) = Rhu0.(6) We use (6) for simplicity in the analysis. Our results may be easily extended to the case where uh(0) = Phu0. Next, we motivate and introduce the method studied here. Suppose that we are interested in the solution uof (1)–(2) at a given time T. At time T, rewriting (1), we have that νAu(T)=−F(u(T)) −d dt u(T). Thus, u(T) can be seen as the solution of an elliptic problem whose right-hand side is not known but can be approximated. The method we propose is as follows. (i) First, integrate (5)–(6) up to Tto obtain the Galerkin approximation uh(T). (ii) Then solve (or, in practice, approximate) the following linear elliptic problem: Find ˜u∈D(A) such that νA˜u=−F(uh(T)) −d dtuh(T).(7) We call this approximation procedure the postprocessed Galerkin method. Under certain assumptions (to be specified in section 2.1) we can state Theorem 1 below, whose proof is deferred until section 2.3 (a version of this result in L∞is presented in section 4.1). In Theorem 1 and in the rest of the paper, the constants µand rare µ=2ifr≥4 and (3) holds, 1 otherwise,r=0ifr= 3 and (3) holds, 1 otherwise. Theorem 1. Let T>0and r≥3.There exists a constant C>0which depends on K(u),defined in (22) below, such that for hsufficiently small the postprocessed Galerkin approximation ˜usolution of (7) satisfies ku(T)−˜ukL2(Ω) +hku(T)−˜ukH1(Ω) ≤Chr+µ|log(h)|r.(8) We notice that the bound (8) is an improvement over the standard Galerkin error bound, ku−uhkL2(Ω) +hku−uhkH1(Ω) ≤Chr, which is optimal even in the case of linear problems. We remark that in practice, of course, ˜ucannot be computed exactly since in general it does not belong to a finite-dimensional space, as opposed to uh. However, one can approximate the solution ˜uof (7) by some ˜uhbelonging to a finite-element space of piecewise polynomials either of degree r+µ−1 over Th, or of degree r−1 but over some finer partition Th0with hr+µ=h0r. It must be noticed that, from the practical point of view, the computational cost of the postprocessing step (7) is a very small fraction of the cost of computing uh(T)by time integration from t=0tot=T. Thus, ˜uhis an approximation to u(T) of higher order (i.e., asymptotically more accurate) than uh(T), but it costs (almost) nothing once uh(T) is available. We will see the implications of this fact in the numerical experiments in section 3. Observe that the combination of the Galerkin method (5)–(6) and the postprocessing step (7) can be seen as a novel two-level or two-grid method for evolution equations. Two-level methods for nonlinear problems have mostly been developed and implemented for steady equations (see, e.g., [5], [31], [47], [48], and the references therein). The basic idea in the steady-state case is to solve the nonlinear problem on a low-order approximation space (or coarse grid) and then obtain a better approximation with one or two Newton iterations on a higher-order approximation space (or Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 473 on a finer grid). In this way, the nonlinearities are treated in the (cheaper) low-order approximation space or coarse grid, and only linear systems are solved in the (more costly) higher-order approximation space. On evolution problems, two-level (or further, multilevel) methods can obviously be used in the nonlinear systems that arise when using implicit time integrators. Other more refined ideas have recently been suggested [3], [9], [33], [34], but, so far, the two levels of discretization are used at every time step along the time evolution. On the contrary, Theorem 1 in this paper states that it is safe to perform the time evolution on the low-level (or coarse grid) approximation space and obtain the higherorder correction only once,att=T, that is, when the time-marching is finished. The price to be paid is that the low-level (or coarse grid) approximation space must be composed of piecewise polynomials with a degree of at least 2 (r≥3), since the postprocessing step (7) does not increase the convergence rate if uh(T) is computed by integrating on time on a space of piecewise-linear polynomials. The postprocessing step (7) is a generalization to finite-element methods of an earlier postprocessing procedure developed for spectral methods in [16], [21], and [22]. To be more specific, suppose that instead of Sh,r, the approximation space is Hmspanned by the eigenfunctions of Aassociated with the lowest meigenvalues or frequencies. Then, the Galerkin approximation (5) would be um:[0,T]→H msuch that dum dt +νAum+PmF(um)=0,(9) where Pmis the L2projection onto Hm. The postprocessing step (7) would be νA˜u=−dum dt (T)−F(um(T)).(10) Projecting (10) onto Hm, since Pmand Acommute, we get νAPm˜u=−PmF(um(T))− dum(t) dt t=T, and hence, Pm˜u=um(T), the Galerkin approximation at time T. Thus, one has only to compute the high-frequency modes ˜q=(I−P m )˜u. Applying (I−Pm) to (10) we obtain νA˜q=−(I−Pm)F(um),(11) which is the postprocessing step proposed in [21], [22]. Observe also that ˜q=Φ(u m )= (νA)−1(I−Pm)F(um). The formulation (11) of the postprocessing step is more natural in the case of spectral methods. However, when developing this paper, (7) (or (10)) was seen to be more meaningful and useful. This point of view was also reached independently in [17], [18]. In fact, in the case of spectral methods, the postprocessed Galerkin method was developed the opposite way. In the context of inertial manifolds [14] and approximate inertial manifolds, it was shown in [13] (see also [15]) that given the low modes Pmuof the solution uof (1)–(2), Φ(Pmu)=(νA)−1(I−Pm)F(Pmu) is a high-order approximation to the high-frequency modes (I−Pm)uof u. Based on this fact, the so-called nonlinear Galerkin methods (NLG) [12], [28], [32] were developed (see also [10], [11], [15], [29], [42], [43], and the many references therein). In these NLG methods, nonlinear terms are typically evaluated in the whole spectrum (or some generous truncation of it), for example, F(ym+Φ(ym)), where ym(t)∈Hmis the low-frequency Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 474 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI component of the NLG approximation ym+Φ(y m ). This is computationally costly, and (to add to the cost) it is done along the time evolution in NLG methods. Currently, there is enough available evidence [20], [21], [35], [46] that NLG methods are not competitive from the point of view of efficiency (or at least, not necessarily competitive). For this reason, the postprocessed Galerkin (spectral) method (9)–(10) was developed. The aim was to exploit the approximation capabilities of Φ but without paying for the computational cost of NLG methods (notice that in (9)–(10) nonlinear terms are applied only to elements of Hm, and Φ is computed only once). The postprocessed Galerkin method was seen to be computationally more efficient than classical Galerkin spectral methods in [16] and [21]. This is also the case of Galerkin finite-element methods, as is shown in section 3 and in [30]. The computational gain is even more significant in the maximum norm. The last section of this paper is devoted to the error analysis in L∞. We end this section by giving an idea of the proof of Theorem 1. Evaluate the PDE (1) at time t=Tand subtract it from (7). After rearranging terms and applying A−1one sees that u(T)−˜u=1 νA−1d dtuh(T)−d dtu(T)+1 νA−1(F(uh(T)) −F(u(T))). Thus, at first sight, u(T)−˜udepends on uh(T)−u(T) (and (duh/dt)−(du/dt)|t=T) which is only O(hr). However, in due course we will show that ° °A−1(F(uh(T)) − F(u(T)))° °L2(Ω) can be bounded in terms of ° °A−1(u(T)−uh(T))° °L2(Ω) (A−1/2instead of A−1when Fis given by (4)). We write u−uh=(u−R h u)+(R hu−u h). Since it is well known (see, e.g., [44]) that ° °A−1(I−Rh)u(T)° °L2(Ω) =O(hmin(2r−2,r+2)) (and the same bound for A−1(I−Rh)ut(T)), we have that ku(T)−˜ukL2(Ω) ≤C νð °A−1(uh(T)−Rhu(T))° °L2(Ω) +° ° ° ° A−1d dtuh(T)−Rh d dtu(T)° ° ° °L2(Ω)!+O(hmin(2r−2,r+2)).(12) The O(hmin(2r−2,r+2)) term above clearly shows that the postprocessed method will have a strictly higher convergence rate than the Galerkin method on which it is based, only if r≥3, that is, uh(T) is a piecewise polynomial with a degree of at least 2. We will duly show that the other two terms on the right-hand side of (12) are O(hr+µ|log(h)|r). The argument above shows that there is nothing especially nonlinear about the idea behind the postprocessed Galerkin method, and the results would be the same if Fwere linear. Notice that even in the linear case the method presented here allows an O(hr+2) (or O(hr+1)) error at the “same” price of a cheaper O(hr) method. In any case, since most of the interesting phenomena arise in nonlinear problems, for completeness we consider Fas in (3)–(4). The rest of this paper is organized as follows. Section 2 is devoted to proving Theorem 1. Numerical experiments are presented in section 3. Finally, we carry the L∞analysis in section 4. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 475 2. The postprocessed Galerkin method. 2.1. Preliminaries. We now present the assumptions under which Theorem 1 holds. We restrict ourselves to quasi-uniform meshes Th, so that the following inverse estimate holds for any element τ∈T hand vh∈Sh,r (see, e.g., [39]): kvhkWm,p(τ)≤Chl−m−d(1 q−1 p)kvhkWl,q (τ),0≤l≤m≤2,1≤q≤p≤∞.(13) Also, we assume the following restriction between the order of approximation r and the dimension d: d≤2r−1.(14) This is required by the presence of the nonlinear term F, but it can be loosened if Fis linear. Since finite elements are most advantageous when dealing with complex geometries, we consider only d≥2. Following the ideas presented here, the onedimensional case admits a much simpler treatment. In this situation, for 1 ≤p≤∞, and for any v∈D(A)∩Wr,p(Ω), an interpolant Ih(v)∈Sh,r exists such that kv−Ih(v)kLp(Ω) +hkv−Ih(v)kW1,p(Ω∩Ωh)≤C¡hr+δ(h)kvkWr,p (Ω) ,(15) (we extend by zero Ih(v)inΩ\Ω h ), where δ(h) = maxx∈Ωhdist(x, ∂Ω). For x∈Ω∩Ωh, this bound follows from the standard theory of interpolation and the Bramble–Hilbert lemma. For x∈Ω\Ωh,v(x) can be bounded using of the mean-value theorem (see, e.g., (19) below). However, for the purpose of analysis, we may assume (see, e.g., [39]) that Ωh⊆Ω, since one may consider a domain Ωδwith a smooth boundary such that Ω ∪Ωh⊂Ωδ, and max {dist(x, ∂(Ω ∪Ωh)) |x∈∂Ωδ}≤c 0 δ(h) for some c0>0. Then, for each t, replace u(t) by the solution uδ(t)of−∆u δ (t)=f u (t), where fu(t) is a suitable extension of −∆u(t)toΩ δ(see [39] for details). The sets Ωδmust be built so that for δsufficiently small we have that kvδkWl+2,p(Ωδ)≤CpkfvkWl,p(Ωδ)=Cpk∆vkWl,p(Ω) ,0≤l≤r, for every p∈[2,∞) and every v∈Wr+2(Ω). Here Cp>0 constant that depends on pbut not on δ. This can be readily done if Ω is bounded and r+ 2 smooth, or uniformly r+ 2 smooth [1, p. 67]. We will also assume that the solution of (1) satisfies M3= max 0≤t≤T° °A3/2u(·,t) ° °L 2(Ω) <+∞.(16) To ensure the estimate (16) for every solution of (3), (4), with initial data in D(A3/2), it is sufficient to require, for example, that both gand bsatisfy g(0) = 0 and b(0)=0. If, in addition, the right-hand side of (1) has a smooth forcing term f=f(x, t) instead of zero, then it is sufficient that f(x, t) has a zero trace at the boundary of Ω. In the understanding of the postprocessed Galerkin method, the elliptic projection Rhuof uwill play a prominent role. To simplify notation we write rh=Rhu. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 476 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI Under these conditions the following bounds hold for 1 ≤l≤r: ku−PhukL2(Ω) ku−rhkL2(Ω) +hku−rhkH1(Ω) ≤ChlkukHl(Ω) ,(17) and, for 0 ≤s≤2, ° °A−s/2(u−Phu)° °L2(Ω) ° °A−s/2(u−rh)° °L2(Ω) )≤Chl+skukHl(Ω) , s = min(s, r −2).(18) For these bounds to hold, care must be taken to approximate the boundary with the necessary accuracy (see, e.g., [44, section 5.1]). Henceforth, we will assume that for some constant cδ>0, δ(h) = max x∈Ωh dist(x, ∂Ω) ≤cδhr+2. We also remark that (16) plays a role in order to bound k∇(u−rh)kL2(Ω) in terms of hl−1kukHl(Ω) in (17), due to the skin-layer Ω\Ωh. To be more precise, let D⊂Ω be such that max{dist(x, ∂Ω) |x∈D}≤cδ(h). For any v∈D(A), since v=0on ∂Ω, one can check that kvkL2(Ω\D)≤Cδ(h)k∇vkL2(Ω) .(19) Hence, taking D=Ω h , we have that k∇vkL2(Ω\Ωh)≤Ckvk1/2 L2(Ω\Ωh)kvk1/2 H2(Ω\Ωh)≤Cδ(h)1/2kvkH2(Ω) . For v∈D(A3/2) better estimations can be obtained as follows. Consider Ω1⊂⊂ Ωhwith a sufficiently smooth boundary such that max{dist(x, ∂Ωh)|x∈Ω1}≤ c 0 δ(h) for some c0>0. Then k∇vkL2(Ω\Ωh)≤k∇vk L 2 (Ω\Ω1). Using the divergence theorem, trace theorems, and the smoothness of Ω, one can show that k∇vk2 L2(Ω\Ω1)≤ k∆vkL2(Ω\Ω1)(kvkL2(Ω\Ω1)+k∇vkL2(Ω\Ω1)). Now, recalling (19), we have that ° °∇v° °L2(Ω\Ωh)≤Cδ(h)kvkH3(Ω) . For the L2projection Phwe will use stronger estimates than those in (17). More precisely, there exists a constant C>0 such that kv−PhvkLq(Ω) ≤CkvkLq(Ω) ,1≤q≤∞,(20) ∀v∈Lq(Ω) (see, e.g., [44, p. 39]). Similarly, for the elliptic projection kRhk∞≤C(21) (see [39]), where, here and in the rest of the paper, k·kpdenotes the operator norm in L(Lp(Ω)) for 1 ≤p≤∞. We remark that for piecewise-linear finite elements (i.e., r= 2) the bound (21) does not hold (Cshould be replaced by C(1+|log(h)|)). In this paper, however, we restrict ourselves to r≥3. The error estimates in sections 2.2 and 2.3, as well as the constant Cin Theorem 1 depend on K(u, t)=ku(·,t)kHr(Ω) +kut(·,t)kHr(Ω) ,K(u) = max 0≤t≤TK(u, t).(22) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 477 Notice then, that we are assuming sufficient smoothness of the solution as well as of the initial condition. It must be stated that if (1) is autonomous and gand bin the nonlinear terms Fin (3)–(4) are analytic, due to the analyticity in time of the solutions of (1)–(2), for most initial conditions u0, it is possible to bound kut(·,t)kHr(Ω) in terms of maxz∈Sku(·,z)kHr(Ω) for zin an adequate neighborhood Sof [0,T] in the complex plane (see, e.g., [7, pp. 104–109]). We will frequently use the following version of Sobolev’s lemma: for p∈[1,∞), there exists a constant C=C(Ω,p) such that kvkLp0(Ω) ≤CkvkWs,p(Ω) for 1 p≥1 p0≥1 p−s dand v∈Ws,p(Ω) whenever p0<∞.Forp 0=∞the above inequality holds for 1 p≥1 p0>1 p−s d, and, furthermore, in this case vis also a continuous function. Also the following inequality, which is easily obtained from H¨older’s inequality, will be frequently used: |(f,vw))|≤kfk L 2 (Ω) kvkLp(Ω) kwkLq(Ω) ,1 p+1 q=1 2.(23) 2.2. Analysis of the Galerkin method. In this section we obtain bounds for uh−rh. The main result is Theorem 2 at the end of the section. Its proof will be obtained as a consequence of several previous lemmas. In the following, we need the residual or truncation error Thof the elliptic projection rh=rh(t)=R h u(t). Notice that rhsatisfies d dtrh+νAhrh+PhF(rh)=T h ,(24) where the truncation error This Th=Ph(F(rh)−F(u)) + d dtrh−Phut.(25) In general, kThkL2(Ω) is not smaller than kut−drh/dtkL2(Ω) =O(hr). In fact, the above is definitely true in the linear case when F= 0. On the other hand, notice that in Lemma 4 below, we bound not only kThkL2(Ω) but also an integral of e−ν(t−s)AhTh(s) with respect to s. It is the smoothing effect of the semigroup e−νtAh that gives the extra power hµ|log(h)|rin Theorem 2 with respect to the O(hr) decay of kThkL2(Ω). Lemma 1. Let vh:[0,T]→S h,r satisfy the threshold condition, krh(t)−vh(t)kL2(Ω) ≤cthhr,(26) for t∈[0,T]. Set µ=0if the nonlinearity Fis given by (3); otherwise set µ=1 when Fis given by (4). Then there exists a constant C, which depends on cth and max0≤t≤Tku(t)kHr(Ω), such that for t∈[0,T] ° °A −µ/2 hPh(F(vh(t)) −F(rh(t)))° °L2(Ω) ≤Ckvh(t)−rh(t)kL2(Ω) .(27) Furthermore, for Fgiven by (4) we have that kPh(F(vh(t)) −F(rh(t)))kL2(Ω) ≤C¡k∇(vh(t)−rh(t))kL2(Ω) +h−1kvh(t)−rh(t)kL2(Ω).(28) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 478 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI Proof. Let us take t∈[0,T], and for simplicity, let us drop the explicit dependence on t. We first notice that both vhand rhare bounded in L∞in terms of u, since kvhkL∞(Ω) ≤kv h−r h k L ∞ (Ω) +krhkL∞(Ω), which, thanks to the inverse estimate (13) and the threshold condition (26), is bounded by cthhr−d/2+krhkL∞(Ω). Then in view of the fact that Rhis bounded in L∞(recall (21)), Sobolev’s lemma, and the restriction (14) between the dimension dand rwe have that kvhkL∞(Ω),krhkL∞(Ω) ≤CkukHr(Ω). As a consequence, for any smooth function f, f(vh),f(r h)∈L p(Ω) for 1 ≤p≤∞.(29) Thus, for Fas in (3), the statement (27) is a direct consequence of the mean-value theorem. For Fas in (4) we first prove (27). Notice that A−µ/2 hPh=A1−µ/2 hA−1 hPh= A1−µ/2 hRhA−1. In addition, for χh∈Sh,r, we have that kA1/2 hχhkL2(Ω) =k∇χhkL2(Ω), and kA−1·kH1(Ω) ≤CkA−1/2·kL 2(Ω). Then, using (17) we have that ° °A−1/2 hPh(F(vh)−F(rh))° °L2(Ω)≤C° °A−1/2(F(vh)−F(rh))° °L2(Ω). Let us estimate the right-hand side above by duality. Since A−1/2is bounded, we need only to study the convection term. Let us denote eh=vh−rh. We have that b(vh)·∇v h−b(r h )·∇r h=b(v h )·∇e h+(b(v h )−b(r h )) ·∇r h .(30) Let φ∈C∞ 0(Ω) be a test function. Taking the inner product in (30) with φand integrating by parts the first term on the right-hand side above becomes ¡b(vh)·∇e h ,A −1/2φ =¡e h,div¡b(vh)A−1/2φ =¡eh,b 0(v h)·∇v h A −1/2φ +¡e h ,b(v h)·∇A −1/2φ .(31) We now take psuch that 1 p=1/4ifd=2 1/2−1/d if d>2. Then, by Sobolev’s lemma, we have that k·kL p(Ω) ≤Ck·kH 1(Ω), and, in particular, kA−1/2φkLp(Ω) ≤CkφkL2(Ω). We also take qsuch that 1/q +1/p =1/2. For this value of q, thanks to (14), it is immediate to check that kukW1,q(Ω) ≤CkukHr(Ω) .(32) Applying (23) to the right-hand side of (31), we have that ¡eh,b 0(v h)·∇v h A −1/2φ  ≤Cke h kL 2(Ω) kb0(vh)kL∞(Ω) k∇vhkLq(Ω) kφkL2(Ω) , ¡eh,b(v h)·∇A −1/2φ  ≤Cke h kL 2(Ω) kb(vh)kL∞(Ω) kφkL2(Ω) . Also, for the second term on the right-hand side of (30), thanks to the mean-value theorem and (29), arguing as above we have that ¡(b(vh)−b(rh))∇rh,A −1/2φ  ≤Cke hkL 2(Ω) k∇rhkLq(Ω) kφkL2(Ω) . Thus, (27) is proved, provided we show that k∇vhkLq(Ω) ,k∇rhkLq(Ω) ≤CkukHr(Ω) (33) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 485 Proof of Theorem 1. Subtracting (7) from (1) and applying νA−1we have that u(T)−˜u=ν−1A−1(F(uh(T)) −F(u(T))) + A−1duh dt (T)−du dt (T).(59) Since k∇·kL2(Ω) =° °A1/2·° °L2(Ω), then by Lemma 3 and an argument similar to that leading from (50) to (47), we obtain the following bounds: ku−˜ukL2(Ω) ≤C° °A−1(uh−u)° °L2(Ω) +° °A−1(uh−u)t° °L2(Ω) +2,(60) ku−˜ukH1(Ω) ≤C° °A−1/2(uh−u)° °L2(Ω) +° °A−1/2(uh−u)t° °L2(Ω) +1,(61) where for µ=1,2, µ=ku−uhkLq(Ω)ku−uhkHi(Ω), with qas indicated in Lemma 3 and i= 1 when Fis given by (4); otherwise i=0. Let us first find an upper bound for µ. We write u−rh=(u−r h )+(r h−u h ).(62) For rh−uh, using the inverse estimate (13) and Theorem 2, we have that krh−uhkLq(Ω) ≤Chd q−d 2krh−uhkL2(Ω) ≤CK(u)hr+µ+d q−d 2|log(h)|r. But d/q−d/2=−d/p, which, in view of (38), allows us to write r+µ−d/p ≥µ+µ−1, where we have used that r≥3 in the case d/2<µ. This and (52) lead to ku−uhkLq(Ω) ≤CK(u)hmin(µ+µ−1,µ)|log(h)|r. Using again (62), applying (17), (13), and Theorem 2, one can show that ku−uhkHi(Ω) ≤CK(u)hr−i, where i= 1 when Fis given by (4); otherwise i= 0. Hence, it follows that 2≤CK(u)hr+µ|log(h)|r, 1 ≤CK(u)hr+1−i|log(h)|r. Let us now turn to the other terms in (60)–(61). Using (62) and taking into account that A−µ/2is bounded, we have that ° °A−µ/2(u−uh)° °L2(Ω) ≤° °A−µ/2(u−rh)° °L2(Ω) +Ckrh−uhkL2(Ω) , ° °A−µ/2(u−uh)t° °L2(Ω) ≤° °A−µ/2(u−rh)t° °L2(Ω) +° °A−µ/2(rh−uh)t° °L2(Ω). Then, (17)–(18) and Lemma 6, together with (54), yield the bound (8). 3. Numerical experiments. We consider the reaction-diffusion system ∂u1 ∂t =1−4u 1+u 2 1 u 2+ν∆u 1 ∂u2 ∂t =3u 1 −u 2 1 u 2 +ν∆u 2  known as the Brusselator (see, e.g., [27]), in Ω ×[0,T], where Ω = [0,1]2and T= 10, subject to the Neumann boundary condition ∂u1 ∂n =∂u2 ∂n =0 on ∂Ω. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 486 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI 0 5 10 15 20 2.5 3 3.5 4 4.5 || u || L 2 t Fig. 1.Evolution of kukL2(Ω)2. 00.5 1 0 0.5 1 1 2 3 4 x y u 2 Fig. 2.Component u2of the solution at T=10. Although for simplicity we have developed the analysis for Dirichelet boundary conditions, our analysis can be modified to cover the Neumann boundary conditions case as well. We set ν=0.002. We complement this system with the initial condition u1(x, y, 0)=1−1 4(2y−1)((2y−1)2−3), u2(x, y, 0)=3−(2x−1)((2x−1)2−3). The corresponding solution becomes nearly periodic in time, and develops moderately large gradients. Figure 1 shows the evolution of kukL2(Ω)2up to time T= 20, and a periodic behavior can be observed. In Figure 2, the component u2of the solution, at time T= 10, is portrayed; the gradient takes values above 50 in some points of the domain. In our calculations we take the so-called regular pattern triangulations of Ω, which are induced by the set of nodes (j/N,k/N), 0 ≤j, k ≤N, where N=|Ω|/h is an integer. For the Galerkin method (5) we used Lagrange quadratic elements on these triangulations (i.e., r= 3) and for the postprocessing step, since the solution u= (u1,u 2) Tis smooth, we choose Lagrange cubic elements on the same triangulations. Since the Laplacian is noninvertible when subject to Neumann boundary conditions, for the postprocessing step we solve ν˜ Ah˜u1,h +4˜u 1,h =1+ ˜ P h (u 2 1,hu2,h)−d dtu1,h,(63) ν˜ Ah˜u2,h +˜ Ph¡˜u2 1,h ˜u2,h)=3˜u 1,h −d dtu2,h,(64) solving first (63) (notice that ν˜ Ah˜ +4Idis invertible) and then (64). Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 487 For the time integration, a variable-formula, variable-coefficient implementation of the backward differentiation formulae (BDF) was used with orders up to six (see, e.g., [26]). This method is widely used in the numerical integration of stiff ordinary differential equations (ODEs), like those arising from the spatial discretizations of dissipative PDEs. This BDF code was previously used in [16],[19], [20]. In [19], this code is shown to be more efficient than other widely used integrators in dissipative PDEs (including semi-implicit methods). Without an efficient time integrator, as shown in [19], [20], it may take an unaffordable amount of computing time to reach the levels of accuracy described below. For each hused in the triangulations of Ω, every experiment was carried out with different values of the tolerance TOL (an input value to the time integrator) below which local time discretization errors are desired. The smallest error obtained was selected for the plots. There was always a point, depending on h, at which further reducing the tolerance did not reduce the errors any more. This means that the error arising from time discretization is much smaller than the error arising from the finite-element spatial discretization. When this happens, one can be sure that the errors observed are due only to the spatial discretization and do not depend on the time integrator used. In view of [19] we were particularly careful to ensure that the dominant error in all the computations presented here was the spatial discretization error in order to avoid wrong conclusions from our numerical experiments. Also, for the computational cost in Figures 4–6 below, the largest tolerance among those with which the spatial discretization error was dominant was used. We refer to this value of the tolerance as the optimal tolerance. A smaller tolerance would result in a larger cost but no smaller error. A larger one, on the other hand, would imply a larger error, but then the same error could be obtained at less cost with the same tolerance and larger h. In Figure 3 we present a convergence diagram showing the errors committed by the methods when used with h=|Ω|/N ,N=12,24,48. By errors we mean the difference between the approximations (u1,h,u 2,h)or(˜u 1,h,˜u2,h) and the exact solution (at the end of the section we explain how it was obtained). In this and in the rest of the plots, the computational results of the classical quadratic elements are joined by continuous lines, and those of the postprocessed method by discontinuous (dashed) lines. With this, we get an idea of what the result would have been for intermediate values of h. It can be seen that as his decreased from left to right, the errors in both methods decrease from top to bottom, but they decrease faster in the postprocessed method. Measures of the slopes of the plots confirm the rates predicted by the theory (i.e., the errors in the plot decrease like Nslope = const.h−slope). The improved convergence rate of the postprocessed method observed in Figure 3 also implies improved efficiency. This can be checked in Figure 4, where the same errors as in Figure 3 are plotted against the smallest amount of CPU time that the methods needed to achieve them, that is, the cost of integrating the method with the optimal tolerance (notice that, for a given h, the optimal tolerance of the postprocessed quadratics is, in general, smaller than that of the quadratics). The reason for this improvement is that the error of the Galerkin method is reduced when the postprocessing is done, but this is done at very little cost (less than 10% of the whole computation, including the extra cost of integrating in time with smaller tolerances). In fact, the cost of the postprocessing step is equivalent to the cost of roughly 10 time steps, and the number of time steps in the results shown in Figure 4 varies between 109 and 344. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 488 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI 101102 105 104 103 102 N=|Ω|/h L2 relative errors Brusselator with nu=0.002: T=10 slope=3.0563 slope=3.8084 Quadratics Postprocessed Quadratics N=48 N=12 Fig. 3.Convergence diagram. 101102103 105 104 103 102 CPU time L2 relative errors Brusselator with nu=0.002: T=10 Quadratics Postprocessed Quadratics N=48 N=24 N=12 Fig. 4.Efficiency diagram. The postprocessed quadratic elements are not only more efficient than the classical quadratics but also more efficient than the classical cubic elements, which have the same convergence rate as the postprocessed quadratics. In Figure 5 we have plotted, besides the results of Figure 4, those of the classical cubic elements, that is, the standard Galerkin method (5)–(6) but with r= 4. These are represented by asterisks joined by dotted lines. Notice how the discontinuous line (postprocessed method) is always on the left (less CPU time) of the dotted line of the cubics. In this example, Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 489 101102103 105 104 103 102 CPU time L2 relative errors Brusselator with nu=0.002: T=10 Quadratics Postprocessed Quadratics Cubics N=6 N=12 N=12 N=24 N=48 N=48 N=48 Fig. 5.Efficiency diagram. 101102103 104 103 102 101Brusselator with nu=0.002: T=10 L∞ relative errors CPU time Cubics Postprocessed Quadratics N=6 N=12 N=12 N=24 N=24 N=48 N=48 Fig. 6.Efficiency diagram in L∞. the postprocessed quadratics save between 15% and 17% of the cost of the cubics. More significant is the gain in the maximum norm in Figure 6, where the postprocessed quadratics require between 21% and 25% less computing time than the classical cubics. A significant detail can be observed in Figures 5 and 6: The postprocessed method is more efficient than classical methods for quite large errors (already above 10−3). This is in contrast to our previous experience with spectral methods, where, typically, postprocessing brings a substantial improvement only for quite small errors (below 10−5) [16], [21]. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 490 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI 0 2 4 6 8 10 0 5 10 15 t iterations CPU time = 1082 s CPU time = 739 s CPU time = 780 s Fig. 7.Evolution of PCG iterations for h=√2/48 and TOL =0.001: −·−MAXIPCG =8,−− MAXIPCG =11,−− MAXIPCG = 18818. We end this section by commenting on the more relevant details of the methods used. Recall that the time integration is done with the BDF. Since these are fully implicit methods, then, at every time step, a system of nonlinear equations of the form αvh+(δt)(νAhvh+PhF(vh))=0,(65) has to be solved, where δt is the time steplength. This was done by Newton iteration. The exact Jacobian matrices were replaced by the Jacobian matrix of the linear terms, that is, αMh+(δt)νSh,(66) where Mhand Share the mass and stiffness matrices (see, e.g., [45]). The corresponding linear systems were solved by the preconditioned conjugate gradient method (PCG) with incomplete Cholesky (ICH) factorization of (66) as preconditioner (see, e.g., [23]). Both PCG and ICH algorithms were taken from the SLAP library. A single ICH factorization of (66) proved to be extremely costly (between 5% to 10% of the whole computation from t=0toT= 10). Thus, it was computed at the first step and reused in later steps until the PCG did not converge, when it was recalculated again. By convergence in the PCG method we mean achieving a relative error below a prescribed value TOLPCG in less than a fixed number of iterations MAXIPCG. In our present case, we chose TOLPCG =10 −3for the linear systems in the Newton iteration of (65), and 10−9for those in the postprocessing (63)–(64). To set the value of MAXIPCG, we performed a practical study for various values of hand settle it to MAXIPCG = 11, which usually forced the computation of two ICH decompositions. In Figure 7, the number of PCG iterations (the maximum for the two or three linear systems solved per step) for different values of MAXIPCG are shown, together with the CPU times of the corresponding runs. It can be seen how if MAXIPCG is generous (dashed line) the code computes only one ICH factorization at t= 0 but in the end performs a crippling number of iterations. If MAXIPCG is too tight (dash-dotted line) the second ICH decomposition is computed too soon, at t=1.6662; the number of PCG iterations is drastically reduced from then onward, but the overall number of iterations (and hence the cost) is larger than for MAXIPCG = 11 (continuous line). Notice how in this last case, the number of iterations is drastically reduced at t= 2.5005 when the second ICH factorization is computed, and remains the lowest from then onward. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 491 For the linear systems in the postprocessing step (63)–(64), it proved far more effective to use the diagonal of the corresponding matrices as a preconditioner. A far greater number of PCG iterations were taken (in the example of Figure 7, 70 for ˜u1,h and 308 for ˜u2,h for systems of 21025 unknowns), but it compensated for the cost of a new ICH factorization (notice that the matrices in the postprocessing step are of larger dimension than those in (66)). The initial guess for Newton’s method in (65) was computed by standard techniques of extrapolation in the BDF code (see, e.g., [26]), and for the PCG method in the corresponding linear systems the initial iterates were set to 0 (notice that the solutions of this systems are the increments of Newton iterations that are supposed to be small). For the linear system in the postprocessing step (63)–(64), the initial guess was the Galerkin approximation uhexpressed in the nodal basis of the cubic elements. It may be asked if other linear algebra techniques would have changed the results shown here. This would be especially the case if a faster (if existing) linear algebra would have favored the cubic elements but not the quadratics (and hence, neither the postprocessed method), although it is difficult to devise why it should be so. Obviously, it is out of the scope of this paper to test all possible linear system solvers, especially when there does not seem to be a clear-cut criterion for when to apply the many existing techniques competing for attention. The PCG used here seems to us a very reasonable choice, especially if we take into account the very low number of iterations that it took to solve the systems in our examples (recall Figure 7). Furthermore, a significant improvement could have been made to the postprocessing step with multilevel techniques (see, e.g., [4]), which would have favored yet further the postprocessed Galerkin method. We tested direct (sparse) solvers, but they were remarkably more costly than the PCG. We also checked that explicit time integrators were much less efficient than the BDF (i.e., the system of ODEs (5) is genuinely stiff). In all the results shown here, the initial condition (6) was replaced by the interpolant Ih(u0), which is much cheaper to compute than the elliptic projection in (6). We also performed the experiments with (6), and observed that the differences with those shown here were insignificant. The theoretical (“exact”) solution uwas computed with cubic elements and h=√2/96, using values of the tolerance ranging from 10−9to 10−12, so that the computed solution taken as exact is reasonably more accurate than those shown in the experiments. All experiments were carried out on a SUN Ultra-1, Mod. 140 workstation with 64 MB RAM, under Solaris 2.5.1. All programs were written in FORTRAN, used double precision arithmetic, and were compiled with the SparcWorks 4.0 compiler with the -fast option. 4. L∞analysis. 4.1. Preliminaries and main result. In comparison to section 2, we further restrict ourselves here in the following sense. 1. The results here apply only to reaction-diffusion equations with nonlinear terms as in (3). 2. We consider only problems in Rd, with d=2,3. (As with the L2analysis, the case d= 1 is particular and easier to study). The error estimates in this section depend on K∞(u, t)=ku(·,t)kWr,∞(Ω) +kut(·,t)kWr,∞(Ω) ,K ∞ (u) = max 0≤t≤TK∞(u, t).(67) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 492 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI Under these conditions, the error bounds of the postprocessed Galerkin method are stated in the following theorem whose proof is given in section 4.3. Theorem 3. There exist constants C=C(K∞(u)) >0and h0>0such that for every h∈(0,h 0]the solution ˜uof (7) satisfies the following bound: ku−˜ukL∞(Ω) ≤Chr+2−d/2|log(h)|4−d/2if r=3.(68) ku−˜ukL∞(Ω) ≤Chr+2 |log(h)|3if r≥4.(69) We now state the equivalent to (17)–(18) in the L∞norm. Instead of (17) we have ku−rhkL∞(Ω) +hku−rhkW1,∞(Ω∩Ωh)≤ChlkukWl,∞(Ω) ,1≤l≤r.(70) This estimate follows from Theorem 5.1 in [39], (15), and (13). We are unaware of an estimate like (18) in L∞, although similar results in Lpcan be found in the literature. Thus, we prove the following result. Lemma 7. Assume that d≥2,r≥3,and δ(h)=O(h min(r+2,2(r−1))). Then there exists a constant C>0such that for v∈Wr,∞(Ω) the following bound holds: ° °A−1(v−Rhv)° °L∞(Ω) ≤CkvkWr,∞(Ω) hmin(r+2,2(r−1)) |log(h)|.(71) Proof. Following [39], we take χh=Ih¡A−1(I−Rh)v) the interpolant of A−1(I− Rh)v. Then ° °A−1(I−Rh)v° °L∞(Ω) ≤° °A−1(I−Rh)v−χh° °L∞(Ω) +kχhkL∞(Ω) . Let us take pwith 2 ≤p≤∞. For the first term on the right-hand side above we have that (cf. [6, Thm. 3.1.6]) ° °A−1(I−Rh)v−χh° °L∞(Ω) ≤Ch2−d/p° °A−1(I−Rh)v° °W2,p(Ω). Let us recall the following inequality (see, e.g., [39] and [49]) which, for linear elliptic equations, illustrates the rate at which the Agmon–Douglis–Nirenberg regularity and well-posedness estimate [2] deteriorates as p→∞. Namely, there exists a constant C>0, which is independent of p, such that ° °A−1w° °W2,p(Ω) ≤CpkwkLp(Ω) ,2≤p<∞.(72) Thus, it follows that ° °A−1(I−Rh)v−χh° °L∞(Ω) ≤Cph2−d/p k(I−Rh)vkLp(Ω) ≤Cphr+2−d/p kvkWr,p (Ω) . Using (13) we have that kχhkL∞(Ω) ≤Ch−d/pkχhkLp(Ω), and thus kχhkL∞(Ω) ≤Ch−d/p° °χh−A−1(I−Rh)v° °Lp(Ω) +° °A−1(I−Rh)v° °Lp(Ω).(73) In view of the interpolation error bound (15), (70) and (72) we obtain the following estimate for the first term on the right-hand side of (73): Ch−d/p° °χh−A−1(I−Rh)v° °Lp(Ω) ≤Cphr+2−d/p kvkWr,p (Ω) . Also, it is well known (see, e.g., [40, Example 3.2]) that Ch−d/p° °A−1(I−Rh)v° °Lp(Ω) ≤Cphmin(r+2,2(r−1))−d/p kvkWr,p (Ω) ,(74) with Cindependent of pand h. Taking p=|log(h)|, (71) follows easily. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php POSTPROCESSING FINITE-ELEMENT METHODS 493 4.2. L∞analysis of the Galerkin method. The following resolvent estimate, which is due to Palencia [38], will be needed later for our L∞error estimates. Lemma 8. Fix α∈(0,π/2) and h0∈(0,e −1].Letβ(h)denote β(h)=k(P h− R h )A −1 k ∞(= O(h2|log(h)|)).Then there exist positive constants C∞and σ<1) such that the bounds ° °(zI −Ah)−1° °∞≤C∞|z|−1if |z|≤σβ(h)−1, C∞h−d/2|z|−1if |z|≥0 hold for z∈Cwith π≥|arg(z)|≥αand for 0<h≤h 0 . Taking for example a=(d+3)/(2σcos(α)), it easily follows that there exists a constant C0>0 such that for hsufficiently small ° °e−tAhχh° °L∞(Ωh)≤C0kχhkL∞(Ωh),χ h ∈S h ,t≥r 0 =aβ(h)|log(h)|.(75) We will also use the following bound from Nitsche and Wheeler [36], ° °e−tAhχh° °L∞(Ωh)≤C1kχhkL∞(Ωh),χ h ∈S h ,t≥0,r≥4.(76) The lack of a bound like (76) ∀t≥0, if r= 3, is the reason that the convergence rate depends also on the dimension d. This will be seen in Lemma 11 below. Lemma 9. Fix cth >0and T>0. Then there exists a constant S>0such that for vh:[0,T]→S h,r satisfying vh(0) = rh(0) and the threshold condition krh(t)−vh(t)kL∞(Ωh)≤cth,(77) for t∈[0,t 1]where t1≤T, the following bound holds for hsufficiently small: max 0≤t≤t1kvh(t)−rh(t)kL∞(Ωh)≤Smax 0≤t≤t1° ° ° °Zt 0 e−ν(t−s)Ah(ˆ Th(s)−Th(s)) ds.° ° ° °L∞(Ωh) , (78) where ˆ Tis as in (35). Proof. We notice that the threshold condition (77) and the fact that Rhis bounded ensures that both kvhkL∞(Ω) and krhkL∞(Ω) are bounded by some constant C= C(cth,max0≤t≤Tku(t)kL∞(Ω)). Since Fis smooth, it follows that kF(vh)−F(rh)kL∞(Ωh)≤Ckvh−rhkL∞(Ωh).(79) Now take norms in (36). For r≥4, thanks to (76), the fact that Phis bounded in L∞(recall (20)) and (79) above, standard arguments with the Gronwall lemma lead to (78). For r= 3, similar arguments (with a generalized Gronwall lemma as in Lemma 2) will also lead to (78) as long as we prove ° ° ° °Zt 0 e−ν(t−s)AhPh(F(vh)−F(rh)) ds° ° ° °L∞(Ωh)≤CZt 0 kvh(s)−rh(s)kL∞(Ωh) (t−s)γds, (80) for some C>0 and γ<1. To do this, we consider the value r0in (75) and let us call t= max(0,t−r 0/ν), so that for s∈[0, t], ν(t−s)≥r0. Then, thanks to (75), (79), and the fact that Phis bounded, we have that ° ° ° °Zt t e−ν(t−s)AhPh(F(vh)−F(rh)) ds° ° ° °L∞(Ωh)≤CZt tkvh(s)−rh(s)kL∞(Ωh)ds. (81) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 494 BOSCO GARC´ IA-ARCHILLA AND EDRISS S. TITI On the other hand, take p=dand κ=d/2−1 so that 1/p =1/2−κ/d. Then, for s∈[t, t] and wh∈Sh,r, thanks to the inverse estimate (13) and Sobolev’s lemma, we have that ° ° °e−ν(t−s)Ahwh° ° °L∞(Ωh)≤C h° ° °e−ν(t−s)Ahwh° ° °Lp(Ωh)≤C h° ° °e−ν(t−s)Ahwh° ° °Hκ(Ωh) ≤C h° ° °e−ν(t−s)Ahwh° ° ° 1−κ L2(Ωh)° ° °A1/2 he−ν(t−s)Ahwh° ° ° κ L2(Ωh) ≤C h¡ν(t−s)κ/2kwhkL2(Ωh).(82) Now observe that for s∈[t, t], since ν(t−s)≤r0≤Ch3/2, it follows that h−1≤ C(ν(t−s)−2/3). Thus, taking into account that κ 2+2 3=d 4−1 2+2 3=d 4+1 6≤3 4+1 6=11 12 =γ<1, we have that ° ° °e−ν(t−s)Ahwh° ° °L∞(Ωh)≤C ¡ν(t−s)11 12 kwhkL2(Ωh)≤C|Ωh|1/2 ¡ν(t−s)11 12 kwhkL∞(Ωh). This (replacing whby Ph(F(vh)−F(rh))) together with (79) and (81) leads to (80)with γ=11/12. Having established stability in the previous lemma, we now face the task of estimating the right-hand side of (78) when vh=uh. Recall that in the L2case, this was straightforward thanks to the very favorable L2-bounds of Ahe−µ(t−s)Ah, which can be obtained by means of the spectrum of Ah. In the present case, i.e., the L∞case, the lack of such bounds for e−µ(t−s)Ah(let alone Ahe−µ(t−s)Ah) makes the task much harder. Since the difficulty lies when t−sis small (recall for example the restriction in (75)), we will decompose the integrals in (78) into two parts corresponding to t−s large (Lemma 10 below) and small (Lemma 11 below). Lemma 10. Fix α∈(0,π/2),let σand β(h)be as in Lemma 8, and set τ1= dβ(h)|log(h)|/(2νσcos(α)).Then there exists a constant C>0such that for t≥τ1 the following bound holds for hsufficiently small: ° ° ° °Zt−τ1 0 eν(t−s)AhTh(s)ds° ° ° °L∞(Ωh)≤C ν log t h max 0≤s≤t−τ1° °A−1 hTh(s)° °L∞(Ωh).(83) Proof. By writing eν(t−s)AhTh(s)=A h e ν(t−s)A hA −1 hT h (s), we only need to bound Aheν(t−s)Ah. Let Γ be the boundary of the sector {z||arg(z)|≤α}, endowed with the positive orientation. Then Ahe−ν(t−s)Ah=1 2πi ZΓ ze−ν(t−s)z(zI −Ah)−1dz. We take ρ1=σβ(h)−1, and divide Γ in the two parts Γ ∩{|z|≤ρ 1 }and the two half-lines Γ ∩{|z|≥ρ 1 }. Taking into account the resolvent estimate in Lemma 8, we can write ° °Ahe−ν(t−s)Ah° °∞≤2C∞Zρ1 0 e−ν(t−s) cos(α)ρdρ +2C ∞ h d/2Z∞ ρ1 e−ν(t−s) cos(α)ρdρ, Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php