Jérôme Lemoine
email: jerome.lemoine@uca.fr

Arnaud Münch
email: arnaud.munch@uca.fr

Resolution of the implicit Euler scheme for the Navier-Stokes equation through a least-squares method

Keywords: Steady Navier-Stokes equation, Implicit time scheme, Leastsquares approach, Damped Newton method

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 construct 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 equation from any initial guess. The method turns out to be related to the globally convergent damped Newton approach applied to the Navier-Stokes operator. 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.

Introduction -Motivation

Let Ω ⊂ R d , d = 2 or d = 3 be a bounded connected open set whose boundary ∂Ω is Lipschitz and T > 0. We endow H 1 0 (Ω) with the scalar product v, w H 1 0 (Ω) = Ω ∇v • ∇w and the associated norm and we endow the dual H -1 (Ω) of H 1 0 (Ω) with the dual norm of H 1 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 0 (Ω) d endowed with the norm of H 1 0 (Ω) d . 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 :

     u t -ν∆u + (u • ∇)u + ∇p = F, ∇ • u = 0 in Ω × (0, T ), u = 0 on ∂Ω × (0, T ), u(•, 0) = u 0 in Ω, (1)
where u is the velocity of the fluid, p its pressure and ν is the viscosity constant assumed smaller than one. We refer to [START_REF] Temam | Theory and numerical analysis[END_REF]. This work is concerned with the approximation of (1) through the time marching fully implicit Euler scheme

               y 0 = u 0 in Ω, y n+1 -y n δt -ν∆y n+1 + (y n+1 • ∇)y n+1 + ∇π n+1 = 1 δt tn+1 tn F (•, s)ds, n ≥ 0, ∇ • y n+1 = 0 in Ω, n ≥ 0, y n+1 = 0 on ∂Ω, n ≥ 0, (2)
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) as δt goes to zero. It achieves a first order convergence with respect to δt. We also refer to [START_REF] Tone | On the long-time stability of the implicit Euler scheme for the two-dimensional Navier-Stokes equations[END_REF] for a stability analysis of the scheme in long time. We refer to [START_REF] Smith | Implicit algorithms and their linearization for the transient incompressible Navier-Stokes equations[END_REF] 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 α y -ν∆y + (y • ∇)y + ∇π = f + α g, ∇ • y = 0 in Ω, y = 0 on ∂Ω,

with

α := 1 δt > 0, f := 1 δt tn+1 tn F (•, s)ds, g = y n . (4)
Recall that for any f ∈ H -1 (Ω) d and g ∈ L 2 (Ω) d , there exists one solution (y, π)

∈ V × L 2 0 (Ω) of (3), unique if g 2 L 2 (Ω) d + α -1 ν -1 f 2 H -1
(Ω) d is small enough (see Proposition 1 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 (3) solves the formulation F (y, w) = 0 for all w ∈ V where F is defined by

F (y, z) := Ω α y • z + ν∇y • ∇z + (y • ∇)y • z -f, z H -1 (Ω) d ×H 1 0 (Ω) d -α Ω g • z = 0, ∀z ∈ V .
(5) If D y F is invertible, one may approximate a weak solution through the iterative Newton method: a sequence {y k } k∈N ∈ V is constructed as follows

y 0 ∈ V , D y F (y k , w) • (y k+1 -y k ) = -F (y k , w), ∀w ∈ V , k ≥ 0. ( 6
)
If the initial guess y 0 is close enough to a weak solution of (3), i.e. a solution satisfying F (y, w) = 0 for all w, then the sequence {y k } k∈N converges. We refer to [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]Section 10.3], [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]Chapter 6]) and for some numerical aspects to [START_REF] Dong | Newton's method for the Navier-Stokes equations with finite-element initial guess of Stokes equations[END_REF].

Alternatively, we may also employ least-squares methods which consist in minimizing a quadratic functional, which measures 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

E(y) := 1 2 Ω α|v| 2 + ν|∇v| 2 (7)
where the corrector v, together with the pressure, is the unique solution in V × L 2 0 (Ω) of the linear boundary value problem:

αv -ν∆v + ∇π + αy -ν∆y + (y • ∇)y -f -αg = 0, ∇ • v = 0 in Ω, v = 0 on ∂Ω. (8)
E(y) vanishes if and only if y ∈ V is a weak solution of (3), equivalently a zero of F (y, w) = 0 for all w ∈ V . As a matter of fact, the infimum is reached. Leastsquares 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 [START_REF] Bochev | Least-squares finite element methods[END_REF], [START_REF] Glowinski | Variational methods for the numerical solution of nonlinear elliptic problems[END_REF]. The minimization of the functional E over V leads to a so-called weak least squares method. Precisely, the equality 2E(y) = sup w∈V ,w =0 F (y,w) |||w||| V where |||w||| V is defined in [START_REF] Jiang | Least-squares finite element method for fluid dynamics[END_REF] -shows that E is equivalent to the V norm of the Navier-Stokes equation (see Remark 1). The terminology "H -1 -least-squares method" is employed in [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF] where the minimization of E has been introduced and numerically implemented to approximate solutions of (1) through the scheme [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF]. We also mention [4, Chapter 4, Section 6] which studied 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 [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF] and show that one may construct minimizing sequences in V for E that converge strongly toward a solution of [START_REF] Deuflhard | Newton methods for nonlinear problems[END_REF]. The second objective is to justify the use of that least-squares method to solve iteratively a weak formulation of the scheme [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF], leading to an approximation of the solution of [START_REF] Bochev | Least-squares finite element methods[END_REF]. 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 ( 7)-( 8) associated to weak solutions of [START_REF] Deuflhard | Newton methods for nonlinear problems[END_REF]. We show that E is differentiable over V and that any critical point for E in the ball B := {y ∈ V , τ d (y) < 1} (see Definition 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 (y) • Y 1 is proportional to E(y). Then, assuming that there exists at least one solution of (3) in B, we show that any minimizing sequence {y k } (k∈N) for E in B strongly converges to a solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. Such limit belongs to B and is actually the unique solution. Eventually, we construct a minimizing sequence (defined in (30)) based on the element Y 1 and initialized with g assumed in V . If α is large enough, we show that this particular sequence belongs to B and converges (quadratically after a finite number of iterates related to the values of ν and α) strongly to the solution of (3) (see Theorem 1). This specific sequence coincides with the one obtained from the damped Newton method, a globally convergent generalization of [START_REF] Glowinski | Finite element methods for incompressible viscous flow[END_REF]. Then, in Section (3), as an application, we consider the least-squares approach to solve iteratively the backward Euler scheme (see (51)), weak formulation of (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, assuming u 0 L 2 (Ω) d + F 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 4). 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), as mentioned in [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]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 [START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF]. We notably exhibit for small values of the viscosity constant the robustness of the damped Newton method (compared to the Newton one).

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 (3) assuming α > 0: we extend [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF] where the particular case α = 0 is addressed.

Technical preliminary results

In the sequel • 2 stands for the norm • L 2 (Ω) d . We shall also use the following notations

|||y||| 2 V := α y 2 2 + ν ∇y 2 2 , ∀y ∈ V (9)
and y, z

V := α Ω y • z + ν Ω ∇y • ∇z so that y, z V ≤ |||y||| V |||z||| V for any y, z ∈ V .
In the sequel, we repeatedly use the following classical estimates (see [START_REF] Temam | Theory and numerical analysis[END_REF]).

Lemma 1 Let any u, v ∈ V . If d = 2, then - Ω u • ∇u • v = Ω u • ∇v • u ≤ √ 2 u 2 ∇v 2 ∇u 2 . ( 10
)
If d = 3, then there exists a constant c = c(Ω) such that Ω u • ∇v • u ≤ c u 1/2 2 ∇v 2 ∇u 3/2 2 . ( 11
)
Definition 1 For any y ∈ V , α > 0 and ν > 0, we define

τ d (y) :=          y H 1 0 (Ω) 2 √ 2αν , if d = 2, M y H 1 0 (Ω) 3 (αν 3 ) 1/4 , if d = 3
with M := 3 3/4 4 c and c from [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF].

We shall also repeatedly use the following Young type inequalities.

Lemma 2 For any u, v ∈ V , the following inequalities hold true :

√ 2 u 2 ∇v 2 ∇u 2 ≤ τ 2 (v)|||u||| 2 V ( 12
)
if d = 2 and c u 1/2 2 ∇v 2 ∇u 3/2 2 ≤ τ 3 (v)|||u||| 2 V ( 13
) if d = 3. Let f ∈ H -1 (Ω) d , g ∈ L 2 (Ω) d and α ∈ R + .
The weak formulation of (3) reads as follows: find y ∈ V solution of

Ω α y•w+ν∇y•∇w+y•∇y•w = f, w H -1 (Ω) d ×H 1 0 (Ω) d +α Ω g•w, ∀w ∈ V . ( 14
)
The following result holds true (we refer to [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]).

Proposition 1 Assume Ω ⊂ R d is bounded and Lipschitz. There exists at least one solution y of ( 14) satisfying

|||y||| 2 V ≤ 1 ν f 2 H -1 (Ω) d + α g 2 2 . ( 15
)
If moreover Ω is C 2 and f ∈ L 2 (Ω) d , then any solution y ∈ V of ( 14) belongs to H 2 (Ω) d .

Lemma 3 Assume that a solution y ∈ V of ( 14) satisfies τ d (y) < 1. Then, such solution is the unique solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF].

Proof Let y 1 ∈ V and y 2 ∈ V be two solutions of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. Set Y = y 1 -y 2 . Then,

α Ω Y • w + ν Ω ∇Y • ∇w + Ω y 2 • ∇Y • w + Ω Y • ∇y 1 • w = 0 ∀w ∈ V .
We now take w = Y and use that

Ω y 2 • ∇Y • Y = 0. If d = 2
, we use ( 10) and ( 12) to get

|||Y ||| 2 V = - Ω 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 ( 14) is unique. In view of [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF], this holds if the data satisfy [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF] and ( 13) to obtain

ν g 2 2 + 1 α f 2 H -1 (Ω) d < 2ν 2 . If d = 3, we use
|||Y ||| 2 V = - Ω Y • ∇y 1 • Y ≤ c Y 1 2 2 ∇Y 3 2 2 ∇y 1 2 ≤ τ 3 (y 1 )|||Y ||| 2 V leading to 1 -τ 3 (y 1 ) Y 2
2 ≤ 0 and to the uniqueness if τ 3 (y 1 ) < 1. In view of [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF], this holds if the data satisfy ν g 2 2

+ 1 α f 2 H -1 (Ω) d < M -2 ν 7/2 α -1/2 .
We now introduce our least-squares functional E : V → R + as follows

E(y) := 1 2 Ω (α|v| 2 + ν|∇v| 2 ) = 1 2 |||v||| 2 V ( 16
)
where the corrector v ∈ V is the unique solution of the linear formulation

α Ω v • w + ν Ω ∇v • ∇w = -α Ω y • w -ν Ω ∇y • ∇w - Ω y • ∇y • w + f, w H -1 (Ω) d ×H 1 0 (Ω) d + α Ω g • w, ∀w ∈ V . (17)
In particular, for d = 2, the corrector v satisfies the estimate:

|||v||| V ≤ |||y||| V 1 + |||y||| V 2 √ αν + f 2 H -1 (Ω) 2 ν + α g 2 2 . (18)
Conversely, we also have

|||y||| V ≤ |||v||| V + f 2 H -1 (Ω) 2 ν + α g 2 2 . ( 19
)
The infimum of E is equal to zero and is reached by a solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. 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 [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF].

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 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 [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF], one cannot expect this property to hold for the functional E in [START_REF] Saramito | A damped Newton algorithm for computing viscoplastic fluid flows[END_REF]. 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 critical points for E correspond to solutions of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. In such a case, the search for an element y solution of ( 14) 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

α Ω Y 1 •w+ν Ω ∇Y 1 •∇w+ Ω (y•∇Y 1 +Y 1 •∇y)•w = -α Ω v•w-ν Ω ∇v•∇w, (20)
for all w ∈ V where v ∈ V is the corrector (associated to y) solution of [START_REF] Smith | Implicit algorithms and their linearization for the transient incompressible Navier-Stokes equations[END_REF]. Y 1 enjoys the following property.

