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

Bosco García-Archilla Departamento de Matemática Aplicada II, Universidad de
Sevilla, Sevilla, Spain. Research is supported by Spanish MCINYU under grants
PGC2018-096265-B-I00 and PID2019-104141GB-I00<EMAIL_ADDRESS>Julia Novo
Departamento de Matemáticas, Universidad Autónoma de Madrid, Spain. Research
is supported by Spanish MINECO under grants PID2019-104141GB-I00 and VA169P20
<EMAIL_ADDRESS>
###### 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 allow us to prove
error bounds in which the constants are independent of inverse powers of the
viscosity. Error bounds of order $r$ in space are obtained for the $L^{2}$
error of the velocity using piecewise polynomials of degree $r$ to 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;
Implicit-explicit 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, 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 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 $\Omega\subset{\mathbb{R}}^{d}$, $d\in\\{2,3\\}$, be a bounded domain with
polyhedral and Lipschitz boundary $\partial\Omega$. The incompressible
Navier–Stokes equations model the conservation of linear momentum and the
conservation of mass (continuity equation) by

$\displaystyle\partial_{t}\boldsymbol{u}-\nu\Delta\boldsymbol{u}+(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}+\nabla
p$ $\displaystyle=\boldsymbol{f}$ $\displaystyle\text{in }\
(0,T)\times\Omega,$ $\displaystyle\nabla\cdot\boldsymbol{u}$ $\displaystyle=0$
$\displaystyle\text{in }\ (0,T)\times\Omega,$ (1)
$\displaystyle\boldsymbol{u}(0,\cdot)$
$\displaystyle=\boldsymbol{u}_{0}(\cdot)$ $\displaystyle\text{in }\ \Omega,$

where $\boldsymbol{u}$ is the velocity field, $p$ the kinematic pressure,
$\nu>0$ the kinematic viscosity coefficient, $\boldsymbol{u}_{0}$ a given
initial velocity, and $\boldsymbol{f}$ represents the external body
accelerations acting on the fluid. The Navier–Stokes equations (1) are
equipped with homogeneous Dirichlet boundary conditions
$\boldsymbol{u}=\boldsymbol{0}$ on $\partial\Omega$.

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 $r$
are obtained for the $L^{2}$ error 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 $\nu$. In particular, a CFL-type condition is needed for
the time step depending on $\nu^{-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 on
$\nu^{-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
$(\Delta t)h^{-1}$ has to be bounded ($h$ being the mesh size) because we
required $(\Delta t)h^{-2}$ to 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 $(\Delta t)h^{-1}\nu^{-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 on inverse powers of the viscosity.

Concerning the BDF2 method with variable time step, we did not found in the
literature any reference with the error analysis for the Navier-Stokes
equations. In [2] a second order backward difference method with variable
steps is analyzed for a linear parabolic problem. The error analysis for
semilinear parabolic problems can be found in [10]. 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].

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 $W^{s,p}(D)$ the Sobolev space of real
valued functions defined on the domain $D\subset\mathbb{R}^{d}$ with
distributional derivatives of order up to $s$ in $L^{p}(D)$, endowed with the
usual norm which is denoted by $\|\cdot\|_{W^{s,p}(D)}$. If $s$ is not a
positive integer, $W^{s,p}(D)$ is defined by interpolation [1]. If $s=0$ we
understand that $W^{0,p}(D)=L^{p}(D)$. As it is standard, $W^{s,p}(D)^{d}$
will be endowed with the product norm that, if no confusion can arise, will be
denoted again by $\|\cdot\|_{W^{s,p}(D)}$. We will distinguish the case $p=2$
using $H^{s}(D)$ to denote the space $W^{s,2}(D)$. We will make use of the
space $H_{0}^{1}(D)$, the closure in $H^{1}(D)$ of the set of infinitely
differentiable functions with compact support in $D$. For simplicity, we use
$\|\cdot\|_{s}$ (resp. $|\cdot|_{s}$) to denote the norm (resp. seminorm) both
in $H^{s}(\Omega)$ or $H^{s}(\Omega)^{d}$. The exact meaning will be clear by
the context. The inner product of $L^{2}(\Omega)$ or $L^{2}(\Omega)^{d}$ will
be denoted by $(\cdot,\cdot)$ and the corresponding norm by $\|\cdot\|_{0}$.
The norm of the space of essentially bounded functions $L^{\infty}(\Omega)$
will be denoted by $\|\cdot\|_{\infty}$. For vector valued function we will
use the same conventions as before. We represent by $\|\cdot\|_{-1}$ the norm
of the dual space of $H^{1}_{0}(\Omega)$ which is denoted by $H^{-1}(\Omega)$.
As usual, we always identify $L^{2}(\Omega)$ with its dual so we have
$H^{1}_{0}(\Omega)\subset L^{2}(\Omega)\subset H^{-1}(\Omega)$ with compact
injection.

Using the function spaces $V=H_{0}^{1}(\Omega)^{d}$, and

$Q=L_{0}^{2}(\Omega)=\left\\{q\in L^{2}(\Omega):(q,1)=0\right\\},$

the weak formulation of problem (1) is: Find
$(\boldsymbol{u},p):(0,T]\rightarrow V\times Q$ such that for all
$(\boldsymbol{v},q)\in V\times Q$,

$(\partial_{t}\boldsymbol{u},\boldsymbol{v})+\nu(\nabla\boldsymbol{u},\nabla\boldsymbol{v})+((\boldsymbol{u}\cdot\nabla)\boldsymbol{u},\boldsymbol{v})-(\nabla\cdot\boldsymbol{v},p)+(\nabla\cdot\boldsymbol{u},q)=(\boldsymbol{f},\boldsymbol{v}).$
(2)

The Hilbert space

$H^{\rm div}=\\{\boldsymbol{u}\in L^{2}(\Omega)^{d}\ \mid\
\nabla\cdot\boldsymbol{u}=0,\,\boldsymbol{u}\cdot\mathbf{n}|_{\partial\Omega}=0\\}$

will be endowed with the inner product of $L^{2}(\Omega)^{d}$ and the space

$V^{\rm div}=\\{\boldsymbol{u}\in V\ \mid\ \nabla\cdot\boldsymbol{u}=0\\}$

with the inner product of $V$.

The following Sobolev’s embedding [1] will be used in the analysis: For $1\leq
p<d/s$ let $q$ be such that $\frac{1}{q}=\frac{1}{p}-\frac{s}{d}$. There
exists a positive constant $C$, independent of $s$, such that

$\|v\|_{L^{q^{\prime}}(\Omega)}\leq
C\|v\|_{W^{s,p}(\Omega)},\qquad\frac{1}{q^{\prime}}\geq\frac{1}{q},\quad v\in
W^{s,p}(\Omega).$ (3)

If $p>d/s$ the above relation is valid for $q^{\prime}=\infty$. A similar
embedding inequality holds trivially for vector valued functions.

Let $V_{h}\subset V$ and $Q_{h}\subset Q$ be two families of finite element
spaces composed of piecewise polynomials of degrees at most $k$ and $l$,
respectively, that correspond to a family of partitions $\mathcal{T}_{h}$ of
$\Omega$ into mesh cells with maximal diameter $h$. In the following, we
consider pairs of finite element spaces that satisfy the inf-sup condition

$\inf_{q_{h}\in Q_{h}}\sup_{\boldsymbol{v}_{h}\in
V_{h}}\frac{(\nabla\cdot\boldsymbol{v}_{h},q_{h})}{\|\nabla\boldsymbol{v}_{h}\|_{0}\|q_{h}\|_{0}}\geq\beta_{0},$
(4)

with $\beta_{0}$ a 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 $v_{h}\in V_{h}$, see e.g., [5, Theorem
3.2.6],

$\|\boldsymbol{v}_{h}\|_{W^{m,p}(K)}\leq
C_{\mathrm{inv}}h_{K}^{n-m-d\left(\frac{1}{q}-\frac{1}{p}\right)}\|\boldsymbol{v}_{h}\|_{W^{n,q}(K)},$
(5)

where $0\leq n\leq m\leq 1$, $1\leq q\leq p\leq\infty$ and $h_{K}$ is the size
(diameter) of the mesh cell $K\in\mathcal{T}_{h}$.

The space of discrete divergence-free functions is denoted by

$V_{h}^{\rm div}=\left\\{\boldsymbol{v}_{h}\in V_{h}\ \mid\
(\nabla\cdot\boldsymbol{v}_{h},q_{h})=0\quad\forall q_{h}\in Q_{h}\right\\}.$

Denoting by $\pi_{h}$ the $L^{2}(\Omega)$ projection of the pressure $p$ onto
$Q_{h}$, we have that for $0\leq m\leq 1$

$\|p-\pi_{h}\|_{m}\leq Ch^{l+1-m}\|p\|_{l+1}\quad\forall p\in
H^{l+1}(\Omega).$ (6)

In the error analysis, the Poincaré–Friedrichs inequality

$\|\boldsymbol{v}\|_{0}\leq
C\left|\Omega\right|^{1/d}\|\nabla\boldsymbol{v}\|_{0},\quad\forall\boldsymbol{v}\in
H_{0}^{1}(\Omega)^{d},$ (7)

will be used.

In the analysis, the Stokes problem

$\displaystyle-\nu\Delta\boldsymbol{u}+\nabla p$
$\displaystyle={\boldsymbol{g}}\quad\mbox{in }\Omega,$
$\displaystyle\boldsymbol{u}$ $\displaystyle=\boldsymbol{0}\quad\mbox{on
}\partial\Omega,$ (8) $\displaystyle\nabla\cdot\boldsymbol{u}$
$\displaystyle=0\quad\mbox{in }\Omega,$

will be considered. Let us denote by $(\boldsymbol{u}_{h},p_{h})\in
V_{h}\times Q_{h}$ the mixed finite element approximation to (2), given by

$\begin{split}\nu(\nabla\boldsymbol{u}_{h},\nabla\boldsymbol{v}_{h})-(\nabla\cdot\boldsymbol{v}_{h},p_{h})&=({\boldsymbol{g}},\boldsymbol{v}_{h})\quad\forall\boldsymbol{v}_{h}\in
V_{h}\\\ (\nabla\cdot\boldsymbol{u}_{h},q_{h})&=0\qquad\quad\,\,\,\forall
q_{h}\in Q_{h}.\end{split}$ (9)

Following [13], one gets the estimates

$\displaystyle
h^{-1}\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0}+\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{1}$
$\displaystyle\leq C\left(\inf_{\boldsymbol{v}_{h}\in
V_{h}}\|\boldsymbol{u}-\boldsymbol{v}_{h}\|_{1}+\nu^{-1}\inf_{q_{h}\in
Q_{h}}\|p-q_{h}\|_{0}\right),$ (10) $\displaystyle\|p-p_{h}\|_{0}$
$\displaystyle\leq C\left(\nu\inf_{\boldsymbol{v}_{h}\in
V_{h}}\|\boldsymbol{u}-\boldsymbol{v}_{h}\|_{1}+\inf_{q_{h}\in
Q_{h}}\|p-q_{h}\|_{0}\right).$ (11)

It can be observed that the error bounds for the velocity depend on negative
powers of $\nu$.

For the analysis, it will be advantageous to use a projection of
$(\boldsymbol{u},p)$ into $V_{h}\times Q_{h}$ with uniform in $\nu$, optimal,
bounds for the velocity. In [6] a projection with this property was
introduced. Let $(\boldsymbol{u},p)$ be the solution of the Navier–Stokes
equations (1) with $\boldsymbol{u}\in V\cap H^{k+1}(\Omega)^{d}$, $p\in Q\cap
H^{k}(\Omega)$, $k\geq 1$, and observe that $(\boldsymbol{u},0)$ is the
solution of the Stokes problem (2) with right-hand side

${\boldsymbol{g}}={\boldsymbol{f}}-\partial_{t}\boldsymbol{u}-(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}-\nabla
p.$ (12)

Denoting the corresponding Galerkin approximation in $V_{h}\times Q_{h}$ by
$({\boldsymbol{s}_{h}},l_{h})$, one obtains from (10)-(11)

$\displaystyle\|\boldsymbol{u}-{\boldsymbol{s}}_{h}\|_{0}+h\|\boldsymbol{u}-{\boldsymbol{s}}_{h}\|_{1}$
$\displaystyle\leq Ch^{k+1}\|\boldsymbol{u}\|_{k+1},$ (13)
$\displaystyle\|l_{h}\|_{0}$ $\displaystyle\leq C\nu
h^{k}\|\boldsymbol{u}\|_{k+1},$ (14)

where the constant $C$ does not depend on $\nu$.

Remark 1 Assuming the necessary smoothness in time and considering (2) with

${\boldsymbol{g}}=\partial_{t}\left({\boldsymbol{f}}-\partial_{t}\boldsymbol{u}-(\boldsymbol{u}\cdot\nabla)\boldsymbol{u}-\nabla
p\right),$

one can derive an error bound of the form (13) also for
$\partial_{t}(\boldsymbol{u}-{\boldsymbol{s}}_{h})$. One can proceed similarly
for higher order derivatives in time.

Following [3], one can also obtain the following bounds for
${\boldsymbol{s}}_{h}$

$\displaystyle\|\boldsymbol{u}-{\boldsymbol{s}}_{h}\|_{\infty}$
$\displaystyle\leq
C_{\infty}h\left(\log(\left|\Omega\right|^{1/d}/h)\right)^{\bar{r}}\|\nabla\boldsymbol{u}\|_{\infty},$
(15) $\displaystyle\|\nabla(\boldsymbol{u}-{\boldsymbol{s}}_{h})\|_{\infty}$
$\displaystyle\leq C_{\infty}\|\nabla\boldsymbol{u}\|_{\infty},$ (16)

where $C_{\infty}$ does not depend on $\nu$ and $\bar{r}=0$ for linear
elements and $\bar{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 $(\boldsymbol{u}_{h},p_{h}):(0,T]\rightarrow V_{h}\times
Q_{h}$ such that for all $(\boldsymbol{v}_{h},q_{h})\in V_{h}\times Q_{h}$ one
has

$\begin{split}(\partial_{t}\boldsymbol{u}_{h},\boldsymbol{v}_{h})+\nu(\nabla\boldsymbol{u}_{h},\nabla\boldsymbol{v}_{h})&+b(\boldsymbol{u}_{h},\boldsymbol{u}_{h},\boldsymbol{v}_{h})-(p_{h},\nabla\cdot\boldsymbol{v}_{h})\\\
&+(\nabla\cdot\boldsymbol{u}_{h},q_{h})+\mu(\nabla\cdot\boldsymbol{u}_{h},\nabla\cdot\boldsymbol{v}_{h})=(\boldsymbol{f},\boldsymbol{v}_{h}),\end{split}$
(17)

with $\boldsymbol{u}_{h}(0)$ given. Here, and in the rest of the paper,

$b(\boldsymbol{u},\boldsymbol{v},\boldsymbol{w})=(B(\boldsymbol{u},\boldsymbol{v}),\boldsymbol{w})\quad\forall\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}\in
H_{0}^{1}(\Omega)^{d},$

where,

$B(\boldsymbol{u},\boldsymbol{v})=(\boldsymbol{u}\cdot\nabla)\boldsymbol{v}+\frac{1}{2}(\nabla\cdot\boldsymbol{u})\boldsymbol{v}\quad\forall\boldsymbol{u},\boldsymbol{v}\in
H_{0}^{1}(\Omega)^{d}$

Notice the well-known property

$b(\boldsymbol{u},\boldsymbol{v},\boldsymbol{w})=-b(\boldsymbol{u},\boldsymbol{w},\boldsymbol{v})\quad\forall\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}\in
V,$ (18)

