Robust error bounds for the Navier-Stokes equations using implicit-explicit second-order BDF method with variable steps

Robust error bounds for the Navier-Stokes equations using implicit-explicit second order BDF method with variable steps Bosco Garc´ıa-Archilla∗Julia Novo† August 24, 2022 Abstract This paper studies fully discrete finite element approximations to the Navier-Stokes equations using inf-sup stable elements and grad-div stabilization. For the time integration two implicit-explicit second order backward differentiation formulae (BDF2) schemes are applied. In both the Laplacian is implicit while the nonlinear term is explicit, in the first one, and semi-implicit, in the second one. The grad-div stabilization allows us to prove error bounds in which the constants are independent of inverse powers of the viscosity. Error bounds of order rin space are obtained for the L2error of the velocity using piecewise polynomials of degree rto approximate the velocity together with second order bounds in time, both for fixed time step methods and for methods with variable time steps. A CFL-type condition is needed for the method in which the nonlinear term is explicit relating time step and spatial mesh sizes parameters. Keywords: Incompressible Navier-Stokes equations; Variable step BDF2; Implicitexplicit methods; Grad-div stabilization; Robust error Bounds 1 Introduction The numerical simulation of the Navier-Stokes equations is still a challenge in which there are several questions that deserve some research. In the present paper we try to contribute to the following two aspects. On the one hand, we consider the possibility of optimize the cost of the temporal integration. On the other hand, we consider spatial approximations that allow to prove error bounds with constants independent of inverse powers of the viscosity. Concerning the first question, following the recent reference [9], we consider implicit-explicit methods with explicit or semi-implicit treatment of the nonlinear term, which reduces the cost of every single time step. Also, we include the error analysis of the variable step case. As stated in [18]: “Adaptive time stepping is an important tool in Computational Fluid Dynamics for controlling the accuracy of simulations and for enhancing their efficiency”. Concerning the second question, we add grad-div stabilization in the spatial discretization of the method. As in [7], adding this stabilization we are able to get bounds with constants independent of inverse powers of the viscosity. This is one of the methods considered in [12] in which methods with robust bounds (error bounds with constants independent of the ∗Departamento de Matem´atica Aplicada II, Universidad de Sevilla, Sevilla, Spain. Research is supported by Spanish MCINYU under grants PGC2018-096265-B-I00 and PID2019-104141GB-I00 (b[email protected]) †Departamento de Matem´aticas, Universidad Aut´onoma de Madrid, Spain. Research is supported by Spanish MINECO under grants PID2019-104141GB-I00 and VA169P20 (julia.no[email protected]) 1 Reynolds number are studied). Methods for whom robust estimates can be derived enable stable flow simulations for small viscosity coefficients on comparatively coarse grids. Let Ω ⊂Rd,d∈ {2,3}, be a bounded domain with polyhedral and Lipschitz boundary ∂Ω. The incompressible Navier–Stokes equations model the conservation of linear momentum and the conservation of mass (continuity equation) by ∂tu−ν∆u+ (u·∇)u+∇p=fin (0, T)×Ω, ∇·u= 0 in (0, T)×Ω,(1) u(0,·) = u0(·) in Ω, u(t, x) = 0 in (0, T )×∂Ω, where uis the velocity field, pthe kinematic pressure, ν > 0 the kinematic viscosity coefficient, u0a given initial velocity, and frepresents the external body accelerations acting on the fluid. In [7] a stabilized method using only grad-div stabilization to approach the evolutionary Navier-Stokes equations is considered and analyzed. Error bounds with constants independent of inverse powers of the viscosity are proved. Using piecewise polynomials of degree r for the velocity bounds of order rare obtained for the L2error of the velocity. Both the continuous-in-time case and the fully discrete scheme with the fully implicit backward Euler method as time integrator are analyzed. Recently, in [9], implicit-explicit (IMEX) variable stepsize methods for the evolutionary Navier-Stokes equations are analyzed. In these algorithms the nonlinear term is treated fully explicitly while the remaining terms are treated implicitly. Stability and convergence is proved for a variable stepsize first order method. However, the error bounds in [9] depend strongly on inverse powers of ν. In particular, a CFL-type condition is needed for the time step depending on ν−1. The authors state that the analysis including grad-div stabilization remains being an open problem. In the present paper, as in [7], we add grad-did stabilization to the mixed finite element approximation to be able to get bounds independent of ν−1. We analyze a second order BDF2 method in time, with explicit treatment of the nonlinear term. We consider both the case of fixed time step and the case of variable time step. A CFL-type condition is needed in our error analysis. This condition is stronger than the usual one in which (∆t)h−1has to be bounded (hbeing the mesh size) because we required (∆t)h−2to be bounded (a similar CFL condition is obtained for the variable step case). However, we improve the condition in [9] in which the quantity that has to be bounded is (∆t)h−1ν−1. Then, in principle, the method we consider could be applied in case of high Reynolds numbers without the need of reducing the size of the time step. Due to the stronger CFL condition of the method with explicit treatment of the nonlinear term, depending on the problem, it could be worth to use instead a semi-implicit form for the nonlinear term. In this paper, we also include the error analysis of such a method in which the nonlinear term is semi-implicit. For this method the error analysis is obtained with the same tools as the previous analysis but it is much simpler. Also, no CFL condition is needed. The same kind of bounds as for the previous method are proved in which the constants are also independent of inverse powers of the viscosity. Concerning the BDF2 method with variable time step, we did not find in the literature any reference with the error analysis for the Navier-Stokes equations for large Reynolds numbers. For low Reynolds numbers and Fourier pseudospectral methods, the analysis can be found in [24], where they treat the two dimensional case in the vorticity-streamfunction formulation. Contrary to the present paper, the error estimates in [24] depend on ν−1 (indeed, on exp(ν−1)) and the stepsizes have to be smaller than ν, so that they only apply when νis not too small (i.e., small Reynolds number). Dealing with other equations, in [2] a second order backward difference method with variable steps is analyzed for a linear parabolic problem. The error analysis for semilinear 2 parabolic problems can be found in [10] (see also [23]). The convergence of the variable two-step BDF time discretisation of nonlinear evolution problems governed by a monotone potential operator is studied in [11] while the analysis for the Cahn-Hilliard equation appears in [4]. Also the case of parabolic integro-differential equation arising in finance discretized with finite differences can be found in [22]. The outline of the paper is as follows. In Section 2 we introduce some notation and state some preliminaries. In Section 3 we carry out the error analysis of the fully discrete methods, both for the fixed step and the variable time step cases. Some technical lemmas are stated and proved along the section. Although the error analysis of the fixed step can be obtained as a particular case of the variable step we decided to start with the analysis of the simpler case since most of the ideas are common and in this way, in our opinion, the skeleton of the analysis can be better understood. We prove optimal convergence of second order in time for the fully discrete methods (both for the fixed and variable time step size cases and both for explicit and semi-implicit treatment of the nonlinear term) and the same rate of convergence in space as in [7]. In Section 4 we show some numerical experiments while some conclusions are stated in the last section. 2 Preliminaries and notation Throughout the paper we will denote by Ws,p(D) the Sobolev space of real valued functions defined on the domain D⊂Rdwith distributional derivatives of order up to sin Lp(D), endowed with the usual norm which is denoted by ∥·∥Ws,p(D). If sis not a positive integer, Ws,p(D) is defined by interpolation [1]. If s= 0 we understand that W0,p(D) = Lp(D). As it is standard, Ws,p(D)dwill be endowed with the product norm that, if no confusion can arise, will be denoted again by ∥·∥Ws,p(D). We will distinguish the case p= 2 using Hs(D) to denote the space Ws,2(D). We will make use of the space H1 0(D), the closure in H1(D) of the set of infinitely differentiable functions with compact support in D. For simplicity, we use ∥·∥s(resp. |·|s) to denote the norm (resp. seminorm) both in Hs(Ω) or Hs(Ω)d. The exact meaning will be clear by the context. The inner product of L2(Ω) or L2(Ω)dwill be denoted by (·,·) and the corresponding norm by ∥·∥0. The norm of the space of essentially bounded functions L∞(Ω) will be denoted by ∥·∥∞. For vector valued function we will use the same conventions as before. We represent by ∥·∥−1the norm of the dual space of H1 0(Ω) which is denoted by H−1(Ω). As usual, we always identify L2(Ω) with its dual so we have H1 0(Ω) ⊂L2(Ω) ⊂H−1(Ω) with compact injection. Using the function spaces V=H1 0(Ω)d, and Q=L2 0(Ω) = q∈L2(Ω) : (q, 1) = 0, the weak formulation of problem (1) is: Find (u, p) : (0, T]→V×Qsuch that for all (v, q)∈V×Q, (∂tu,v) + ν(∇u,∇v) + ((u·∇)u,v)−(∇·v, p) + (∇·u, q) = (f,v).(2) The Hilbert space Hdiv ={u∈L2(Ω)d,∇·u∈L2(Ω),| ∇·u= 0,u·n|∂Ω= 0} will be endowed with the inner product of L2(Ω)dand the space Vdiv ={u∈V| ∇·u= 0} with the inner product of V. The following Sobolev’s embedding [1] will be used in the analysis: For 1 ≤p < d/s let qbe such that 1 q=1 p−s d. There exists a positive constant C, independent of s, such that ∥v∥Lq′(Ω) ≤C∥v∥Ws,p(Ω),1 q′≥1 q, v ∈Ws,p(Ω).(3) 3 If p > d/s the above relation is valid for q′=∞. A similar embedding inequality holds trivially for vector valued functions. Let Vh⊂Vand Qh⊂Qbe two families of finite element spaces composed of piecewise polynomials of degrees at most kand l, respectively, that correspond to a family of partitions Thof Ω into mesh cells with maximal diameter h. In the following, we consider pairs of finite element spaces that satisfy the inf-sup condition inf qh∈Qh sup vh∈Vh (∇·vh, qh) ∥∇vh∥0∥qh∥0≥β0,(4) with β0>0a constant independent of the mesh size h. It will be assumed that the meshes are quasi-uniform and that the following inverse inequality holds for each vh∈Vh, see e.g., [5, Theorem 3.2.6], ∥vh∥Wm,p(K)≤cinvhn−m−d1 q−1 p K∥vh∥Wn,q (K),(5) where 0 ≤n≤m≤1, 1 ≤q≤p≤ ∞ and hKis the size (diameter) of the mesh cell K∈ Th. The space of discrete divergence-free functions is denoted by Vdiv h={vh∈Vh|(∇·vh, qh) = 0 ∀qh∈Qh}. Denoting by πhthe L2(Ω) projection of the pressure ponto Qh, we have that for 0 ≤ m≤1 ∥p−πh∥m≤Chl+1−m∥p∥l+1,∀p∈Hl+1(Ω).(6) In the error analysis, the Poincar´e–Friedrichs inequality ∥v∥0≤C|Ω|1/d ∥∇v∥0,∀v∈H1 0(Ω)d,(7) will be used. In the analysis, the Stokes problem −ν∆u+∇p=gin Ω, u=0on ∂Ω,(8) ∇·u= 0 in Ω, will be considered. Let us denote by (uh, ph)∈Vh×Qhthe mixed finite element approximation to (8), given by ν(∇uh,∇vh)−(∇·vh, ph) = (g,vh)∀vh∈Vh (∇·uh, qh) = 0 ∀qh∈Qh.(9) Following [13], one gets the estimates h−1∥u−uh∥0+∥u−uh∥1≤Cinf vh∈Vh∥u−vh∥1+ν−1inf qh∈Qh∥p−qh∥0,(10) ∥p−ph∥0≤Cνinf vh∈Vh∥u−vh∥1+ inf qh∈Qh∥p−qh∥0.(11) It can be observed that the error bounds for the velocity depend on negative powers of ν. For the analysis, it will be advantageous to use a projection of (u, p) into Vh×Qh with uniform in ν, optimal, bounds for the velocity. In [6] a projection with this property was introduced. Let (u, p) be the solution of the Navier–Stokes equations (1) with u∈ V∩Hk+1(Ω)d,p∈Q∩Hk(Ω), k≥1, and observe that (u,0) is the solution of the Stokes problem (8) with right-hand side g=f−∂tu−(u·∇)u−∇p. (12) 4 Denoting the corresponding Galerkin approximation in Vh×Qhby (sh, lh), one obtains from (10)-(11) ∥u−sh∥0+h∥u−sh∥1≤Chk+1∥u∥k+1,(13) ∥lh∥0≤Cνhk∥u∥k+1,(14) where the constant Cdoes not depend on ν. Remark 1 Assuming the necessary smoothness in time and considering (8) with g=∂t(f−∂tu−(u·∇)u−∇p), one can derive an error bound of the form (13) also for ∂t(u−sh). One can proceed similarly for higher order derivatives in time. Following [3], one can also obtain the following bounds for sh ∥u−sh∥∞≤C∞hlog(|Ω|1/d /h)¯r∥∇u∥∞,(15) ∥∇(u−sh)∥∞≤C∞∥∇u∥∞,(16) where C∞does not depend on νand ¯r= 0 for linear elements and ¯r= 1 otherwise. 3 The IMEX and semi-implicit BDF2 methods in time with grad-div stabilization in space In the sequel we well call IMEX to the method with explicit treatment of the nonlinear term and semi-implicit to the method with semi-implicit nonlinear term. In this section we state and analyze a fixed step and a variable step versions of both an IMEX and a semi-implicit BDF2 methods with grad-div stabilization in space. We first introduce the grad-div stabilization and a couple of Gronwall lemmas that will be applied for the error analysis of the methods. Then, we introduce the fully discrete schemes. Some lemmas are then proved that will be used both for the fixed and variable step cases analyzed in next two subsections. In Section 3.1 we carry out the error analysis for the case in which the step size is fixed. The variable step size case is analyzed in Section 3.2. The ideas in both sections are similar, the second one having more technical difficulties. For this reason, we think it is convenient to show first the analysis for the fixed step size, although this case can be obviously obtained as a particular case of the error analysis of Section 3.2. The spatial discretization that will be studied for the approximation of the solution of the Navier–Stokes equations (1) is obtained by adding to the Galerkin equations a control of the divergence constraint (grad-div stabilization). More precisely, the following grad-div method will be considered: Find (uh, ph) : (0, T]→Vh×Qhsuch that for all (vh, qh)∈Vh×Qhone has (∂tuh,vh) + ν(∇uh,∇vh) + b(uh,uh,vh)−(ph,∇·vh) + (∇·uh, qh) + µ(∇·uh,∇·vh) = (f,vh),(17) with uh(0) given. Here, and in the rest of the paper, b(u,v,w) = (B(u,v),w)∀u,v,w∈H1 0(Ω)d, where, B(u,v) = (u·∇)v+1 2(∇·u)v∀u,v∈H1 0(Ω)d. Notice the well-known property b(u,v,w) = −b(u,w,v)∀u,v,w∈V, (18) such that, in particular, b(u,w,w) = 0 ∀u,w∈V. We will use the following discrete Gronwall inequality whose proof can be found in [15]. 5 Lemma 1 Let k, B, aj, bj, cj, γjbe nonnegative numbers such that an+k n X j=0 bj≤k n X j=0 γjaj+k n X j=0 cj+B, for n ≥0. Suppose that kγj<1, for all j, and set σj= (1 −kγj)−1. Then an+k n X j=0 bj≤exp k n X j=0 σjγj!(k n X j=0 cj+B), for n ≥0. We will also apply the following Gronwall lemma. Lemma 2 Let B, aj, bj, cj, γj, δjbe nonnegative numbers such that an+ n X j=0 bj≤γnan+ n−1 X j=0 (γj+δj)aj+ n X j=0 cj+B, for n ≥0. Suppose that γj<1, for all j, and set σj= (1 −γj)−1. Then an+ n X j=0 bj≤exp σnγn+ n−1 X j=0 (σjγj+δj)!( n X j=0 cj+B), for n ≥0. Proof We follow the proof of [15, Lemma 5.1]. Let dnbe defined by the relation an+ n X j=0 bj+dn=γnan+ n−1 X j=0 (γj+δj)aj+ n X j=0 cj+B, for n ≥0, and let us denote by Snthe left hand side above. Notice then Sn−Sn−1=γnan+δn−1an−1+cn≤γnSn+δn−1Sn−1+cn, n ≥1. For n= 0, we have S0≤(1 −γ0)−1(c0+B) = (1 + σ0γ0) (c0+B)≤exp(σ0γ0)(c0+B). Assume that Sk≤expσkγk+ k−1 X j=0 (σjγj+δj) k X j=0 cj+B, k = 0,...,n, and let us see that the inequality also holds for k=n+ 1, which will finish the proof. We have Sn+1 ≤(1 −γn+1)−1((1 + δn)Sn+cn+1)≤(1 + σn+1γn+1) (exp(δn)Sn+cn+1) ≤exp (σn+1γn+1)expn X k=0 (σkγk+δk) n X k=0 ck+B+cn+1. □ Now, we consider a fully discrete approximation to the solution of (17). Let us consider time levels 0 = t0< t1< . . . < tN=T, and let us denote ∆tn=tn+1 −tn, n = 0,1...,N −1, and ωn= ∆tn/∆tn−1, n = 1,...,N −1. 6 In the sequel we will assume that for γ∈(0,1) and Γ >1, we have γ≤ωn≤Γ, n = 0,1,...,N −1.(19) We wil also use the following notation, Dun=un−un−1,ˆ un=un+ωnDun, and Dun=Dun+ωn−1 1 + ωn−1un−ˆ un−1. Observe that ˆ unis the extrapolated value at tn+1 of the linear interpolant taking values unand un−1at tnand tn−1, respectively. Notice also that for fixed step size, that is when ∆tn=T/N,n= 0,1,...,N −1, we have ˆ un= 2un−un−1and Dun= (D+1 2D2)un=3 2un−2un−1+1 2un−2. Now, we consider the following implicit-explicit (IMEX) approximation to (17) based on the BDF2 method. Given u0 h,u1 hsolve, for n≥1 and for all (vh, qh)∈Vh×Qh 1 ∆tnDun+1 h,vh+ν(∇un+1 h,∇vh) + b(ˆ un h,ˆ un h,vh)−(pn+1 h,∇·vh) (20) +(∇·un+1 h, qh) + µ(∇·un+1 h,∇·vh) = (fn+1,vh), where, here and in the sequel, the notation fn+1 means f(tn+1). We also consider the semi-implicit method obtained when the third term on the left-hand side of (20) is replaced by b(ˆ un h,un+1 h,vh) (21) that is, 1 ∆tnDun+1 h,vh+ν(∇un+1 h,∇vh) + b(ˆ un h,un+1 h,vh)−(pn+1 h,∇·vh) (22) +(∇·un+1 h, qh) + µ(∇·un+1 h,∇·vh) = (fn+1,vh). Taking vh∈Vdiv hin (20) we get 1 ∆tnDun+1 h,vh+ν(∇un+1 h,∇vh) + b(ˆ un h,ˆ un h,vh) + µ(∇·un+1 h,∇·vh) = (fn+1,vh).(23) To get u1 hwe apply the IMEX Euler method with the same structure as BDF2 (explicit non-linear term) for both the IMEX and the semi-implicit methods. In fact, formulae (20) and (23) are valid for n= 0 if we set ˆ u0=u0,Du1=Du1.(24) We will compare un hwith sn h. For simplicity we will assume u0 h=s0 h. It is easy to prove that for n≥1 and vh∈Vdiv h 1 ∆tnDsn+1 h,vh+ν(∇sn+1 h,∇vh) + b(ˆ sn h,ˆ sn h,vh) + µ(∇·sn+1 h,∇·vh) (25) = (fn+1,vh) + (pn+1 −πn+1 h,∇·vh) + µ(∇·(sn+1 h−un+1),∇·vh) +1 ∆tnDsn+1 h−un+1 t,vh+b(ˆ sn h,ˆ sn h,vh)−b(un+1,un+1,vh). Then, denoting by en h=sn h−un h. 7 and subtracting (23) from (25) we get 1 ∆tnDen+1 h,vh+ν(∇en+1 h,∇vh) + b(ˆ sn h,ˆ sn h,vh)−b(ˆ un h,ˆ un h,vh) +µ(∇·en+1 h,∇·vh) = (τn+1 1+τn+1 3,vh) + (τn+1 2,∇·vh),(26) where τn+1 1=1 ∆tnDsn+1 h−un+1 t, τn+1 2= (pn+1 −πn+1 h) + µ∇·(sn+1 h−un+1),(27) (τn+1 3,vh) = b(ˆ sn h,ˆ sn h,vh)−b(un+1,un+1,vh). Similarly, for the semi-implicit method we get (23) with the third and forth terms replaced by b(ˆ sn h,sn+1 h,vh)−b(ˆ un h,un+1 h,vh), and τ3replaced by τ4defined by (τn+1 4,vh) = b(ˆ sn h,sn+1 h,vh)−b(un+1,un+1,vh).(28) In the following lemma we bound the difference of the nonlinear terms for the IMEX method. Lemma 3 Considering the fixed stepsize case ∆tn= ∆t=T/N, the following bound holds |b(ˆ sn h,ˆ sn h,en+1 h)−b(ˆ un h,ˆ un h,en+1 h)| ≤ ∥∇ˆ sn h∥∞ 2+9 2µ∥ˆ sn h∥2 ∞∥en+1 h∥2 0 +3∥∇ˆ sn h∥∞∥en h∥2 0+3 2∥∇ˆ sn h∥∞∥en−1 h∥2 0+µ 3∥∇·en h∥2 0+µ 6∥∇·en−1 h∥2 0(29) +2∥ˆ un h∥∞cinvh−1∥en+1 h−ˆ en h∥0∥en+1 h∥0. In view of the upper bound in (19), for the variable stepsize case the bound for the difference of the nonlinear terms is |b(ˆ sn h,ˆ sn h,en+1 h)−b(ˆ un h,ˆ un h,en+1 h)| ≤ ∥∇ˆ sn h∥∞ 2+9α 2µ∥ˆ sn h∥2 ∞∥en+1 h∥2 0 +(1 + 2Γ)(1 + Γ) 2∥∇ˆ sn h∥∞∥en h∥2 0+Γ(1 + 2Γ) 2∥∇ˆ sn h∥∞∥en−1 h∥2 0 +µ(1 + 2Γ)(1 + Γ) 18α∥∇·en h∥2 0+µΓ(1 + 2Γ) 18α∥∇·en−1 h∥2 0(30) +2∥ˆ un h∥∞cinvh−1∥en+1 h−ˆ en h∥0∥en+1 h∥0, where αis a positive constant. Proof We start by proving (29). Adding ±b(ˆ un h,ˆ sn h,en+1 h), and noticing that b(ˆ un h,en+1 h,en+1 h) = 0, we have |b(ˆ sn h,ˆ sn h,en+1 h)−b(ˆ un h,ˆ un h,en+1 h)|=|b(ˆ en h,ˆ sn h,en+1 h)−b(ˆ un h,en+1 h−ˆ en h,en+1 h)|.(31) For the first term on the right-hand side above we write |b(ˆ en h,ˆ sn h,en+1 h)|≤∥ˆ en h∥0∥∇ˆ sn h∥∞∥en+1 h∥0+∥ˆ sh∥∞∥∇· ˆ en h∥0∥en+1 h∥0(32) ≤∥∇ˆ sn h∥∞ 2+9 2µ∥ˆ sn h∥2 ∞∥en+1 h∥2 0+∥∇ˆ sn h∥∞ 1 2∥ˆ en h∥2 0+µ 18∥∇· ˆ en h∥2 0. Let us now observe that, ∥ˆ en h∥2 0≤6∥en h∥2 0+ 3∥en−1 h∥2 0, 8 and a similar expression for ∥∇· ˆ en h∥2 0, so that, |b(ˆ en h,ˆ sn h,en+1 h)| ≤ ∥∇ˆ sn h∥∞ 2+9 2µ∥ˆ sn h∥2 ∞∥en+1 h∥2 0+µ 3∥∇·en h∥2 0+µ 6∥∇·en−1 h∥2 0 +3∥∇ˆ sn h∥∞∥en h∥2 0+3 2∥∇ˆ sn h∥∞∥en−1 h∥2 0. Finally, from the definition of the nonlinear term and using inverse inequality (5) it is easy to check that |b(ˆ un h,en+1 h−ˆ en h,en+1 h)| ≤ 2∥ˆ un h∥∞cinvh−1∥en+1 h−ˆ en h∥0∥en+1 h∥0. The proof of (30) can be obtained arguing in the same way with only two differences. The first one is that we include a parameter αto bound the second term on the first line of the right-hand side of (32). The second one is that we apply (19) to bound both ∥ˆ en h∥2 0and ∥∇· ˆ en h∥2 0.□ For the semi-implicit method we have the following result. Lemma 4 The following bound holds in the fixed stepsize case ∆tn= ∆t=T/N, |b(ˆ sn h,sn+1 h,en+1 h)−b(ˆ un h,un+1 h,en+1 h)| ≤ ∥∇sn+1 h∥∞ 2+9 2µ∥sn+1 h∥2 ∞∥en+1 h∥2 0 +3∥∇sn+1 h∥∞∥en h∥2 0+3 2∥∇sn+1 h∥∞∥en−1 h∥2 0+µ 3∥∇·en h∥2 0+µ 6∥∇·en−1 h∥2 0,(33) and, in the variable stepsize case, for Γthe constant in (19) and αany positive constant, the following bound holds: |b(ˆ sn h,sn+1 h,en+1 h)−b(ˆ un h,un+1 h,en+1 h)| ≤ ∥∇sn+1 h∥∞ 2+9α 2µ∥sn+1 h∥2 ∞∥en+1 h∥2 0 +(1 + 2Γ)(1 + Γ) 2∥∇sn+1 h∥∞∥en h∥2 0+Γ(1 + 2Γ) 2∥∇sn+1 h∥∞∥en−1 h∥2 0 +µ(1 + 2Γ)(1 + Γ) 18α∥∇·en h∥2 0+µΓ(1 + 2Γ) 18α∥∇·en−1 h∥2 0.(34) Proof By adding ±b(ˆ un h,sn+1 h,en+1 h) and taking into account that b(ˆ un h,en+1 h,en+1 h) = 0, we have |b(ˆ sn h,sn+1 h,en+1 h)−b(ˆ un h,un+1 h,en+1 h)|=|b(ˆ en h,sn+1 h,en+1 h)|, which can be bounded as the term |b(ˆ en h,ˆ sn h,en+1 h)|in the previous lemma. □ Remark 2 We notice the similarity between the bounds in Lemmas 3 and 4. Indeed, bounds (33) and (34) differ from (29) and (30) in that ˆ sn hin (29) and (30) is replaced by sn+1 h in (33) and (34) and, more importantly, the last term, 2∥ˆ un h∥∞cinvh−1∥en+1 h−ˆ en h∥0∥en+1 h∥0 on the right-hand side of (29) and (30) is not present in (33) and (34). We will see that this term gives rise to a CFL type condition that affects the IMEX (20) method but not the semi-implicit method (22). We now estimate the truncation errors. For τn+1 2, in view of (6) and (13), it easily follows that  τn+1 2 0≤Cµ un+1 k+1 hk+ pn+1 l+1 hl+1.(35) For τn+1 1,τn+1 3and τn+1 4we treat separately the cases n≥1 and n= 0. 9 Remark 4 From Theorems 1 and 2, using triangle inequality and applying (13), we reach max 1≤n≤N∥un h−un∥0=O((∆t)2+hk+hl+1), which means that the rate of convergence is of optimal order 2 in time. Choosing for example Hood-Taylor elements with l=k−1 we get a rate of convergence of order kin space, which matches previously results in [7], as stated in the introduction. 3.2 Variable stepsize In this section we carry out the error analysis for the variable stepsize case. We start with a technical lemma that is needed to prove the main result. Lemma 7 Assuming condition (19) and denoting by Gn=ωn−1 2(1 + ωn−1)∥en h∥2 0+1 4 en h+ (en h−en−1 h)  2 0, n = 1,2,...,N, (70) the following inequality holds, for a constant K0≤75/14, (Den+1 h,en+1 h)≥Gn+1 −Gn+ωn 2(1 + ωn) en+1 h−ˆ en h  2 0(71) −K0 γ(1 + Γ)2(|ωn−1−1|+|ωn−1|)Gn. Proof We first rewrite adequately the term (Den+1 h,en+1 h). Using the identity (v−w,v) = 1 2∥v∥2 0−∥w∥2 0+∥v−w∥2 0),(72) we have that (Den+1 h,en+1 h) =1 2 en+1 h  2 0−∥en h∥2 0+ en+1 h−en h  2 0(73) +ωn 2(1 + ωn) en+1 h  2 0−∥ˆ en h∥2 0+ en+1 h−ˆ en h  2 0. For the first term on the right-hand side above, using the identity (72), we write 1 2 en+1 h  2 0−∥en h∥2 0+ en+1 h−en h  2 0=1 4 en+1 h  2 0−∥en h∥2 0+ en+1 h−en h  2 0 +1 2(en+1 h,en+1 h−en h) =1 4 en+1 h+ (en+1 h−en h)  2 0−1 4∥en h∥2 0. Then (Den+1 h,en+1 h) = Gn+1 −1 4∥en h∥2 0+ωn 2(1 + ωn)∥ˆ en h∥2 0+ωn 2(1 + ωn) en+1 h−ˆ en h  2 0 =Gn+1 −Gn+Rn+ωn 2(1 + ωn) en+1 h−ˆ en h  2 0,(74) where Rn=Gn−1 4∥en h∥2 0+ωn 2(1 + ωn)∥ˆ en h∥2 0. Noticing that  en h+ωn(en h−en−1 h)  2 0= (1 + ωn)2∥en h∥2 0−2ωn(1 + ωn)(en h,en−1 h) + ω2 n en−1 h  2 0, 16 a straightforward computation shows Rn=c1,n ∥en h∥2 0−c2,n(en h,en−1 h) + c3,n  en−1 h  2 0,(75) where c1,n =ωn−1 2(1 + ωn−1)−1 4+ 1 −ωn(1 + ωn) 2,(76) c2,n = 1 −ω2 n,(77) c3,n =1 4−ω3 n 2(1 + ωn).(78) We now express these three coefficients in terms of ωn−1 and ωn−1−1. For c1,n, we have c1,n =ωn−1−1 4(1 + ωn−1)+2−ωn−ω2 n 2=ωn−1−1 4(1 + ωn−1)+2 + ωn 2(1 −ωn). For c3,n, we have c3,n =1 + ωn−2ω3 n 4(1 + ωn)=1 + 2ωn+ 2ω2 n 4(1 + ωn)(1 −ωn). From (19) it is easy to check that |c1,n| ≤ 1 4|ωn−1−1|+ (1 + Γ/2) |ωn−1|, |c2,n| ≤ (1 + Γ) |ωn−1|, |c3,n| ≤ 1 4+1 + Γ 2|ωn−1|, so that max 1≤j≤3|cj,n| ≤ 1 4|ωn−1−1|+ (1 + Γ) |ωn−1|.(79) Furthermore, since Gn=1 + ωn−1 2(1 + ωn−1)∥en h∥2 0−(en h,en−1 h) + 1 4 en−1 h  2 0, and the smallest eigenvalue of matrix 1 + x−1/2 −1/2 1/4(80) is λ=(5 + 4x)−p(5 + 4x)2−16x/8, which can be seen to be λ≥0.14x, for 0 ≤x≤ 1/2 (see Fig. 1). Then, if follows that Gn≥0.07 (1 + Γ)ωn−1 en−1 h  2 0+∥en h∥2 0.(81) Thus, if (19) holds, from (79) and (81) it follows that |Rn| ≤ K0 1 + Γ 4γ(|ωn−1−1|+ 4(1 + Γ) |ωn−1|)Gn ≤K0 γ(1 + Γ)2(|ωn−1−1|+|ωn−1|)Gn,(82) with K0≤(3/8)100/7 = 75/14. Inserting (82) into (74) we finally obtain (71). □ 17 0 0.1 0.2 0.3 0.4 0.5 0 0.02 0.04 0.06 0.08 Figure 1: Graphs of smallest eigenvalue of matrix (80) and y= 0.14x. In the sequel, we define ∆t= max 1≤n≤N∆tn. Also, for α∈[1,10], we consider the polynomial x(1 + 2x) αx2+(1 + 2x)(1 + x) αx−17.(83) One can check that it has two complex roots and two real ones, one negative and one positive. Let Γ∗(α) be this positive real root. It is possible to check that Γ∗(α) is an increasing function of αwhose values for α= 1,5,10 are larger than, 1.25, 2.12 and 2.61, respectively. Fixed α∈[1,10], we restrict ourselves to meshes satisfying Γ≤Γ∗(α).(84) We also assume N−1 X n=2 (|ωn−1−1|+|ωn−1|)≤Λ,(85) for some positive Λ. We now state and prove the convergence result for the IMEX method (20). Theorem 3 Fix κ > 0, and σ > 1, and let ∆tnand hsatisfy, the following CFL-type condition 3c2 inv∥u∥2 L∞(L∞)(1 + 2Γ)2(1 + Γ) γ ∆tn h2≤κ T, n ≥1,(86) and 1−∆tn (1 + Γ)Ln h 0.07ωn >1 σ, n ≥1,(87) where Ln his defined in (95). Let Gnbe as defined in (70) for the error en hof method (20). Then, the following bound holds with mkdefined in (99) Gn+1 +ν n+1 X k=2 ∆tk−1∥∇ek h∥2 0+µ n+1 X k=2 ∆tkmk∥∇·ek h∥2 0≤(88) Cn(C2 0+TC2 3)(∆t)4+TC2 1h2k+C2 2h2l+2, where Cn= exp σn+1∆tn (1 + Γ)Ln h 0.07ωn + n X k=2 σk∆tk−1 (1 + Γ)Lk−1 h 0.07ωk−1 +fk!!,(89) and the constants fkand σkare defined in (97),(102), respectively, and the constants C2 i are defined in (59) and (63). 18 Proof Taking vh=en+1 hin (26) and applying (71) from Lemma 7 we reach Gn+1 −Gn+ωn 2(1 + ωn) en+1 h−ˆ en h  2 0+ν∆tn∥∇en+1 h∥2 0+µ∆tn∥∇·en+1 h∥2 0≤(90) ∆tn|b(ˆ sn h,ˆ sn h,en+1 h)−b(ˆ un h,ˆ un h,en+1 h)|+ ∆tn(τn+1 1+τn+1 3,en+1 h) + ∆tn(τn+1 2,∇·en+1 h) +K0 γ(1 + Γ)2(|ωn−1−1|+|ωn−1|)Gn. As before, we now apply Lemma 3 to estimate the nonlinear term and, further, we bound 2∥ˆ un h∥∞cinvh−1∥en+1 h−ˆen h∥0∥en+1 h∥0≤2c2 inv∥ˆ un h∥2 ∞ (1 + ωn) ωn ∆tn h2∥en+1 h∥2 0(91) +ωn 2(1 + ωn)∥en+1 h−ˆ en h∥2 0 ∆tn . Instead of (86) we will now assume 2c2 inv∥ˆ un h∥2 ∞ (1 + ωn) ωn ∆tn h2≤κ T.(92) As before, we will see at the end of the proof that (92) holds if (86) holds and his sufficiently small. Thus, from (30) in Lemma 3, (91) and assuming (92) we obtain ∆tn|b(ˆ sn h,ˆ sn h,en+1 h)−b(ˆ un h,ˆ un h,en+1 h)| ≤ ∆tn∥∇ˆ sn h∥∞ 2+9α 2µ∥ˆ sn h∥2 ∞+κ T∥en+1 h∥2 0 +∆tn (1 + 2Γ)(1 + Γ) 2∥∇ˆ sn h∥∞∥en h∥2 0+ ∆tn Γ(1 + 2Γ) 2∥∇ˆ sn h∥∞∥en−1 h∥2 0 +∆tn µ(1 + 2Γ)(1 + Γ) 18α∥∇·en h∥2 0+ ∆tn µΓ(1 + 2Γ) 18α∥∇·en−1 h∥2 0(93) +ωn 2(1 + ωn)∥en+1 h−ˆ en h∥2 0. Inserting (93) into (90) and, arguing as before with the terms involving τn+1 1+τn+1 3and τn+1 2, we get Gn+1 −Gn+ν∆tn∥∇en+1 h∥2 0+17 18µ∆tn∥∇·en+1 h∥2 0≤ ∆tn∥∇ˆ sn h∥∞ 2+9α 2µ∥ˆ sn h∥2 ∞+κ+ 1 T∥en+1 h∥2 0+(1 + 2Γ)(1 + Γ) 2∆tn∥∇ˆ sn h∥∞∥en h∥2 0 +Γ(1 + 2Γ) 2∆tn∥∇ˆ sn h∥∞∥en−1 h∥2 0+ ∆tn µ(1 + 2Γ)(1 + Γ) 18α∥∇·en h∥2 0 +∆tn µΓ(1 + 2Γ) 18α∥∇·en−1 h∥2 0+ ∆tn T 2∥τn+1 1∥2 0 +∆tn T 2∥τn+1 3∥2 0+9∆tn 2µ∥τn+1 2∥2 0+K0 γ(1 + Γ)2(|ωn−1−1|+|ωn−1|)Gn.(94) Let us denote by Ln h=∥∇ˆ sn h∥∞ 2+9α 2µ∥ˆ sn h∥2 ∞+κ+ 1 T,(95) so that, using (81), the first term on the right-hand side of (94) can be written as ∆tnLn h∥en+1 h∥2 0≤∆tn (1 + Γ)Ln h 0.07ωn Gn+1. Similarly, for the terms involving ∥en h∥2 0and ∥en−1 h∥2 0on the right-hand side in (94), we have (1 + 2Γ)(1 + Γ) 2∆tn∥∇ˆ sn h∥∞∥en h∥2 0+Γ(1 + 2Γ) 2∆tn∥∇ˆ sn h∥∞∥en−1 h∥2 0≤ (1 + Γ)(1 + 2Γ)2 0.14ωn−1 ∆tn∥∇ˆ sn h∥∞Gn. 19 Thus, going back to (94) we reach Gn+1 −Gn+ν∆tn∥∇en+1 h∥2 0+17 18µ∆tn∥∇·en+1 h∥2 0≤∆tn (1 + Γ)Ln h 0.07ωn Gn+1 +(1 + Γ)(1 + 2Γ)2 0.14ωn−1 ∆tn∥∇ˆ sn h∥∞+K0 γ(1 + Γ)2(|ωn−1−1|+|ωn−1|)Gn +∆tn µ(1 + 2Γ)(1 + Γ) 18α∥∇·en h∥2 0+ ∆tn µΓ(1 + 2Γ) 18α∥∇·en−1 h∥2 0(96) +∆tn 2T(∥τn+1 1∥2 0+∥τn+1 3∥2 0) + 9 µ∥τn+1 2∥2 0. Denoting by fn=(1 + Γ)(1 + 2Γ)2 0.14ωn−1 ∆tn∥∇ˆ sn h∥∞+K0 γ(1 + Γ)2(|ωn−1−1|+|ωn−1|),(97) and adding terms we get Gn+1 +ν n+1 X k=2 ∆tk−1∥∇ek h∥2 0 +µ n+1 X k=2 ∆tk17 18 1 ωk−(1 + 2Γ)(1 + Γ) 18α−Γ(1 + 2Γ)ωk+1 18α∥∇·ek h∥2 0≤ G1+µ(1 + 2Γ)(1 + Γ)∆t1 18α+Γ(1 + 2Γ)∆t2 18α∥∇·e1 h∥2 0+µΓ(1 + 2Γ)∆t1 18α∥∇·e0 h∥2 0 + n+1 X k=2 ∆tk−1 (1 + Γ)Lk−1 h 0.07ωk−1 Gk+ n+1 X k=2 fk−1Gk−1 +1 2 n+1 X k=2 ∆tk−1T(∥τk 1∥2 0+∥τk 3∥2 0) + 9 µ∥τk 2∥2 0.(98) As mentioned before, for any fixed value of α, the number Γ∗(α) in (84) is chosen to be the the positive real root of (83). Since ωk≤Γ∗,k= 1,...,N −1, it follows that mk:= 17 18 1 ωk−(1 + 2Γ)(1 + Γ) 18α−Γ(1 + 2Γ)ωk+1 18α≥0.(99) Let us also denote by Ev 0= (1 + f1)G1+µ(1 + 2Γ)(1 + Γ)∆t1 18α+Γ(1 + 2Γ)∆t2 18α∥∇·e1 h∥2 0 +µ(1 + 2Γ)(1 + Γ)∆t1 18α∥∇·e0 h∥2 0.(100) Taking into account (100), (98) can be written as Gn+1 +ν n+1 X k=2 ∆tk−1∥∇ek h∥2 0+µ n+1 X k=2 ∆tkmk∥∇·ek h∥2 0≤(101) ∆tn (1 + Γ)Ln h 0.07ωn Gn+1 + n X k=2 ∆tk−1 (1 + Γ)Lk−1 h 0.07ωk−1 +fk!Gk +1 2 n+1 X k=2 ∆tk−1T(∥τk 1∥2 0+∥τk 3∥2 0) + 9 µ∥τk 2∥2 0+Ev 0. 20 Now assumption (87) implies that σn+1 defined by σn+1 =1−∆tn (1 + Γ)Ln h 0.07ωn−1 ,(102) is positive. Thus, applying Lemma 2 we have Gn+1 +ν n+1 X k=2 ∆tk−1∥∇ek h∥2 0+µ n+1 X k=2 ∆tkmk∥∇·ek h∥2 0≤(103) CnEv 0+1 2 n+1 X k=2 ∆tk−1T(∥τk 1∥2 0+∥τk 3∥2 0) + 9 µ∥τk 2∥2 0, where Cnis defined in (89). To conclude we bound the initial and truncation errors. For the initial errors, we argue as in the fixed stepsize case. Since we assume u0 h=s0 hso that e0 h= 0. Consequently, G1=ω0 2(1 + ω0) e1 h  2 0+1 4 2e1 h  2 0≤3 2 e1 h  2 0, and, in view of (100), Ev 0≤(1 + f1)3 2 e1 h  2 0+ω0(1 + 2Γ)(1 + Γ) 18α+Γ(1 + 2Γ)ω1 18αµ∆t0 ∇·e1 h  2 0.(104) The first step is carried out with the implicit-explicit Euler method and then, in view of (104), and arguing as we did with (64) we get Ev 0≤C1 + ∆t∥u∥L∞(W1,∞)C2 0(∆t)4+TC2 1h2k+C2 2h2l+2,(105) where Cdepends on Γ, 1/γ,ω0and ω1. For the truncation errors we also argue as in the fixed stepsize case. Then, arguing exactly as in (58) we get 1 2 n+1 X k=2 ∆tk−1T(∥τk 1∥2 0+∥τk 3∥2 0) + 9 µ∥τk 2∥2 0≤TC2 1h2k+C2 2h2l+2 +C2 3(∆t)4,(106) with the constants C2 idefined in (59) and avoiding writing the explicit dependency on ωn. Inserting (105) and (106) into (103) we finally reach (88). To conclude the proof, as in the fixed step size, we now show that taking hsufficiently small, the CFL type condition (86) implies (92). Arguing as in the fixed step-size, it is possible to choose h0,1such that ∥sn h−un∥∞≤ ∥u∥L∞(L∞)/12 for h≤h0,1, which implies ∥sn h∥∞≤13 12 ∥u∥L∞(L∞),(107) and, since ∥ˆ sn h∥ ≤ (1 + 2Γ) ∥sh∥L∞(L∞), by writing ˆ un h= (ˆ un h−ˆ sn h) + ˆ sn h, it also holds ∥ˆ un h∥∞≤cinvh−d/2(1 + Γ) ∥sn h−un h∥0+ Γ  sn−1 h−un−1 h 0+ (1 + 2Γ)13 12 ∥u∥L∞(L∞). (108) If an h0,2>0 exists such that for h≤h0,2both the right-hand side of (62) and the righthand side of (88) multiplied by 0.07(1 + Γ)/γ are smaller than hd∥u∥2 L∞(L∞)/(12cinv)2it follows that ∥sk h−uk h∥0≤hd/2 12cinv ∥u∥L∞(L∞), k = 1,··· , n. 21 Thus, in view of (108), it follows that ∥ˆ un h∥∞≤(1 + 2Γ)7 6∥u∥L∞(L∞), so that (92) holds as a consequence of assumption (86). We now argue why h0,2exists, and this will finish the proof. The exponents of the powers of hon right-hand sides of (62) and (88) are larger than d, so that they will decay with hfaster than hdif the exponential term in (89) remains bounded as h→0. To show that this is the case, we first notice that, in view of the expression of Ln hin (95) and recalling (16) and (107), we have Ln h≤(1 + 2Γ) (1 + C∞) 2∥u∥L∞(W1,∞)+9α 2µ13 122 ∥u∥2 L∞(L∞)!+κ+ 1 T and in view of the expression of fkin (97), we also have fn≤(1 + Γ)(1 + 2Γ)2 0.14γ∆tn(1 + C∞)∥u∥L∞(W1,∞)+K0 γ(1 + Γ)2(|ωn−1−1|+|ωn−1|). Finally, due to assumption (87), σndefined in (102) satisfies σn< σ and noticing that the sum of the increments ∆tnis always bounded by Tand recalling assumption (85), it is clear that the constants Cnin (89) remain bounded as h→0. □ Remark 5 In view of (86) we notice that the ratios ∆tn/h2must remain bounded. As stated in Remark 3, the size of the constant κin the CFL condition (86) reflects in the error bound (88) through an exponential factor. More precisely, considering also the size of the constant σin (87), taking into account σn< σ (with σndefined in (102)) and in view of the expression of Ln hin (95) the influence of the constants κand σaffect the constant Cn in (103) by a factor of size expκσ 1 + Γ 0.07T n X k=1 ∆tk ωk= expκσ 1 + Γ 0.07T n X k=1 ∆tk−1≤expκσ 1 + Γ 0.07 . On the other hand, the value of Γ∗(α) in (84) increases with the value of the free parameter α. However, arguing as above, we can observe that the effect of increasing α(to allow for larger increases in step length) also affects the constant Cnin (89) since the term 9α 2µ∥ˆ sn h∥2 ∞ is present in the definition of Ln hin (95). Theorem 4 Fix σ > 1, and let ∆tnand hsatisfy 1−∆tn (1 + Γ)Ln h 0.07ωn >1 σ, n ≥1, where Ln h=∥∇sn+1 h∥∞ 2+9α 2µ∥sn+1 h∥2 ∞+1 T. Let Gnbe as defined in (70) for the error en hof the semi-implicit method (22). Then, the following bound holds with mkdefined in (99) Gn+1 +ν n+1 X k=2 ∆tk−1∥∇ek h∥2 0+µ n+1 X k=2 ∆tkmk∥∇·ek h∥2 0≤(109) Cn(C2 0+TC2 3)(∆t)4+TC2 1h2k+C2 2h2l+2, where Cn= exp σn+1∆tn (1 + Γ)Ln h 0.07ωn + n X k=2 σk∆tk−1 (1 + Γ)Lk−1 h 0.07ωk−1 +ωkfk!!, 22 and the constants fkand σkare defined in (97),(102), respectively, and the constants C2 i are defined in (59) and (63). Proof The proof follows that of Theorem 3 with the changes that we now comment on. Indeed, it is easy to see that (90) also holds for the semi-implicit method (22) with the first term on the right-hand side of (90) replaced by |b(ˆ sn h,sn+1 h,en+1 h)−b(ˆ un h,un+1 h,en+1 h)| and τ3by τ4. Then, recalling Lemma 4 and Remark 2, we have that (94) holds with κ= 0, ˆ sn hreplaced by sn+1 hand τ3by τ4. The proof is concluded following the steps of the proof of Theorem 3, but without having to prove that the CFL condition (92) holds. □ Remark 6 Analogous comments to those in Remark 4 apply in this case taking into account (88), (109) and the definition of Gnin (70). Remark 7 Concerning condition (85), stability of multistep methods with variable step size under similar assumptions can be found in the literature ([19], [20]). However, condition (85) is stronger than some less restrictive requirements that can be found in more recent literature. For example, in [2], it is only required that PN−2 n=2 [ωn−ωn+2]+to be bounded, where [·]+denotes the positive part, which allows for graded meshes. However, bounds obtained under these milder restrictions require to use as test functions in the error equations (26) a linear combination ˜ en+1 h=αnen+1 h+βnen h, with αnβn= 0, n≥2. Then, for the semi-implicit method, when estimating the nonlinear terms in Lemma 4, instead of term b(ˆ un h,en+1 h,en+1 h), which is zero, we would have b(ˆ un h,en+1 h,˜ en+1 h), which cannot be assumed to be zero. Obtaining error bounds independent of ν−1, which is our aim in this paper, requires estimating this term in terms of ∥ˆ un+1 h∥2 ∞∥en+1 h∥2 0and ∥en h∥2 0 which is impossible unless the inverse inequality (5) is used, which, in turn, would imply a factor 1/h in the exponent when applying Gronwall Lemma, preventing convergence. On the other hand, one could try to find out if using ˜ en+1 has test function in the error equations in the IMEX method may allow for a milder condition than 85. For this case, taking into account that αn= 0, the CLF-type condition (51) would be still required. In practice, as it will be seen in Section 4 below, condition (51) makes difficult to find the IMEX method preferable to semi-implicit from the practical point of view. For this reason, and in order to keep a single analysis for both IMEX and the semi-implicit methods valid for high Reynolds numbers we do not pursue the issue of taking ˜ en+1 has test function in the IMEX method. Remark 8 Condition (85) can be slightly improved by studying the eigenvalues of matrix Bn=c1,n c2,n/2 c2,n/2c3,n , where c1,n,c2,n,c3,n are the coefficients in the expression of Rnin (75). In fact, it is easy to check that det(Bn) = 0 when ωn= 1 or when ωn−1= (1 + 4ωn)/(1 + 4ω2 n), and that it is positive semidefinite only when ωn−1≥(1 + 4ωn)/(1 + 4ω2 n) and ωn≤1, so that the corresponding values of nwhere this happens can be removed from the sum in condition (85). In order not to increase the already extensive notation, we have preferred not to reflect this fact in condition (85). 4 Numerical Experiments We test methods (20) and (22) with variable step size. We now describe the algorithm used, which follows standard procedures in the numerical integration of ordinary differential equations (ODEs). For the implementation of the method, we used a variable formula strategy as described in [14, §III.5] for the fully nonlinear BDF formulae, but with explicit 23 treatment of the convective term as specified in (20). In particular we notice that the local error estimate of the method of order kis given by ESTn=∆tn tn+1 −tn−k  Un,k h  0, where Un,k h= k−1 Y i=0 (tn+1 −tn−i)uh[tn+1,...,tn−k], where uh[tn+1,...,tn−k] is the standard (k+1)-th divided difference based on tn+1,...,tn−k and un+1 h,...,un−k h. This corresponds to the fully nonlinear BDF, but, nevertheless, it was used with the methods (20) and (22). At each time level tn, the estimation ESTnis compared with the quantity TOLn=TOLrmax  un+1 h 0,∥un h∥0+ 0.001, where TOLris a given tolerance. If ESTn>TOLn,un+1 his not considered sufficiently accurate and is rejected. It is then recomputed from tnwith a new step length given by ∆tnew n= 0.9∆tnk+1 rTOLn ESTn .(110) On the contrary, if ESTn<TOLn, then, un+1 his accepted and the algorithm proceeds to compute un+2 hwith a step length ∆tn+1 given by the right-hand side of (110). The algorithm starts with k= 1 and ∆t0= ∆t1=√T OLr/100. At t1the first local error estimation EST1is computed. If EST1>TOL1the computation is restarted again from t0with the step size changed according to (110). From t2onwards, estimates corresponding to the two-step BDF can be estimated. When these estimates are smaller than those corresponding to k= 1, then kis switched to k= 2. In the experiments below, no more than eight or ten steps were performed with k= 1. We now comment on the solution of the linear systems to be solved at each step, where, again, we follow standard procedures in the numerical integration of ODEs. If we denote by ynthe vector containing the coefficients of the velocity un+1 hand the pressure pn+1 hin the corresponding nodal basis of the finite element spaces, then, since (20) and (22) are linear in un+1 hand pn+1 h, to find ynone has to solve a linear system Anyn=bn.(111) This was solved by iterative refinement. Starting with an approximation yn,[0] =ˆ yn, extrapolated from the previous kvalues (k= 1,2), the approximation is updated as yn,[j+1] = yn,[j]+en,[j], where en,[j]is obtained by solving the system Amen,[j]=bn−Anyn,[j], where, we notice, Amhas been factored on time level tm≤tn. The iteration is stopped either when the velocity in two consecutive iterations satisfies   un,[j+1] h−un,[j] h  0≤min 10−8,TOLr 100  un−1 h 0+ 0.001,(112) in which case we set un n=un,[j+1] n, or if five iterations are performed without ever (112) being satisfied, in which case, Anis factored, replaces Amand system (111) is solved. As we will see below, the number of factorizations never exceeded 5% of the total number of steps. The code was programmed in Matlab and for the factorizations in linear systems we used command decomposition with no extra arguments. The set of equations obtained from (20) when vh= 0 was rewriten as −(∇·un+1 h, qh) = 0, so that matrices Anwere symmetric (although indefinite). 24 10-2 10-1 10-6 10-5 10-4 10-3 10-2 10-1 Figure 2: Velocity errors (115) for T= 4 and µ= 0.05. Continuous lines correspond to method (20) and discontinuous lines to method (22). 4.1 Problem with known solution We consider the Navier-Stokes equations in the domain Ω = [0,1]2with T= 4, and the forcing term fchosen so that the solution uand pare given by u(x, y, t) = 6 + 4 cos(4t) 10 8 sin2(πx)2y(1 −y)(1 −2y) −8πsin(2πx)(y(1 −y))2,(113) p(x, y, t) = 6 + 4 cos(4t) 10 sin(πx) cos(πy).(114) We present results with the P2/P1pair of mixed finite-elements on a regular triangulation with SW-NE diagonals. The meshes have N= 6, 12, 24 and 48 subdivisions in each coordinate directions, and the tolerances for the BDF2 on these meshes are, respectively, 10−4, 10−5, 10−6and 10−7. We checked that they are sufficiently small so that the error arising from time discretization is negligible when compared to that arising form spatial discretization. In the experiments in this section, the value of the grad-div parameter is µ= 0.05. In Fig. 2, for the final time T= 4, we show the L2errors in velocity ∥uN h−Ih(uN)∥0,(115) where Ihis the standard Lagrange interpolant. In Fig. 2, as in the rest of the pictures in this section, results of IMEX method (20) are joined by continuous lines, while those of the semiimplict method (22) by discontinuous lines. It can be seen that continuous and discontinuous lines are super imposed, reflecting the fact that the errors shown in Fig. 2 correspond to spatial discretization and are independent of the method chosen for time integration. It can be seen that, as hdecreases, the errors also decrease, but with a higher rate for the two largest values of the viscosity, ν= 10−2and ν= 10−4. The results corresponding to ν= 10−6,ν= 10−8and ν= 10−10 coincide and are on top of one another. We also show the slopes of a least squares fit to the results corresponding to ν= 10−2(in blue) and ν= 10−10, confirming the analysis that the rate of decay is O(h2) independent of νin the convection dominated regime, i.e., for νsmall enough. For ν= 10−2, and the last three results of ν= 10−4, the slope close to 4 reflects that the numerical approximation uhis supraconvergent, in the sense that the error (115) decays faster than ∥un−Ih(un)∥0, which is O(h3). 25 •We prove that using an explicit form for the nonlinear term a CFL type condition is required. This condition is stronger than the usual one in which (∆t)h−1has to be bounded (we required (∆t)h−2to be bounded). Our numerical experiments confirm experimentally that the CFL-type restriction cannot be weakened. Using a semi-implicit form of the nonlinear term, however, the CFL-type condition is neither required in the analysis nor in practice. •Methods for which robust estimates can be derived enable stable flow simulations for small viscosity coefficients on comparatively coarse grids. By robust error estimates we mean estimates with error constants independent of the Reynolds number or of ν−1. Adding grad-div stabilization to the standard mixed finite element formulation we are able to prove robust error bounds for the two methods we analyzed. •Adaptive time stepping is an important tool in Computational Fluid Dynamics. In the present paper, we include the error analysis of a second order backward differentiation formulae (BDF2) scheme with variable time step both for an explicit and semi-implicit form of the nonlinear term. Concerning BDF2 method with variable time step this is the first time where robust error bounds are obtained for the Navier-Stokes equations. •The error bounds in the present work for variable stepsizes are obtained under the condition (86). It remains as an open problem if robust error bounds can be obtained under milder requirements. References [1] R. A. Adams and J. J. F. Fournier. Sobolev spaces, volume 140 of Pure and Applied Mathematics (Amsterdam). Elsevier/Academic Press, Amsterdam, second edition, 2003. [2] J. Becker. A second order backward difference method with variable steps for a parabolic problem. BIT, 38(4):644–662, 1998. [3] H. Chen. Pointwise error estimates for finite element solutions of the Stokes problem. SIAM J. Numer. Anal., 44(1):1–28, 2006. [4] W. Chen, X. Wang, Y. Yan, and Z. Zhang. A second order BDF numerical scheme with variable steps for the Cahn-Hilliard equation. SIAM J. Numer. Anal., 57(1):495–525, 2019. [5] P. G. Ciarlet. The finite element method for elliptic problems. North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978. Studies in Mathematics and its Applications, Vol. 4. [6] J. de Frutos, B. Garc´ıa-Archilla, V. John, and J. Novo. Grad-div stabilization for the evolutionary Oseen problem with inf-sup stable finite elements. J. Sci. Comput., 66(3):991–1024, 2016. [7] J. de Frutos, B. Garc´ıa-Archilla, V. John, and J. Novo. Analysis of the grad-div stabilization for the time-dependent Navier-Stokes equations with inf-sup stable finite elements. Adv. Comput. Math., 44(1):195–225, 2018. [8] J. de Frutos, B. Garc´ıa-Archilla, and J. Novo. Fully discrete approximations to the time-dependent Navier-Stokes equations with a projection method in time and graddiv stabilization. J. Sci. Comput., 80(2):1330–1368, 2019. [9] V. DeCaria and M. Schneier. An embedded variable step IMEX scheme for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., 376:113661, 26, 2021. [10] E. Emmrich. Stability and error of the variable two-step BDF for semilinear parabolic problems. J. Appl. Math. Comput., 19(1-2):33–55, 2005. 32 [11] E. Emmrich. Convergence of the variable two-step BDF time discretisation of nonlinear evolution problems governed by a monotone potential operator. BIT, 49(2):297–323, 2009. [12] B. Garc´ıa-Archilla, V. John, and J. Novo. On the convergence order of the finite element error in the kinetic energy for high Reynolds number incompressible flows. Comput. Methods Appl. Mech. Engrg., 385:Paper No. 114032, 54, 2021. [13] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations, volume 5 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986. Theory and algorithms. [14] E. Hairer and G. Wanner. Solving ordinary differential equations. II, volume 14 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1991. Stiff and differential-algebraic problems. [15] J. G. Heywood and R. Rannacher. Finite-element approximation of the nonstationary Navier-Stokes problem. IV. Error analysis for second-order time discretization. SIAM J. Numer. Anal., 27(2):353–384, 1990. [16] V. John. Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder. International Journal for Numerical Methods in Fluids, 44(7):777– 788, 2004. [17] V. John. Finite element methods for incompressible flow problems, volume 51 of Springer Series in Computational Mathematics. Springer, Cham, 2016. [18] V. John and J. Rang. Adaptive time step control for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., 199(9-12):514–524, 2010. [19] M.-N. Le Roux. Variable step size multistep methods for parabolic problems. SIAM J. Numer. Anal., 19(4):725–741, 1982. [20] C. Palencia and B. Garc´ıa-Archilla. Stability of linear multistep methods for sectorial operators in Banach spaces. Appl. Numer. Math., 12(6):503–520, 1993. [21] M. Sch¨afer, S. Turek, F. Durst, E. Krause, and R. Rannacher. Benchmark Computations of Laminar Flow Around a Cylinder, pages 547–566. Vieweg+Teubner Verlag, Wiesbaden, 1996. [22] W. Wang, Y. Chen, and H. Fang. On the variable two-step imex bdf method for parabolic integro-differential equations with nonsmooth initial data arising in finance. SIAM Journal on Numerical Analysis, 57(3):1289–1317, 2019. [23] W. Wang, M. Mao, and Z. Wang. Stability and error estimates for the variable stepsize BDF2 method for linear and semilinear parabolic equations. Adv. Comput. Math., 47(1):Paper No. 8, 28, 2021. [24] W. Wang, Z. Wang, and M. Mao. Linearly implicit variable step-size BDF schemes with Fourier pseudospectral approximation for incompressible Navier-Stokes equations. Appl. Numer. Math., 172:393–412, 2022. 33