Proposition 2 For all y ∈ V satisfying τ d (y) < 1, there exists a unique solution Y 1 of (20) associated to y. Moreover, this solution satisfies

(1 -τ d (y))|||Y 1 ||| V ≤ 2E(y). ( 21
)
Proof We define the bilinear and continuous form a :

V × V → R by a(Y, w) = α Ω Y • w + ν Ω ∇Y • ∇w + Ω (y • ∇Y + Y • ∇y) • w so that a(Y, Y ) = |||Y ||| 2 V + Ω Y • ∇y • Y . If d = 2, using (12), we obtain a(Y, Y ) ≥ (1 -τ 2 (y))|||Y ||| 2 V for all Y ∈ V . Lax-Milgram lemma leads to the existence and uniqueness of Y 1 assuming that τ 2 (y) < 1. Then, putting w = Y 1 in (20) implies a(Y 1 , Y 1 ) ≤ -α Ω v • Y 1 -ν Ω ∇v • ∇Y 1 ≤ |||Y 1 ||| V |||v||| V = |||Y 1 ||| V 2E(y) leading to (21). If d = 3, using (13), we obtain a(Y, Y ) ≥ (1 -τ 3 (y))|||Y ||| 2
V for all Y ∈ V and we conclude as before.

We now check the differentiability of the least-squares functional.

Proposition 3 For all y ∈ V , the map Y → E(y + Y ) is a differentiable function on the Hilbert space V and for any Y ∈ V , we have

E (y) • Y = Ω α v • V + ν∇v • ∇V (22)
where V ∈ V is the unique solution of

Ω α V •w+ν∇V •∇w = -α Ω Y •w-ν Ω ∇Y •∇w- Ω (y•∇Y +Y •∇y)•w, ∀w ∈ V . ( 23
) Proof Let y ∈ V and Y ∈ V . We have E(y + Y ) = 1 2 V 2 V where V ∈ V is the unique solution of α Ω V • w + ν Ω ∇V • ∇w + α Ω (y + Y ) • w + ν Ω ∇(y + Y ) • ∇w+ Ω (y + Y ) • ∇(y + Y ) • w -f, w H -1 (Ω) d ×H 1 0 (Ω) d -α Ω g • w = 0, ∀w ∈ V . If v ∈ V is the solution of (17) associated to y, v ∈ V is the unique solution of α Ω v • w + ν Ω ∇v • ∇w + Ω Y • ∇Y • w = 0, ∀w ∈ V (24)
and V ∈ V is the unique solution of (23), then it is straightforward to check that

V -v -v -V ∈ V is solution of α Ω (V -v -v -V ) • w + ν Ω ∇(V -v -v -V ) • ∇w = 0, ∀w ∈ V and therefore V -v -v -V = 0. Thus E(y + Y ) = 1 2 |||v + v + V ||| 2 V = 1 2 |||v||| 2 V + 1 2 |||v ||| 2 V + 1 2 |||V ||| 2 V + V, v V + V, v V + v, v V .
(25) Assume d = 2. Then, writing (23) with w = V and using [START_REF] Dong | Newton's method for the Navier-Stokes equations with finite-element initial guess of Stokes equations[END_REF], we obtain

|||V ||| 2 V ≤ |||V ||| V |||Y ||| V + √ 2( y 2 ∇Y 2 + Y 2 ∇y 2 ) ∇V 2 ≤ |||V ||| V |||Y ||| V + √ 2 √ αν |||y||| V |||Y ||| V ∇V 2 leading to |||V ||| V ≤ |||Y ||| V (1 + √ 2 √ αν |||y||| V ). Similarly, using (24), we obtain |||v ||| V ≤ 1 √ 2αν |||Y ||| 2 V . It follows that 1 2 |||v ||| 2 V + 1 2 |||V ||| 2 V + V, v V + v, v V = o(|||Y ||| V ) and from (25) that E(y + Y ) = E(y) + v, V V + o(|||Y ||| V ). Eventually, the estimate | v, V V | ≤ |||v||| V |||V ||| V ≤ √ 2(1+ √ 2 √ αν |||y||| V ) E(y)|||Y ||| V
gives the continuity of the linear map Y → v, V V . The case d = 3 is similar.

We are now in position to prove the following result which indicates that, in the ball B, any critical point for E is also a zero of E.

Proposition 4 For all y ∈ V satisfying τ d (y) < 1, (1 -τ d (y)) 2E(y) ≤ 1 √ ν E (y) H -1 (Ω) d . Proof For any Y ∈ V , E (y) • Y = Ω α v • V + ν∇v • ∇V where V ∈ V is
the unique solution of (23). In particular, taking Y = Y 1 defined by [START_REF] Tone | On the long-time stability of the implicit Euler scheme for the two-dimensional Navier-Stokes equations[END_REF], we obtain an element

V 1 ∈ V solution of Ω α V 1 •w+ν∇V 1 •∇w = -α Ω Y 1 •w-ν Ω ∇Y 1 •∇w- Ω (y•∇Y 1 +Y 1 •∇y)•w, (26)
for all w ∈ V . Summing [START_REF] Tone | On the long-time stability of the implicit Euler scheme for the two-dimensional Navier-Stokes equations[END_REF] and (26), we obtain that v -

V 1 ∈ V solves α Ω (v -V 1 ) • w + ν Ω (∇v -∇V 1 ) • w = 0, ∀w ∈ V .
This implies that v and V 1 coincide and then that

E (y) • Y 1 = Ω α|v| 2 + ν|∇v| 2 = 2E(y), ∀y ∈ V . (27)
It follows that

2E(y) = E (y) • Y 1 ≤ E (y) H -1 (Ω) d Y 1 H 1 0 (Ω) d ≤ E (y) H -1 (Ω) d |||Y 1 ||| V √ ν .
Proposition 2 allows to conclude.

Eventually, we prove the following coercivity type inequality for the functional E.

Proposition 5 Assume that a solution y ∈ V of (14) satisfies τ d (y) < 1.

Then, for all y ∈ V ,

|||y -y||| V ≤ 1 -τ d (y) -1 2E(y). ( 28
)
Proof For any y ∈ V , let v be the corresponding corrector and let Y = y -y.

We have

α Ω Y •w+ν Ω ∇Y •∇w+ Ω y•∇Y •w+ Ω Y •∇y•w = -α Ω v•w-ν Ω ∇v•∇w (29) for all w ∈ V . For w = Y , this equality rewrites |||Y ||| 2 V = - Ω Y • ∇y • Y -α Ω v • Y -ν Ω ∇v • ∇Y.
Repeating the arguments of the proof of Proposition 2, the result follows.

Assuming the existence of a solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF] in the ball B = {y ∈ V , τ d (y) < 1}, Proposition 4 and Proposition 5 imply that any minimizing sequence {y k } (k∈N) for E in B strongly converges to a solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. Remark that, from Lemma 3, such solution is unique. In the next section, assuming the parameter α large enough, we construct such sequence {y k } (k∈N) .

Remark 1 In order to simplify notations, we have introduced the corrector variable v leading to the functional E. Instead, we may consider the functional E : V → R defined by

E(y) := 1 2 αy + νB 1 (y) + B(y, y) -f -αg 2 V with B 1 : V → L 2 (Ω) d and B : V × V → L 2 (Ω) d defined
E(y) = 1 2 |||v||| 2 V ≤ 1 2ν αy + νB 1 (y) + B(y, y) -f -αg 2 V = 1 ν E(y), ∀y ∈ V .
Conversely,

αy + νB 1 (y) + B(y, y) -f -αg V = sup w∈V ,w =0 Ω (αv • w + ν∇v • ∇w) w H 1 0 (Ω) d ≤|||v||| V sup w∈V ,w =0 |||w||| V w H 1 0 (Ω) d ≤ √ α + ν|||v||| V
so that E(y) ≤ (α + ν)E(y) for all y ∈ V .

A strongly convergent minimizing sequence for E

We define in this section a sequence converging strongly to a solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF] for which E vanishes. According to Proposition 4, it suffices to define a minimizing sequence for E included in the ball B. In this respect, the equality (27) shows that -Y 1 (see [START_REF] Tone | On the long-time stability of the implicit Euler scheme for the two-dimensional Navier-Stokes equations[END_REF]) is a descent direction for the functional E. Therefore, we can define at least formally, for any m ≥ 1, the minimizing sequence {y k } (k≥0) :

     y 0 ∈ V given, y k+1 = y k -λ k Y 1,k , k ≥ 0, λ k = argmin λ∈[0,m] E(y k -λY 1,k ) (30) with Y 1,k ∈ V the solution of the formulation α Ω Y 1,k • w + ν Ω ∇Y 1,k • ∇w + Ω (y k • ∇Y 1,k + Y 1,k • ∇y k ) • w = -α Ω v k • w -ν Ω ∇v k • ∇w, ∀w ∈ V (31)
and v k ∈ V the corrector (associated to y k ) solution of [START_REF] Smith | Implicit algorithms and their linearization for the transient incompressible Navier-Stokes equations[END_REF]. The algorithm (30) can be expanded as follows:

                     y 0 ∈ V given, Ω α y k+1 • w + ν∇y k+1 • ∇w + y k • ∇y k+1 • w + y k+1 • ∇y k • w = (1 -λ k ) Ω α y k • w + ν∇y k • ∇w + y k • ∇y k • w + λ k Ω α g • w + f, w H -1 (Ω)×H 1 0 (Ω) + Ω y k • ∇y k • w k ≥ 0, ∀w ∈ V .
(32) From ( 19), the sequence {y k } k>0 is bounded. However, we insist that, in order to justify the existence of the element Y 1,k , y k should satisfy τ d (y k ) < 1 for all k. We proceed in two steps: assuming that the sequence {y k } (k>0) defined by (30) satisfies τ d (y k ) ≤ c 1 < 1 for all k, we show that E(y k ) → 0 and that {y k } k∈N converges strongly in V to a solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. Then, we determine sufficient conditions on the initial guess y 0 ∈ V so that τ d (y k ) < 1 for all k ∈ N.

We start with the following lemma which provides the main property of the sequence {E(y k )} (k≥0) .

Lemma 4 Assume that the sequence {y k } (k≥0) defined by (30) satisfies

τ d (y k ) < 1 for all k ≥ 0. Then, for all λ ∈ R, E(y k -λY 1,k ) ≤ E(y k ) |1 -λ| + λ 2 (1 -τ 2 (y k )) -2 √ αν E(y k ) 2 (33) if d = 2 and E(y k -λY 1,k ) ≤ E(y k ) |1-λ|+λ 2 √ 2 √ ν M (αν 3 ) 1/4 (1-τ 3 (y k )) -2 E(y k ) 2 (34) if d = 3.
Proof For any real λ and any y k , w k ∈ V we get the expansion

E(y k -λw k ) = E(y k ) -λ v k , v k V + λ 2 2 v k , v k V + 2 v k , v k V -λ 3 v k , v k V + λ 4 2 v k , v k V (35) where v k , v k ∈ V and v k ∈ V solves respectively α Ω v k • w + ν Ω ∇v k • ∇w + α Ω y k • w + ν Ω ∇y k • ∇w + Ω y k • ∇y k • w = f, w H -1 (Ω) d ×H 1 0 (Ω) d + α Ω g • w, ∀w ∈ V , (36)
α

Ω v k • w + ν Ω ∇v k • ∇w + α Ω w k • w + ν Ω ∇w k • ∇w + Ω w k • ∇y k • w + y k • ∇w k • w = 0, ∀w ∈ V , (37) and α
Ω v k • w + ν Ω ∇v k • ∇w + Ω w k • ∇w k • w = 0, ∀w ∈ V . ( 38
)
Since the corrector v k associated to Y 1,k coincides with the corrector v k associated to y k (see proof of Proposition ( 4)), expansion (35) reduces to