such that, in particular, $b(\boldsymbol{u},\boldsymbol{w},\boldsymbol{w})=0$
$\forall\boldsymbol{u},\boldsymbol{w}\in V$. We will use the following
discrete Gronwall inequality whose proof can be found in [15].

###### Lemma 1

Let $k,B,a_{j},b_{j},c_{j},\gamma_{j}$ be nonnegative numbers such that

$a_{n}+k\sum_{j=0}^{n}b_{j}\leq
k\sum_{j=0}^{n}\gamma_{j}a_{j}+k\sum_{j=0}^{n}c_{j}+B,\quad for\quad n\geq 0.$

Suppose that $k\gamma_{j}<1$, for all j, and set
$\sigma_{j}=(1-k\gamma_{j})^{-1}$. Then

$a_{n}+k\sum_{j=0}^{n}b_{j}\leq\exp\left(k\sum_{j=0}^{n}\sigma_{j}\gamma_{j}\right)\left\\{k\sum_{j=0}^{n}c_{j}+B\right\\},\quad
for\quad n\geq 0.$

We will also apply the following Gronwall lemma.

###### Lemma 2

Let $B,a_{j},b_{j},c_{j},\gamma_{j},\delta_{j}$ be nonnegative numbers such
that

$a_{n}+\sum_{j=0}^{n}b_{j}\leq\gamma_{n}a_{n}+\sum_{j=0}^{n-1}(\gamma_{j}+\delta_{j})a_{j}+\sum_{j=0}^{n}c_{j}+B,\quad
for\quad n\geq 0.$

Suppose that $\gamma_{j}<1$, for all j, and set
$\sigma_{j}=(1-\gamma_{j})^{-1}$. Then

$a_{n}+\sum_{j=0}^{n}b_{j}\leq\exp\left(\sigma_{n}\gamma_{n}+\sum_{j=0}^{n-1}(\sigma_{j}\gamma_{j}+\delta_{j})\right)\left\\{\sum_{j=0}^{n}c_{j}+B\right\\},\quad
for\quad n\geq 0.$

###### Proof

We follow the proof of [15, Lemma 5.1]. Let $d_{n}$ be defined by the relation

$a_{n}+\sum_{j=0}^{n}b_{j}+d_{n}=\gamma_{n}a_{n}+\sum_{j=0}^{n-1}(\gamma_{j}+\delta_{j})a_{j}+\sum_{j=0}^{n}c_{j}+B,\quad
for\quad n\geq 0,$

and let us denote by $S_{n}$ the left hand side above. Notice then

$S_{n}-S_{n-1}=\gamma_{n}a_{n}+\delta_{n-1}a_{n-1}+c_{n}\leq\gamma_{n}S_{n}+\delta_{n-1}S_{n-1}+c_{n},\qquad
n\geq 1.$

For $n=0$, we have

$S_{0}\leq\left(1+\gamma_{0}\right)^{-1}(c_{0}+B)=\left(1+\sigma_{0}\gamma_{0}\right)(c_{0}+B)\leq\exp(\sigma_{0}\gamma_{0})(c_{0}+B).$

Assume that

$S_{k}\leq\exp\biggl{(}\sigma_{k}\gamma_{k}+\sum_{j=0}^{k-1}(\sigma_{j}\gamma_{j}+\delta_{j})\biggr{)}\biggl{(}\sum_{j=0}^{k}c_{j}+B\biggr{)},\qquad
k=0,\ldots,n.$

and let us see that the inequality also holds for $k=n+1$, which will finish
the proof. We have

$\displaystyle S_{n+1}$
$\displaystyle\leq(1-\gamma_{n+1})^{-1}\left((1+\delta_{n})S_{n}+c_{n+1}\right)\leq(1+\sigma_{n+1}\gamma_{n+1})\left(\exp(\delta_{n})S_{n}+c_{n+1}\right)$
$\displaystyle\leq\exp\left(\sigma_{n+1}\gamma_{n+1}\right)\biggl{(}\exp\biggl{(}\sum_{k=0}^{n}(\sigma_{k}\gamma_{k}+\delta_{k})\biggr{)}\biggl{(}\sum_{k=0}^{n}c_{k}+B\biggr{)}+c_{n+1}\biggr{)}.$

$\Box$

Now, we consider a fully discrete approximation to the solution of (17). Let
us consider time levels $0=t_{0}<t_{1}<\ldots<t_{N}=T$, and let us denote

$\Delta t_{n}=t_{n+1}-t_{n},\qquad n=0,1\ldots,N-1,$

and

$\omega_{n}=\Delta t_{n}/\Delta t_{n-1},\qquad n=1,\ldots,N-1.$

In the sequel we will assume that for $\gamma\in(0,1)$ and $\Gamma>1$, we have

$\gamma\leq\omega_{n}\leq\Gamma,\qquad n=0,1,\ldots,N-1.$ (19)

We wil also use the following notation,

$D\boldsymbol{u}^{n}=\boldsymbol{u}^{n}-\boldsymbol{u}^{n-1},\quad\hat{\boldsymbol{u}}^{n}=\boldsymbol{u}^{n}+\omega_{n}D\boldsymbol{u}^{n},$

and

${\cal
D}\boldsymbol{u}^{n}=D\boldsymbol{u}^{n}+\frac{\omega_{n-1}}{1+\omega_{n-1}}\left(\boldsymbol{u}^{n}-\hat{\boldsymbol{u}}^{n-1}\right).$

Observe that $\hat{\boldsymbol{u}}^{n}$ is the extrapolated value at $t_{n+1}$
of the linear interpolant taking values $\boldsymbol{u}^{n}$ and
$\boldsymbol{u}^{n-1}$ at $t_{n}$ and $t_{n-1}$, respectively. Notice also
that for fixed step size, that is when $\Delta t_{n}=T/N$, $n=0,1,\ldots,N-1$,
we have

$\hat{\boldsymbol{u}}^{n}=2\boldsymbol{u}^{n}-\boldsymbol{u}^{n-1}\quad\hbox{\rm
and}\quad{\cal
D}\boldsymbol{u}^{n}=(D+\frac{1}{2}D^{2})\boldsymbol{u}^{n}=\frac{3}{2}\boldsymbol{u}^{n}-2\boldsymbol{u}^{n-1}+\frac{1}{2}\boldsymbol{u}^{n-2}.$

Now, we consider the following implicit-explicit (IMEX) approximation to (17)
based on the BDF2 method. Given
$\boldsymbol{u}_{h}^{0},\boldsymbol{u}_{h}^{1}$ solve, for $n\geq 1$ and for
all $(\boldsymbol{v}_{h},q_{h})\in V_{h}\times Q_{h}$

$\displaystyle\left(\frac{1}{\Delta t_{n}}{\cal
D}\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{u}_{h}^{n+1},\nabla\boldsymbol{v}_{h})+b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{v}_{h})-(p_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})\quad$
$\displaystyle\quad\quad+(\nabla\cdot\boldsymbol{u}_{h}^{n+1},q_{h})+\mu(\nabla\cdot\boldsymbol{u}_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})=(\boldsymbol{f}^{n+1},\boldsymbol{v}_{h}),$

where, here and in the sequel, the notation $\boldsymbol{f}^{n+1}$ means
$\boldsymbol{f}(t_{n+1})$. We also consider the semi-implicit method obtained
when the third term on the left-hand side of (3) is replaced by

$b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h})$
(21)

that is,

$\displaystyle\left(\frac{1}{\Delta t_{n}}{\cal
D}\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{u}_{h}^{n+1},\nabla\boldsymbol{v}_{h})+b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h})-(p_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})\quad$
$\displaystyle\quad\quad+(\nabla\cdot\boldsymbol{u}_{h}^{n+1},q_{h})+\mu(\nabla\cdot\boldsymbol{u}_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})=(\boldsymbol{f}^{n+1},\boldsymbol{v}_{h}),$

Taking $\boldsymbol{v}_{h}\in V_{h}^{\rm div}$ in (3) we get

$\displaystyle\left(\frac{1}{\Delta t_{n}}{\cal
D}\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{u}_{h}^{n+1},\nabla\boldsymbol{v}_{h})+b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{v}_{h})+\mu(\nabla\cdot\boldsymbol{u}_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})=(\boldsymbol{f}^{n+1},\boldsymbol{v}_{h}).$
(23)

To get $\boldsymbol{u}_{h}^{1}$ we 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 (3) and (23) are valid for $n=0$ if we set

$\hat{\boldsymbol{u}}^{0}=\boldsymbol{u}^{0},\quad{\cal
D}\boldsymbol{u}^{1}=D\boldsymbol{u}^{1}.$ (24)

We will compare $\boldsymbol{u}_{h}^{n}$ with $\boldsymbol{s}_{h}^{n}$. For
simplicity we will assume $\boldsymbol{u}_{h}^{0}=\boldsymbol{s}_{h}^{0}$. It
is easy to prove that for $n\geq 1$ and $\boldsymbol{v}_{h}\in V_{h}^{\rm
div}$

$\displaystyle\left(\frac{1}{\Delta t_{n}}{\cal
D}\boldsymbol{s}_{h}^{n+1},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{s}_{h}^{n+1},\nabla\boldsymbol{v}_{h})+b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{v}_{h})+\mu(\nabla\cdot\boldsymbol{s}_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})\quad$
$\displaystyle\quad=(\boldsymbol{f}^{n+1},\boldsymbol{v}_{h})+(p^{n+1}-\pi_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})+\mu(\nabla\cdot(\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}^{n+1}),\nabla\cdot\boldsymbol{v}_{h})$
$\displaystyle\quad+\left(\frac{1}{\Delta t_{n}}{\cal
D}\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}_{t}^{n+1},\boldsymbol{v}_{h}\right)+b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{v}_{h})-b(\boldsymbol{u}^{n+1},\boldsymbol{u}^{n+1},\boldsymbol{v}_{h}).$

Then, denoting by

$\boldsymbol{e}_{h}^{n}=\boldsymbol{s}_{h}^{n}-\boldsymbol{u}_{h}^{n}.$

and subtracting (23) from (3) we get

$\displaystyle\left(\frac{1}{\Delta t_{n}}{\cal
D}\boldsymbol{e}_{h}^{n+1},\boldsymbol{v}_{h}\right)+\nu(\nabla\boldsymbol{e}_{h}^{n+1},\nabla\boldsymbol{v}_{h})+b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{v}_{h})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{v}_{h})\quad$
(26)
$\displaystyle\quad+\mu(\nabla\cdot\boldsymbol{e}_{h}^{n+1},\nabla\cdot\boldsymbol{v}_{h})=(\tau_{1}^{n+1}+\tau_{3}^{n+1},\boldsymbol{v}_{h})+(\tau_{2}^{n+1},\nabla\cdot\boldsymbol{v}_{h}).$

where

$\displaystyle\tau_{1}^{n+1}$ $\displaystyle=$ $\displaystyle\frac{1}{\Delta
t_{n}}{\cal D}\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}_{t}^{n+1},$
$\displaystyle\tau_{2}^{n+1}$ $\displaystyle=$
$\displaystyle(p^{n+1}-\pi_{h}^{n+1})+\mu\nabla\cdot(\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}^{n+1}),$
(27) $\displaystyle(\tau_{3}^{n+1},\boldsymbol{v}_{h})$ $\displaystyle=$
$\displaystyle
b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{v}_{h})-b(\boldsymbol{u}^{n+1},\boldsymbol{u}^{n+1},\boldsymbol{v}_{h}).$

Similarly, for the semi-implicit method we get (23) with the third and foruth
terms replaced by
$b(\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{v}_{h})-b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{v}_{h})$,
and $\tau_{3}$ replaced by $\tau_{4}$ defined by

$(\tau_{4}^{n+1},\boldsymbol{v}_{h})=b(\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{v}_{h})-b(\boldsymbol{u}^{n+1},\boldsymbol{u}^{n+1},\boldsymbol{v}_{h}).$
(28)

In the following lemma we bound the difference of the nonlinear terms for the
IMEX method.

###### Lemma 3

Considering the fixed stepsize case $\Delta t_{n}=\Delta t=T/N$, the following
bound holds

$\displaystyle|b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|\leq\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+3\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}+\frac{\mu}{3}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\mu}{6}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
(29) $\displaystyle\quad+2\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}\|\boldsymbol{e}_{h}^{n+1}\|_{0}.$

In view of the upper bound in (19), for the variable stepsize case the bound
for the difference of the nonlinear terms is

$\displaystyle|b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|\leq\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9\alpha}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+\frac{(1+2\Gamma)(1+\Gamma)}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\Gamma(1+2\Gamma)}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
$\displaystyle\quad+\frac{\mu(1+2\Gamma)(1+\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}\quad+\frac{\mu\Gamma(1+2\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
(30) $\displaystyle\quad+2\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}\|\boldsymbol{e}_{h}^{n+1}\|_{0},$

where $\alpha$ is a positive constant.

###### Proof

We start by proving (3). Adding $\pm
b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})$,
and noticing that
$b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})=0$,
we have

$|b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|=|b(\hat{\boldsymbol{e}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|.$
(31)

For the first term on the right-hand side above we write

$\displaystyle|b(\hat{\boldsymbol{e}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|$
$\displaystyle\leq$
$\displaystyle\|\hat{\boldsymbol{e}}_{h}^{n}\|_{0}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n+1}\|_{0}+\|\hat{\boldsymbol{s}}_{h}\|_{\infty}\|\nabla\cdot\hat{\boldsymbol{e}}_{h}^{n}\|_{0}\|\boldsymbol{e}_{h}^{n+1}\|_{0}$
$\displaystyle\leq$
$\displaystyle\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\frac{1}{2}\|\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}+\frac{\mu}{18}\|\nabla\cdot\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}.$

Let us now observe that,

$\|\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}\leq
6\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+3\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2},$

and a similar expression for
$\|\nabla\cdot\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}$, so that,

