RESOLUTION OF THE IMPLICIT EULER SCHEME FOR THE
NAVIER-STOKES EQUATION THROUGH A LEAST-SQUARES METHOD
J
´
ER
ˆ
OME LEMOINE AND ARNAUD M
¨
UNCH
Abstract. This work analyzes a least-squares method in order to solve implicit time schemes
associated to the 2D and 3D Navier-Stokes system, introduced in 1979 by Bristeau, Glowinksi,
Periaux, Perrier and Pironneau. Implicit time schemes reduce the numerical resolution of the
Navier-Stokes system to multiple resolutions of steady Navier-Stokes equations. We first con-
struct a minimizing sequence (by a gradient type method) for the least-squares functional
which converges strongly and quadratically toward a solution of a steady Navier-Stokes equa-
tion from any initial guess. The method turns out to be related to the globally convergent
damped Newton approach applied to the Navier-Stokes operator, in contrast to standard
Newton method used to solve the weak formulation of the Navier-Stokes system. Then, we
apply iteratively the analysis on the fully implicit Euler scheme and show the convergence of
the method uniformly with respect to the time discretization. Numerical experiments for 2D
examples support our analysis.
Key Words. Steady Navier-Stokes equation, Implicit time scheme, Least-squares approach,
Damped Newton method.
AMS subject classifications. 35Q30, 93E24.
1. Introduction - Motivation
Let Ω ⊂ R
d
, d = 2 or d = 3 be a bounded connected open set whose boundary ∂Ω is Lipschitz
and T > 0. We denote by V = {v ∈ D(Ω)
d
, ∇·v = 0}, H the closure of V in L
2
(Ω)
d
and V the
closure of V in H
1
(Ω)
d
, endowed with the norm kvk
V
= k∇vk
2
.
The Navier-Stokes system describes a viscous incompressible fluid flow in the bounded domain
Ω during the time interval (0, T ) submitted to the external force F . It reads as follows :
(1.1)
u
t
− ν∆u + (u · ∇)u + ∇p = F, ∇ · u = 0 in Ω × (0, T ),
u = 0 on ∂Ω × (0, T ),
u(·, 0) = u
0
in Ω,
where u is the velocity of the fluid, p its pressure and ν is the viscosity constant assumed smaller
than one. We refer to [19]. This work is concerned with the numerical approximation of (1.1)
through the time marching fully implicit Euler scheme
(1.2)
y
0
= u
0
in Ω,
y
n+1
− y
n
δt
− ν∆y
n+1
+ (y
n+1
· ∇)y
n+1
+ ∇π
n+1
=
1
δt
Z
t
n+1
t
n
F (·, s)ds, n ≥ 0,
∇ · y
n+1
= 0 in Ω, n ≥ 0,
y
n+1
= 0 on ∂Ω, n ≥ 0,
Date: July 8, 2020.
Laboratoire de Math´ematiques, Universit´e Clermont Auvergne, UMR CNRS 6620, Campus des C´ezeaux,
63177 Aubi`ere, France. e-mail: jerome.lemoine@uca.fr.
Laboratoire de Math´ematiques, Universit´e Clermont Auvergne, UMR CNRS 6620, Campus des C´ezeaux,
63177 Aubi`ere, France. e-mail: arnaud.munch@uca.fr.
1
2 J
´
ER
ˆ
OME LEMOINE AND ARNAUD M
¨
UNCH
where {t
n
}
n=0...N
, for a given N ∈ N, is a uniform discretization of the time interval (0, T ).
δt = T/N is the time discretization step. This also-called backward Euler scheme is studied for
instance in [19, chapter 3, section 4]. It is proved there that the piecewise linear interpolation
(in time) of {y
n
}
n∈[0,N]
weakly converges in L
2
(0, T, V ) toward a solution u of (1.1) as δt goes
to zero. It achieves a first order convergence with respect to δt. We also refer to [20] for a
stability analysis of the scheme in long time. We refer to [17] and the references therein for
Crank-Nicolson schemes achieving second order convergence.
The determination of y
n+1
from y
n
requires the resolution of a nonlinear partial differential
equation. Precisely y
n+1
together with the pressure π
n+1
solve the following problem: find y ∈ V
and π ∈ L
2
0
(Ω), solution of
(1.3)
(
α y − ν∆y + (y · ∇)y + ∇π = f + α g, ∇ · y = 0 in Ω,
y = 0 on ∂Ω,
with
(1.4) α :=
1
δt
> 0, f :=
1
δt
Z
t
n+1
t
n
F (·, s)ds, g = y
n
.
Recall that for any f ∈ H
−1
(Ω)
d
and g ∈ L
2
(Ω)
d
, there exists one solution (y, π) ∈ V ×L
2
0
(Ω) of
(1.3), unique if kgk
2
L
2
(Ω)
d
+ α
−1
ν
−1
kfk
2
H
−1
(Ω)
d
is small enough (see Proposition 2.3 for a precise
statement). L
2
0
(Ω) stands for the space of functions in L
2
(Ω)
d
with zero means.
A weak solution y ∈ V of (1.3) solves the formulation F (y, w) = 0 for all w ∈ V where F is
defined by
(1.5)
F (y, z) :=
Z
Ω
α y · z + ν∇y · ∇z + (y · ∇)y · z
− < f, z >
H
−1
(Ω)
d
×H
1
0
(Ω)
d
−α
Z
Ω
g · z = 0, ∀z ∈ V .
If D
y
F is invertible, one may approximate a weak solution through the iterative Newton method:
a sequence {y
k
}
k
∈ V is constructed as follows
(1.6)
(
y
0
∈ V ,
D
y
F (y
k
, w) · (y
k+1
− y
k
) = −F (y
k
, w), ∀w ∈ V , k ≥ 0.
If the initial guess y
0
is close enough to a weak solution of (1.3), i.e. a solution satisfying
F (y, w) = 0 for all w, then the sequence {y
k
}
k
converges. We refer to [15, Section 10.3] and [4,
Chapter 6]) and for some numerical aspects to [10].
Alternatively, we may also employ least-squares methods which consists in minimizing a qua-
dratic functional, which measure how an element y is close to the solution. For instance, we may
introduce the extremal problem : inf
y∈V
E(y) with E : V → R
+
defined by
(1.7) E(y) :=
1
2
Z
Ω
α|v|
2
+ ν|∇v|
2
where the corrector v, together with the pressure, is the unique solution in V × L
2
0
(Ω) of the
linear boundary value problem:
(1.8)
(
αv − ν∆v + ∇π +
αy − ν∆y + (y · ∇)y − f − αg
= 0, ∇ · v = 0 in Ω,
v = 0 on ∂Ω.
Remark that E(y) vanishes if and only if y ∈ V is a weak solution of (1.3), equivalently a zero of
F (y, w) = 0 for all w ∈ V . As a matter of fact, the infimum is reached. Least-squares methods
to solve nonlinear boundary value problems have been the subject of intensive developments
in the last decades, as they present several advantages, notably on computational and stability
viewpoints. We refer to the books [1, 7]. The minimization of the functional E over V leads
WEAK LEAST-SQUARES METHODS FOR NAVIER-STOKES 3
to a so-called weak least squares method. Precisely, the equality
p
2E(y) = sup
w∈V ,w6=0
F (y,w)
kwk
V
shows that E is equivalent to the V
0
norm of the Navier-Stokes equation (we refer to Remark
2.9). The terminology “H
−1
-least-squares method” is employed in [2] where the minimization
of E has been introduced and numerically implemented to approximate the solutions of (1.1)
through the scheme (1.2). We also mention [4, Chapter 4, Section 6] which studied later the use
of a least-squares strategy to solve a steady Navier-Stokes equation without incompressibility
constraint.
The first objective of the present work is to analyze rigorously the method introduced in [2]
and show that one may construct minimizing sequences in V for E that converge strongly toward
a solution of (1.3). The second objective is to justify the use of that weak least-squares method
to solve iteratively a weak formulation of the scheme (1.2), leading to an approximation of the
solution of (1.1). This requires to show some convergence properties of a minimizing sequence
for E, uniformly with respect to the parameter n related to the time discretization. As we shall
see, this requires smallness assumptions on the data u
0
and F .
The rest of the paper is organized as follows. Section 2 is devoted to the analyze the least-
squares method (1.7)-1.8 associated to weak solutions of (1.3). We show that E is differentiable
over V and that any critical point for E in the ball B defined in Definition 2.1 is also a zero of E.
This is done by introducing a descent direction Y
1
for E at any point y ∈ V for which E
0
(y) ·Y
1
is proportional to E(y). Then, assuming that there exists at least one solution of (1.3) in the ball
B, we show that any minimizing sequence {y
k
}
(k∈N)
for E uniformly in B strongly converges to
a solution of (2.5). Such limit belongs to B and is actually the unique solution. Eventually, we
construct a minimizing sequence (defined in (2.22)) based on the element Y
1
and initialized with
g assume in V . We show that, if α is large enough, then this particular sequence is uniformly in
B and converges (quadratically after a finite number of iterates related to the values of ν and α)
strongly to the solution of (1.3) (see Theorem 2.15). This specific sequence coincides with the
one obtained from the damped Newton method, a globally convergent generalization of (1.6).
Then, in Section (3), as an application, we consider the least-squares approach to solve it-
eratively the backward Euler scheme (see (3.1)), weak formulation of (1.2). For each n > 0,
in order to approximate y
n+1
, we define a minimizing sequence {y
n+1
k
}
k≥0
based on Y
n+1
1
and
initialized with y
n
. Adapting the global convergence result of Section 2, we then show, as-
suming ku
0
k
L
2
(Ω)
d
+ kF k
L
2
(0,T ;H
−1
(Ω)
d
)
small enough, the strong convergence of the minimizing
sequences, uniformly with respect to the time discretization parameter n (see Theorem 3.7). The
analysis is performed for d = 2 for both weak and regular solutions and for d = 3 for regular
solutions. Our analysis justifies the use of Newton type methods to solve implicit time schemes
for (1.1), as mentioned in [15, Section 10.3]. To the best of our knowledge, such analysis of
convergence is original.
In Section 4, we discuss numerical experiments based on finite element approximations in
space for two 2D geometries: the celebrated example of the channel with a backward facing step
and the semi-circular driven cavity introduced in [5]. We notably exhibit for small values of the
viscosity constant the robustness of the damped Newton method (compared to the Newton one).
2. Analysis of a Least-squares method for a steady Navier-Stokes equation
We analyse in this section a least-squares method to solve the steady Navier-Stokes equation
(1.3) assuming α > 0: we extend [11] where the particular case α = 0 is addressed.
2.1. Technical preliminary results. We endow the space V with the norm kyk
V
:= k∇yk
2
,
for all y ∈ V . k·k
2
stands for the L
2
norm k·k
L
2
(Ω)
d
. We shall also use the following notations
|||y|||
2
V
:= αkyk
2
2
+ νk∇yk
2
2
, ∀y ∈ V
and < y, z >
V
:= α
R
Ω
y · z + ν
R
Ω
∇y · ∇z so that < y, z >
V
≤ |||y|||
V
|||z|||
V
for any y, z ∈ V .
4 J
´
ER
ˆ
OME LEMOINE AND ARNAUD M
¨
UNCH
In the sequel, we repeatedly use the following classical estimates (see [19]).
Lemma 2.1. Let any u, v ∈ V . If d = 2, then
(2.1) −
Z
Ω
u · ∇u · v =
Z
Ω
u · ∇v · u ≤
√
2kuk
2
k∇vk
2
k∇uk
2
.
If d = 3, then there exists a constant c = c(Ω) such that
(2.2)
Z
Ω
u · ∇v · u ≤ ckuk
1/2
2
k∇vk
2
k∇uk
3/2
2
.
Definition 2.1. For any y ∈ V , we define
τ
d
(y) :=
kyk
V
√
2αν
, if d = 2,
Mkyk
V
(αν
3
)
1/4
, if d = 3
with M :=
3
3/4
4
c and c from (2.2).
We shall also repeatedly use the following Young type inequalities.
Lemma 2.2. For any u, v ∈ V , the following inequalities hold true :
(2.3)
√
2kuk
2
k∇vk
2
k∇uk
2
≤ τ
2
(v)|||u|||
2
V
if d = 2 and
(2.4) ckuk
1/2
2
k∇vk
2
k∇uk
3/2
2
≤ τ
3
(v)|||u|||
2
V
if d = 3.
Let f ∈ H
−1
(Ω)
d
, g ∈ L
2
(Ω)
d
and α ∈ R
?
+
. The weak formulation of (1.3) reads as follows:
find y ∈ V solution of
(2.5) α
Z
Ω
y ·w + ν
Z
Ω
∇y ·∇w +
Z
Ω
y ·∇y ·w =< f, w >
H
−1
(Ω)
d
×H
1
0
(Ω)
d
+α
Z
Ω
g ·w, ∀w ∈ V .
The following result holds true (we refer to [13]).
Proposition 2.3. Assume Ω ⊂ R
d
is bounded and Lipschitz. There exists at least one solution
y of (2.5) satisfying
(2.6) |||y|||
2
V
≤ c
0
νkfk
2
H
−1
(Ω)
d
+ αkgk
2
2
where c
0
> 0, only connected to the Poincar´e constant, depends on Ω. If moreover Ω is C
2
and
f ∈ L
2
(Ω)
d
, then any solution y ∈ V of (2.5) belongs to H
2
(Ω)
d
.
Lemma 2.4. Assume that a solution y ∈ V of (2.5) satisfies τ
d
(y) < 1. Then, such solution is
the unique solution of (2.5).
Proof. Let y
1
∈ V and y
2
∈ V be two solutions of (2.5). Set Y = y
1
− y
2
. Then,
α
Z
Ω
Y · w + ν
Z
Ω
∇Y · ∇w +
Z
Ω
y
2
· ∇Y · w +
Z
Ω
Y · ∇y
1
· w = 0 ∀w ∈ V .
We now take w = Y and use that
R
Ω
y
2
· ∇Y · Y = 0. If d = 2, we use (2.1) and (2.3) to get
|||Y |||
2
V
= −
Z
Ω
Y · ∇y
1
· Y ≤ τ
2
(y
1
)|||Y |||
2
V
leading to (1−τ
2
(y
1
))|||Y |||
2
V
≤ 0. Consequently, if τ
2
(y
1
) < 1 then Y = 0 and the solution of (2.5)
is unique. In particular, in view of (2.6), this holds if the data satisfy νkgk
2
2
+
c
0
α
kfk
2
H
−1
(Ω)
d
< 2ν
3
.
WEAK LEAST-SQUARES METHODS FOR NAVIER-STOKES 5
If d = 3, we use (2.2) and (2.4) to obtain
|||Y |||
2
V
= −
Z
Ω
Y · ∇y
1
· Y ≤ ckY k
1
2
2
k∇Y k
3
2
2
k∇y
1
k
2
≤ τ
3
(y
1
)|||Y |||
2
V
leading to
1−τ
3
(y
1
)
kY k
2
2
≤ 0 and to the uniqueness of the solution if τ
3
(y
1
) < 1. In particular,
in view of (2.6), this holds if the data satisfy νkgk
2
2
+
c
0
α
kfk
2
H
−1
(Ω)
d
< M
−2
ν
7/2
α
−1/2
.
We now introduce our least-squares functional E : V → R
+
as follows
(2.7) E(y) :=
1
2
Z
Ω
(α|v|
2
+ ν|∇v|
2
) =
1
2
|||v|||
2
V
where the corrector v ∈ V is the unique solution of the linear formulation
(2.8)
α
Z
Ω
v · w + ν
Z
Ω
∇v · ∇w = −α
Z
Ω
y · w − ν
Z
Ω
∇y · ∇w −
Z
Ω
y · ∇y · w
+ < f, w >
H
−1
(Ω)
d
×H
1
0
(Ω)
d
+α
Z
Ω
g · w, ∀w ∈ V .
In particular, for d = 2, the corrector v satisfies the estimate:
(2.9) |||v|||
V
≤ |||y|||
V
1 +
|||y|||
V
2
√
αν
+
s
c
0
kfk
2
H
−1
ν
+ αkgk
2
2
.
Conversely, we also have
(2.10) |||y|||
V
≤ |||v|||
V
+
s
c
0
kfk
2
H
−1
ν
+ αkgk
2
2
.
The infimum of E is equal to zero and is reached by a solution of (2.5). In this sense, the
functional E is a so-called error functional which measures, through the corrector variable v, the
deviation of the pair y from being a solution of the underlying equation (2.5).
A practical way of taking a functional to its minimum is through some (clever) use of descent
directions, i.e. the use of its derivative. In doing so, the presence of local minima is always
something that may dramatically spoil the whole scheme. The unique structural property that
discards this possibility is the strict convexity of the functional. However, for non-linear equations
like (2.5), one cannot expect this property to hold for the functional E in (2.7). Nevertheless,
we insist in that for a descent strategy applied to the extremal problem min
y∈V
E(y) numerical
procedures cannot converge except to a global minimizer leading E down to zero.
Indeed, we would like to show that the only critical points for E correspond to solutions of
(2.5). In such a case, the search for an element y solution of (2.5) is reduced to the minimization
of E.
For any y ∈ V , we now look for an element Y
1
∈ V solution of the following formulation
(2.11) α
Z
Ω
Y
1
·w+ν
Z
Ω
∇Y
1
·∇w+
Z
Ω
(y·∇Y
1
+Y
1
·∇y)·w = −α
Z
Ω
v·w −ν
Z
Ω
∇v·∇w, ∀w ∈ V
where v ∈ V is the corrector (associated to y) solution of (2.8). Y
1
enjoys the following property.
Proposition 2.5. For all y ∈ V satisfying τ
d
(y) < 1, there exists a unique solution Y
1
of (2.11)
associated to y. Moreover, this solution satisfies
(2.12) (1 − τ
d
(y))|||Y
1
|||
V
≤
p
2E(y).
Proof. The proof uses the arguments of Lemma 2.4. We define the bilinear and continuous form
a : V × V → R by
(2.13) a(Y, w) = α
Z
Ω
Y · w + ν
Z
Ω
∇Y · ∇w +
Z
Ω
(y · ∇Y + Y · ∇y) · w