E(y k -λY 1,k ) = (1 -λ) 2 E(y k ) + λ 2 (1 -λ) v k , v k V + λ 4 2 v k , v k V ≤ |1 -λ| E(y k ) + λ 2 √ 2 v k V 2 . ( 39
) If d = 2, then (38) leads to v k V ≤ |||Y 1,k ||| 2 V √ 2αν ≤ √ 2(1 -τ 2 (y k )) -2 E(y k ) √ αν and then to (33). If d = 3, then v k V ≤ 1 √ ν M (αν 3 ) 1/4 |||Y 1,k ||| 2 V ≤ 1 √ ν M (αν 3 ) 1/4 2(1 -τ 3 (y k )) -2 E(y k ) leading to (34).
We are now in position to prove the convergence of the sequence {E(y k )} (k≥0) .

Proposition 6 Let {y k } k≥0 be the sequence defined by (30). Assume that there exists a constant c 1 ∈ (0, 1) such that τ d (y k ) ≤ c 1 for all k. Then E(y k ) → 0 as k → ∞. Moreover, there exists k 0 ∈ N such that the sequence {E(y k )} (k≥k0) decays quadratically.

Proof Consider the case d = 2. The inequality τ 2 (y k ) ≤ c 1 and (33) imply that

E(y k -λY 1,k ) ≤ E(y k ) |1 -λ| + λ 2 c α,ν E(y k ) 2 , c α,ν := (1 -c 1 ) -2 √ αν .
Let us denote the function

p k (λ) = |1-λ|+λ 2 c α,ν E(y k ) for all λ ∈ [0, m]. We can write E(y k+1 ) = min λ∈[0,m] E(y k -λY 1,k ) ≤ p k ( λ k ) E(y k ). with p k ( λ k ) := min λ∈[0,m] p k (λ).
Suppose first that c α,ν E(y 0 ) ≥ 1 and prove that the set

I := {k ∈ N, c α,ν E(y k ) ≥ 1} is a finite subset of N. For all k ∈ I, we get min λ∈[0,m] p k (λ) = min λ∈[0,1] p k (λ) = p k 1 2c α,ν E(y k ) = 1 - 1 4c α,ν E(y k )
and thus, for all k ∈ I,

c α,ν E(y k+1 ) ≤ 1 - 1 4c α,ν E(y k ) c α,ν E(y k ) = c α,ν E(y k ) - 1 4 .
Consequently, the sequence {c α,ν E(y k )} k∈I strictly decreases and thus, there exists k 0 ∈ N such that for all k ≥ k 0 , c α,ν E(y k ) < 1. Thus I is a finite subset of N. Then, for all k ≥ k 0 , we get that

p k ( λ k ) ≤ p k (1) = c α,ν E(y k )
and thus, for all k ≥ k 0 ,

c α,ν E(y k+1 ) ≤ c α,ν E(y k ) 2 (40)
implying that c α,ν E(y k ) → 0 as k → ∞ with a quadratic rate.

On the other hand, if c α,ν E(y 0 ) < 1 (and thus c α,ν E(y k ) < 1 for all k ∈ N, since by construction the sequence {E(y k )} k decreases), then (40) holds true for all k ≥ 0.

In both cases, remark that p k ( λ k ) decreases with respect to k. The case d = 3 is similar with

c α,ν = √ 2 √ ν M (αν 3 ) 1/4 (1 -c 1 ) -2 .
Lemma 5 Assume that the sequence {y k } (k≥0) defined by (30) satisfies τ d (y k ) ≤ c 1 for all k and some c 1 ∈ (0, 1).

Then λ k → 1 as k → ∞.
Proof In view of (39), we have, as long as

E(y k ) > 0, (1 -λ k ) 2 = E(y k+1 ) E(y k ) -λ 2 k (1 -λ k ) v k , v k V E(y k ) -λ 4 k v k 2 V 2E(y k ) . From the proof of Lemma 4, v k ,v k V E(y k ) ≤ C(α, ν)(1 -c 1 ) -2 E(y k ) while |||vk||| 2 V E(y k ) ≤ C(α, ν) 2 (1-c 1 ) -4 E(y k ). Consequently, since λ k ∈ [0, m] and E(y k+1 ) E(y k ) → 0, we deduce that (1 -λ k ) 2 → 0 as k → ∞.
Proposition 7 Let {y k } (k≥0) be the sequence defined by (30). Assume that there exists a constant c 1 ∈ (0, 1) such that τ d (y k ) ≤ c 1 for all k. Then, y k → y in V where y ∈ V is the unique solution of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF].

Proof Remark that we can not use Proposition 5 since we do not know yet that there exists a solution, say z, of ( 14) satisfying τ d (z) < 1. In view of

y k+1 = y 0 - k n=0 λ n Y 1,n , we write k n=0 |λ n ||||Y 1,n ||| V ≤ m k n=0 |||Y 1,n ||| V ≤ m √ 2 k n=0 E(y n ) 1 -τ d (y n ) ≤ m √ 2 1 -c 1 k n=0 E(y n ).
Using that p n ( λ n ) ≤ p 0 ( λ 0 ) for all n ≥ 0, we obtain for all n > 0,

E(y n ) ≤ p n-1 ( λ n-1 ) E(y n-1 ) ≤ p 0 ( λ 0 ) E(y n-1 ) ≤ p 0 ( λ 0 ) n E(y 0 ).
Recalling that p 0 ( λ 0 ) = min λ∈[0,1] p 0 (λ) < 1 since p 0 (0) = 1 and p 0 (0) = -1, we finally obtain

k n=0 |λ n ||||Y 1,n ||| V ≤ m √ 2 1 -c 1 E(y 0 ) 1 -p 0 ( λ 0 ) from which we deduce that the series k≥0 λ k Y 1,k converges in V . Then, y k converges in V to y := y 0 + k≥0 λ k Y 1,k .
Eventually, the convergence of E(y k ) to 0 implies the convergence of the corrector v k to 0 in V ; taking the limit in the corrector equation (36) shows that y solves [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. Since τ d (y) ≤ c 1 < 1, Lemma 3 shows that this solution is unique.

As mentioned earlier, the remaining and crucial point is to show that the sequence {y k } (k≥0) satisfies the uniform property τ d (y k ) ≤ c 1 for some c 1 < 1.

Lemma 6 Assume that y 0 = g ∈ V . For all c 1 ∈ (0, 1) there exists α 0 > 0, such that, for any α ≥ α 0 , the unique sequence defined by (30) satisfies τ d (y k ) ≤ c 1 for all k ≥ 0.

Proof Let c 1 ∈ (0, 1) and assume that y 0 belongs to V . Since τ d (y 0 ) → 0 as α → ∞, there exists α 1 > 0 such that for all α ≥ α 1 τ d (y 0 ) ≤ c1 2 . Moreover, in view of the above computation and using that v

H 1 0 (Ω) d ≤ 1 √ ν |||v||| V for all v ∈ V and α > 0, we obtain, for all k ∈ N y k+1 H 1 0 (Ω) d ≤ y 0 H 1 0 (Ω) d + m √ 2 √ ν(1 -c 1 ) E(y 0 ) 1 -p 0 ( λ 0 )
where

E(y 0 ) 1 -p 0 ( λ 0 ) ≤      E(y 0 ) 1 -c α,ν E(y 0 ) , if c α,ν E(y 0 ) < 1, 4c α,ν E(y 0 ), if c α,ν E(y 0 ) ≥ 1.
Assume d = 2. From [START_REF] Smith | Implicit algorithms and their linearization for the transient incompressible Navier-Stokes equations[END_REF], we obtain for all y ∈ V that

|||v||| 2 V ≤ α g -y 2 2 + 1 ν ν ∇y 2 + √ 2 y 2 ∇y 2 + f H -1 (Ω) 2 2 .
In particular, taking y = g allows to remove the α term and gives

E(g) ≤ 1 2ν g H 1 0 (Ω) d (ν + √ 2 g 2 ) + f H -1 (Ω) 2 2 := 1 2ν c 2 (f, g). ( 41
)
If c α1,ν E(g) ≥ 1 then for all α ≥ α 1 such that c α,ν E(g) ≥ 1 and for all k ∈ N :

y k+1 H 1 0 (Ω) d ≤ g H 1 0 (Ω) d + m √ 2 √ ν(1 -c 1 ) E(g) 1 -p 0 ( λ 0 ) ≤ g H 1 0 (Ω) d + 2m √ 2 ν 5/2 √ α(1 -c 1 ) 3 c 2 (f, g). (42)
If c α1,ν E(g) < 1 then there exists 0 < K < 1 such that for all α ≥ α 1 we have c α,ν E(g) ≤ K. We therefore have for all α ≥ α 1 E(g)

1 -p 0 ( λ 0 ) ≤ E(g) 1 -K
and thus for all k ∈ N :

y k+1 H 1 0 (Ω) d ≤ g H 1 0 (Ω) d + m √ 2 √ ν(1 -c 1 ) E(g) 1 -p 0 ( λ 0 ) ≤ g H 1 0 (Ω) d + m ν(1 -c 1 )(1 -K) c 2 (f, g). (43)
On the other hand, there exists α 0 ≥ α 1 such that, for all α ≥ α 0 we have

2m √ 2 ν 5/2 √ α(1 -c 1 ) 3 c 2 (f, g) ≤ c 1 2 √ 2αν and m ν(1 -c 1 )(1 -K) c 2 (f, g) ≤ c 1 2 √ 2αν.
We then deduce from (42) and ( 43) that for all α ≥ α 0 and for all k ∈ N :

y k+1 H 1 0 (Ω) d ≤ c 1 2 √ 2αν + c 1 2 √ 2αν = c 1 √ 2αν that is τ 2 (y k+1 ) ≤ c 1 .
Assume d = 3. We argue as in the case d = 2 and deduce from [START_REF] Smith | Implicit algorithms and their linearization for the transient incompressible Navier-Stokes equations[END_REF], since

y 0 = g, that E(g) ≤ 1 2ν g H 1 0 (Ω) 3 (ν + c y 0 1/2 2 g 1/2 H 1 0 (Ω) 3 ) + f H -1 (Ω) 3 2 := 1 2ν c 3 (f, g) (44
) and thus, if c α1,ν E(g) ≥ 1, then for all α ≥ α 1 such that c α,ν E(g) ≥ 1 and for all k ∈ N :

y k+1 H 1 0 (Ω) 3 ≤ g H 1 0 (Ω) 3 + m √ 2 √ ν(1 -c 1 ) E(g) 1 -p 0 ( λ 0 ) ≤ g H 1 0 (Ω) 3 + 4mM ν 2 (αν 3 ) 1/4 (1 -c 1 ) 3 c 3 (f, g). ( 45
)
If c α1,ν E(g) < 1 then there exists 0 < K < 1 such that for all α ≥ α 1 we have c α,ν E(g) ≤ K. We therefore have for all α ≥ α 1 E(g)

1 -p 0 ( λ 0 ) ≤ E(y 0 ) 1 -K
and thus for all k ∈ N :

y k+1 H 1 0 (Ω) 3 ≤ g H 1 0 (Ω) 3 + m √ 2 √ ν(1 -c 1 ) E(g) 1 -p 0 ( λ 0 ) ≤ g H 1 0 (Ω) 3 + m ν(1 -c 1 )(1 -K) c 3 (f, g). (46)
On the other hand, there exists α 0 ≥ α 1 such that, for all α ≥ α 0 we have

4mM ν 2 (αν 3 ) 1/4 (1 -c 1 ) 3 c 3 (f, g) ≤ c 1 2 (αν 3 ) 1/4 M and m ν(1 -c 1 )(1 -K) c 3 (f, g) ≤ c 1 2 (αν 3 ) 1/4 M .
We then deduce from ( 45) and ( 46) that for all α ≥ α 0 and for all k ∈ N :

y k+1 H 1 0 (Ω) d ≤ c 1 2 (αν 3 ) 1/4 M + c 1 2 (αν 3 ) 1/4 M = c 1 (αν 3 ) 1/4 M that is τ 3 (y k+1 ) ≤ c 1 .
Gathering the previous lemmas and propositions, we deduce the strong convergence of the sequence {y k } k≥0 defined by (30), initialized by y 0 = g. Theorem 1 Let c 1 ∈ (0, 1). Assume that y 0 = g ∈ V and α is large enough so that

c 2 (f, g) ≤ max 1 -c 1 2 , c 1 (1 -K) 2 m c 1 4m ν 2 (1 -c 1 ) 2 2αν, if d = 2, c 3 (f, g) ≤ max 1 -c 1 2 , c 1 (1 -K) 2 m c 1 4mM 2 ν 2 (1 -c 1 ) 2 (αν 3 ) 1/2 , if d = 3,
(47) where c 2 (f, g) and c 3 (f, g) are defined in (41) and ( 44) respectively. The sequence {y k } (k∈N) defined by (30) strongly converges to the unique solution y of (14). Moreover, there exists k 0 ∈ N such that the sequence {y k } k≥k0 converges quadratically to y. Moreover, this solution satisfies τ d (y) < 1.

Additional comments

1) Estimate ( 15) is usually used to obtain a sufficient condition on the data f, g to ensure the uniqueness of the solution of (14) (i.e. τ d (y) < 1): it leads to