$\displaystyle|b(\hat{\boldsymbol{e}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|$
$\displaystyle\leq$
$\displaystyle\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}{}+\frac{\mu}{3}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\mu}{6}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
$\displaystyle+3\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}.$

Finally, from the definition of the nonlinear term and using inverse
inequality (5) it is easy to check that

$|b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|\leq
2\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}\|\boldsymbol{e}_{h}^{n+1}\|_{0}.$

The proof of (3) can be obtained arguing in the same way with only two
differences. The first one is that we include a parameter $\alpha$ to bound
the second term on the first line of the right-hand side of (Proof). The
second one is that we apply (19) to bound both
$\|\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}$ and
$\|\nabla\cdot\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}$. $\Box$

For the semi-implicit method we have the following result.

###### Lemma 4

The following bound holds in the fixed stepsize case $\Delta t_{n}=\Delta
t=T/N$,

$\displaystyle|b(\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})|\leq\left(\frac{\|\nabla\boldsymbol{s}_{h}^{n+1}\|_{\infty}}{2}+\frac{9}{2\mu}\|\boldsymbol{s}_{h}^{n+1}\|_{\infty}^{2}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+3\|\nabla\boldsymbol{s}_{h}^{n+1}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{3}{2}\|\nabla\boldsymbol{s}_{h}^{n+1}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}+\frac{\mu}{3}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\mu}{6}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2},$
(33)

and, in the variable stepsize case, for $\Gamma$ the constant in (19) and
$\alpha$ any positive constant, the following bound holds:

$\displaystyle|b(\hat{\boldsymbol{s}}_{h}^{n}\boldsymbol{s}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})|\leq\left(\frac{\|\nabla\boldsymbol{s}_{h}^{n+1}\|_{\infty}}{2}+\frac{9\alpha}{2\mu}\|\boldsymbol{s}_{h}^{n+1}\|_{\infty}^{2}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+\frac{(1+2\Gamma)(1+\Gamma)}{2}\|\nabla\boldsymbol{s}_{h}^{n+1}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\Gamma(1+2\Gamma)}{2}\|\nabla\boldsymbol{s}_{h}^{n+1}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
$\displaystyle\quad+\frac{\mu(1+2\Gamma)(1+\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}\quad+\frac{\mu\Gamma(1+2\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}.$
(34)

###### Proof

By adding $\pm
b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})$
and taking into account that
$b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})=0$,
we have

$|b(\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})|=|b(\hat{\boldsymbol{e}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})|,$

which can be bounded as the term
$|b(\hat{\boldsymbol{e}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|$
in the previous lemma.

$\Box$

Remark 2 We notice the similarity between the bounds in Lemmas 3 and 4.
Indeed, bounds (4) and (4) differ from (3) and (3) in that
$\hat{\boldsymbol{s}}_{h}^{n}$ in (3) and (3) is replaced by
$\boldsymbol{s}_{h}^{n+1}$ in (4) and (4) and, more importantly, the last
term, $2\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}\|\boldsymbol{e}_{h}^{n+1}\|_{0}$
on the right-hand side of (3) and (3) is not present in (4) and (4). We will
see that this term gives rise to a CFL type condition that affects the IMEX
(3) method but not the semi-implicit method (3).

We now estimate the truncation errors. For $\tau_{2}^{n+1}$, in view of (6)
and (13), it easily follows that

$\left\|\tau_{2}^{n+1}\right\|_{0}\leq
C\bigl{(}\mu\left\|\boldsymbol{u}^{n+1}\right\|_{k+1}h^{k}+\left\|p^{n+1}\right\|_{l+1}h^{l+1}\bigr{)}.$
(35)

For $\tau_{1}^{n+1}$, $\tau_{3}^{n+1}$ and $\tau_{4}^{n+1}$ we treat
separately the cases $n\geq 1$ and $n=0$.

###### Lemma 5

The following bounds hold for $n\geq 1$.

$\displaystyle\left\|\tau_{1}^{n+1}\right\|_{0}^{2}$
$\displaystyle\leq\frac{4h^{2k}}{\Delta
t_{n}}\int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{u}_{t}(t)\right\|_{k}^{2}\,dt+\frac{2h^{2k}}{\Delta
t_{n-1}}\int_{t_{n-1}}^{t_{n}}\left\|\boldsymbol{u}_{t}(t)\right\|_{k}^{2}\,dt$
(36) $\displaystyle{}+\frac{(1+\omega_{n})^{2}}{10}(\Delta
t_{n})^{3}\int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{u}_{ttt}(t)\right\|_{0}^{2}\,ds+\frac{\omega_{n}^{2}}{10}(\Delta
t_{n}+\Delta
t_{n-1})^{3}\int_{t_{n-1}}^{t_{n+1}}\left\|\boldsymbol{u}_{ttt}(t)\right\|_{0}^{2}\,dt.$
$\displaystyle\left|(\tau_{3}^{n+1},\boldsymbol{v}_{h})\right|\leq$
$\displaystyle
C\Bigl{(}(1+\omega_{n})\left|\Omega\right|^{1/d}\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}h^{k}\left\|\boldsymbol{u}\right\|_{L^{\infty}(H^{k+1})}$
$\displaystyle{}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}\left\|\hat{\boldsymbol{u}}^{n}-\boldsymbol{u}^{n+1}\right\|_{0}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}\left\|\hat{\boldsymbol{u}}^{n}-\boldsymbol{u}^{n+1}\right\|_{1}\Bigr{)}\left\|\boldsymbol{v}_{h}\right\|_{0},$
(37) $\displaystyle\left|(\tau_{4}^{n+1},\boldsymbol{v}_{h})\right|\leq$
$\displaystyle
C\Bigl{(}(1+\omega_{n})\left|\Omega\right|^{1/d}\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}h^{k}\left\|\boldsymbol{u}\right\|_{L^{\infty}(H^{k+1})}$
$\displaystyle{}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}\left\|\hat{\boldsymbol{u}}^{n}-\boldsymbol{u}^{n+1}\right\|_{0}\Bigr{)}\left\|\boldsymbol{v}_{h}\right\|_{0},$
(38)

where

$\displaystyle\left\|\boldsymbol{u}^{n+1}-\hat{\boldsymbol{u}}^{n}\right\|_{l}\leq$
$\displaystyle\frac{1+\omega_{n}}{\sqrt{3}}(\Delta
t_{n})^{3/2}\biggl{(}\int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{u}_{tt}(t)\right\|_{l}^{2}\,dt\biggr{)}^{1/2}$
$\displaystyle{}+\frac{\omega_{n}}{\sqrt{3}}(\Delta t_{n}+\Delta
t_{n-1})^{3/2}\biggl{(}\int_{t_{n-1}}^{t_{n+1}}\left\|\boldsymbol{u}_{tt}(t)\right\|_{l}^{2}\,dt\biggr{)}^{1/2},\qquad
l=0,1.$ (39)

###### Proof

We first write $\tau_{1}=(\Delta t_{n})^{-1}{\cal
D}(\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}^{n+1})+(\Delta t_{n})^{-1}({\cal
D}\boldsymbol{u}^{n+1}-\boldsymbol{u}_{t}^{n+1})$ and, further, we notice that
we can write for the first term

$\displaystyle\frac{1}{\Delta t_{n}}{\cal
D}(\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}^{n+1})=$
$\displaystyle\Bigl{(}1+\frac{\omega_{n}}{1+\omega_{n}}\Bigr{)}\frac{1}{\Delta
t_{n}}D(\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}^{n+1})-\frac{\omega_{n}^{2}}{1+\omega_{n}}\frac{1}{\Delta
t_{n}}D(\boldsymbol{s}_{h}^{n}-\boldsymbol{u}^{n})$ $\displaystyle{}=$
$\displaystyle\Bigl{(}1+\frac{\omega_{n}}{1+\omega_{n}}\Bigr{)}\frac{1}{\Delta
t_{n}}\int_{t_{n}}^{t_{n+1}}\frac{d}{dt}(\boldsymbol{s}_{h}(t)-\boldsymbol{u}(t))\,dt$
$\displaystyle{}-\frac{\omega_{n}}{1+\omega_{n}}\frac{1}{\Delta
t_{n-1}}\int_{t_{n-1}}^{t_{n}}\frac{d}{dt}(\boldsymbol{s}_{h}(t)-\boldsymbol{u}(t))\,dt.$

Furthermore, Taylor expansion with integral remainder shows that the second
term can be bounded by

$\displaystyle\frac{1}{\Delta t_{n}}({\cal
D}\boldsymbol{u}^{n+1}-\boldsymbol{u}_{t}^{n+1})=$
$\displaystyle(1+\omega_{n})\frac{1}{2\Delta
t_{n}}\int_{t_{n}}^{t_{n+1}}(t-t_{n})^{2}\boldsymbol{u}_{ttt}(t)\,dt$
$\displaystyle{}-\frac{\omega_{n}^{2}}{2(1+\omega_{n})}\frac{1}{\Delta
t_{n}}\int_{t_{n-1}}^{t_{n+1}}(t-t_{n-1})^{2}\boldsymbol{u}_{ttt}(t)\,dt.$

We also notice that

$\frac{\omega_{n}^{2}}{2(1+\omega_{n})}\frac{1}{\Delta
t_{n}}=\frac{\omega_{n}}{2(1+\omega_{n})}\frac{1}{\Delta
t_{n-1}}=\frac{\omega_{n}}{2(\Delta t_{n}+\Delta t_{n-1})}.$

Thus, applying (13) and Hölder’s inequality an easy calculation shows that
(36) holds.

For $\tau_{3}^{n+1}$, adding $\pm
b(\hat{\boldsymbol{u}}^{n},\hat{\boldsymbol{s}}^{n}_{h},\boldsymbol{v}_{h})\pm
b(\hat{\boldsymbol{u}}^{n},\hat{\boldsymbol{u}}^{n},\boldsymbol{v}_{h})\pm
b(\hat{\boldsymbol{u}}^{n},\boldsymbol{u}^{n+1},\boldsymbol{v}_{h})$, we can
write

$\displaystyle(\tau_{3}^{n+1},\boldsymbol{v}_{h})=$ $\displaystyle
b(\hat{\boldsymbol{s}}_{h}^{n}-\hat{\boldsymbol{u}}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{v}_{h})+b(\hat{\boldsymbol{u}}^{n},\hat{\boldsymbol{s}}_{h}^{n}-\hat{\boldsymbol{u}}^{n},\boldsymbol{v}_{h})+b(\hat{\boldsymbol{u}}^{n},\hat{\boldsymbol{u}}^{n}-\boldsymbol{u}^{n+1},\boldsymbol{v}_{h})$
$\displaystyle{}+b(\hat{\boldsymbol{u}}^{n}-\boldsymbol{u}^{n+1},\boldsymbol{u}^{n+1},\boldsymbol{v}_{h}).$
(40)

Thus, by noticing that
$\left\|\hat{\boldsymbol{s}}_{h}^{n}-\hat{\boldsymbol{u}}^{n}\right\|_{l}\leq
2(1+\omega_{n})\left\|\boldsymbol{s}_{h}-\boldsymbol{u}\right\|_{L^{\infty}(H^{l})}$,
$l=0,1$, applying (13), Hölder’s inequality and bounding
$h\leq\left|\Omega\right|^{1/d}$, when necessary, we get (5). For $\tau_{4}$,
adding $\pm
b(\hat{\boldsymbol{u}}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{v}_{h})\pm
b(\hat{\boldsymbol{u}}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{v}_{h})$ one
gets

$\displaystyle(\tau_{4}^{n+1},\boldsymbol{v}_{h})=$ $\displaystyle
b(\hat{\boldsymbol{s}}_{h}^{n}-\hat{\boldsymbol{u}}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{v}_{h})+b(\hat{\boldsymbol{u}}^{n}-\boldsymbol{u}^{n+1},\boldsymbol{s}_{h}^{n+1},\boldsymbol{v}_{h})+b(\boldsymbol{u}^{n+1},\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}^{n+1},\boldsymbol{v}_{h}),$
(41)

from where, arguing as with $\tau_{3}^{n+1}$, the bound (5) follows.

Also, by means of Taylor expansion with integral remainder and applying
Hölder’s inequality, it is easy to show (5). $\Box$

We now deal with the truncation errors in the first step. Recall that

$\tau_{1}^{1}=\frac{1}{\Delta
t_{0}}{D}\boldsymbol{s}_{h}^{1}-\boldsymbol{u}_{t}^{1}$

and $\tau_{3}^{1}$ are defined as in (3) but with
$\hat{\boldsymbol{s}}_{h}^{0}=\boldsymbol{s}_{h}^{0}$.

###### Lemma 6

We have the following bounds,

$\displaystyle\|\tau_{1}^{1}\|_{0}^{2}\leq$ $\displaystyle
C\Bigl{(}\frac{1}{\Delta
t_{0}}h^{2k}\|\boldsymbol{u}_{t}\|_{L^{2}(H^{k})}+(\Delta
t_{0})^{2}\|\boldsymbol{u}_{tt}\|_{L^{\infty}(L^{2})}^{2}\Bigr{)},$ (42)
$\displaystyle\left\|\tau_{3}^{1}\right\|_{0}^{2}\leq$ $\displaystyle
C\Bigl{(}\left|\Omega\right|^{2/d}\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}^{2}h^{2k}$
$\displaystyle\quad{}+(\Delta
t_{0})^{2}\bigl{(}\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}^{2}\left\|\boldsymbol{u}_{t}\right\|_{L^{\infty}(L^{2})}^{2}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}\left\|\boldsymbol{u}_{t}\right\|_{L^{\infty}(H^{1})}^{2}\bigr{)}\Bigr{)}.$
(43)

###### Proof

For $\tau_{1}^{1}$, by arguing as in (36) we have

$\|\tau_{1}^{1}\|_{0}^{2}\leq\frac{C}{\Delta
t_{0}}h^{2k}\int_{t_{0}}^{t_{1}}\|\boldsymbol{u}_{t}\|_{k}^{2}\,dt+C\Delta
t_{0}\int_{t_{0}}^{t_{1}}\|\boldsymbol{u}_{tt}\|_{0}^{2},\,dt,$

from where (42) follows. To prove (6) we go back to (Proof), and repeat the
argument there, but now taking into account that for $n=0$,
$\hat{\boldsymbol{s}}_{h}^{0}=\boldsymbol{s}_{0}$ and
$\hat{\boldsymbol{u}}^{0}=\boldsymbol{u}^{0}$, so that
$\left\|\hat{\boldsymbol{s}}_{h}^{0}-\boldsymbol{u}^{0}\right\|_{l}\leq
Ch^{k+1-l}\left\|\boldsymbol{u}^{0}\right\|_{k+1}$ and

$\left\|\hat{\boldsymbol{u}}^{0}-\boldsymbol{u}^{1}\right\|_{l}=\biggl{\|}\int_{t_{0}}^{t_{1}}\boldsymbol{u}_{t}\,dt\biggr{\|}_{l}\leq\Delta
t_{0}\left\|\boldsymbol{u}_{t}\right\|_{L^{\infty}(H^{l})}.$

Then (6) follows easily. $\Box$

###  3.1 Fixed stepsize

We consider now the case where $\Delta t_{n}=\Delta t=T/N$,
$n=0,1,\ldots,N-1$.

###### Theorem 1

Fix $\kappa>0$ and let $\Delta t$ and $h$ satisfy the following CFL-type
condition

$\displaystyle 64c_{\rm
inv}^{2}\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}\frac{\Delta
t}{h^{2}}\leq\frac{\kappa}{T},$ (44)

and

$\displaystyle\Delta
t\Bigl{(}4\frac{\kappa+1}{T}+M^{k}\Bigr{)}\leq\frac{1}{2},\quad
k=2,\ldots,n+1,$ (45)

with $M^{k}$ defined in (54). Then, there exist a positive constant $h_{0}$,
such that for $h<h_{0}$ the error
$\boldsymbol{e}_{h}^{n+1}=\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}_{h}^{n+1}$
of the IMEX method (3), satisfies the following bound for $1\leq n\leq
N-1=T/\Delta t-1$,

$\displaystyle\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\sum_{k=2}^{n+1}\nu\Delta
t\|\nabla\boldsymbol{e}_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$ $\displaystyle\quad
e^{8(\kappa+1)+2T\max_{j}M^{j}}\left(\left(C_{0}^{2}+TC_{3}^{2}\right)(\Delta
t)^{4}+T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}\right)\right),$ (46)

where the constants $C_{i}^{2}$ are defined in (59) and (63).

###### Proof

Taking $\boldsymbol{v}_{h}=\boldsymbol{e}_{h}^{n+1}$ in (26) we obtain

$\displaystyle\frac{1}{4}\frac{1}{\Delta
t}\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{1}{4}\frac{1}{\Delta
t}\|\hat{\boldsymbol{e}}_{h}^{n+1}\|_{0}^{2}+\frac{1}{4}\frac{1}{\Delta
t}\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}$ (47)
$\displaystyle\quad-\frac{1}{4}\frac{1}{\Delta
t}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}-\frac{1}{4}\frac{1}{\Delta
t}\|\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}+\nu\|\nabla\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+\mu\|\nabla\cdot\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}\leq|b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|+\frac{T}{2}\|\tau_{1}^{n+1}\|_{0}^{2}+\frac{T}{2}\|\tau_{3}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+\frac{1}{T}\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{9}{2\mu}\|\tau_{2}^{n+1}\|_{0}^{2}+\frac{\mu}{18}\|\nabla\cdot\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}.$
(48)

We now apply (3) from Lemma 3 and, further, we bound

$2\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h}^{n+1}-\hat{e}_{h}^{n}\|_{0}\|\boldsymbol{e}_{h}^{n+1}\|_{0}\leq
4c_{\rm inv}^{2}\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}^{2}\frac{\Delta
t}{h^{2}}\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}}{4\Delta
t}.$ (49)

Then

$\displaystyle|b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h})|$
$\displaystyle\leq$
$\displaystyle\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}+4c_{\rm
inv}^{2}\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}^{2}\frac{\Delta
t}{h^{2}}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+3\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
$\displaystyle\quad+\frac{\mu}{3}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\mu}{6}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}+\frac{\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}}{4\Delta
t}.$

The last term on the right-hand side of (Proof) will be absorbed with the
left-hand side of (47). In the sequel, we assume that, for $0\leq n\leq N-1$,

$\displaystyle 4c_{\rm
inv}^{2}\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}^{2}\frac{\Delta
t}{h^{2}}\leq\frac{\kappa}{T}.$ (51)

At the end of the proof we will show that (51) always holds for $h$
sufficiently small.

Inserting (Proof) into (47) we reach

$\displaystyle\frac{1}{4}\frac{1}{\Delta
t}\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{1}{4}\frac{1}{\Delta
t}\|\hat{\boldsymbol{e}}_{h}^{n+1}\|_{0}^{2}-\frac{1}{4}\frac{1}{\Delta
t}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}-\frac{1}{4}\frac{1}{\Delta
t}\|\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}$
$\displaystyle\quad+\nu\|\nabla\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{17}{18}\mu\|\nabla\cdot\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}\leq$
(52)
$\displaystyle\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}+\frac{\kappa+1}{T}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+3\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
$\displaystyle\quad+\frac{\mu}{3}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\mu}{6}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}+\frac{T}{2}\|\tau_{1}^{n+1}\|_{0}^{2}+\frac{T}{2}\|\tau_{3}^{n+1}\|_{0}^{2}+\frac{9}{2\mu}\|\tau_{2}^{n+1}\|_{0}^{2}.$

Taking the sum of terms in (Proof) we get

$\displaystyle\frac{1}{4}\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{1}{4}\|\hat{\boldsymbol{e}}_{h}^{n+1}\|_{0}^{2}+\sum_{k=2}^{n+1}\nu\Delta
t\|\nabla\boldsymbol{e}_{h}^{k}\|_{0}^{2}+\frac{4\mu}{9}\sum_{k=2}^{n+1}\Delta
t\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq\frac{1}{4}\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\frac{1}{4}\|\hat{\boldsymbol{e}}_{h}^{1}\|_{0}^{2}$
$\displaystyle\quad+\frac{\mu}{2}(\Delta
t)\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\frac{\mu}{6}(\Delta
t)\|\nabla\cdot\boldsymbol{e}_{h}^{0}\|_{0}^{2}+\frac{3}{2}(\Delta
t)\|\nabla\hat{\boldsymbol{s}}_{h}^{1}\|_{\infty}\|\boldsymbol{e}_{h}^{0}\|_{0}^{2}$
(53) $\displaystyle\quad+(\Delta
t)\left(\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{2}\|_{\infty}+3\|\nabla\hat{\boldsymbol{s}}_{h}^{1}\|_{\infty}\right)\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\sum_{k=2}^{n+1}\Delta
t\left(\frac{T}{2}\|\tau_{1}^{k}\|_{0}^{2}+\frac{T}{2}\|\tau_{3}^{k}\|_{0}^{2}+\frac{9}{2\mu}\|\tau_{2}^{k}\|_{0}^{2}\right)$
$\displaystyle\quad+\sum_{k=2}^{n+1}\Delta
t\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{k-1}\|_{\infty}}{2}+\frac{9}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{k-1}\|_{\infty}^{2}+\frac{\kappa+1}{T}+3\|\nabla\hat{\boldsymbol{s}}_{h}^{k}\|_{\infty}+\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{k+1}\|_{\infty}\right)\|\boldsymbol{e}_{h}^{k}\|_{0}^{2},$

where in the last term of (Proof) to simplify writing we assume
$\|\nabla\hat{\boldsymbol{s}}_{h}^{k+1}\|_{\infty}=0$ for $k=n+1$. In the
sequel, to simplify notation we define

$M^{k}=4\biggl{(}\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{k-1}\|_{\infty}}{2}+\frac{9}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{k-1}\|_{\infty}^{2}+3\|\nabla\hat{\boldsymbol{s}}_{h}^{k}\|_{\infty}+\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{k+1}\|_{\infty}\biggr{)}$
(54)

and

$\displaystyle E^{0}$ $\displaystyle=$ $\displaystyle
4\biggl{(}\frac{1}{4}\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\frac{1}{4}\|\hat{\boldsymbol{e}}_{h}^{1}\|_{0}^{2}+\frac{\mu}{2}(\Delta
t)\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\frac{\mu}{6}(\Delta
t)\|\nabla\cdot\boldsymbol{e}_{h}^{0}\|_{0}^{2}+\frac{3}{2}(\Delta
t)\|\nabla\hat{\boldsymbol{s}}_{h}^{1}\|_{\infty}\|\boldsymbol{e}_{h}^{0}\|_{0}^{2}$
(55) $\displaystyle\quad+(\Delta
t)\left(\frac{3}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{2}\|_{\infty}+3\|\nabla\hat{\boldsymbol{s}}_{h}^{1}\|_{\infty}\right)\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}\biggr{)}.$

Let us observe that all the terms in (54) can be bounded in terms of
$\|\boldsymbol{u}\|_{L^{\infty}(L^{\infty})}$ and
$\|\nabla\boldsymbol{u}\|_{L^{\infty}(L^{\infty})}$ applying (15) and (16).

With the above simplifications we can write equation (Proof) in the form

$\displaystyle\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\sum_{k=2}^{n+1}\nu\Delta
t\|\nabla\boldsymbol{e}_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}$ $\displaystyle\leq
E^{0}+\sum_{k=2}^{n+1}\Delta
t\Bigl{(}4\frac{\kappa+1}{T}+M^{k}\Bigr{)}\|\boldsymbol{e}_{h}^{k}\|_{0}^{2}$
(56) $\displaystyle{}+2\sum_{k=2}^{n+1}\Delta
t\left(T(\|\tau_{1}^{k}\|_{0}^{2}+\|\tau_{3}^{k}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{k}\|_{0}^{2}\right).$

Assuming (45) we can apply Gronwall Lemma 1 to obtain

$\displaystyle\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\sum_{k=2}^{n+1}\nu\Delta
t\|\nabla\boldsymbol{e}_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$ (57) $\displaystyle\quad
e^{8(\kappa+1)+2T\max_{j}M^{j}}\left(E^{0}+2\sum_{k=2}^{n+1}\Delta
t\left(T(\|\tau_{1}^{k}\|_{0}^{2}+\|\tau_{3}^{k}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{k}\|_{0}^{2}\right)\right).$

In view of (35) and (36)-(5)-(5) from Lemma 5 and taking into account that
$\omega_{n}=1$, $n=1,\ldots,N=T/\Delta t$, we have

$\Delta t\sum_{k=2}^{n}\left\|\tau_{1}^{k}\right\|_{0}^{2}\leq
C\Bigl{(}\left\|\boldsymbol{u}_{t}\right\|_{L^{2}(H^{k})}^{2}h^{2k}+\left\|\boldsymbol{u}_{ttt}\right\|_{L^{2}(L^{2})}^{2}(\Delta
t)^{4}\Bigr{)},$ $\Delta
t\sum_{k=2}^{n}\left\|\tau_{2}^{k}\right\|_{0}^{2}\leq
CT\left(\mu^{2}\left\|\boldsymbol{u}\right\|_{L^{\infty}(H^{k+1})}^{2}h^{2k}+\left\|p\right\|_{L^{\infty}(H^{l+1})}^{2}h^{2(l+1)}\right),$

and

$\displaystyle\Delta t\sum_{k=2}^{n}\left\|\tau_{3}^{k}\right\|_{0}^{2}\leq$
$\displaystyle
C\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}^{2}\left(T\left|\Omega\right|^{2/d}\left\|\boldsymbol{u}\right\|_{L^{\infty}(H^{k+1})}^{2}h^{2k}+\left\|\boldsymbol{u}_{tt}\right\|_{L^{2}(L^{2})}^{2}(\Delta
t)^{4}\right)$
$\displaystyle+C\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}\left\|\boldsymbol{u}_{tt}\right\|_{L^{2}(H^{1})}^{2}(\Delta
t)^{4}.$

Consecuently,

$\sum_{k=2}^{n+1}\Delta
t\left(T(\|\tau_{1}^{k}\|_{0}^{2}+\|\tau_{3}^{k}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{k}\|_{0}^{2}\right)\leq
T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}+C_{3}^{2}(\Delta t)^{4}\right),$ (58)

where

$\displaystyle C_{1}^{2}$
$\displaystyle=C\left(\left\|\boldsymbol{u}_{t}\right\|_{L^{2}(H^{k})}^{2}+\left(\mu+T\left|\Omega\right|^{2/d}\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}^{2}\right)\left\|\boldsymbol{u}\right\|_{L^{\infty}(H^{k+1})}^{2}\right),$
$\displaystyle C_{2}^{2}$
$\displaystyle=C\frac{1}{\mu}\left\|p\right\|_{L^{\infty}(H^{l+1})}^{2},$ (59)
$\displaystyle C_{3}^{2}$
$\displaystyle=C\left(\left\|\boldsymbol{u}_{ttt}\right\|_{L^{2}(L^{2})}^{2}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}^{2}\left\|\boldsymbol{u}_{tt}\right\|_{L^{2}(L^{2})}^{2}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}\left\|\boldsymbol{u}_{tt}\right\|_{L^{2}(H^{1})}^{2}\right).$

Finally, since we have assumed $\boldsymbol{u}_{h}^{0}=\boldsymbol{s}_{h}^{0}$
to bound $E_{0}$ we only need to bound $\boldsymbol{e}_{h}^{1}$. To this end,
let us observe that after one step using Euler method with the explicit form
of the nonlinear term, and arguing as before, it is easy to get

$\frac{1}{\Delta
t_{0}}\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\nu\|\nabla\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\mu\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2}=(\tau_{1}^{1}+\tau_{3}^{1},\boldsymbol{e}_{h}^{1})+(\tau_{2}^{1},\nabla\cdot\boldsymbol{e}_{h}^{1}).$
(60)

From (60) we get

$\frac{1}{2}\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\frac{\mu}{2}\Delta
t_{0}\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2}\leq\frac{(\Delta
t_{0})^{2}}{2}\|\tau_{1}^{1}+\tau_{3}^{1}\|_{0}^{2}+\frac{\Delta
t_{0}}{2\mu}\|\tau_{2}^{1}\|_{0}^{2}.$

Now, in view of (42) and taking into account that $\Delta t\leq T$, we have

$(\Delta t)^{2}\|\tau_{1}^{1}\|_{0}^{2}\leq
C\|\boldsymbol{u}_{tt}\|_{L^{\infty}(L^{2})}^{2}(\Delta
t)^{4}+TC_{1}^{2}h^{2k},$ (61)

where $C_{1}$ is the constant in (59). Also in view of (6) we have

$(\Delta t)^{2}\left\|\tau_{3}^{1}\right\|_{0}^{2}\leq C(\Delta
t)^{4}\bigl{(}\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}^{2}\left\|\boldsymbol{u}_{t}\right\|_{L^{\infty}(L^{2})}^{2}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}\left\|\boldsymbol{u}_{t}\right\|_{L^{\infty}(H^{1})}^{2}\bigr{)}+TC_{1}h^{2k}.$

Thus, we finally obtain

$\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\mu\Delta
t\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2}\leq C_{0}^{2}(\Delta
t)^{4}+T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}\right),$ (62)

where

$C_{0}^{2}=C\bigl{(}\|\boldsymbol{u}_{tt}\|_{L^{\infty}(L^{2})}^{2}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}^{2}\left\|\boldsymbol{u}_{t}\right\|_{L^{\infty}(L^{2})}^{2}+\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}\left\|\boldsymbol{u}_{t}\right\|_{L^{\infty}(H^{1})}^{2}\bigr{)}.$
(63)