α g 2 2 + 1 ν f 2 H -1 (Ω) 2 ≤ 2αν 2 , if d = 2, α g 2 2 + 1 ν f 2 H -1 (Ω) 3 ≤ ν(αν 3 ) 1/2 M 2 , if d = 3.
We emphasize that such (sufficient) conditions are more restrictive than (47), as they impose smallness properties on g:

precisely g 2 2 ≤ 2ν 2 if d = 2 and g 2 2 ≤ ν 5/2 M 2 α 1/2 if d = 3
. This latter yields a restrictive condition for α large contrary to (47).

2) Let F : V → V the application be defined as F(y) = αy + νB 1 (y) + B(y, y) -f -αg. The sequence {y k } (k>0) associated to the Newton method to find the zero of F is formally defined as follows:

y 0 ∈ V , DF(y k ) • (y k+1 -y k ) = -F(y k ), k ≥ 0. ( 48
)
We check that this sequence coincides with the sequence obtained from (30) if λ k is fixed equal to one. The algorithm (30) which consists in optimizing the parameter λ k ∈ [0, m], m ≥ 1, in order to minimize E(y k ), equivalently F(y k ) V , corresponds to the so-called damped Newton method for the application F (see [START_REF] Deuflhard | Newton methods for nonlinear problems[END_REF]). As the iterates increase, the optimal parameter λ k converges to one (according to Lemma 5), this globally convergent method behaves like the standard Newton method (for which λ k is fixed equal to one): this explains the quadratic rate of convergence after a finite number of iterates. To the best of our knowledge, this is the first analysis of the damped Newton method for a partiel differential equation. Among the few numerical works devoted to the damped Newton method for partial differential equations, we mention [START_REF] Saramito | A damped Newton algorithm for computing viscoplastic fluid flows[END_REF] for computing viscoplastic fluid flows.

3) Section 6, chapter 6 of the book [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] introduces a least-squares method in order to solve an Oseen type equation (without incompressibility constraint). The convergence of any minimizing sequence toward a solution y is proved under the a priori assumption that the operator DF(y) :

V → V DF(y) • w = α w -ν∆w + [(w • ∇)y + (y • ∇)w] (49)
(for some α > 0) is an isomorphism. y is then said to be a nonsingular point. According to Proposition 2, a sufficient condition for y to be a nonsingular point is τ d (y) < 1. Recall that τ d depends on α. As far as we know, determining a weaker condition ensuring that DF(y) is an isomorphism is an open question. Moreover, according to Lemma 3, it turns out that this condition is also a sufficient condition for the uniqueness of [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF]. Theorem 1 asserts that, if α is large enough, then the sequence {y k } (k∈N) defined in (30), initialized with y 0 = g, is a convergent sequence of nonsingular points.

4)

We may also define a minimizing sequence for E using the derivative E (see ( 22)):

     y 0 ∈ V given, y k+1 = y k -λ k g k , k ≥ 0, λ k = argmin λ∈[0,m] E(y k -λg k ) (50) with g k ∈ V such that (g k , w) V = E (y k ), w V ×V for all w ∈ V . In partic- ular, g k V = E (y k ) V .
Using the expansion (25) with w k = g k , we can prove the linear decrease of the sequence {E(y k )} k>0 to zero assuming however that E(y 0 ) is small enough, of the order of ν 2 , independently of the value of α (we refer to [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF]Lemma 4.1] in a similar context).

3 Application to the backward Euler scheme

We now use the results of the previous section to discuss the resolution of the backward Euler scheme (2) through a least-squares method. The weak formulation of this scheme reads as follows: given y 0 = u 0 ∈ V , the sequence {y n } n>0 in V is defined by recurrence as follows:

Ω y n+1 -y n δt •w+ν Ω ∇y n+1 •∇w+ Ω y n+1 •∇y n+1 •w = f n , w H -1 (Ω) d ×H 1 0 (Ω) d
(51) with f n defined by (4) in term of the external force of the Navier-Stokes model [START_REF] Bochev | Least-squares finite element methods[END_REF]. We recall that a piecewise linear interpolation in time of {y n } n≥0 weakly converges in L 2 (0, T ; V ) toward a solution of (1).

As done in [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF], one may use the least-squares method (analyzed in Section 2) to solve iteratively (51). Precisely, in order to approximate y n+1 from y n , one may consider the following extremal problem

inf y∈V E n (y), E n (y) = 1 2 |||v||| 2 V ( 52
)
where the corrector v ∈ V solves

α Ω v • w + ν Ω ∇v • ∇w = -α Ω y • w -ν Ω ∇y • ∇w - Ω y • ∇y • w + f n , w H -1 (Ω) d ×H 1 0 (Ω) d + α Ω y n • w, ∀w ∈ V (53
) with α and f n given by (4). For any n ≥ 0, a minimizing sequence {y n k } (k≥0) for E n is defined as follows :

       y n+1 0 = y n , y n+1 k+1 = y n+1 k -λ k Y n+1 1,k , k ≥ 0, λ k = argmin λ∈[0,m] E n (y n+1 k -λY n+1 1,k ) (54)
where Y n 1,k ∈ V solves (31) for y k = y n+1 k . For each n, algorithm (54) can be expanded as follows

                     y n+1 0 = y n , Ω α y n+1 k+1 • w + ν∇y n+1 k+1 • ∇w + y n+1 k • ∇y n+1 k+1 • w + y n+1 k+1 • ∇y n+1 k • w = (1 -λ k ) Ω α y n+1 k • w + ν∇y n+1 k • ∇w + y n+1 k • ∇y n+1 k • w + λ k Ω α y n • w + f n , w H -1 (Ω)×H 1 0 (Ω) + Ω y n+1 k • ∇y n+1 k • w k ≥ 0, (55)
for all w ∈ V . In view of Theorem 1, the first element of the minimizing sequence is chosen equal to y n , i.e. the minimizer of E n-1 .

The main goal of this section is to prove that for all n ∈ N, the minimizing sequence {y n+1 k } k∈N converges to a solution y n+1 of (51). Arguing as in Lemma 6, we have to prove the existence of a constant c 1 ∈ (0, 1), such that τ d (y n k ) ≤ c 1 for all n and k in N. Remark that the initialization y n+1 0 is fixed as the minimizer of the functional E n-1 , obtained at the previous iterate. Consequently, the uniform property τ d (y n k ) ≤ c 1 is related to the initial guess y 0 0 equal to the initial position u 0 , to the external force F (see [START_REF] Bochev | Least-squares finite element methods[END_REF]) and to the value of α. u 0 and F are given a priori. On the other hand, the parameter α, related to the discretization parameter δt, can be chosen as large as necessary. As we shall see, this uniform property, which is essential to set up the least-squares procedure, requires smallness properties on u 0 and F . We start with the following result analogue to Proposition 1.

Proposition 8 Let (f n ) n∈N be a sequence in H -1 (Ω) d , α > 0 and y 0 = u 0 ∈ H. For any n ∈ N, there exists a solution

y n+1 ∈ V of α Ω (y n+1 -y n )•w+ν Ω ∇y n+1 •∇w+ Ω y n+1 •∇y n+1 •w = f n , w H -1 (Ω) d ×H 1 0 (Ω) d
(56) for all w ∈ V . Moreover, for all n ∈ N, y n+1 satisfies

y n+1 2 V ≤ 1 ν f n 2 H -1 (Ω) d + α y n 2 2 . (57)
Moreover, for all n ∈ N :

y n 2 2 + ν α n k=1 ∇y k 2 2 ≤ 1 ν 1 α n-1 k=0 f k 2 H -1 (Ω) d + ν u 0 2 2 . ( 58
)
Proof The existence of y n+1 is given in Proposition 1. ( 58) is obtained by summing (57).

Remark 2 Arguing as in Lemma 3, if there exists a solution y n+1 in V of (53) satisfying τ d (y n+1 ) < 1, then such solution is unique. In view of Proposition 8, this holds true if the quantity M(f, α, ν) defined as follows

M(f, α, ν) =            1 ν 3 1 α n k=0 f k 2 H -1 (Ω) 2 + ν u 0 2 2 , if d = 2, α 1/2 ν 7/2 1 α n k=0 f k 2 H -1 (Ω) 3 + ν u 0 2 2 , if d = 3 (59)
is small enough.

We now distinguish the case d = 2 from the case d = 3 and consider weak and regular solutions.

Two dimensional case

We have the following convergence for weak solutions of (56).

Theorem 2 Suppose F ∈ L 2 (0, T ; H -1 (Ω) 2 ), u 0 ∈ V . Let α be large enough and f n be given by ( 4) for all n ∈ {0, • • • , N -1}. Let c(u 0 , F ) be defined as follows :

c(u 0 , F ) := max 1 α u 0 2 H 1 0 (Ω) d (ν + √ 2 u 0 2 ) 2 + F 2 L 2 (0,T ;H -1 (Ω) 2 ) , 2 F 2 L 2 (0,T ;H -1 (Ω) 2 ) + 2ν u 0 2 2 .
There exists a constant c > 0 such that if

c(u 0 , F ) ≤ cν 3 (60)
then for all n ∈ {0, • • • , N -1} the solution y n+1 ∈ V of (56) is unique and the minimizing sequence {y n+1 k } k∈N defined by (54) strongly converges in V to y n+1 .

Proof According to Proposition 7, we have to prove the existence of a constant c 1 ∈ (0, 1) such that, for all n ∈ {0, • • • , N -1} and all k ∈ N, τ 2 (y n k ) ≤ c 1 . For n = 0, as in the previous section, it suffices to take α large enough to ensure the conditions (47) with g = y 0 0 = u 0 leading to the property τ 2 (y 0 k ) < c 1 for all k ∈ N and therefore τ 2 (y 1 ) ≤ c 1 .

For the next minimizing sequences, we recall (see Lemma 6) that for all

n ∈ {0, • • • , N -1} and all k ∈ N y n+1 k H 1 0 (Ω) 2 ≤ y n H 1 0 (Ω) 2 + m √ 2 √ ν(1 -c 1 ) E n (y n ) 1 -p n,0 ( λ n,0 )
where p n,0 ( λ n,0 ) is defined as in the proof of Proposition 4.

First, since for all

n ∈ {0, • • • , N -1}, f n 2 H -1 (Ω) 2 ≤ α F 2 L 2 (0,T ;H -1 (Ω) 2 ) , we can write E 0 (y 0 ) = E 0 (u 0 ) ≤ 1 2ν u 0 H 1 0 (Ω) d (ν + √ 2 u 0 2 ) + f 0 H -1 (Ω) 2 2 ≤ 1 ν u 0 2 H 1 0 (Ω) d (ν + √ 2 u 0 2 ) 2 + f 0 2 H -1 (Ω) 2 ≤ α ν 1 α u 0 2 H 1 0 (Ω) d (ν + √ 2 u 0 2 ) 2 + F 2 L 2 (0,T ;H -1 (Ω) 2 ) .
Since y n is solution of (56), it follows from (53) and (58), that for all n in {1, • • • , N -1}:

E n (y n ) ≤ 1 2ν f n -f n-1 2 H -1 (Ω) 2 + α 2 y n -y n-1 2 2 ≤ 1 ν f n 2 H -1 (Ω) 2 + f n-1 2 H -1 (Ω) 2 + α y n 2 2 + y n-1 2 2 ≤ 1 ν f n 2 H -1 (Ω) 2 + f n-1 2 H -1 (Ω) 2 + 1 ν n-1 k=0 f k 2 H -1 (Ω) 2 + n-2 k=0 f k 2 H -1 (Ω) 2 ≤ 2 ν n k=0 f k 2 H -1 (Ω) 2 + 2α u 0 2 2 ≤ α ν 2 F 2 L 2 (0,T ;H -1 (Ω) 2 ) + 2ν u 0 2 2 .
Therefore, for all n ∈ {0, • • • , N -1}, E n (y n ) ≤ α ν c(u 0 , F ). Let c 1 ∈ (0, 1) and suppose that c(u 0 , F ) < (1 -c 1 ) 4 ν 3 . Then, there exists K ∈ (0, 1) and α 0 > 0 such that, for all α ≥ α 0 , c α,ν E n (y n ) ≤ K < 1. We therefore have (see Lemma 6), for all α ≥ α 0 , all n ∈ {0, • • • , N -1} and all k ∈ N :

y n+1 k H 1 0 (Ω) 2 ≤ y n H 1 0 (Ω) 2 + m √ 2 √ ν(1 -c 1 ) E n (y n ) 1 -c α,ν E n (y n ) ≤ y n H 1 0 (Ω) 2 + m √ 2 √ ν(1 -c 1 ) E n (y n ) 1 -K ≤ y n H 1 0 (Ω) 2 + m √ 2α ν(1 -c 1 )(1 -K) c(u 0 , F ). (61)
From (58) we then obtain, for all n ∈ {0, • • • , N -1},

y n H 1 0 (Ω) 2 ≤ √ α ν 1 α n-1 k=0 f k 2 H -1 (Ω) 2 + ν u 0 2 2 and since 1 α n-1 k=0 f k 2 H -1 (Ω) 2 ≤ F 2 L 2 (0,T ;H -1 (Ω) 2 )
, we deduce that if

F 2 L 2 (0,T ;H -1 (Ω) 2 ) + ν u 0 2 2 ≤ c 2 1 2 ν 3 then y n H 1 0 (Ω) 2 ≤ c1 2 √ 2αν. Moreover, assuming c(u 0 , F ) ≤ c 2 1 (1-c1) 2 (1-K) 2 4m 2
ν 3 , we deduce from (61), for all n ∈ {0, • • • , N -1} and for all k ∈ N :

y n k H 1 0 (Ω) 2 ≤ c 1 2 √ 2αν + c 1 2 √ 2αν = c 1 √ 2αν that is τ 2 (y n k ) ≤ c 1 .
The result follows from Proposition 7.

We emphasize that, for each n ∈ N, the limit y n+1 of the sequence {y n+1 k } k∈N satisfies τ 2 (y n+1 ) < 1 and is therefore the unique solution of (56). Moreover, for α large enough, the condition (60) reads as the following smallness property on the data:

F 2 L 2 (0,T ;H -1 (Ω) 2 ) + ν u 0 2 2 ≤ cν 3 .
In contrast with the static case of Section (2) where the unique condition (47) on the data g is fulfilled as soon as α is large, the iterated case requires a condition on the data u 0 and F , whatever be the amplitude of α. Again, this smallness property is introduced in order to guarantees the condition τ 2 (y n ) < 1 for all n. In view of (58), this condition implies notably that y n 2 ≤ cν 3/2 for all n > 0. For regular solutions of (56) which we now consider, we may slightly improve the results, notably based on the control of two consecutive elements of the corresponding sequence {y n } n∈N for the L 2 (Ω) norm. 2 and that u 0 ∈ V . Then, for all n ∈ N, any solution y n+1 ∈ V of (56) belongs to H 2 (Ω) 2 .

Proposition 9 Assume that Ω is C 2 , that (f n ) n is a sequence in L 2 (Ω)
If moreover, there exists C > 0 such that

1 α n k=0 f k 2 H -1 (Ω) 2 + ν y 0 2 2 < Cν 3 , (62)
then y n+1 satisfies

Ω |∇y n+1 | 2 + ν 2α n+1 k=1 Ω |P ∆y k | 2 ≤ 1 ν 1 α n k=0 f k 2 2 + ν ∇u 0 2 2 ( 63
)
where P is the operator of projection from L 2 (Ω) d into H.

Proof From Proposition 1, we know that for all n ∈ N * , y n ∈ H 2 (Ω) 2 ∩ V . Thus, integrating by part (56) we obtain, using density argument:

α Ω (y n+1 -y n )•w-ν Ω ∆y n+1 •w+ Ω y n+1 •∇y n+1 •w = f n , w H -1 (Ω) 2 ×H 1 0 (Ω) 2
for all w ∈ H. Then, taking w = P ∆y n+1 and integrating by part leads to :

α Ω |∇y n+1 | 2 +ν Ω |P ∆y n+1 | 2 = - Ω f n P ∆y n+1 + Ω y n+1 • ∇y n+1 • P ∆y n+1 + α Ω ∇y n • ∇y n+1 . (64)
Recall that Ω y n+1 • ∇y n+1 • P ∆y n+1 ≤ y n+1 ∞ ∇y n+1 2 P ∆y n+1 2 .

We now use (see [18, chapter 25]) that there exist three constants c 1 , c 2 and c 3 such that

∆y n+1 2 ≤ c 1 P ∆y n+1 2 , y n+1 ∞ ≤ c 2 y
α 2 Ω |∇y n+1 | 2 + ν 2 -c y n+1 2 Ω |P ∆y n+1 | 2 ≤ 1 2ν f n 2 2 + α 2 Ω |∇y n | 2 .
But, from estimate (58), the assumption (62) implies that y n+1 2 ≤ ν 4c and

Ω |∇y n+1 | 2 + ν 2α Ω |P ∆y n+1 | 2 ≤ 1 να f n 2 2 + Ω |∇y n | 2 .
Summing then implies (63) for all n ∈ N.

Remark 3 Under the hypothesis of Proposition 9, suppose that

B α,ν := (αν 5 ) -1 α -1 n k=0 f k 2 H -1 (Ω) 2 +ν y 0 2 2 α -1 n-1 k=0 f k 2 2 +ν ∇y 0 2 2
is small (which is satisfied as soon as α is large enough). Then, the solution of (56) is unique. Indeed, let n ∈ N and let y n+1 1 , y n+1 2 ∈ V be two solutions of (56). Then

Y := y n+1 1 -y n+1 2 satisfies α Ω Y •w +ν Ω ∇Y •∇w + Ω y n+1 2 •∇Y •w + Ω Y •∇y n+1 1 •w = 0 ∀w ∈ V and in particular, for w = Y (since Ω y n+1 2 • ∇Y • Y = 0) α Ω |Y | 2 + ν Ω |∇Y | 2 = - Ω Y • ∇y n+1 1 • Y = Ω Y • ∇Y • y n+1 1 ≤ c y n+1 1 ∞ ∇Y 2 Y 2 ≤ c y n+1 1 1/2 2 P ∆y n+1 1 1/2 2 ∇Y 2 Y 2 ≤ α Y 2 2 + c α y n+1 1 2 P ∆y n+1 1 2 ∇Y 2 2 leading to ν - c α y n+1 1 2 P ∆y n+1 1 2 ∇Y 2 2 ≤ 0. If y n+1 1 2 P ∆y n+1 1 2 < να c , (65)
then Y = 0 and the solution is unique. But, from ( 58) and (63),

y n+1 1 2 2 P ∆y n+1 1 2 2 ≤ 4α ν 3 1 α n k=0 f k 2 H -1 (Ω) 2 +ν y 0 2 2 1 α n k=0 f k 2 2 +ν ∇y 0 2 2 .
Therefore, there exists a constant C such that if B α,ν < C, then (65) holds true.

Proposition 9 then allows to obtain the following estimate of y n+1 -y n 2 in term of the parameter α.

Theorem 3 We assume that Ω is C 2 , that {f n } n∈N is a sequence in L 2 (Ω) 2 and that α -1 +∞ k=0 f k 2 2 < +∞, that u 0 ∈ V and that for all n ∈ N, y n+1 ∈ H 2 (Ω) 2 ∩ V is a solution of (56) satisfying y n+1 2 ≤ ν 4c . There exists C 1 > 0 such that y n+1 -y n 2 2 ≤ C 1 α ν 3/2 . ( 66
)
Proof For all n ∈ N, w = y n+1 -y n in (56) gives :

α y n+1 -y n 2 2 + ν ∇y n+1 2 2 ≤ Ω y n+1 .∇y n+1 .(y n+1 -y n ) + Ω f n .(y n+1 -y n ) + ν Ω ∇y n .∇y n+1 .
Moreover,

Ω y n+1 .∇y n+1 .(y n+1 -y n ) ≤ c ∇y n+1 2 2 ( ∇y n+1 2 + ∇y n ) 2 ). Ω y n+1 .∇y n+1 .(y n+1 -y n ) ≤ c ∇y n+1 2 2 ∇(y n+1 -y n ) 2 ≤ c ∇y n+1 2 2 ( ∇y n+1 2 + ∇y n ) 2 ).
Therefore,

α y n+1 -y n 2 2 +ν ∇y n+1 2 2 ≤ c ∇y n+1 2 2 ( ∇y n+1 2 + ∇y n 2 )+ 1 α f n 2 2 +ν ∇y n 2 2 .
But, (63) implies that for all n ∈ N

Ω |∇y n+1 | 2 ≤ 1 ν 1 α +∞ k=0 f k 2 2 + ν ∇y 0 2 2 := C ν and thus, since ν < 1 α y n+1 -y n 2 2 + ν ∇y n+1 2 2 ≤ 2cC 3/2 ν 3/2 + 2C ≤ C 1 ν 3/2
leading to y n+1 -y n 2 2 = O( 1 αν 3/2 ) as announced. This result asserts that two consecutive elements of the sequence {y n } n≥0 defined by recurrence from the scheme (2) are close to each other as soon as δt, the time step discretization, is small enough. In particular, this justifies the choice of the initial term y n+1 0 = y n of the minimizing sequence in order to approximate y n+1 .

We end this section devoted to the case d = 2 with the analogue of Theorem 2 for regular data.

Theorem 4 Suppose F ∈ L 2 (0, T ; L 2 (Ω) 2 ), u 0 ∈ V , for all n ∈ {0, • • • , N - 1},
α and f n are given by ( 4) and y n+1 ∈ V solution of (56

). If C(u 0 , F ) := F 2 L 2 (0,T ;L 2 (Ω) 2 ) +ν u 0 2 H 1 0 (Ω) 2 ≤ Cν 2
for some C and α is large enough, then, for any n ≥ 0, the minimizing sequence {y n+1 k } k∈N defined by (54) strongly converges to the unique of solution of (56).

Proof As for Theorem 2, it suffices to prove that there exists c 1 ∈ (0, 1) such that, for all n ∈ {0,

• • • , N -1} and all k ∈ N, τ 2 (y n k ) ≤ c 1 . Let us recall that for all n ∈ {0, • • • , N -1} and all k ∈ N y n+1 k+1 H 1 0 (Ω) 2 ≤ y n H 1 0 (Ω) 2 + m √ 2 ν(1 -c 1 ) E n (y n ) 1 -p n,0 ( λ n,0 )
where p n,0 ( λ n,0 ) is defined as in the proof of Proposition 4. From (53), since for all n ∈ {0,

• • • , N -1}, f n 2 2 ≤ α F 2 L 2 (0,T ;L 2 (Ω) 2 ) : E 0 (y 0 ) = E 0 (u 0 ) ≤ 1 2ν u 0 H 1 0 (Ω) 2 (ν + u 0 2 ) + ν α f 1 2 2 ≤ 1 ν u 0 2 H 1 0 (Ω) 2 (ν + u 0 2 ) 2 + F 2 L 2 (0,T ;L 2 (Ω) 2 )
and, since y n is solution of (56), then for all n ∈ {1, • • • , N -1} :

E n (y n ) ≤ 1 α f n -f n-1 2 2 + α y n -y n-1 2 2 ≤ 2 F 2 L 2 (0,T ;L 2 (Ω) 2 ) + α y n -y n-1 2 2 .
From the proof of Theorem 3, we deduce that for all n ∈ {0, • • • , N -1} :

α y n+1 -y n 2 2 ≤ 2cC(u 0 , F ) 3/2 ν 3/2 + 2C(u 0 , F ) and thus, for all n ∈ {1, • • • , N -1} E n (y n ) ≤ 2cC(u 0 , F ) 3/2 ν 3/2 + 4C(u 0 , F ).
Moreover, from (63), for all n ∈ {0, • • • , N -1} :