Going back to the definition of $E_{0}$, (55), we reach

$E_{0}\leq\left(5+\Delta
t\left[6\|\nabla\hat{\boldsymbol{s}}_{h}^{2}\|_{\infty}+12\|\nabla\hat{\boldsymbol{s}}_{h}^{1}\|_{\infty}\right]\right)\|\boldsymbol{e}_{h}^{1}\|_{0}^{2}+2\mu(\Delta
t)\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2},$

so that applying (62) and (16) we finally obtain

$E_{0}\leq C\left(1+\Delta
t\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}\right)\left(C_{0}^{2}(\Delta
t)^{4}+T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}\right)\right).$ (64)

Inserting (58), (64) into (57) we conclude (1) as long as condition (51)
holds. Thus, in order to finish the proof, we have to check that this is the
case if $h$ and $\Delta t$ are taken sufficiently small. To this end we devote
the rest of the proof.

We will show that for $h$ sufficiently small it holds

$\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}\leq
4\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})},$ (65)

for $0\leq n\leq N-1$, so that (51) will be a consequence of (44). For this
purpose, we start by taking $h_{0,1}$ such that the right-hand side of (15) is
bounded by $\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}/6$, so that
by writing
$\boldsymbol{s}_{h}=\boldsymbol{s}_{h}-\boldsymbol{u}+\boldsymbol{u}$, and in
view of (15), we have that, if $h\leq h_{0,1}$, then,
$\left\|\boldsymbol{s}_{h}\right\|_{L^{\infty}(L^{\infty})}\leq(7/6)\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}$,
and, consequently,

$\left\|\hat{\boldsymbol{s}}_{h}^{n}\right\|\leq(7/2)\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})},\qquad
0\leq t_{n}\leq T.$ (66)

For $\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}$, adding and subtracting
$\hat{\boldsymbol{s}}_{h}^{n}$ and using inverse inequality (5) we get

$\displaystyle\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}$
$\displaystyle\leq\|\hat{\boldsymbol{u}}_{h}^{n}-\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}+\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\leq
c_{\rm
inv}h^{-d/2}\|\hat{\boldsymbol{u}}_{h}^{n}-\hat{\boldsymbol{s}}_{h}^{n}\|_{0}+\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}$
$\displaystyle\leq c_{\rm
inv}h^{-d/2}2\|\boldsymbol{u}_{h}^{n}-\boldsymbol{s}_{h}^{n}\|_{0}+c_{\rm
inv}h^{-d/2}\|\boldsymbol{u}_{h}^{n-1}-\boldsymbol{s}_{h}^{n-1}\|_{0}+(7/2)\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})},$
(67)

where, in the last inequality, we have also applied (66). We now take
$h_{0,2}\leq h_{0,1}$ such that if $h\leq h_{0,2}$ the right-hand sides of (1)
and (62) are smaller that

$r(h)=\frac{h^{d}}{36c_{\rm inv}^{2}}.$ (68)

Notice that this is possible since, on the one hand, due to (44), $\Delta t$
can be bounded in terms of $h^{2}$, and, on the other hand, the right-hand
sides of (1) and (62), being at least $O(h^{4})$, decay faster with $h$ than
$r(h)$ in (68), which is, at most, $O(h^{3})$.

For $h\leq h_{0,2}$, and since we are taking
$\boldsymbol{u}_{h}^{0}=\boldsymbol{s}_{h}^{0}$, in view of (Proof), we have
that (65), holds for $n=0,1$. Assuming that (65) holds for $n\leq m$, we will
now show that it also holds for $n=m+1$, and this will finish the proof.
Indeed, if (65) holds for $n\leq m$, then, as argued above, (1) holds for
$n=m$. But, since $h\leq h_{0,2}$, the right-hand side of (1) is smaller that
$r(h)$ in (68), and, consequently, for $n=m+1$, the right-hand side of (Proof)
is smaller that
$(7/2+1/2)\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}$, that is,
(65) also holds for $n=m+1$. $\Box$

Remark 3 In view of (44), we observe that as $\Delta t,h\rightarrow 0$ the
ratio $\Delta t/h^{2}$ must be bounded. One can allow the upper bound to be
larger or smaller, but its size reflects in the error bound (1) through the
factor $\exp(8(\kappa+1))$ on the right-hand side of the error bound (1).

As pointed out in the introduction, arguing as in [9], it is possible to get a
different CFL condition in which the quantity that has to be bounded is
$(\Delta t)h^{-1}\nu^{-1}$ instead of $(\Delta t)h^{-2}$ and that involves a
weaker norm for $\boldsymbol{u}$ in (44). However, since we focus on getting
bounds that hold for high Reynolds numbers we will develop this line of
research in future works.

For the semi-implicit method (3) we have the following result.

###### Theorem 2

Let $\Delta t$ satisfy the following bound for $k=2,\ldots,n+1$,

$\Delta t\Bigl{(}\frac{1}{T}+N^{k}\Bigr{)}\leq\frac{1}{2},\quad
k=2,\ldots,n+1,$

where

$N_{k}=4\biggl{(}\frac{\|\nabla\boldsymbol{s}_{h}^{k}\|_{\infty}}{2}+\frac{9}{2\mu}\|\boldsymbol{s}_{h}^{k}\|_{\infty}^{2}+3\|\nabla\boldsymbol{s}_{h}^{k+1}\|_{\infty}+\frac{3}{2}\|\nabla\boldsymbol{s}_{h}^{k+2}\|_{\infty}\biggr{)}.$

Then, there exist a positive constant $h_{0}$, such that for $h<h_{0}$ the
error
$\boldsymbol{e}_{h}^{n+1}=\boldsymbol{s}_{h}^{n+1}-\boldsymbol{u}_{h}^{n+1}$
of the semi-implicit method (3), satisfies the following bound for $1\leq
n\leq N-1=T/\Delta t-1$,

$\displaystyle\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\sum_{k=2}^{n+1}\nu\Delta
t\|\nabla\boldsymbol{e}_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$ $\displaystyle\quad
e^{8+2T\max_{j}N^{j}}\left(\left(C_{0}^{2}+TC_{3}^{2}\right)(\Delta
t)^{4}+T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}\right)\right),$

where the constants $C_{i}^{2}$ are defined in (59) and (63).

###### Proof

Arguing as in the case of method (3), we have that (47) also holds for the
semi-implicit method but with the first term on the right-hand side of (47)
replaced by

$|b(\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})|$

and $\tau_{3}^{n}$ replaced by $\tau_{4}^{n}$. Then, recalling Lemma 4 and
taking into account Remark Proof, we have that (Proof) also holds true for
method (3) with $\kappa=0$, $\hat{\boldsymbol{s}}_{h}^{n}$ replaced by
$\boldsymbol{s}_{h}^{n+1}$ and $\tau_{3}^{n}$ replaced by $\tau_{4}^{n}$. Then
the proof follows the steps of the proof of Theorem 1. In the present case,
there is no need to argue about the boundedness of $\boldsymbol{u}_{h}^{n}$
since it is not necessary that condition (51) is satisfied. $\Box$

Remark 4 From Theorems 1 and 2, using triangle inequality and applying (13),
we reach

$\max_{1\leq n\leq
N}\|\boldsymbol{u}_{h}^{n}-\boldsymbol{u}^{n}\|_{0}=O((\Delta
t)^{2}+h^{k}+h^{l+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 $k$ in 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

$\displaystyle G_{n}$
$\displaystyle=\frac{\omega_{n-1}}{2(1+\omega_{n-1})}\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}+\frac{1}{4}\left\|\boldsymbol{e}_{h}^{n}+(\boldsymbol{e}_{h}^{n}-\boldsymbol{e}_{h}^{n-1})\right\|_{0}^{2},\qquad
n=1,2,\ldots,N,$ (69)

the following inequality holds, for a constant $K_{0}\leq 75/14$,

$\displaystyle\Delta t_{n}({\cal
D}\boldsymbol{e}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})$ $\displaystyle\geq
G_{n+1}-G_{n}+\frac{\omega_{n}}{2(1+\omega_{n})}\left\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2}$
(70)
$\displaystyle\quad-\frac{K_{0}}{\gamma}(1+\Gamma)^{2}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right)G_{n}.$

###### Proof

We first rewrite adequately the term $({\cal
D}\boldsymbol{e}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})$. Using the identity

$(\boldsymbol{v}-\boldsymbol{w},\boldsymbol{v})=\frac{1}{2}\bigl{(}\left\|\boldsymbol{v}\right\|_{0}^{2}-\left\|\boldsymbol{w}\right\|_{0}^{2}+\left\|\boldsymbol{v}-\boldsymbol{w}\right\|_{0}^{2}),$
(71)

we have that

$\displaystyle\Delta t_{n}({\cal
D}\boldsymbol{e}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})=$
$\displaystyle\frac{1}{2}\bigl{(}\left\|\boldsymbol{e}_{h}^{n+1}\right\|_{0}^{2}-\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}+\left\|\boldsymbol{e}_{h}^{n+1}-\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}\bigr{)}$
(72)
$\displaystyle{}+\frac{\omega_{n}}{2(1+\omega_{n})}\bigl{(}\left\|\boldsymbol{e}_{h}^{n+1}\right\|_{0}^{2}-\left\|\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2}+\left\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2}\bigr{)}.$

For the first term on the right-hand side above, using the identity (71), we
write

$\displaystyle\frac{1}{2}\bigl{(}\left\|\boldsymbol{e}_{h}^{n+1}\right\|_{0}^{2}-\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}+\left\|\boldsymbol{e}_{h}^{n+1}-\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}\bigr{)}=$
$\displaystyle\frac{1}{4}\bigl{(}\left\|\boldsymbol{e}_{h}^{n+1}\right\|_{0}^{2}-\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}+\left\|\boldsymbol{e}_{h}^{n+1}-\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}\bigr{)}$
$\displaystyle{}+\frac{1}{2}(\boldsymbol{e}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1}-\boldsymbol{e}_{h}^{n})$
$\displaystyle=$
$\displaystyle\frac{1}{4}\left\|\boldsymbol{e}_{h}^{n+1}+(\boldsymbol{e}_{h}^{n+1}-\boldsymbol{e}_{h}^{n})\right\|_{0}^{2}-\frac{1}{4}\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}.$

Then

$\displaystyle\Delta t_{n}({\cal
D}\boldsymbol{e}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})$
$\displaystyle=G_{n+1}-\Bigl{(}\frac{1}{4}\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}+\frac{\omega_{n}}{2(1+\omega_{n})}\left\|\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2}\Bigr{)}+\frac{\omega_{n}}{2(1+\omega_{n})}\left\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2}$
$\displaystyle=G_{n+1}-G_{n}+{\cal
R}_{n}+\frac{\omega_{n}}{2(1+\omega_{n})}\left\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2},$
(73)

where

${\cal
R}_{n}=G_{n}-\Bigl{(}\frac{1}{4}\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}+\frac{\omega_{n}}{2(1+\omega_{n})}\left\|\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2}\Bigr{)}.$

Noticing that

$\left\|\boldsymbol{e}_{h}^{n}+\omega_{n}(\boldsymbol{e}_{h}^{n}-\boldsymbol{e}_{h}^{n-1})\right\|_{0}^{2}=(1+\omega_{n})^{2}\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}-2\omega_{n}(1+\omega_{n})(\boldsymbol{e}_{h}^{n},\boldsymbol{e}_{h}^{n-1})+\omega_{n}^{2}\left\|\boldsymbol{e}_{h}^{n-1}\right\|_{0}^{2},$

a straightforward computation shows

${\cal
R}_{n}=c_{1,n}\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}-c_{2,n}(\boldsymbol{e}_{h}^{n},\boldsymbol{e}_{h}^{n-1})+c_{3,n}\left\|\boldsymbol{e}_{h}^{n-1}\right\|_{0}^{2},$

where

$\displaystyle c_{1,n}$
$\displaystyle=\frac{\omega_{n-1}}{2(1+\omega_{n-1})}-\frac{1}{4}+1-\frac{\omega_{n}(1+\omega_{n})}{2},$
(74) $\displaystyle c_{2,n}$ $\displaystyle=1-\omega_{n}^{2},$ (75)
$\displaystyle c_{3,n}$
$\displaystyle=\frac{1}{4}-\frac{\omega_{n}^{3}}{2(1+\omega_{n})}.$ (76)

We now express these three coefficients in terms of $\omega_{n}-1$ and
$\omega_{n-1}-1$. For $c_{1,n}$, we have

$c_{1,n}=\frac{\omega_{n-1}-1}{4(1+\omega_{n-1})}+\frac{2-\omega_{n}-\omega_{n}^{2}}{2}=\frac{\omega_{n-1}-1}{4(1+\omega_{n-1})}+\frac{2+\omega_{n}}{2}(1-\omega_{n}).$

For $c_{3,n}$, we have

$c_{3,n}=\frac{1+\omega_{n}-2\omega_{n}^{3}}{4(1+\omega_{n})}=\frac{1+2\omega_{n}+2\omega_{n}^{2}}{4(1+\omega_{n})}(1-\omega_{n}).$

From (19) it is easy to check that

$\displaystyle\left|c_{1,n}\right|$
$\displaystyle\leq\frac{1}{4}\left|\omega_{n-1}-1\right|+(1+\Gamma/2)\left|\omega_{n}-1\right|,$
$\displaystyle\left|c_{2,n}\right|$
$\displaystyle\leq(1+\Gamma)\left|\omega_{n}-1\right|,$ (77)
$\displaystyle\left|c_{3,n}\right|$
$\displaystyle\leq\Bigl{(}\frac{1}{4}+\frac{1+\Gamma}{2}\Bigr{)}\left|\omega_{n}-1\right|,$
(78)

so that

$\max_{1\leq j\leq
3}\left|c_{j,n}\right|\leq\frac{1}{4}\left|\omega_{n-1}-1\right|+(1+\Gamma)\left|\omega_{n}-1\right|.$
(79)

Furthermore, since

$G_{n}=\Bigl{(}1+\frac{\omega_{n-1}}{2(1+\omega_{n-1})}\Bigr{)}\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}-(\boldsymbol{e}_{h}^{n},\boldsymbol{e}_{h}^{n-1})+\frac{1}{4}\left\|\boldsymbol{e}_{h}^{n-1}\right\|_{0}^{2},$

and the smallest eigenvalue of matrix

$\left[\begin{array}[]{cc}1+x&-1/2\\\ -1/2&1/4\end{array}\right]$ (80)

is $\lambda=\left((5+4x)-\sqrt{(5+4x)^{2}-16x}\right)/8$, which can be seen to
be $\lambda\geq 0.14x$, for $0\leq x\leq 1/2$ (see Fig. 1). Then, if follows
that