y n 2 H 1 0 (Ω) 2 ≤ 1 ν 1 α n k=0 f k 2 2 +ν u 0 2 H 1 0 (Ω) 2 ≤ 1 ν F 2 L 2 (Q T ) 2 +ν u 0 2 H 1 0 (Ω) 2 .
Eventually, let c 1 ∈ (0, 1). Then there exists α 0 > 0 such that, for all α ≥ α 0 c α,ν E n (y n ) ≤ K < 1. We therefore have (see Theorem 2), for all α ≥ α 0 , all n ∈ {0, • • • , N -1} and all k ∈ N :

y n+1 k+1 H 1 0 (Ω) 2 ≤ y n H 1 0 (Ω) 2 + m √ 2 ν(1 -c 1 ) E n (y n ) 1 -K which gives a bound of y n+1 k+1 H 1 0 (Ω) 2 independent of α ≥ α 0 . Taking α 1 ≥ α 0 large enough, we deduce that, for all α ≥ α 1 , all n ∈ {0, • • • , N -1} and all k ∈ N, y n k H 1 0 (Ω) 2 ≤ c 1 √ 2αν, that is τ 2 (y n k ) ≤ c 1 .
The announced convergence follows from Proposition 7.

Three dimensional case

We now consider regular solutions for the case d = 3. The following intermediate regularity result holds true. 3 and that u 0 ∈ V . Then any solution y n+1 ∈ V of (56) belongs to H 2 (Ω) 3 .

Proposition 10 Assume that Ω is C 2 , that (f n ) n is a sequence in L 2 (Ω)
If moreover, there exists C > 0 such that

1 α n k=0 f k 2 2 + ν Ω |∇u 0 | 2 ≤ Cν 3 , (67)
then the inequality (63) holds true.

Proof From Proposition 1, we know that for all n ∈ N * , y n ∈ H 2 (Ω) 3 ∩ V . Let now P be the operator of projection from L 2 (Ω) 3 into H. Taking w = P ∆y n+1 in (56) leads to :

α Ω |∇y n+1 | 2 + ν Ω |P ∆y n+1 | 2 = - Ω f n • P ∆y n+1 + Ω y n+1 • ∇y n+1 • P ∆y n+1 + α Ω ∇y n • ∇y n+1 . ( 68
)
In view of the inequality

Ω y n+1 • ∇y n+1 • P ∆y n+1 ≤ y n+1 3 ∇y n+1 6 P ∆y n+1 2 ,
we use again that there exist constants c 1 , c 2 > 0 such that

∆y n+1 2 ≤ c 1 P ∆y n+1 2 , ∇y n+1 6 ≤ c 2 ∆y n+1 2 ≤ c 1 c 2 P ∆y n+1 2 so that, for c = c 1 c 2 , we obtain Ω y n+1 • ∇y n+1 • P ∆y n+1 ≤ c y n+1 3 P ∆y n+1 2 2 .
It results from (68) that

α 2 Ω |∇y n+1 | 2 + ν 2 -c y n+1 3 Ω |P ∆y n+1 | 2 ≤ 1 2ν f n 2 2 + α 2 Ω |∇y n | 2 .
(69) Assume that, for all n ∈ N * , we have constructed by recurrence an element y n solution of (56

) such that ν 4 -c y n 3 > 0. (70)
Then, for all n ∈ N

Ω |∇y n+1 | 2 + ν 2α Ω |P ∆y n+1 | 2 ≤ 1 να f n 2 2 + Ω |∇y n | 2 (71)
and recursively, for all n ∈ N * , we get (63).

It remains to construct a sequence {y n } n∈N * solution of (56) and satisfying for all n ∈ N * the property (70). Let us first remark that the hypothesis (67) implies that y 0 satisfies (70). Let then n ∈ N fixed. Assume now, that we have constructed, for k ∈ {0, • • • , n} a solution y k satisfying (56

) if k ≥ 1 and ν 4 -c y k 3 > 0 for c = c 1 c 2 introduced above. Let y 1 ∈ V and let y 2 ∈ H 2 (Ω) 3 ∩ V be the unique solution of α Ω (y 2 -y n )•w+ν Ω ∇y 2 •∇w+ Ω y 1 •∇y 2 •w = f n , w H -1 (Ω) d ×H 1 0 (Ω) d , ∀w ∈ V .
If y 1 satisfies y 1 3 ≤ ν 4c , then in view of (69),

α 2 Ω |∇y 2 | 2 + ν 2 -c y 1 3 Ω |P ∆y 2 | 2 ≤ 1 2ν f n 2 2 + α 2 Ω |∇y n | 2
and consequently

α 2 Ω |∇y 2 | 2 + ν 4 Ω |P ∆y 2 | 2 ≤ 1 2ν f n 2 2 + α 2 Ω |∇y n | 2 .
(71) then implies

Ω |∇y 2 | 2 + ν 2α Ω |P ∆y 2 | 2 ≤ 1 ν 1 α n k=0 f k 2 2 + ν Ω |∇y 0 | 2 . (72)
We now use that there exists a constant c 3 > 0 such that, for all n ∈ N y 2 3 ≤ c 3 ∇y 2 2 to obtain

y 2 2 3 ≤ c 2 3 ν 1 α n k=0 f k 2 2 + ν Ω |∇u 0 | 2 .
Invoking assumption (67), we conclude that y 2 3 ≤ ν 4c . Eventually, we introduce the application T : C → C, y 1 → y 2 where C is the closed convex set of V defined by

C := {y ∈ V , ν 4c ≥ y 3 }. Let us check that T is continuous. Let y 1 ∈ C et z 1 ∈ C, y 2 = T (y 1 ) et z 2 = T (z 1 ) so that α Ω (z 2 -y 2 )•w+ν Ω ∇(z 2 -y 2 )•∇w+ Ω y 1 •∇(y 2 -z 2 )•w+ Ω (y 1 -z 1 )•∇z 2 •w = 0 for all w ∈ V and then for w = z 2 -y 2 α Ω |z 2 -y 2 | 2 + ν Ω |∇(z 2 -y 2 )| 2 ≤ Ω (y 1 -z 1 ) • ∇z 2 .(z 2 -y 2 ) ≤ c ∇(y 1 -z 1 ) 2 ∇z 2 2 z 2 -y 2 3
≤ c ∇(y 1 -z 1 ) 2 using (72); this implies the continuity of T . On the other hand, since T (C) is a bounded set of H 2 (Ω) 3 , T is relatively compact. The Schauder Theorem allows to affirm that T has a fixed point y ∈ C, that is, a solution y n+1 ∈ C of (56).

Remark 4 Under the hypothesis of Proposition 10, suppose moreover that

C α,ν := ν -5/2 α -1/2 α -1 n-1 k=0 f k 2 2 + ν ∇y 0 2 2
is small enough, then, the solution of (56) is unique. Indeed, let n ∈ N and let y n+1

1 , y n+1 2 ∈ V be two solutions of (56). Let Y := y n+1 1 -y n+1 2 . Then, α Ω Y •w +ν Ω ∇Y •∇w + Ω y n+1 2 •∇Y •w + Ω Y •∇y n+1 1 •w = 0 ∀w ∈ V and in particular, for w = Y (since Ω y n+1 2 • ∇Y • Y = 0) α Ω |Y | 2 + ν Ω |∇Y | 2 = - Ω Y • ∇y n+1 1 • Y = Ω Y • ∇Y • y n+1 1 ≤ c y n+1 1 ∞ ∇Y 2 Y 2 ≤ c ∇y n+1 1 1/2 2 P ∆y n+1 1 1/2 2 ∇Y 2 Y 2 ≤ α Y 2 2 + c α ∇y n+1 1 2 P ∆y n+1 1 2 ∇Y 2 2
and therefore (ν -c α ∇y n+1

1 2 P ∆y n+1 1
2 ) ∇Y 2 2 ≤ 0. Moreover, from (63),

∇y n+1 1 2 P ∆y n+1 1 2 ≤ 2α 1/2 ν 3/2 1 α n k=0 f k 2 2 + ν ∇y 0 2 2 .
Therefore, there exists a constant c > 0 such that if C α,ν < c, then, arguing as in the 2D case, ∇Y 2 2 ≤ 0 and Y = 0.

As in the 2D case, Proposition 10 then allows, following the proof of Theorem 3, to obtain an estimate of y n+1 -y n 2 in term of the parameter α.

Theorem 5 Assume that Ω is C 2 , that {f n } n∈N is a sequence in L 2 (Ω) 3 satisfying α -1 +∞ k=0 f k 2 2 < +∞.
Assume moreover that y 0 ∈ V and that for all n ∈ N, y n+1 ∈ H 2 (Ω) 3 ∩ V is a solution of (56) satisfying (63). Then, the sequence (y n ) n satisfies (66).

Eventually, adapting the proof of Theorem 4, we get the following convergence result.

Theorem 6 Suppose F ∈ L 2 (0, T ; L 2 (Ω) 3 ), y 0 ∈ V , for all n ∈ {0, • • • , N -1}, α and f n are given by ( 4) and y n+1 ∈ V solution of (56

). If C(y 0 , F ) := F 2 L 2 (0,T ;L 2 (Ω) 3 ) + ν y 0 2 H 1 0 (Ω) 3 ≤ Cν 3
for some C > 0 and α is large enough, then for any n ≥ 0, the minimizing sequence {y n+1 k } k∈N defined by (54) strongly converges to the unique of solution of (56).

Remark 5 We have considered regular solutions in the case d = 3 in order to be able to prove the uniform property τ 3 (y n k ) ≤ c 1 < 1 for some c 1 independent of k and n, i.e. y n k H 1 0 (Ω) 3 ≤ c 1 M -1 (αν 3 ) 1/4 . Actually, for regular solutions, Proposition 10 implies that y n k H 1 0 (Ω) 3 ≤ C for some C independent of α, which is sufficient for α large enough. By considering weak solutions, we can only prove that y n k H 1 0 (Ω) 3 ≤ Cα 1/2 for some C > 0 (see (58)) which does not imply τ 3 (y n k ) ≤ c 1 < 1.

Numerical illustrations

We discuss in this section numerical experiments based on finite element approximations in space for two geometries of R 2 : the celebrated channel with a backward facing step and the semi-circular driven cavity introduced in [START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF].

In both cases, the velocity of the fluid is imposed on the boundary. We first start with the case α = 0 in (14) (discussed in [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF]) allowing, first to get the solution of (1) as time becomes large and secondly, to enhance the gain of the optimization of the descent step parameter λ k in (30). Then, for the semi driven cavity, we consider the cases α = 0 and α > 0 applied to the resolution of the backward Euler scheme (51). In a final part, we briefly compare the computational cost of this least-squares approach with standard explicit and semi-explicit scheme.

The numerical simulations are performed with the FreeFem++ package (see [START_REF] Hecht | New development in freefem++[END_REF]). Regular triangular meshes are used together with the P 2 /P 1 Taylor-Hood finite element, satisfying the Ladyzenskaia-Babushka-Brezzi condition of stability.

Steady case : Two dimensional channel with a backward facing step

We consider in the steady situation the test problem of a two-dimensional channel with a backward facing step, described for instance in Section 45 of [START_REF] Glowinski | Finite element methods for incompressible viscous flow[END_REF] (see also [START_REF] Jiang | Least-squares finite element method for fluid dynamics[END_REF]). The geometry is depicted in Figure 1. Dirichlet conditions of the Poiseuille type are imposed on the inflow and outflow sides Γ 1 and Γ 2 of the channel: we impose y = (4(H -x 2 )(x 2 -h)/(H -h) 2 , 0) on Γ 1 and y = (4(H -h)x 2 (H -x 2 )/H 2 , 0) on Γ 2 , with h = 1, H = 3, l = 3 and L = 30. On the remaining part ∂Ω \ (Γ 1 ∪ Γ 2 ), the fluid flow is imposed to zero. The external force f is zero.

Γ 1 Γ 2 x 2 x 1 (l, 0) (L, 0) (0, h) (0, H)

Ω

Fig. 1 A two-dimensional channel with a step.

We consider the extremal problem [START_REF] Saramito | A damped Newton algorithm for computing viscoplastic fluid flows[END_REF] to solve the steady Navier-Stokes equation ( 3) with here α = 0. We compare the descent algorithm (30) based on the descent direction Y 1 with the conjugate gradient (CG) algorithm used in [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF]. In both cases, the initial guess is defined as the solution of the corresponding Stokes problem. Moreover, the scalar extremal problem with respect to λ k in (30) is performed with the Newton-Rasphon method.

We start with a large value of ν = 1/150. Table 1 reports the evolution of the quantity y k+1 -y k H 1 0 (Ω) 2 / y k H 1 0 (Ω) 2 with respect to the iterate k associated to the algorithms (30), (30) with fixed step λ k = 1 and CG respectively. We also consider the so-called by analogy damped quasi newton method

       y 0 ∈ V given, y k+1 = y k -λ k Y 1,k , k ≥ 0, λ k = argmin λ∈[0,m] E(y k -λ Y 1,k ) (73) with Y 1,k ∈ V the solution of the formulation α Ω Y 1,k • w + ν Ω ∇ Y 1,k • ∇w + Ω (y 0 • ∇ Y 1,k + Y 1,k • ∇y 0 ) • w = -α Ω v k • w -ν Ω ∇v k • ∇w, ∀w ∈ V
and v k ∈ V the corrector (associated to y k ) solution of [START_REF] Smith | Implicit algorithms and their linearization for the transient incompressible Navier-Stokes equations[END_REF].

A regular mesh composed of 14 143 triangles and 7 360 vertices is used. Table 2 reports the evolution of the norm of the corrector v k H 1 0 (Ω) 2 , an upper bound of y -y k H 1 0 (Ω) 2 , according to Proposition 5. As expected in view of the results in Section 2.2, the descent algorithm (30) based on Y 1,k is much faster than the CG algorithm. Moreover, the optimal values for the optimal step λ k are close to one, so that the Newton method provides a similar speed of convergence. As the norm of Y 1,k goes to zero with k, the term factor of λ 2 in (33) gets small, and the optimal λ k gets close to one. Remark as well that the algorithm (73) offers an excellent speed of convergence. In term of CPU times, algorithms (30) and (73) require about 53 and 108 seconds respectively and leads to the same approximation. We have notably For smaller values of ν, the results are qualitatively differents. Table 3 reports some norms with respect to k for ν = 1/700. We observe, from the last column, that the Newton method for which λ k is fixed to one does not converge anymore. Actually, the Newton method, when initialized with the solution of the corresponding Stokes problem, diverges for ν ≤ 1/250. On the other hand, the optimization of the step λ k produces a very fast convergence of the sequence {y k } (k>0) . Observe here that the values for the optimal λ k are not close to one, during the first iterations. We obtain notably ∇ • y L 2 (Ω) /|Ω| ≈ 5.78 × 10 -2 . In agreement with Theorem 1, we observe from Table 3 that the decrease of E(y k ) to zero is first linear and then becomes quadratic.

The algorithm (73) is a bit more robust than the Newton one as it converges for all ν satisfying ν ≥ 1/290 approximately. Finally, as discussed in [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF], the CG algorithm converges and produces similar numerical values: the convergence is however much slower since about 350 iterates are needed to achieve 2E(y k ) of the order 10 -3 .

The algorithm (30) requires however the initial guess to be close enough to the solution. Initialized with the solution of the corresponding Stokes problem, it diverges for ν ≤ 1/720. A continuation method with respect to ν is then necessary in that case. Algorithm (30) is also robust with respect to the mesh size: with a twice finer mesh composed of 84 707 triangles and 43 069 vertices, the convergence y k+1 -y k H 1 0 (Ω) 2 ≤ 10 -12 y k H 1 0 (Ω) 2 is observed after k = 18 iterates (instead of 14 for the coarser mesh) leading notably to

∇ • y L 2 (Ω) /|Ω| ≈ 3.91 × 10 -2 . iterate k y k+1 -y k H 1 0 (Ω) 2 y k H 1 0 (Ω) 2 2E(y k ) λ k 2E(y k ) with λ k = 1 1

Steady case: 2D semi-circular cavity

We now consider the test discussed in [START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF]. The geometry is a semi-disk

Ω = {(x 1 , x 2 ) ∈ R 2 , x 2 1 + x 2 2 < 1/4, x 2 ≤ 0} depicted on Figure 2. The velocity is imposed to y = (g, 0) on Γ 0 = {(x 1 , 0) ∈ R 2 , |x 1 | < 1/2} with g vanishing at
x 1 = ±1/2 and close to one elsewhere: we take g(x 1 ) = (1 -e 100(x1-1/2) )(1e -100(x1+1/2) ). On the rest Γ

1 = {(x 1 , x 2 ) ∈ R 2 , x 2 < 0, x 2 1 + x 2 2 = 1/4} of the boundary the velocity is fixed to zero. (-1 2 , 0) ( 1 2 , 0) Γ 0 : y = (1, 0) Γ 1 : y = (0, 0) Fig. 2 Semi-disk geometry.
For a regular triangular mesh, composed of 79 628 triangles and 40 205 vertices, leading to a mesh size h ≈ 6.23 × 10 -3 , the Newton method (λ k = 1) initialized with the corresponding Stokes solution, converges up to ν -1 ≈ 500. On the other hand, the algorithm (30) still converges up to ν -1 ≈ 910. Figures 5 depicts the streamlines of the steady state solution corresponding to ν -1 = 500 and to ν -1 = i × 10 3 for i = 1, • • • , 9. The values used to plot the stream function are given in Table 6. The figures are in very good agreements with those depicted in [START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF]. The solution corresponding to ν -1 = 500 is obtained from the sequence given (30) initialized with the Stokes solution. Eight iterates are necessary to achieve 2E(y) ≈ 3.4 × 10 -17 . The stopping criterion is y k+1 -y k H 1 0 (Ω) 2 ≤ 10 -12 y k H 1 0 (Ω) 2 . Tables 4 and5 collect some values for ν = 1/500 and ν = 1/700. Then, the other solutions are obtained by a continuation method with respect to ν taking δν -1 = 500. For instance, the solution corresponding to ν -1 = 5000 is obtained from the algorithm (30) initialized with the steady solution corresponding to ν -1 = 4500. Table 9 reports the history of the continuation method and highlights the efficiency of the algorithm (30): up to ν -1 = 9500, few iterations achieve the convergence of the minimizing sequence {y k } k∈N . From ν -1 = 10 4 , with a finer mesh (for which the mesh size is h ≈ 4.37 × 10 -3 ), δν is reduced to δν -1 = 100 and leads to convergence beyond ν -1 = 15000. Table also reports the minimal values of the streamline function ψ which compare very well with those of [START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF].

The case α > 0 leads to similar results, in full agreement with the theoretical Section 2. For ν = 1/1000, Table 7 reports results of the algorithm (30) for α ∈ {10 -1 , 1., 10., 100.}. As expected, the gain of coercivity of the functional E involves a notable robustness and speed up of the algorithm. Recall that for α = 0. and ν = 1/1000, algorithm (30) does not converge. For α = 10 -1 , we observe the convergence after 8 iterates. For a fixed value of ν, this number of iterates decreases as α gets larger. Actually, we observe that when α -0.07, -0.0675, -0.065, -.05, -0.04, -0.03, -0.02, -0.01, ±10 -4 , ±10 -5 ±10 -7 , -10 -10 , 0., 10 -8 , 10 -6 , 5 × 10 -4 , 10 -3 , 2 × 10 -3 , 3 × 10 -3 , 4 × 10 -3 5 × 10 -3 , 6 × 10 -3 , 7 × 10 -3 , 8 × 10 -3 , 9 × 10 -3 , 0.01 Table 6 Values used to plot the contours of the stream function.