$G_{n}\geq\frac{0.07}{(1+\Gamma)}{\omega_{n-1}}\bigl{(}\left\|\boldsymbol{e}_{h}^{n-1}\right\|_{0}^{2}+\left\|\boldsymbol{e}_{h}^{n}\right\|_{0}^{2}\bigr{)}.$
(81)

Thus, if (19) holds, from (79) and (81) it follows that

$\displaystyle\left|{\cal R}_{n}\right|$ $\displaystyle\leq$ $\displaystyle
K_{0}\frac{1+\Gamma}{4\gamma}\left(\left|\omega_{n-1}-1\right|+4(1+\Gamma)\left|\omega_{n}-1\right|\right)G_{n}$
(82) $\displaystyle\leq$
$\displaystyle\frac{K_{0}}{\gamma}(1+\Gamma)^{2}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right)G_{n},$

with $K_{0}\leq(3/8)100/7=75/14$.

Figure 1: Graphs of smallest eigenvalue of matrix (80) and $y=0.14x$.

Inserting (82) into (Proof) we finally obtain (70). $\Box$

In the sequel, we define

$\Delta t=\max_{1\leq n\leq N}\Delta t_{n}.$

Also, for $\alpha\in[1,10]$, we consider the polynomial

$\frac{x(1+2x)}{\alpha}x^{2}+\frac{(1+2x)(1+x)}{\alpha}x-17.$ (83)

One can check that it has two complex roots and two real ones, one negative
and one positive. Let $\Gamma^{*}(\alpha)$ be this positive real root. It is
possible to check that $\Gamma^{*}(\alpha)$ is an increasing function of
$\alpha$ whose vales for $\alpha=1,5,10$ are larger than, $1.25$, $2.12$ and
$2.61$, respectively. Fixed $\alpha\in[1,10]$, we restrict ourselves to meshes
satisfying

$\Gamma\leq\Gamma^{*}(\alpha).$ (84)

We also assume

$\sum_{n=2}^{N-1}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right)\leq\Lambda,$
(85)

for some positive $\Lambda$.

We now state and prove the convergence result for the IMEX method (3).

###### Theorem 3

Fix $\kappa>0$, and $\sigma>1$, and let $\Delta t_{n}$ and $h$ satisfy, the
following CFL-type condition

$\displaystyle 3c_{\rm
inv}^{2}\|\boldsymbol{u}\|_{L^{\infty}(L^{\infty})}^{2}(1+2\Gamma)^{2}\frac{(1+\Gamma)}{\gamma}\frac{\Delta
t_{n}}{h^{2}}\leq\frac{\kappa}{T},\quad n\geq 1,$ (86)

and

$\displaystyle 1-\Delta t_{n}\frac{(1+\Gamma)L_{h}^{n}}{0.07\omega_{n}}$
$\displaystyle>$ $\displaystyle\frac{1}{\sigma},\quad n\geq 1.$ (87)

where $L_{h}^{n}$ is defined in (95). Let $G_{n}$ be as defined in (69) for
the error $\boldsymbol{e}_{h}^{n}$ of method (3). Then, the following bound
holds with $m_{k}$ defined in (99)

$\displaystyle G_{n+1}+\nu\sum_{k=2}^{n+1}\Delta t_{k-1}\|\nabla
e_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t_{k}m_{k}\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$ (88)
$\displaystyle C_{n}\biggl{(}(C_{0}^{2}+TC_{3}^{2})(\Delta
t)^{4}+T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}\right)\biggr{)},$

where

$C_{n}=\exp\left(\sigma_{n+1}\Delta
t_{n}\frac{(1+\Gamma)L_{h}^{n}}{0.07\omega_{n}}+\sum_{k=2}^{n}\left(\sigma_{k}\Delta
t_{k-1}\frac{(1+\Gamma)L_{h}^{k-1}}{0.07\omega_{k-1}}+\omega_{k}f_{k}\right)\right),$
(89)

and the constants $f_{k}$ and $\sigma_{k}$ are defined in (97), (102),
respectively, and the constants $C_{i}^{2}$ are defined in (59) and (63).

###### Proof

Taking $\boldsymbol{v}_{h}=\boldsymbol{e}_{h}^{n+1}$ in (26) and applying (70)
from Lemma 7 we reach

$\displaystyle
G_{n+1}-G_{n}+\frac{\omega_{n}}{2(1+\omega_{n})}\left\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\right\|_{0}^{2}+\nu\Delta
t_{n}\|\nabla\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\mu\Delta
t_{n}\|\nabla\cdot\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}\leq\qquad$ (90)
$\displaystyle\Delta
t_{n}|b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h}^{n+1})|+\Delta
t_{n}(\tau_{1}^{n+1}+\tau_{3}^{n+1},\boldsymbol{e}_{h}^{n+1})+\Delta
t_{n}(\tau_{2}^{n+1},\nabla\cdot\boldsymbol{e}_{h}^{n+1})$
$\displaystyle+\frac{K_{0}}{\gamma}(1+\Gamma)^{2}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right)G_{n}.$

As before, we now apply Lemma 3 to estimate the nonlinear term and, further,
we bound

$\displaystyle 2\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}c_{\rm
inv}h^{-1}\|\boldsymbol{e}_{h}^{n+1}-\hat{e}_{h}^{n}\|_{0}\|\boldsymbol{e}_{h}^{n+1}\|_{0}\leq$
$\displaystyle 2c_{\rm
inv}^{2}\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}^{2}\frac{(1+\omega_{n})}{\omega_{n}}\frac{\Delta
t_{n}}{h^{2}}\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$ (91)
$\displaystyle{}+\frac{\omega_{n}}{2(1+\omega_{n})}\frac{\|\boldsymbol{e}_{h}^{n+1}-\hat{e}_{h}^{n}\|_{0}^{2}}{\Delta
t_{n}}.$

Instead of (86) we will now assume

$\displaystyle 2c_{\rm
inv}^{2}\|\hat{\boldsymbol{u}}_{h}^{n}\|_{\infty}^{2}\frac{(1+\omega_{n})}{\omega_{n}}\frac{\Delta
t_{n}}{h^{2}}\leq\frac{\kappa}{T}.$ (92)

As before, we will see at the end of the proof that (92) holds if (86) holds
and $h$ is sufficiently small. Thus, from (3) in Lemma 3, (91) and assuming
(92) we obtain

$\displaystyle\Delta
t_{n}|b(\hat{\boldsymbol{s}}_{h}^{n},\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{e}_{h})-b(\hat{\boldsymbol{u}}_{h}^{n},\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{e}_{h})|\leq\Delta
t_{n}\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9\alpha}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}+\frac{\kappa}{T}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}$
$\displaystyle\quad+\Delta
t_{n}\frac{(1+2\Gamma)(1+\Gamma)}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\Delta
t_{n}\frac{\Gamma(1+2\Gamma)}{2}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
$\displaystyle\quad+\Delta
t_{n}\frac{\mu(1+2\Gamma)(1+\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\Delta
t_{n}\frac{\mu\Gamma(1+2\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
(93)
$\displaystyle\quad+\frac{\omega_{n}}{2(1+\omega_{n})}{\|\boldsymbol{e}_{h}^{n+1}-\hat{\boldsymbol{e}}_{h}^{n}\|_{0}^{2}}.$

Inserting (Proof) into (90) and, arguing as before with the terms involving
$\tau_{1}^{n+1}+\tau_{3}^{n+1}$ and $\tau_{2}^{n+1}$, we get

$\displaystyle G_{n+1}-G_{n}+\nu\Delta
t_{n}\|\nabla\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{17}{18}\mu\Delta
t_{n}\|\nabla\cdot\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}\leq$
$\displaystyle\quad\Delta
t_{n}\left(\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9\alpha}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}+\frac{\kappa+1}{T}\right)\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{(1+2\Gamma)(1+\Gamma)}{2}\Delta
t_{n}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}$
$\displaystyle\quad+\frac{\Gamma(1+2\Gamma)}{2}\Delta
t_{n}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}+\Delta
t_{n}\frac{\mu(1+2\Gamma)(1+\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}$
$\displaystyle\quad+\Delta
t_{n}\frac{\mu\Gamma(1+2\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}+\Delta
t_{n}\frac{T}{2}\|\tau_{1}^{n+1}\|_{0}^{2}$ $\displaystyle\quad+\Delta
t_{n}\frac{T}{2}\|\tau_{3}^{n+1}\|_{0}^{2}+\frac{9\Delta
t_{n}}{2\mu}\|\tau_{2}^{n+1}\|_{0}^{2}+\frac{K_{0}}{\gamma}(1+\Gamma)^{2}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right)G_{n}.$
(94)

Let us denote by

$L_{h}^{n}=\frac{\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}}{2}+\frac{9\alpha}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}+\frac{\kappa+1}{T},$
(95)

so that, using (81), the first term on the right-hand side of (Proof) can be
written as

$\Delta t_{n}L_{h}^{n}\|\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}\leq\Delta
t_{n}\frac{(1+\Gamma)L_{h}^{n}}{0.07\omega_{n}}G_{n+1}.$

Similarly, for the terms involving $\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}$ and
$\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$ on the right-hand side in (Proof), we
have

$\displaystyle\frac{(1+2\Gamma)(1+\Gamma)}{2}\Delta
t_{n}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\frac{\Gamma(1+2\Gamma)}{2}\Delta
t_{n}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}\|\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}\leq$
$\displaystyle\frac{(1+\Gamma)(1+2\Gamma)^{2}}{0.14\omega_{n-1}}\Delta
t_{n}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}G_{n}.$

Thus, going back to (Proof) we reach

$\displaystyle G_{n+1}-G_{n}+\nu\Delta
t_{n}\|\nabla\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}+\frac{17}{18}\mu\Delta
t_{n}\|\nabla\cdot\boldsymbol{e}_{h}^{n+1}\|_{0}^{2}\leq\Delta
t_{n}\frac{(1+\Gamma)L_{h}^{n}}{0.07\omega_{n}}G_{n+1}$
$\displaystyle\quad+\left(\frac{(1+\Gamma)(1+2\Gamma)^{2}}{0.14\omega_{n-1}}\Delta
t_{n}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}+\frac{K_{0}}{\gamma}(1+\Gamma)^{2}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right)\right)G_{n}$
$\displaystyle\quad+\Delta
t_{n}\frac{\mu(1+2\Gamma)(1+\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n}\|_{0}^{2}+\Delta
t_{n}\frac{\mu\Gamma(1+2\Gamma)}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{n-1}\|_{0}^{2}$
(96) $\displaystyle\quad+\frac{\Delta
t_{n}}{2}\left(T(\|\tau_{1}^{n+1}\|_{0}^{2}+\|\tau_{3}^{n+1}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{n+1}\|_{0}^{2}\right).$

Denoting by

$f_{n}=\frac{(1+\Gamma)(1+2\Gamma)^{2}}{0.14\omega_{n-1}}\Delta
t_{n}\|\nabla\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}+\frac{K_{0}}{\gamma}(1+\Gamma)^{2}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right),$
(97)

and adding terms we get

$\displaystyle G_{n+1}+\nu\sum_{k=2}^{n+1}\Delta
t_{k-1}\|\nabla\boldsymbol{e}_{h}^{k}\|_{0}^{2}$
$\displaystyle\quad+\mu\sum_{k=2}^{n+1}\Delta
t_{k}\left(\frac{17}{18}\frac{1}{\omega_{k}}-\frac{(1+2\Gamma)(1+\Gamma)}{18\alpha}-\frac{\Gamma(1+2\Gamma)\omega_{k+1}}{18\alpha}\right)\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$
$\displaystyle G_{1}+\mu\left(\frac{(1+2\Gamma)(1+\Gamma)\Delta
t_{1}}{18\alpha}+\frac{\Gamma(1+2\Gamma)\Delta
t_{2}}{18\alpha}\right)\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2}+\mu\frac{\Gamma(1+2\Gamma)\Delta
t_{1}}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{0}\|_{0}^{2}$
$\displaystyle\quad+\sum_{k=2}^{n+1}\Delta
t_{k-1}\frac{(1+\Gamma)L_{h}^{k-1}}{0.07\omega_{k-1}}G_{k}+\sum_{k=2}^{n+1}f_{k-1}G_{k-1}$
$\displaystyle\quad+\frac{1}{2}\sum_{k=2}^{n+1}\Delta
t_{k-1}\left(T(\|\tau_{1}^{k}\|_{0}^{2}+\|\tau_{3}^{k}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{k}\|_{0}^{2}\right).$
(98)

As mentioned before, for any fixed value of $\alpha$, the number
$\Gamma^{*}(\alpha)$ in (84) is chosen to be the the positive real root of
(83). Since $\omega_{k}\leq\Gamma^{*}$, $k=1,\ldots,N-1$, it follows that

$\displaystyle
m_{k}:=\frac{17}{18}\frac{1}{\omega_{k}}-\frac{(1+2\Gamma)(1+\Gamma)}{18\alpha}-\frac{\Gamma(1+2\Gamma)\omega_{k+1}}{18\alpha}$
$\displaystyle\geq$ $\displaystyle 0.$ (99)

Let us also denote by

$\displaystyle E_{0}^{v}$ $\displaystyle=$
$\displaystyle(1+f_{1})G_{1}+\mu\left(\frac{(1+2\Gamma)(1+\Gamma)\Delta
t_{1}}{18\alpha}+\frac{\Gamma(1+2\Gamma)\Delta
t_{2}}{18\alpha}\right)\|\nabla\cdot\boldsymbol{e}_{h}^{1}\|_{0}^{2}$ (100)
$\displaystyle\quad+\mu\frac{(1+2\Gamma)(1+\Gamma)\Delta
t_{1}}{18\alpha}\|\nabla\cdot\boldsymbol{e}_{h}^{0}\|_{0}^{2}.$

Taking into account (100), (Proof) can be written as

$\displaystyle G_{n+1}+\nu\sum_{k=2}^{n+1}\Delta
t_{k-1}\|\nabla\boldsymbol{e}_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t_{k}m_{k}\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$ (101)
$\displaystyle\quad\quad\Delta
t_{n}\frac{(1+\Gamma)L_{h}^{n}}{0.07\omega_{n}}G_{n+1}+\sum_{k=2}^{n}\left(\Delta
t_{k-1}\frac{(1+\Gamma)L_{h}^{k-1}}{0.07\omega_{k-1}}+\omega_{k}f_{k}\right)G_{k}$
$\displaystyle\quad+\frac{1}{2}\sum_{k=2}^{n+1}\Delta
t_{k-1}\left(T(\|\tau_{1}^{k}\|_{0}^{2}+\|\tau_{3}^{k}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{k}\|_{0}^{2}\right)+E_{0}^{v}.$

Now assumption (87) implies that $\sigma_{n+1}$ defined by

$\sigma_{n+1}=\biggl{(}1-\Delta
t_{n}\frac{(1+\Gamma)L_{h}^{n}}{0.07\omega_{n}}\biggr{)}^{-1},$ (102)

is positive. Thus, applying Lemma 2 we have

$\displaystyle G_{n+1}+\nu\sum_{k=2}^{n+1}\Delta t_{k-1}\|\nabla
e_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t_{k}m_{k}\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$ (103)
$\displaystyle C_{n}\biggl{(}E_{0}^{v}+\frac{1}{2}\sum_{k=2}^{n+1}\Delta
t_{k-1}\left(T(\|\tau_{1}^{k}\|_{0}^{2}+\|\tau_{3}^{k}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{k}\|_{0}^{2}\right)\biggr{)},$

where $C_{n}$ is 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
$\boldsymbol{u}_{h}^{0}=\boldsymbol{s}_{h}^{0}$ so that
$\boldsymbol{e}_{h}^{0}=0$. Consequently,

$G_{1}=\frac{\omega_{0}}{2(1+\omega_{0})}\left\|\boldsymbol{e}_{h}^{1}\right\|_{0}^{2}+\frac{1}{4}\left\|2\boldsymbol{e}_{h}^{1}\right\|_{0}^{2}\leq\frac{3}{2}\left\|\boldsymbol{e}_{h}^{1}\right\|_{0}^{2},$

and, in view of (100),

$E_{0}^{v}\leq(1+f_{1})\frac{3}{2}\left\|\boldsymbol{e}_{h}^{1}\right\|_{0}^{2}+\omega_{0}\left(\frac{(1+2\Gamma)(1+\Gamma)}{18\alpha}+\frac{\Gamma(1+2\Gamma)\omega_{1}}{18\alpha}\right)\mu\Delta
t_{0}\left\|\nabla\cdot\boldsymbol{e}_{h}^{1}\right\|_{0}^{2}.$ (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

$E_{0}^{v}\leq C\left(1+\Delta
t\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}\right)\left(C_{0}^{2}(\Delta
t)^{4}+T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}\right)\right),$ (105)

where $C$ depends on $\Gamma$, $1/\gamma$, $\omega_{0}$ and $\omega_{1}$.

For the truncation errors we also argue as in the fixed stepsize case. Then,
arguing exactly as in (58) we get

$\displaystyle\frac{1}{2}\sum_{k=2}^{n+1}\Delta
t_{k-1}\left(T(\|\tau_{1}^{k}\|_{0}^{2}+\|\tau_{3}^{k}\|_{0}^{2})+\frac{9}{\mu}\|\tau_{2}^{k}\|_{0}^{2}\right)\leq
T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}+C_{3}^{2}(\Delta t)^{4}\right),$
(106)