√ ν = O(
to 1 for all k. Moreover, when αν = O(1), the convergence is achieved after one iterate only. This behavior suggests that one may recover the solution of the steady Navier-Stokes system (corresponding to α = 0) by using a continuation procedure with respect to the parameter α decreasing to zero. We easily check, for any α ≥ 0, the estimate ∇(y α -y α=0 ) 2 ≤ c p α ν y α=0 2 where y α solves (5) with g = 0, and c p the Poincaré constant. For ν = 1/5000, Table 8 reports the history of the continuation approach starting from α = 1. to α = 0 with intermediate steps α ∈ {10 -1 , 10 -2 , 10 -3 , 10 -4 }. Figure 4 depicts the evolution of the sequences { 2E(y k )} k and {λ k } k obtained from the algorithm (30) with α = 10 -3 . The algorithm is initialized with the solution corresponding to α = 10 -2 . The behavior of these sequences fully illustrates Theorem 1 and the robustness of the method: as λ k increases, the decay of 2E(y k ), initially low, gets larger and becomes very fast when λ k is close to one. Figure 5

Unsteady case: 2D semi-circular cavity

We now use the least-squares method in order to solve iteratively the implicit Euler scheme [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF]. The parameter α = 1/δt is strictly positive. We remind that for ν approximatively larger than 1/6600, the unsteady solution converges as time evolves to the steady solution (corresponding to α = 0) obtained in the previous section by a continuation technique. Actually, the iterative process due to the time discretization can also be seen as a continuation approach. We consider the value ν = 1/1000. Following [START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF], we take as initial condition u 0 the steady-state solution corresponding to ν = 1/500. Before comparing with standard time marching schemes, let us make a comment on the algorithm (55). For α large, the optimal step λ k in (55) equals one and the convergence with respect to k is achieved after one iterate. The convergent approximation y n+1 := y n+1 1 then simply solves, for each n to solve the fully implicitly Euler scheme associated to the unsteady Navier-Stokes equation. When the time discretization is fine enough, each step of the damped Newton method simply reduced to the Newton one. In such a case, we obtained a proof of convergence of the Newton scheme to solve the unsteady Navier-Stokes. As far as we know, this proof is original. Numerical experiments have highlighted the robustness of the method, including for values of the viscosity coefficient of order 10 -4 . We also emphasize that the least-squares approaches, employed here to treat the Navier-Stokes nonlinearity, can be used to solve other nonlinear equations, as formally done in [START_REF] Münch | Numerical null controllability of the heat equation through a least squares and variational approach[END_REF] for a sublinear heat equation. Eventually, we may solve the unsteady Navier-Stokes system by a fully L 2 (0, T ; H -1 (Ω)) least-squares approach. The underlying corrector solves an unsteady Stokes type equation; we refer to [START_REF] Lemoine | A fully space-time least-squares method for the unsteady Navier-Stokes system[END_REF].

  Fig. 3 Streamlines of the steady state solution for ν -1 = 500, 1000, 2000, 3000, 4000, 5000, 6000, 7000 and ν -1 = 8000.

Fig. 4

 4 Fig. 4 Semi-disk geometry; ν = 1/5000; Evolution of 2E(y k ) and λ k w.r.t. k for α from 10 -2 to 10 -3 .

Fig. 5

 5 Fig. 5 Streamlines of the α-steady state solution for α = 1., 10 -1 , 10 -2 (Top) and α = 10 -3 , 10 -4 and 0. (Bottom).

2 yα -y α=0 H 1 0

 21 (Ω)

  For the value α = 200 corresponding to the time discretization parameter δt = 5 × 10 -3 , we observe the convergence of the sequence {y n+1 k } k>0 after at most three iterations, for each n (except for n = 0 requiring 6 iterations). For the value α = 2000 corresponding to δt = 5 × 10 -4 (used in [5]), we observe the convergence of the sequence after one iterate. At time T = 10, the unsteady state solution is close to the solution of the steady Navier-Stokes equation: we compute that the sequence {y n } n=0,...,2000 satisfies y 2000 -y 1999 L 2 (Ω) / y 2000 L 2 (Ω) ≈ 1.19 × 10 -5 . n = 2000 corresponds to T = 2000 × δt = 10. Figures 6 display the streamlines of the unsteady state solution corresponding to ν = 1/1000 at time 0, 1, 2, 3, 4, 5, 6 and 7 seconds to be compared with the streamlines of the steady solution depicted in Figure 5. These figures are in full agreement with [5].

Fig. 6

 6 Fig. 6 Streamlines of the unsteady state solution for ν -1 = 1000 at time t = i, i = 0, • • • , 8s.

  ∇ • y L 2 (Ω) /|Ω| ≈ 1.83 × 10 -4 .

	iterate k	(30) with λ k = 1		(30)	(73)	CG
	1		4.44 × 10 -1	3.79 × 10 -1 (0.854)	3.79 × 10 -1	5.21 × 10 -2
	2		1.95 × 10 -1	1.81 × 10 -1 (0.957)	1.59 × 10 -1	4.19 × 10 -2
	3		5.60 × 10 -2	4.04 × 10 -2 (0.994)	4.37 × 10 -2	3.27 × 10 -2
	4		3.98 × 10 -3	2.22 × 10 -3 (1.001)	6.05 × 10 -3	2.94 × 10 -2
	5		2.08 × 10 -5	5.71 × 10 -6 (0.999)	6.80 × 10 -3	2.56 × 10 -2
	6		5.91 × 10 -10		4.95 × 10 -11 (1)	9.89 × 10 -4	2.29 × 10 -2
	7		4.88 × 10 -15		3.29 × 10 -15 (1)	9.00 × 10 -4	2.21 × 10 -2
	8		-		-	1.48 × 10 -4	2.02 × 10 -2
	9		-		-	9.55 × 10 -5	1.95 × 10 -2
	10		-		-	2.09 × 10 -5	1.81 × 10 -2
	11		-		-	1.39 × 10 -5	1.76 × 10 -2
	12		-		-	3.17 × 10 -6	1.72 × 10 -2
	13		-		-	1.83 × 10 -6	1.67 × 10 -2
	14		-		-	3.80 × 10 -7	1.65 × 10 -2
	26		-		-	4.32 × 10 -13	1.12 × 10 -2
	200		-		-	-	2.09 × 10 -5
	Table 1 2D channel geometry; ν = 1/150; Evolution of y k+1 -y k H 1 0 (Ω) 2 / y k H 1 0 (Ω) 2
	with respect to k.
	iterate k	(30) with λ k = 1	(30)	(73)	CG
	1		5.46 × 10 -2		5.46 × 10 -2	5.47 × 10 -2	5.46 × 10 -2
	2		2.39 × 10 -2		2.22 × 10 -2	2.22 × 10 -2	3.70 × 10 -2
	3		4.95 × 10 -3		4.60 × 10 -3	5.45 × 10 -3	2.91 × 10 -2
	4		3.20 × 10 -4		1.56 × 10 -4	9.32 × 10 -4	2.49 × 10 -2
	5		1.53 × 10 -6		5.43 × 10 -7	5.19 × 10 -4	2.20 × 10 -2
	6		3.65 × 10 -11		4.22 × 10 -12	1.71 × 10 -4	1.99 × 10 -2
	7		6.54 × 10 -16		2.54 × 10 -16	1.71 × 10 -4	1.84 × 10 -2
	8		-		-	7.85 × 10 -5	1.70 × 10 -2
	9		-		-	2.47 × 10 -5	1.60 × 10 -2
	10		-		-	8.95 × 10 -6	1.51 × 10 -2
	11		-		-	3.42 × 10 -6	1.43 × 10 -2
	12		-		-	1.20 × 10 -6	1.36 × 10 -2
	13		-		-	4.25 × 10 -7	1.30 × 10 -2
	14		-		-	1.36 × 10 -7	1.24 × 10 -2
	26		-		-	1.59 × 10 -14	6.25 × 10 -3
	200		-		-	-	1.55 × 10 -5

Table 2

 2

	2D channel geometry; ν = 1/150; Evolution of v k H 1 0 (Ω) 2 =	2E(y k ) with
	respect to k.

Table 3

 3 2D channel geometry; ν = 1/700; Results for the algorithm (30).

  [START_REF] Bochev | Least-squares finite element methods[END_REF], algorithm (30) converges after few iterates with λ k close

	iterate k 1	y k+1 -y k H 1 0 (Ω) 2 y k H 1 0 (Ω) 2 4.66 × 10 -1	2E(y k ) 8.51 × 10 -3	λ k 0.810	2E(y k ) with λ k = 1 8.51 × 10 -3
	2	2.03 × 10 -1	3.57 × 10 -3	0.716	4.09 × 10 -3
	3	9.54 × 10 -2	1.36 × 10 -3	0.860	2.26 × 10 -3
	4	2.64 × 10 -2	3.38 × 10 -4	0.982	2.51 × 10 -4
	5	1.53 × 10 -3	2.43 × 10 -5	1.002	3.47 × 10 -6
	6	4.21 × 10 -6	6.36 × 10 -8	1.	1.54 × 10 -9
	7	4.20 × 10 -11	6.50 × 10 -13	1.	1.52 × 10 -16
	8	3.25 × 10 -15	3.45 × 10 -17	1.	3.39 × 10 -17

Table 4

 4 Semi-disk geometry ; ν = 1/500; Results for the algorithm (30).

	iterate k 1	y k+1 -y k H 1 0 (Ω) 2 y k H 1 0 (Ω) 2 4.89 × 10 -1	2E(y k ) 8.51 × 10 -3	λ k 0.702	2E(y k ) with λ k = 1 8.51 × 10 -3
	2	2.41 × 10 -1	4.43 × 10 -3	0.583	5.75 × 10 -3
	3	1.48 × 10 -1	2.15 × 10 -3	0.510	5.27 × 10 -3
	4	1.01 × 10 -1	1.15 × 10 -3	0.568	1.23 × 10 -2
	5	6.83 × 10 -2	6.50 × 10 -4	0.931	4.40 × 10 -3
	6	1.19 × 10 -2	1.76 × 10 -4	1.018	1.52 × 10 -2
	7	3.35 × 10 -4	3.42 × 10 -6	0.999	5.06 × 10 -3
	8	3.19 × 10 -7	4.22 × 10 -9	1.	1.45 × 10 -2
	9	2.97 × 10 -13	2.30 × 10 -15	1.	1.85 × 10 -2
	20	-	-	-	5.39

Table 5 Semi

 5

-disk geometry ; ν = 1/700; Results for the algorithm (30).

Table 7

 7 reports the streamlines of the solution for various values of α. 10 -3 (0.65) 3.11 × 10 -3 (0.87) 1.04 × 10 -3 (0.99) 2.99 × 10 -4 (0.99) 2 3.54 × 10 -3 (0.63) 7.59 × 10 -4 (0.98) 3.68 × 10 -5 (0.99) 7.81 × 10 -7 (1.) 3 1.55 × 10 -3 (0.64) 5.94 × 10 -5 (1.00) 5.15 × 10 -8 (1.) 6.92 × 10 -12 (1.) 4 7.05 × 10 -4 (0.74) 4.84 × 10 -7 (1.) 1.20 × 10 -13 5.50 × 10 -17 5 2.76 × 10 -4 (1.01) 2.63 × 10 -11 (1.) Semi-disk geometry ; ν = 1/1000; Results for the algorithm (30); 2E(y k ) and λ k with respect to α.

	k	α = 0.1	α = 1.	α = 10.	α = 100.
	1	6.59 × -	-
	6	2.95 × 10 -5 (1.)	1.03 × 10 -17	-	-
	7	1.55 × 10 -7 (1.)	-	-	-
	8	1.99 × 10 -11 (1.)	-	-	-
	9	1.60 × 10 -17	-	-	-

  2

	Stokes → 1.	6	0.00684	7.659	4.28
	1. → 10 -1	37	0.011	7.4896	3.32
	10 -1 → 10 -2	49	0.0297	7.8951	2.67
	10 -2 → 10 -3	70	0.0360	8.0166	4.02 × 10 -1
	10 -3 → 10 -4	12	0.0366	8.0278	4.03 × 10 -2
	10 -4 → 0	4	0.0367	8.029	4.03 × 10 -3

Table 8

 8 Semi-disk geometry; Continuation method with respect to α for the solution of the steady Navier-Stokes equation; ν = 1/5000.

	ν -1 Stokes → 500	it. 7	π L 2 (Ω) 4.31 × 10 -2	y 1 H 1 0 (Ω) 2 4.462	y 2 H 1 0 (Ω) 2 2.489	min Ω ψ -0.0766	min Ω ψ [5] -
	500 → 1000	7	4.07 × 10 -2	4.919	2.883	-0.0780	-0.0779
	1000 → 1500	6	3.99 × 10 -2	5.296	3.153	-0.0775	-
	1500 → 2000	6	3.93 × 10 -2	5.612	3.361	-0.0766	-0.0763
	2000 → 2500	5	3.88 × 10 -2	5.884	3.531	-0.0756	-
	2500 → 3000	5	3.83 × 10 -2	6.126	3.675	-0.0744	-0.0742
	3000 → 3500	5	3.79 × 10 -2	6.345	3.801	-0.0733	-
	3500 → 4000	5	3.75 × 10 -2	6.545	3.911	-0.0721	-
	4000 → 4500	6	3.71 × 10 -2	6.731	4.010	-0.0710	-
	4500 → 5000	6	3.67 × 10 -2	6.903	4.099	-0.0699	-0.0700
	5000 → 5500	6	3.64 × 10 -2	7.065	4.181	-0.0689	-
	5500 → 6000	6	3.60 × 10 -2	7.217	4.256	-0.0679	-
	6000 → 6500	6	3.57 × 10 -2	7.362	4.325	-0.0669	-
	6500 → 7000	6	3.53 × 10 -2	7.499	4.389	-0.0660	-
	7000 → 7500	5	3.50 × 10 -2	7.631	4.448	-0.0651	-
	7500 → 8000	5	3.47 × 10 -2	7.756	4.504	-0.0643	-
	8000 → 8500	6	3.44 × 10 -2	7.876	4.557	-0.0634	-
	8500 → 9000	6	3.41 × 10 -2	7.992	4.606	-0.0626	-

Table 9

 9 Semi-disk geometry; Continuation method with respect to ν for the solution of the steady Navier-Stokes equation; α = 0.

the following semi-implicit scheme (mentioned in [15, section 13.4])

For δt of the order 10 -3 , this first-order in time scheme displays, for ν = 1/1000 and the semi-disk geometry of Figure 2, very similar results than the conditionally stable partially explicit scheme (we refer to [START_REF] Temam | Theory and numerical analysis[END_REF]Section 5.1])

for all w ∈ V and than the unconditionally stable scheme

(76) for all w ∈ V . In term of computational cost, the scheme (75) is as expected faster than the scheme (74): the ratio of the computational time to perform 10000 iterates (leading to T = 10) between ( 74) and ( 75) is approximatively equal to 1.65. The regular triangulation used corresponds to a mesh size h of the order of 6.23 × 10 -3 making (75) stable. On the other hand, the computational times of (74) and ( 76) are equivalent: we observe a ratio equal to 1.05. For δt = 10 -2 , scheme (75) is unstable. The convergence with respect to k of (55) (not anymore equivalent to (74)) is observed after two iterates. The ratio of computational time between (55) and (76) raises to 1.89. We observe that the approximation for (55) is much less sensitive to the variation of δt: we observe similar results than with δt = 10 -3 and than [START_REF] Glowinski | Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity[END_REF] where δt = 5 × 10 -4 is used.

Conclusions and perspectives

We have rigorously analyzed a weak least-squares method introduced forty years ago in [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF] allowing to solve a steady nonlinear Navier-Stokes equation, in the incompressible regime. This equation with a zero order term appears after any fully implicit time discretization of the unsteady Navier-Stokes equation. We have constructed a sequence converging strongly to the solution of the steady equation. Using a particular descent direction very appropriate for the analysis, this convergent sequence turns out to coincide with the sequence obtained using the damped Newton method to solve the underlying variational weak formulation. This globally convergent approach enjoys a quadratic rate of convergence after a finite number of iterates and is in particular much faster than the conjugate gradient method used in [START_REF] Bristeau | On the numerical solution of nonlinear problems in fluid dynamics by least squares and finite element methods. I. Least square formulations and conjugate gradie[END_REF], [START_REF] Lemoine | Analysis of continuous H -1 -least-squares approaches for the steady Navier-Stokes system[END_REF]. Then, we have shown the convergence of the method, uniformly with respect to the time discretization,