with the constants $C_{i}^{2}$ defined in (59) and avoiding writing the
explicit dependency on $\omega_{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 $h$
sufficiently small, the CFL type condition (86) implies (92). Arguing as in
the fixed step-size, it is possible to choose $h_{0,1}$ such that
$\left\|\boldsymbol{s}_{h}^{n}-\boldsymbol{u}^{n}\right\|_{\infty}\leq\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}/12$
for $h\leq h_{0,1}$, which implies

$\left\|\boldsymbol{s}_{h}^{n}\right\|_{\infty}\leq\frac{13}{12}\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})},$
(107)

and, since
$\left\|\hat{\boldsymbol{s}}_{h}^{n}\right\|\leq(1+2\Gamma)\left\|\boldsymbol{s}_{h}\right\|_{L^{\infty}(L^{\infty})}$,
by writing
$\hat{\boldsymbol{u}}_{h}^{n}=(\hat{\boldsymbol{u}}_{h}^{n}-\hat{\boldsymbol{s}}_{h}^{n})+\hat{\boldsymbol{s}}_{h}^{n}$,
it also holds

$\left\|\hat{\boldsymbol{u}}_{h}^{n}\right\|_{\infty}\leq c_{\rm
inv}h^{-d/2}\left((1+\Gamma)\left\|\boldsymbol{s}_{h}^{n}-\boldsymbol{u}_{h}^{n}\right\|_{0}+\Gamma\left\|\boldsymbol{s}_{h}^{n-1}-\boldsymbol{u}_{h}^{n-1}\right\|_{0}\right)+(1+2\Gamma)\frac{13}{12}\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}.$
(108)

If an $h_{0,2}>0$ exists such that for $h\leq h_{0,2}$ both the right-hand
side of (62) and the right-hand side of (88) multiplied by
$0.07(1+\Gamma)/\gamma$ are smaller than
$h^{d}\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}/(12c_{\rm
inv})^{2}$ it follows that

$\|\boldsymbol{s}_{h}^{k}-\boldsymbol{u}_{h}^{k}\|_{0}\leq\frac{h^{d/2}}{12c_{\rm
inv}}\|\boldsymbol{u}\|_{L^{\infty}(L^{\infty})},\quad k=1,\cdots,n.$

Thus, in view of (108), it follows that

$\left\|\hat{\boldsymbol{u}}_{h}^{n}\right\|_{\infty}\leq(1+2\Gamma)\frac{7}{6}\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})},$

so that (92) holds as a consequence of assumption (86). We now argue why
$h_{0,2}$ exists, and this will finish the proof. The exponents of the powers
of $h$ on right-hand sides of (62) and (88) are larger than $d$, so that they
will decay with $h$ faster than $h^{d}$ if the exponential term in (89)
remains bounded as $h\rightarrow 0$. To show that this is the case, we first
notice that, in view of the expression of $L_{h}^{n}$ in (95) and recalling
(16) and (107), we have

$L_{h}^{n}\leq(1+2\Gamma)\left(\frac{(1+C_{\infty})}{2}\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}+\frac{9\alpha}{2\mu}\left(\frac{13}{12}\right)^{2}\left\|\boldsymbol{u}\right\|_{L^{\infty}(L^{\infty})}^{2}\right)+\frac{\kappa+1}{T}$

and in view of the expression of $f_{k}$ in (97), we also have

$f_{k}\leq\frac{(1+\Gamma)(1+2\Gamma)^{3}}{0.14\gamma}\Delta
t_{n}(1+C_{\infty})\left\|\boldsymbol{u}\right\|_{L^{\infty}(W^{1,\infty})}+\frac{K_{0}}{\gamma}(1+\Gamma)^{2}\left(\left|\omega_{n-1}-1\right|+\left|\omega_{n}-1\right|\right).$

Finally, due to assumption (87), $\sigma_{n}$ defined in (102) satisfies
$\sigma_{n}<\sigma$ and noticing that the sum of the increments $\Delta t_{n}$
is always bounded by $T$ and recalling assumption (85), it is clear that the
constants $C_{n}$ in (89) remain bounded as $h\rightarrow 0$. $\Box$

Remark 5 In view of (86) we notice that the ratios $\Delta t_{n}/h^{2}$ must
remain bounded. As stated in Remark Proof, the size of the constant $\kappa$
in the CFL condition (86) reflects in the error bound (88) through an
exponential factor. More precisely, considering also the size of the constant
$\sigma$ in (87), taking into account $\sigma_{n}<\sigma$ (with $\sigma_{n}$
defined in (102)) and in view of the expression of $L_{h}^{n}$ in (95) the
influence of the constants $\kappa$ and $\sigma$ affect the constant $C_{n}$
in (103) by a factor of size

$\exp\biggl{(}\kappa\sigma\frac{1+\Gamma}{0.07T}\sum_{k=1}^{n}\frac{\Delta
t_{k}}{\omega_{k}}\biggr{)}=\exp\biggl{(}\kappa\sigma\frac{1+\Gamma}{0.07T}\sum_{k=1}^{n}\Delta
t_{k-1}\biggr{)}\leq\exp\biggl{(}\kappa\sigma\frac{1+\Gamma}{0.07}\biggr{)}.$

On the other hand, the value of $\Gamma^{*}(\alpha)$ in (84) increases with
the value of the free parameter $\alpha$. However, arguing as above, we can
observe that the effect of increasing $\alpha$ (to allow for larger increases
in step length) also affects the constant $C_{n}$ in (89) since the term
$\frac{9\alpha}{2\mu}\|\hat{\boldsymbol{s}}_{h}^{n}\|_{\infty}^{2}$ is present
in the definition of $L_{h}^{n}$ in (95).

###### Theorem 4

Fix $\sigma>1$, and let $\Delta t_{n}$ and $h$ satisfy

$1-\Delta t_{n}\frac{(1+\Gamma){\cal
L}_{h}^{n}}{0.07\omega_{n}}>\frac{1}{\sigma},\quad n\geq 1.$

where

${\cal
L}_{h}^{n}=\frac{\|\nabla\boldsymbol{s}_{h}^{n+1}\|_{\infty}}{2}+\frac{9\alpha}{2\mu}\|\boldsymbol{s}_{h}^{n+1}\|_{\infty}^{2}+\frac{1}{T},$

Let $G_{n}$ be as defined in (69) for the error $\boldsymbol{e}_{h}^{n}$ of
the semi-implicit method (3). Then, the following bound holds with $m_{k}$
defined in (99)

$\displaystyle G_{n+1}+\nu\sum_{k=2}^{n+1}\Delta t_{k-1}\|\nabla
e_{h}^{k}\|_{0}^{2}+\mu\sum_{k=2}^{n+1}\Delta
t_{k}m_{k}\|\nabla\cdot\boldsymbol{e}_{h}^{k}\|_{0}^{2}\leq$ (109)
$\displaystyle C_{n}\biggl{(}(C_{0}^{2}+TC_{3}^{2})(\Delta
t)^{4}+T\left(C_{1}^{2}h^{2k}+C_{2}^{2}h^{2l+2}\right)\biggr{)},$

where

$C_{n}=\exp\left(\sigma_{n+1}\Delta t_{n}\frac{(1+\Gamma){\cal
L}_{h}^{n}}{0.07\omega_{n}}+\sum_{k=2}^{n}\left(\sigma_{k}\Delta
t_{k-1}\frac{(1+\Gamma){\cal
L}_{h}^{k-1}}{0.07\omega_{k-1}}+\omega_{k}f_{k}\right)\right),$

and the constants $f_{k}$ and $\sigma_{k}$ are defined in (97), (102),
respectively, and the constants $C_{i}^{2}$ 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
(3) with the first term on the right-hand side of (90) replaced by

$|b(\hat{\boldsymbol{s}}_{h}^{n},\boldsymbol{s}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})-b(\hat{\boldsymbol{u}}_{h}^{n},\boldsymbol{u}_{h}^{n+1},\boldsymbol{e}_{h}^{n+1})|$

and $\tau_{3}$ by $\tau_{4}$. Then, recalling Lemma 4 and Remark Proof, we
have that (Proof) holds with $\kappa=0$, $\hat{\boldsymbol{s}}_{h}^{n}$
replaced by $\boldsymbol{s}_{h}^{n+1}$ and $\tau_{3}$ by $\tau_{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. $\Box$

Remark 6 Analogous comments to those in Remark Proof apply in this case taking
into account (88), (109) and the definition of $G_{n}$ in (69).

##  4 Numerical Experiments

We test methods (3) and (3) 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 treatment of the
convective term as specified in (3). In particular we notice that the local
error estimate of the method of order $k$ is given by

${\rm EST}_{n}=\frac{\Delta t_{n}}{t_{n+1}-t_{n-k}}\left\|{\cal
U}^{n,k}_{h}\right\|_{0},$

where

${\cal
U}^{n,k}_{h}=\prod_{i=0}^{k-1}(t_{n+1}-t_{n-i})\boldsymbol{u}_{h}[t_{n+1},\ldots,t_{n-k}],$

where $\boldsymbol{u}_{h}[t_{n+1},\ldots,t_{n-k}]$ is the standard $(k+1)$-th
divided difference based on $t_{n+1},\ldots,t_{n-k}$ and
$\boldsymbol{u}_{h}^{n+1},\ldots,\boldsymbol{u}_{h}^{n-k}$. This corresponds
to the fully nonlinear BDF, but, nevertheless, it was used with the methods
(3) and (3). At each time level $t_{n}$, the estimation ${\rm EST}_{n}$ is
compared with the quantity

${\rm
TOL}_{n}=TOL_{r}\left(\max\left\\{\left\|\boldsymbol{u}_{h}^{n+1}\right\|_{0},\left\|\boldsymbol{u}_{h}^{n}\right\|_{0}\right\\}+0.001\right),$

where $TOL_{r}$ is a given tolerance. If ${\rm EST_{n}}>{\rm TOL_{n}}$,
$\boldsymbol{u}_{h}^{n+1}$ is not considered sufficiently accurate and is
rejected. It is then recomputed from $t_{n}$ with a new step length given by

$\Delta t_{n}^{\rm new}=0.9\Delta t_{n}\sqrt[k+1]{\frac{{\rm TOL}_{n}}{{\rm
EST}_{n}}}.$ (110)

On the contrary, if ${\rm EST_{n}}<{\rm TOL_{n}}$, then,
$\boldsymbol{u}_{h}^{n+1}$ is accepted and the algorithm proceeds to compute
$\boldsymbol{u}_{h}^{n+2}$ with a step length $\Delta t_{n+1}$ given by the
right-hand side of (110).

The algorithm starts with $k=1$ and $\Delta t_{0}=\Delta
t_{1}=\sqrt{TOL_{r}}/100$. At $t_{1}$ the first local error estimation ${\rm
EST}_{1}$ is computed. If ${\rm EST}_{1}>{\rm TOL}_{1}$ the computation is
restarted again from $t_{0}$ with the step size changed according to (110).
From $t_{2}$ onwards, estimates corresponding to the two-step BDF can be
estimated. When these estimates are smaller than those corresponding to $k=1$,
then $k$ is 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 $\boldsymbol{y}^{n}$ the vector containing the
coefficients of the velocity $\boldsymbol{u}_{h}^{n+1}$ and the pressure
$p_{h}^{n+1}$ in the corresponding nodal basis of the finite element spaces,
then, since (3) and (3) are linear in $\boldsymbol{u}_{h}^{n+1}$ and
$p_{h}^{n+1}$, to find $\boldsymbol{y}^{n}$ one has to solve a linear system

$A_{n}\boldsymbol{y}^{n}=\boldsymbol{b}^{n}.$ (111)

This was solved by iterative refinement. Starting with an approximation
$\boldsymbol{y}^{n,[0]}=\hat{\boldsymbol{y}}^{n}$, extrapolated from the
previous $k$ values ($k=1,2$), the approximation is updated as
$\boldsymbol{y}^{n,[j+1]}=\boldsymbol{y}^{n,[j]}+\boldsymbol{e}^{n,[j]}$,
where $\boldsymbol{e}^{n,[j]}$ is obtained by solving the system
$A_{m}\boldsymbol{e}^{n,[j]}=\boldsymbol{b}^{n}-A_{n}\boldsymbol{y}^{n,[j]}$,
where, we notice, $A_{m}$ has been factored on time level $t_{m}\leq t_{n}$.
The iteration is stopped either when the velocity in two consecutive
iterations satisfies

$\left\|\boldsymbol{u}_{h}^{n,[j+1]}-\boldsymbol{u}_{h}^{n,[j]}\right\|_{0}\leq\min\left\\{10^{-8},\frac{TOLr}{100}\right\\}\left(\left\|\boldsymbol{u}_{h}^{n-1}\right\|_{0}+0.001\right),$
(112)

in which case we set $\boldsymbol{u}_{n}^{n}=\boldsymbol{u}_{n}^{n,[j+1]}$, or
if five iterations are performed without ever (112) being satisfied, in which
case, $A_{n}$ is factored , replaces $A_{m}$ and 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 (3) when $\boldsymbol{v}_{h}=0$ was rewriten as

$-(\nabla\cdot\boldsymbol{u}_{h}^{n+1},q_{h})=0,$

so that matrices $A_{n}$ were symmetric (although indefinite).

###  4.1 Problem with known solution

We consider the Navier-Stokes equations in the domain $\Omega=[0,1]^{2}$ with
$T=4$, and the forcing term $\boldsymbol{f}$ chosen so that the solution
$\boldsymbol{u}$ and $p$ are given by

$\displaystyle\boldsymbol{u}(x,y,t)$ $\displaystyle=$
$\displaystyle\frac{6+4\cos(4t)}{10}\left[\begin{array}[]{c}8\sin^{2}(\pi
x)(2y(1-y)(1-2y)\\\ -8\pi\sin(2*\pi x)(y(1-y))^{2}\end{array}\right]$ (115)
$\displaystyle p(x,y,t)$ $\displaystyle=$
$\displaystyle\frac{6+4\cos(4t)}{10}\sin(\pi x)\cos(\pi y).$ (116)

We present results with the $P_{2}/P_{1}$ pair 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^{-6}$
and $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 $\mu=0.05$.

In Fig. 2, for the final time $t=4$, we show the $L^{2}$ errors in velocity

$\|\boldsymbol{u}_{h}^{n}-I_{h}(\boldsymbol{u}^{n})\|_{0},$ (117)

where $I_{h}$ is the standard Lagrange interpolant. In Fig. 2, as in the rest
of the pictures in this section, results of IMEX method (3) are joined by
continuous lines, while those of the semi-implict method (3) 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.

Figure 2: Velocity errors (117) for $T=4$ and $\mu=0.05$. Continuous lines
correspond to method ((3) and continuous lines to method (3).

It can be seen that, as $h$ decreases, the errors also decrease, but with a
higher rate for the two largest values of the viscosity, $\nu=10^{-2}$ and
$\nu=10^{-4}$. The results corresponding to $\nu=10^{-6}$, $\nu=10^{-8}$ and
$\nu=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 $\nu=10^{-2}$ (in blue)
and $\nu=10^{-10}$, confirming the analysis that the rate of decay is
$O(h^{2})$ independent of $\nu$ in the convection dominated regime, i.e., for
$\nu$ small enough. For $\nu=10^{-2}$, and the last three results of
$\nu=10^{-4}$, the slope close to 4 reflects that the numerical approximation
$\boldsymbol{u}_{h}$ is supraconvergent, in the sense that the error (117)
decays faster that $\|\boldsymbol{u}^{n}-I_{h}(\boldsymbol{u}^{n})\|_{0}$,
which is $O(h^{3})$.

We check that our implementation uses variable step sizes in Fig. 3, where we
show the step lengths of method (3) (the ones corresponding to method (3)
being hardly distinguishable). They correspond to $\nu=10^{-6}$,
$TOL_{r}=10^{-5}$ and the mesh with $N=12$ subdivisions in each coordinate
directions. After the initial steps corresponding to the first order method,
the step lengths grow quickly, and then vary between $0.001$ and $0.005$
approximately. In Fig. 3 we marked with a red cross the locations where a step
was rejected. We see that the number of rejections is really low. Most of the
results corresponding to different values of $\nu$, $TOL_{r}$ and $N$ showed
similar variations.

Figure 3: Step lengths corresponding to $\nu=10^{-6}$, $TOL_{r}=10^{-5}$ and
$N=12$. Rejections marked with a red cross.

In Fig. 4 we compare the number of factorizations with the number of steps for
the results shown in Fig. 2 (recall each factorization is reused in subsequent
steps for iterative refinement until failure to converge, where it is
updated). We can see that the number of factorizations hardly grows as $h$
(and the corresponding tolerances) are decreased, and that at worst (large
$h$, small $\nu$) they are hardly more than 5% of the number of steps. Some
differences can be seen between methods (3) and (3), for large values of
$\nu$. We do not have an explanation for them. We notice, however, that the
differences between the methods become less significant for smaller values of
$\nu$.

Figure 4: Number of factorizations (left) and steps (right) for results shown
in Fig. 2.

Next we check if the CFL-type condition (86) is necessary in practice for the
IMEX method (3), whereas this is not the case of the semi-implicit method (3).
For that purpose, we repeat the above computations (and including also meshes
with $N=96$ subdivisions in each coordinate direction) but with tolerance
$TOL_{r}=10^{-4}$ for all meshes. Notice that in principle the number of steps
depends on the tolerance, so, if the tolerance is kept constant while the
spatial mesh size $h$ is decreased, the number of steps should not change
much. On the other hand, if stability worsens as $h$ is decreased, errors
arising from time discretization increase with every step, and will become
larger than the tolerance more easily, forcing the algorithm to take smaller
step sizes, this resulting in a larger number of steps. In Fig. 5 we show the
number of steps (left plot). For the IMEX method (continuous lines), we see
that although they remain fairly constant for $\nu=10^{-2}$, this is not the
case of the rest of the values of the viscosity. Indeed, for the three
smallest values of the viscosity, the number of steps roughly doubles from
$h=1/24$ to $h=1/48$, and it is multiplied by three from $h=1/48$ to $h=1/96$,
approaching the factor $4$ predicted by the CFL-type condition (86). On the
other hand, for the semi-implicit method (discontinuous lines) we see that,
contrary to the IMEX method, the number of steps fairly constant as $h$ is
refined. As we have shown, this method does not require the CFL-type condition
(86) to converge. On the contrary, in view of Fig. 5, we may conclude that
method with explicit convection (3) indeed suffers from restriction (86).

Figure 5: Number of steps (left) and factorizations (right) for results for
$TOL=10^{-4}$ for all meshes. Color and symbol codes on the right plot are
those of the left one. Discontinuous lines correspond to the semi-implicit
method (3) and continuous ones for the IMEX method (3).

In Fig 5 we can also see that the number of factorizations remains fairly
independent of the number of steps, as it was the case in the previous example
where tolerances decreased with $h$.

###  4.2 Flow past a cylinder

We now present results for the well-known benchmark problem defined in [19].
The domain is given by

$\Omega=(0,2.2)\times(0,0.41)/\left\\{(x,y)\mid(x-0.2)^{2}+(y-0.2)^{2}\leq
0.0025\right\\}$

and the time interval is $[0,8]$. The velocity is identical on both vertical
sides, and is given by

$\boldsymbol{u}(0,y)=\boldsymbol{u}(2.2,y)=\frac{6}{0.41^{2}}\sin\Bigl{(}\frac{\pi
t}{8}\Bigr{)}\left(\begin{array}[]{c}y(0.41-y)\\\ 0\end{array}\right).$

In the rest of the boundary the velocity is set $\boldsymbol{u}={\bf 0}$. At
$t=0$, the initial velocity is ${\boldsymbol{u}}={\bf 0}$. The viscosity is
set to $\nu=10^{-3}$ and the forcing term is $\boldsymbol{f}={\bf 0}$.

It is well-known that around $t=4$ a vortex sheet develops behind the
cylinder, as it can be seen in Fig. 6 where we show the speed and velocity
fields for $t=5$ to $t=8$. We use the same scale in the four plots of the
velocity fields, and the results we obtain are virtually identical to those in
[16, Fig. 2].

Figure 6: Speed contours (left) and velocity field for times $t=5,6,7,8$ (from
top to bottom).

In Fig. 7 we show the mesh used in all the results below, which has 6624
elements whose diameters range from $5.53\times 10^{-3}$ to $3.38\times
10^{-2}$. The number of degrees of freedom for the velocity on this mesh is
27168, and for for the pressure, 3480. Following experience in [8], the grad-
div parameter for the experiments in this section is set to $\mu=0.01$.

In Fig 8 we show the step sizes used with tolerances $10^{-4}$, $10^{-6}$, and
$10^{-8}$. We notice that, for each tolerance, after the very small steps near
$t=0$, step sizes vary by a factor of at least ten, so that, as in the
previous section, the algorithm uses variable steps sizes. It can also be seen
that the smaller the tolerance, the smaller the step sizes, except for the
IMEX method (3) with tolerances $10^{-4}$ and $10^{-6}$ on the interval
$[1.7,3.7]$ where step sizes are the same for both tolerances and clearly
smaller than those of the semi-implicit method (3) for the same tolerances.
This suggest that, for these two tolerances, the step size in method (3) is
limited by the CFL-type condition (86) rather than by the roughness of the
solution, suggesting, once again, that restriction (86) is not a consequence
of the techniques used in the proof of Theorem 3.

Figure 7: The mesh used.

Figure 8: Flow around a cylinder: step lengths for $TOLr=10^{-4}$, $10^{-6}$,
and $10^{-8}$. Top: IMEX method (3) bottom; semi-implicit method (3).

Figure 9: Flow around a cylinder: drag and lift coefficients for
$TOL_{r}=10^{-4}$, and $10^{-6}$. Reference values for the maximum of the lift
and drag coefficients marked with an asterisk.

From the practical point of view, we may mention that the number of steps for
the three tolerances, $10^{-4}$, $10^{-6}$, and $10^{-8}$, were $11221$,
$21187$ and $93671$, respectively, for the IMEX method, and 4384, 19885 and
92086 for the semi-implicit method, but the number of factorizations were,
respectively, $67$, $53$, and $92$ for the IMEX method and, 72, 84 and 89 for
the semi-implicit method, low figures for both methods, in line with those of
the previous section.

Quantities of practical interest in this example are the drag and lift
coefficients, $c_{d}$ and $c_{l}$, respectively. In Fig. 9, for tolerances
$10^{-4}$ and $10^{-6}$, we show the evolution of the drag and lift
coefficients obtained with the IMEX method, which were computed with the
following formulae, [16], [17]

$\displaystyle c_{d}(t_{n})$ $\displaystyle=-20\left(\left(\frac{1}{\Delta
t_{n-1}}{\cal
D}\boldsymbol{u}_{h}^{n},\boldsymbol{v}_{d}\right)+\nu(\nabla\boldsymbol{u}_{h}^{n},\nabla\boldsymbol{v}_{d})+b(\boldsymbol{u}_{h}^{n},\boldsymbol{u}_{h}^{n},\boldsymbol{v}_{d})-(p_{h}^{n}(t),\nabla\cdot\boldsymbol{v}_{d})\right),$
$\displaystyle c_{l}(t_{n})$ $\displaystyle=-20\left(\left(\frac{1}{\Delta
t_{n-1}}{\cal
D}\boldsymbol{u}_{h}^{n},\boldsymbol{v}_{l}\right)+\nu(\nabla\boldsymbol{u}_{h}^{n},\nabla\boldsymbol{v}_{l})+b(\boldsymbol{u}_{h}^{n},\boldsymbol{u}_{h}^{n},\boldsymbol{v}_{l})-(p_{h}^{n}(t),\nabla\cdot\boldsymbol{v}_{l})\right),$

where $\boldsymbol{v}_{d}$ and $\boldsymbol{v}_{l}$ are piecewise quadratic
functions taking values $\boldsymbol{v}_{d}=[1,0]^{T}$ and
$\boldsymbol{v}_{l}=[0,1]^{T}$ on those nodes on the circumference
$c\equiv(x-0.2)^{2}+(y-0.2)^{2}=0.0025$, and vanishing on triangles without
vertices on $c$.

The results of both tolerances are indistinguishable, suggesting sufficient
accuracy already for $TOL_{r}=10^{-4}$. We also show the reference values of
the maximun values, $c_{d,\text{max}}^{\text{ref}}=2.950921575$ and
$c_{l,\text{max}}^{\text{ref}}=0.47795$, plotted with a red asterisk at their
corresponding times $t_{d,\text{max}}^{\text{ref}}=3.93625$, and
$t_{l,\text{max}}^{\text{ref}}=5.693125$ (all taken from [16]) showing a very
good coincidence. In fact, the relative errors for the times where the maximum
values are reached are in both cases below $0.0006$, and for the drag and lift
coefficients are $0.012$ and $0.043$, respectively. Furthermore, there is good
agreement between the plots in Fig. 9 and those in [16, Fig. 4]. Finally, also
of interest is the pressure difference between the front and the back of the
cylinder, whose reference value for $t=8$ is, according to [16], $\Delta
p_{8}=-0.1116$. For the tolerance $10^{-4}$, the relative error in $\Delta
p_{8}$ is below $0.002$.

The results of these two sections suggest that, as predicted by the theory,
the IMEX method (3) is stable and convergent, but indeed suffers from the CFL-
type restriction (86). On the other hand, the computational cost of the semi-
implicit method (3) is not much larger than the IMEX one, see Figure 5, and
has the advantage of having no step restriction.

##  5 Conclusions

The numerical simulation of the Navier-Stokes equations is still a challenge
in which there are several questions that deserve some research. Important
aspects of the present paper are the following:

  * •

We study the possibility of optimizing the cost of the temporal integration.
To this end we consider two IMEX schemes for the incompressible Navier-Stokes
equations, which reduce the cost of every time step, compared with an implicit
method. In the first method the nonlinear term is explicit and in the second
one is semi-implicit.

  * •

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
$(\Delta t)h^{-1}$ has to be bounded (we required $(\Delta t)h^{-2}$ to be
bounded). Our numerical experiments confirm experimentally that the CFL-type
restriction cannot be weakened.

  * •

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 we did not found in the literature any reference with
the error analysis for the Navier-Stokes equations.

  * •

Methods for whom robust estimates can be derived enable stable flow
simulations for small viscosity coefficients on comparatively coarse grids.
Adding grad-div stabilization to the standard mixed finite element formulation
we are able to prove error bounds with constants independent of the Reynolds
number for the two methods we analized.

  * •

In summary, to our knowledge, this is the first time in the literature
variable step BDF2 methods are analyzed for the Navier-Stokes equations and
also the first time in the literature robust bounds (independent of viscosity
coefficient) are proved for Navier-Stokes with fully discrete methods using
variable time-step.

## 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 grad-div 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.
  * [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., to appear.
  * [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. Schäfer, 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.