G Girault

Y Guevel

J M Cadou

An algorithm for the computation of multiple Hopf bifurcation points based on Pade approximants

Recently, a numerical method was proposed to compute a Hopf bifurcation point in fluid mechanics. This numerical method associates a bifurcation indicator and a Newton method. The former gives initial guesses to the iterative method. These initial values are the minima of the bifurcation indicator. However, sometimes, these minima do not lead to the convergence of the Newton method. Moreover, as only a single initial guess is obtained for each computation of the indicator, the computational time to obtain a Hopf bifurcation point can be quite long. The present algorithm is an enhancement of the previous one. It consists in automatically computing several initial guesses for each indicator curve. The majority of these initial values leads to the convergence of the Newton method. This method is evaluated through the problem of the lid-driven cavity with several aspect ratios in the framework of the finite element analysis of the 2D Navier-Stokes equations. The results prove the efficiency and the robustness of the proposed algorithm.

INTRODUCTION

This paper focuses on the numerical predictions of a Hopf bifurcation for incompressible viscous flows. Such instability occurs when a stationary flow becomes a time-periodic one.

From a mathematical point of view, if ¢ stands for a complex eigenvalue of the Jacobian matrix, then a Hopf bifurcation occurs when the real part of two eigenvalues becomes null for a critical Reynolds number, noted Re c .

The numerical methods used to compute accurately such instabilities are divided into two categories: the direct methods and the indirect ones with the monitoring of an indicator. The latter consists, for example, in computing the eigenvalues of Jacobian matrix and determining when a pair of complex conjugated values crosses the imaginary axis [START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF][START_REF] Fortin | Localization of Hopf bifurcations in fluid flow problems[END_REF]. In [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF], the authors have proposed a bifurcation indicator which has the property to be null at a Hopf bifurcation point. This method gives accurate values of a bifurcation point, but it requires much calculi and, consequently, much central processing unit (CPU) time. Let us remark that a stationary bifurcation indicator also exists and has been successfully used to study the Coanda effect [START_REF] Allery | Application of the asymptotic numerical method to the Coanda effect study[END_REF].

In the direct methods [START_REF] Jepson | Numerical Hopf bifurcation[END_REF][START_REF] Jackson | A finite-element study of the onset of vortex shedding in flow past variously shaped bodies[END_REF], a Newton algorithm is used to iteratively determine the solution of a nonlinear problem which exactly satisfies a Hopf bifurcation point. The drawback of this approach is that the convergence of the iterative method is governed by the quality of the computed initial guesses. Recently, in [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], the authors proposed to link the indicator proposed in [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF] and the Newton method. The principle of this method, which can be denoted as a hybrid method, is that the bifurcation indicator automatically provides an initial value to the iterative method. This initial guess is the minimum value of the indicator which is computed for a fixed value of the Reynolds number and only depends on the frequency. Although the majority of these initial values permits the computing of a Hopf bifurcation point, it does happen that, for some candidates, the Newton method fails to converge. As the bifurcation indicator only gives a single value for the Newton method per computation, a significant amount of computing time is then needed without assuring that this value is a good initial guess.

The purpose of the study is to develop an improved algorithm which allows one to compute Hopf bifurcation points in an automatic and robust way. As in [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], this algorithm combines both indirect and direct methods. Thus, the indirect method helps isolate several admissible initial values, and the direct method efficiently performs a Newton calculation. The initial values are determined as minima of an indicator curve, noted ¥ Re .¨/, after the computation of which the first and second derivatives respect the pulsation ¨of the flow. These indicator curves are determined by using a perturbation method (the asymptotic numerical method [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]) associated to rational fractions called Padé approximants [START_REF] Baker | Padé approximants: Volume 59 of encyclopaedia of mathematics and its applications[END_REF][START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF]. As these indicator curves are analytical functions of the angular frequency, their minima are easily determined by computing the roots of the numerator of the Padé approximants.

The proposed method permits one to give more than one initial value per indicator computation. This lets one increase the number of possible bifurcation points without increasing the computational times. Compared to the previous study [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF] and other existing algorithms, the presented method is characterized by its robustness and reliability.

This paper is divided into three main parts. The first part is devoted to theoretical aspects for incompressible viscous flow and stability analysis. The second part of the paper recalls the fundamentals of asymptotical numerical methods. The proposed algorithm, based on the use of an indicator curve and its derivative functions, is exposed. Finally, some classical results for the academic problem of the lid-driven cavity are presented. They show the efficiency and the accuracy of the studied method.

THEORY

Equations to be solved

The movement of a laminar incompressible flow in two dimensions is described by the Navier-Stokes equations:

8 < : P u i u i,jj C u j u i,j C 1 p ,i D 0 in u i,i D 0 in u D u d on @ u (1)
In Equation ( 1), u i , p, ¡ and stand, respectively, for the velocity in the ith-direction, the pressure, the fluid density and the kinematic viscosity. The notations . / and ( ,i ) stand for the derivatives with respect to time and space, respectively, in the direction i. The term is a control parameter indicating the intensity of the imposed velocity and can be identified with the Reynolds number Re D uL=v using a reference length L. Equation (1) can be written in the following form:

M. P U / C L.U / C Q.U , U / F D 0 in (2)
where U designates a concatenated vector containing the velocity u and the pressure p, and M is a mass matrix. L.U / is a linear operator containing the pressure and diffusion terms, whereas the quadratic operator Q.U ,U / contains the convective terms of the Navier-Stokes equations. The term F stands for a load vector proportional to the imposed velocity on @ u . The reader can find explicit descriptions of these operators in [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF].

The stability of the flow is studied by adding a perturbation term V(t/ to the stationary solution U S of the Navier-Stokes Equation [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF]. For a Hopf bifurcation, this perturbation term is written as follows:

_ .t / D V .x/e i¨t (3)
In expression (3), ¨and V stand, respectively, for the pulsation and the complex mode of the perturbation. The vector x designates the spatial coordinates. Introducing term Equation (3) in the Navier-Stokes Equation (2) and neglecting quadratic and higher order terms in V, the following system is obtained:

8 < : i¨M.V / C L.V / C Q.V , U S / C Q.U S , V / D 0 in V U D 0 on @ U L.U S / C Q.U S , U S / S F D 0 in (4)
The solution of system (4) is the vector ƒ D ¹V a , V b , ¨, U S º, which corresponds to a Hopf bifurcation for the studied stationary flow. The notations V a and V b stand, respectively, for the real and imaginary part of the mode V . The quantity V U represents the velocity part of the mixed unknown vector V .

Newton method

The previous nonlinear system of Equation ( 4) can be written in the following form:

² R.ƒ/ D 0 in V U D 0 on @ U ( 5
)
where ƒ is the unknown vector. This solution is sought under the following expression:

ƒ D ƒ 0 C ƒ (6)
ƒ 0 is the given initial value, and ƒ is the unknown increment. By introducing the previous relation into Equation [START_REF] Fortin | Localization of Hopf bifurcations in fluid flow problems[END_REF] and by neglecting the second order term in ƒ, one obtains the linear system to be solved at each iteration:

² L ƒ 0 t .ƒ/ D R.ƒ 0 / in V U D 0 on @ U (7)
where R(ƒ 0 ) is the residual vector computed at the initial value ƒ 0 (i.e. ƒ 0 is introduced in expressions [START_REF] Fortin | Localization of Hopf bifurcations in fluid flow problems[END_REF]) and L ƒ 0 t .ƒ/ is the tangent operator computed at the initial point ƒ 0 . To obtain the same number of equations and unknowns, two additional conditions are introduced, which give a well posed problem (see [START_REF] Jepson | Numerical Hopf bifurcation[END_REF][START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF][START_REF] Griewank | The calculation of Hopf points by a direct method[END_REF]). Finally, finding a Hopf bifurcation point consists in determining iteratively solution ƒ which verifies:

kR.ƒ/k < " (8)
where k k is the Euclidian norm of vector , and © is a chosen accuracy parameter. In the numerical tests studied, the value of this parameter is set and equal to 10 8 . As in reference [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], choosing as initial value ƒ 0 of the quantities coming from the computation of a bifurcation indicator is proposed.

Bifurcation indicator

In reference [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF], the authors have proposed to determine a Hopf bifurcation by introducing an indicator. This indicator was first proposed in references [START_REF] Bensaadi | Méthode Asymptotique-Numérique pour le calcul de bifurcations de Hopf et de solutions périodiques[END_REF][START_REF] Tri | Méthodes Asymptotiques-Numériques pour les fluides visqueux incompressibles et la détection de la bifurcation de Hopf[END_REF] in a solid mechanics framework. In this method, the steady solution U S is computed independently by means of an asymptotic numerical method ([8, 14, 15]). The determination of pulsation ¨and mode V is then solved by introducing an indicator ¥ in Equation ( 4), where the quantities (U S , oe S / are known. This indicator is a scalar function, which has the properties to be always positive and null at the Hopf bifurcation points [START_REF] Bensaadi | Méthode Asymptotique-Numérique pour le calcul de bifurcations de Hopf et de solutions périodiques[END_REF]. Considering that the steady solution is known and Equation (4) haS been rewritten with the indicator ¥, the following problem has to be solved:

² i¨M.V / C L.V / C Q.V , U s / C Q.U s , V / D ¥f in V U D 0 o n @ U (9)
where f is a known random vector. This problem is solved by using the finite element method. The following discrete problem has to be solved:

8 < : Ä K t ¨M ¨M K t ² V a V b ³ D ² ¥f 0 ³ in V U D 0 on @ U (10)
where matrix K t is the discrete tangent operator which is defined by

L t . / D L. / C Q. , U S / C Q.U S , /.
To have a well-formulated problem, a condition of normalization is imposed and has to be applied to both the real and imaginary parts of mode V. In past studies [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Potier-Ferry | Treatment of strong nonlinearities by the asymptotic numerical method[END_REF], it has been shown that some conditions lead to numeric instabilities, and that quadratic conditions are the most stable and the most adapted conditions to the perturbation method. According to these studies [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF][START_REF] Potier-Ferry | Treatment of strong nonlinearities by the asymptotic numerical method[END_REF], the following quadratic condition is applied:

kV k 2 D kV initial k 2 (11)
where V initial is the solution of Equation ( 9) for ¨D 0.0 and ¥ D 1.0. Finally, expressions [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] and [START_REF] Griewank | The calculation of Hopf points by a direct method[END_REF] permit one to compute the evolution of the indicator ¥ versus the angular frequency for a fixed value of the stationary solution U S . In other words, an indicator curve, ¥ Re .¨/, is computed for each Reynolds number. This indicator can be computed by varying the value of the angular pulsation ¨in Equation [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF]. Thus, each value of the indicator requires the triangulation of the matrix defined in Equation [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF]. Finally, with this method, the computational time needed for the calculus of this indicator can be large. To avoid this problem, the indicator curve is determined by using a perturbation method. Hence, only one matrix triangulation is carried out to compute an analytic part of the indicator curve.

First, the indicator is computed as a polynomial expansion of pulsation ¨D ¨0 C O ¨, where ! 0 is the initial angular frequency and O ¨is the perturbated parameter: 8 < :

¥. O ¨/D ¥ 0 C O ¨¥1 C O ¨2¥ 2 C : : : C O ¨p¥ p V a D V a 0 C O ¨V a 1 C O ¨2V a 2 C : : : C O ¨pV a p V b D V b 0 C O ¨V b 1 C O ¨2V b 2 C : : : C O ¨pV b p (12)
In Equations ( 12), V a 0 , V b 0 and ¥ 0 are known and regular solutions of problem [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF] and p is the order of truncature of the polynomial expansions.

Inserting expansion [START_REF] Bensaadi | Méthode Asymptotique-Numérique pour le calcul de bifurcations de Hopf et de solutions périodiques[END_REF] in Equation [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF], and balancing terms with identical powers of O ¨, the following set of linear equations is obtained: Order 0:

K t .V a 0 / ¨0M.V b 0 / D ¥ 0 f in K t .V b 0 / C ¨0M.V a 0 / D 0 in V a 0 D V b 0 D 0 on @ u kV 0 k 2 D kV initial k 2 (13) Order p: K t .V a p / ¨0M.V b p / D ¥ p f C M.V b p 1 / in K t .V b p / C ¨0M.V a p / D M.V a p 1 / in V a p D V b p D 0 on @ u 2 < V p , V 0 > C p 1 P rD1 < V p r , V r > D 0 (14)
Finally, the previous linear system of equations helps compute all the unknowns at each order p with a single matrix triangulation and p backward and forward substitutions.

Once all the terms of the polynomial expansion are determined, the latter is replaced by equivalent rational approximations, called Padé approximants first introduced in [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF]:

V Padé . O ¨/ V 0 D p 1 X kD1 P .p 1 k/ . O ¨/ Q .p 1/ . O ¨/ O ¨kV k (15)
where P k and Q k are polynoms of degree k. These rational approximations increase the range of validity of the polynomial approximations [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF] and consequently help decrease the number of matrix triangulations required to compute the curve ¥ Re .¨/. The validity range of the Padé approximants is determined by the following expression [START_REF] Elhage-Hussein | A numerical continuation method based on Padé approximants[END_REF]:

ı D V Padé,p . O ¨max / V Padé,p 1 . O ¨max / V Padé,p . O ¨max / V 0 ( 16
)
where ı is a small chosen parameter (10 3 Ä ı Ä10 6 /. With the latter expression, one can define the maximum value of the perturbation parameter, O ¨max , and then compute a new starting point of the perturbation method by introducing O ¨max in the asymptotic expansions [START_REF] Bensaadi | Méthode Asymptotique-Numérique pour le calcul de bifurcations de Hopf et de solutions périodiques[END_REF]. Finally, the indicator curve, ¥ Re .¨/, is entirely computed by applying, in a step by step manner, the previous continuation method.

Additional remarks

The present method can be viewed as an eigenvalue problem. The eigenvalues, corresponding to Hopf's bifurcations, are obtained when the indicator is null.

For that kind of problem, iterative methods such as the inverse power method [START_REF] Golub | Matrix Computations, 3rd ed[END_REF] can be used. In solid mechanics [START_REF] Cadou | A solver combining reduced basis and convergence acceleration with applications to nonlinear elasticity[END_REF], it has been shown that the results obtained by an iterative method and those obtained by a perturbation method are the same. Nevertheless, it has also been shown in the latter study, that the acceleration convergence method has to be used to obtain acceptable efficiency. In the studies [START_REF] Cadou | A solver combining reduced basis and convergence acceleration with applications to nonlinear elasticity[END_REF][START_REF] Damil | Mathematical and numerical connections between polynomial extrapolation and Padé approximants: applications in structural mechanics[END_REF], the convergence acceleration methods are based on the use of Padé approximants or the extrapolation method.

PROPOSED ALGORITHM

Before introducing in details the method used in this study to determine the Hopf bifurcation points, let us first give some specifications.

In reference [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF], the indicator, as previously defined, was used to determine the instability points. Consequently, many computations of the indicator curves were necessary to determine the precise values of the critical Reynolds numbers and angular frequency. In [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], this indicator was associated with the Newton method. This method is denoted 'mixed method'. The indicator is first computed and gives initial values to the iterative method. These initial values are the minima of the indicator in the computed range of angular frequency. Nevertheless, the Newton method is efficient if the initial guess is not far from the solution. Thus, in some cases, the initial value does not lead to the convergence of the direct method, for example, when the Reynolds number, for which the indicator is computed, is far from the critical value. Consequently, much computation is spent in indicator calculus without giving an interesting value for the Newton method. In the results shown in [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], the indicator curves present many extrema and consequently many possible initial values for the Newton method.

Thus, in this study, we propose to determine the indicator curve ¥ Re D ¥ Re .¨/ for a fixed Reynolds number Re and its first and second derivatives with respect to pulsation d¥ Re /d¨and d 2 ¥ Re /d¨2. One has to firstly compute the roots of the first derivative and secondly check that these roots are minima of the indicator curves by introducing them in the second derivative. Therefore, these minima are determined and a set of possible candidates are easily and automatically obtained and used as initial values for the Newton procedure.

Expression (8) allows for obtaining an analytical expression of the indicator function with respect to pulsation ¨. Then, the two derivative functions d¥=d O ¨and d 2 ¥=d O ¨2 with respect to pulsation O can be derived in an analytical way:

d¥ d O ¨. O ¨/ D ¥ 1 C 2¥ 2 O ¨C 3¥ 3 O ¨2 C : : : C p¥ p O ¨p 1 (17) d 2 ¥ d O ¨2 . O ¨/ D 2¥ 2 C 6¥ 3 O ¨C : : : C p.p 1/¥ p O ¨p 2 (18)
To increase the convergence, these derivative functions are also computed with the help of Padé approximants by using the two previous asymptotic expansions [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF]:

d¥ d O ¨ˇP ad K e . O ¨/ ¥ 1 D p 2 X kD1 N P .p 2 k/ . O ¨/ N Q .p 2/ . O ¨/ O ¨k¥ k (19) d 2 ¥ d O ¨2 ˇPad K e . O ¨/ 2¥ 2 D p 3 X kD1 Q P .p 3 k/ . O ¨/ Q Q .p 3/ . O ¨/ O ¨k¥ k (20)
In expressions [START_REF] Cadou | A solver combining reduced basis and convergence acceleration with applications to nonlinear elasticity[END_REF] and [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolation and Padé approximants: applications in structural mechanics[END_REF],

N P k . O ¨/, N Q. O ¨/, Q P k . O ¨/, Q Q. O
¨/ appear as polynomial expressions and are combined into a rational one. All the fractions are noted to have the same respective denominator.

The roots and the poles of Padé expansions are evaluated by a classical method [START_REF] Press | Numerical Recipes in Fortran 77: The Art of Scientific Computing[END_REF]. Among all the possible values, the minima of function ¥.¨/ are extracted. The latter are the initial guesses for the Newton algorithm. The main advantage offered by the present method is to determine automatically a set of possible candidates instead of calculating a unique one which does not ensure the convergence of the Newton calculation as in reference [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF].

NUMERICAL RESULTS

To highlight the performances of the proposed algorithm, we study the Hopf bifurcation points of the lid-driven cavity problem for three geometric configurations, characterized by the aspect ratio values A D 0.8, A D 1.0 and A D 1.5 (Figure 1). The description of the meshes used in the numerical tests, for all the aspect ratio A, is given in Table I. For these configurations, the results issued from the literature are exposed in Table II. A particular emphasis in this section is given for the case A D 1.0. [START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF] 8026.6 0.4497 2268.0 Fortin [START_REF] Fortin | Localization of Hopf bifurcations in fluid flow problems[END_REF] 8000.0 0.4500 2261.9 Abouhamza [START_REF] Abouhamza | A neutral stability curve for incompressible flows in a rectangular driven cavity[END_REF] 8004.5 --Tiesinga [START_REF] Tiesinga | Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method[END_REF] 8375.0 0.4399 2314.8 Poliashenko [START_REF] Poliashenko | A direct method for computation of simple bifurcations[END_REF] 7763.4 0.4500 2195.0 Bruneau [START_REF] Bruneau | The 2D lid-driven cavity problem revisited[END_REF] 8000-8050 --Peng [START_REF] Peng | Transition in 2-D lid-driven cavity flow[END_REF] 7402.0 0.5900 2744.0 Gervais [START_REF] Gervais | Some experiments with stability analysis of discrete incompressible flows in the lid-driven cavity[END_REF] 7960.0 0.4500 2250.6 8040.0 0.4500 2273.3 Cazemier [START_REF] Cazemier | Proper orthogonal decomposition and low-dimensional models for driven cavity flows[END_REF] 7819.0 0.6100 2996.8 7972.0 0.4500 2254.0 Auteri [START_REF] Auteri | Numerical investigation on the stability of singular driven cavity flow[END_REF] 8018.0 0.4496 2265.0 A=1.5 This study 4722.0 0.3280 973.0 Boppana [START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF] 5326.0 0.3400 1137.8 Abouhamza [START_REF] Abouhamza | A neutral stability curve for incompressible flows in a rectangular driven cavity[END_REF] 5674.0 0.4900 1746.9 Poliashenko [START_REF] Poliashenko | A direct method for computation of simple bifurcations[END_REF] 7216.0 0.4400 1994.9

x = 0.0 x = 1.0 y = 1.0 y = 0.0 (u = v = 0.0) (u = v = 0.0) (u = v = 0.0) (u = 1.0, v = 0.0) L = 1.0 A = 1.5 x = 0.0 x = 1.0 y = 1.0 y = 0.0 (u = v = 0.0) (u = v = 0.0) (u = v = 0.0) (u = 1.0, v = 0.0) L = 1.0 A = 1.5
A=0.8 This study 5297.0 0.3460 1152.0 Poliashenko [START_REF] Poliashenko | A direct method for computation of simple bifurcations[END_REF] 5225.0 0.3500 1149.0

Results with the proposed algorithm

In the case of a lid-driven cavity parameter A D 1.0, three initial values of Reynolds number Re i are chosen: Re i D 6332, 7257, 8089 (see Table III). For this example, the indicator curve ¥.¨/ and the derivative functions d¥/d¨and d 2 ¥/d¨2 are plotted in Figures 234for a Reynolds number equal to 6332. In the latter, the curves are exposed in a non-dimensional form using the Strouhal number St D ¨L=.2 u/. These curves are computed by using a truncature order of the asymptotic expansions equal to 20 and a parameter ı equal to 10 6 . The indicator is computed by using the Padé-based continuation method up to an angular frequency equal to 5000. Let us remark that these parameters are the same for all the numerical tests presented in this work. For each of these three chosen values, the present algorithm determines several initial values of pulsation ¨i and vector V i : 13, 14 and 9 initial couples (¨i V i ) are obtained for ReD6332, 7257 and 8089, respectively. For each Reynolds number, at least one critical Reynolds number is obtained, that is to say at least one Hopf bifurcation is determined (see results in Table III). It is considered in this study that the Newton algorithm converges if the norm of the residual vector is lower than 10 4 and 10 8 , respectively, for the stationary and dynamic solution. The maximum number of iterations is set as equal to 20. When the Newton method converges to a critical value, the required number of iterations is given in Table III. Symbol '-' indicates that the iterative method does not find any critical values after 20 iterations.

However, it is noticed that the convergence to a critical number is not necessarily obtained for the minimum value of indicator ¥ i . This is a noticeable fact and constitutes a major feature compared with the previous work [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], where all the calculations are based on the values obtained for the minimum of the indicator ¥ i . In other words, according to the previous algorithm [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF] and by considering the case ReD6332, all the calculations with the retained values ¨i D1085 and ¥ i D0.653 were unnecessary, as convergence was not obtained for this set of values. This is also the case when the considered initial Reynolds number is equal to 7257, and this can be noted as one of the drawbacks of the previously proposed algorithm [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF].

The minimum value of the indicator is then not a sufficient condition to ensure the algorithm convergence. Thus, the present algorithm offers the real improvement of having a set of candidates to calculate Hopf bifurcation points. For example, six, five and five different critical points are found with an initial Reynolds number respectively equal to 6332, 7257 and 8089.

To verify that the bifurcation indicator gives good initial values, we define the following correlation coefficient: where vector V i is the initial velocity (computed from the minimum value of the indicator) and V c is the critical Hopf mode found. For example, if one considers the results obtained in the case Re i D 8089 (Table III), the proposed algorithm gives six initial values for the Newton method. Among these values, from the couple (¨i D 2236, ¥ i D 0.395) four Newton iterations are needed to obtain the first Hopf bifurcation point (i.e. Re c D 7890, ¨c D 2185). Thus, if this initial vector is introduced in the correlation coefficient as V i and the computed first bifurcation vector as V c , this correlation coefficient is then equal to 0.97 (see Table IV), which is very close to 1. In Table IV, for all the computed Hopf points, this correlation coefficient is given and is never lower than 0.6 even for critical Reynolds numbers far from the initial value. For example, for the critical Reynolds number equal to Re c D11418 and ¨c D5637, this correlation coefficient is equal to 0.67. The results obtained from the proposed algorithm are very close to the critical values issued from the literature for the first bifurcation point (Re c D 7890, ¨c D 2185). The latter are summarized in Table II. Indeed, for a geometric aspect ratio AD1.0, the critical Reynolds number is between 7400, and 8050 and the Strouhal number is close to 0.44.

C.V i , V c / D < V i , V c > kV i k kV c k (21)
The values exposed in Table III are directly issued from the calculations. Surprisingly, it can be observed that the critical points are not found in function with the increasing Reynolds number. The algorithm seems to compute the solutions in function with the increasing Strouhal numbers.

The present method is designed to offer an automatic algorithm to compute Hopf bifurcation points. The results exposed in Table V show that the algorithm is able to find, automatically, 10 critical values declared as Hopf bifurcation without any user intervention. This is an advantageous characteristic for parametric study, for which many calculations are required. It can be noticed that the fundamental Hopf bifurcation is also obtained. These 10 critical values are obtained by computing the indicator curves for some Reynolds numbers in the range [2000-16 000]. The streamline real and imaginary vectors of these Hopf bifurcation points are plotted respectively in Figures 5 and6. For these plots, we only consider the first five critical values. For the first Hopf bifurcation, Re c D7890 and St c D0.44, one can remark that these streamlines are quite similar to those shown in references [START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF][START_REF] Fortin | Localization of Hopf bifurcations in fluid flow problems[END_REF][START_REF] Abouhamza | A neutral stability curve for incompressible flows in a rectangular driven cavity[END_REF].

Thus in this range, the performance of the presented algorithm is evaluated through Figure 7 and8. Figure 7 presents the evolution of the number of roots and bifurcation points in function with the Reynolds number in the case of the lid-driven cavity with AD1.0. It can be seen that the number of possible candidates is comprised between 0 and 16 for the Reynolds numbers tested. Moreover, Figure 8 exposes the rank of the critical points in function with the initial Reynolds number Re i . This rank is given in Table V. For example, it can be observed that the fundamental Hopf Bifurcation (point n ı 1) is computed from 13 previously isolated initial Reynolds numbers Re i . It can be noted that the convergence range is significant: the fundamental Hopf bifurcation point is calculated for initial Reynolds number Re i comprised between [ 3500-9500]. For information, it can be noticed that a 10 th critical point can be computed for high values of initial Reynolds numbers. Finally, it is observed that the Hopf bifurcation n˚5, corresponding to a minimum Strouhal number (St c D 0.33) is found 22 times in a large range of Reynolds number. Hence, the proposed algorithm gives many possible Hopf bifurcation points. In the literature, few references give critical parameters about the second or third instability point. Nevertheless, in reference [START_REF] Tiesinga | Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method[END_REF], the authors, by computing the eigenvalues of the Jacobian matrix, give an approximation of the first critical points. These authors have found that the Hopf bifurcation points appear for the following Reynolds number Re c1 2[8000-8375], Re c2 2[8375-8875], Re c3 2[8875-9375], Re c4 2[9875-10 375] and respectively for the following Strouhal numbers: 0.44, 0.6, 0.52, 0.7 and 0.77. A quick comparison between the latter results and ours shows that the Strouhal numbers are the same. Indeed, results in Table III show that we found exactly these five Strouhal numbers among the first six bifurcation points. Nevertheless, these numbers are not found for the same critical Reynolds numbers in both studies. For example, in reference [START_REF] Tiesinga | Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method[END_REF], they have found the first bifurcation for a Reynolds number between [8000-8375] (Re c1 D 7890 in this study). Moreover, in this reference, the second Hopf bifurcation point is obtained for a Reynolds number included in [8375-8875] and for a Strouhal number equal to 0.6. It seems that this point does not correspond to the second, but to the third bifurcation point we found (see Table V). The second and third Hopf bifurcation points are then reversed compared with those computed in reference [START_REF] Tiesinga | Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method[END_REF]. This dissimilarity can be explained by the fact that both studies use different numerical schemes. For example, the space discretization is performed with a finite volume method in reference [START_REF] Tiesinga | Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method[END_REF], whereas in this study, a finite element method is used.

Comparison with other numerical methods

In this section, we compare the critical values obtained with four numerical methods. The first one is the initial bifurcation indicator presented in [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF], the second one is the hybrid method [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], the third one is the proposed method and the last one consists in computing the eigenvalues of the Jacobian matrix for a chosen Reynolds number and in checking if a pair of eigenvalues crosses the imaginary axis.

The computation of the eigenvalues of the Jacobian matrix allows the bifurcation points to be determined. Hence, a Hopf bifurcation is characterized by two eigenvalues, which are purely imaginary. However, the Jacobian matrix is non symmetric, and specific algorithms have to be used. In this paper, we propose comparing the critical Reynolds numbers obtained with the proposed algorithm and those calculated with the ARPACK [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF] package based on the Arnoldi method [START_REF] Lehoucq | Arpack User's Guide: Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]. The results obtained from the eigenvalue calculations for different Reynolds numbers are exposed in small. This indicates that for the numerical scheme used in this study, the critical parameters found with the proposed method seem to be correct. We now compare the computational times required for the four considered methods to determine Hopf bifurcation. It is quite difficult to seriously compare these four methods. Indeed, two methods give accurate values of the Hopf bifurcation points (the Hybrid and the proposed method), and the two others only give approximate values of the critical parameters (the bifurcation indicator and ARPACK [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF]). Nevertheless, we give hereafter some elements, which permit us to exhibit the advantages and drawbacks of each method. Thus, in Table VII, the basic CPU times are given for three methods to compute Hopf bifurcation points in the case of the lid-driven problem with an aspect ratio A D 1.0 and a number of D.o.F. equal to 13222.

The first method is the indicator presented in [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF]. With this method, 66 steps are necessary to obtain the indicator curve in the range ¨2 OE0, 5000. Each step requires a computational time of nearly 68 s (according to the results presented in reference [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF]). Finally, for each initial Reynolds number, the indicator curve is determined with nearly 4500 s. In Table VII, the computational time for a single Newton iteration is also given. This time is equal to 60 s. Finally, the computational time required with ARPACK [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF] to give a part of the spectrum is also indicated in this table. This computational time is approximatively 1440 s. Indicator [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF] for ¨2[0, 5000] 66 steps 68 s = 4500 s Newton's iteration 60 s Arnoldi computation [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF] 1440 s

With these basic times, we are able to give an estimation of the CPU times needed to compute Hopf bifurcation points with the four considered methods.

Therefore, in Table VIII, these estimated CPU times are given to obtain the Hopf bifurcation point with the Hybrid method [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], the bifurcation indicator [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF], an eigenvalue computation with ARPACK [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF] and the proposed method. For example, the bifurcation indicator [START_REF] Cadou | A numerical method for the computation of bifurcation points in fluid mechanics[END_REF] requires at least three computations to check if a Hopf bifurcation occurs. By taking into account the basic CPU times given in Table VII, this leads to an approximated computational time equal to 13 500 s. Nevertheless, these three indicator curves are not sufficient to give an accurate value of the critical point. For this, additional computations are necessary and consequently the CPU times will increase. That is also the case with ARPACK [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF], where more than one eigenvalue computation is required to give accurate critical parameters. Nevertheless, a single one helps certify if a Hopf bifurcation occurs by checking if a pair of conjugated eigenvalues crosses the imaginary axis. For these reasons, we have indicated in Table VIII that with the two previous methods, only estimated critical points can be found. Contrarily, with the two other methods, the Hybrid [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF] and the proposed method, accurate critical parameters can be determined. Hence, the Hybrid method requires approximately 4740 s to determine a Hopf bifurcation point with an initial Reynolds number equal to 8089. This computational time has been evaluated by taking into account an indicator calculus (4500 s) and four iterations of Newton's method (according to the results presented in reference [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF]). The proposed method requires a computational time equal to 6660 s which is determined by also considering an indicator calculus (4500 s) and 36 iterations of Newton's method. These 36 iterations are issued from the results presented in Table III by only taking into account the initial guesses leading to a bifurcation point. This means that, at the end of the computations, six Hopf bifurcation points are accurately determined with the proposed method instead of a single one with the Hybrid method [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF]. Finally, one can conclude that with the proposed method and for the considered numerical examples, one can determine six bifurcation points with a computational time, which is approximately equal to five eigenvalue calculi with ARPACK [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF]. Let us remark that to compute the same number of critical values with ARPACK [START_REF] Nayar | Computation of selected eigenvalues of generalized eigenvalue problems[END_REF], more than five eigenvalue computations would certainly be necessary.

The computational times of the hybrid and the proposed methods show that it is the indicator phase, which is the step that consumes the most amount of time. This topic will be dealt with in further works, which means a decrease of the CPU times required to obtain the initial guesses of Newton's method with the indicator.

Supplementary aspect ratios

One considers now the lid-driven cavity with two other aspect ratios: AD0.8 and AD1.5. For these two examples, we give, in Table II, the first Hopf bifurcation point determined by the proposed method and the ones reported in the literature. For the case AD0.8, with the proposed method, the first bifurcation point appears for a Reynolds number and a Strouhal number respectively equal to 5297 and 0.346. These numbers are relatively close to the ones reported by Poliashenko and Aidun [START_REF] Poliashenko | A direct method for computation of simple bifurcations[END_REF], which are respectively 5225 and 0.35. Let us remark that these authors also use a Newton method to detect instabilities, but by taking as initial guesses the quantities resulting from a time-dependent Navier-Stokes computation. Concerning the geometric ratio A D 1.5, the results determined with the proposed method are close to the results reported by Boppana and Gajjar [START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF] essentially for the Strouhal number (0.328 for the proposed method and 0.34 in reference [START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF]). The critical Reynolds numbers are quite different (4722 for this study and 5326 for reference [START_REF] Boppana | Global flow instability in a lid-driven cavity[END_REF]). Nevertheless, comparison with the results given in references [START_REF] Abouhamza | A neutral stability curve for incompressible flows in a rectangular driven cavity[END_REF] and [START_REF] Poliashenko | A direct method for computation of simple bifurcations[END_REF] shows that for this example, the critical numbers for which the first Hopf bifurcation appears is still an open question.

As with the case A D 1.0, the proposed method can give more than one bifurcation point. Thus, in Tables IX and X, we give the first seven critical values respectively for A D 0.8 and A D 1.5. These values are also obtained by computing a single indicator curve and by determining along this curve the minima of the indicator. These values are then introduced as previously performed in the Newton scheme.

CONCLUSION

In this article, an extension of a previous study is presented [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF]. The exposed algorithm allows for the automatic determination of bifurcation points by computing an indicator curve and its first two derivatives with respect to the pulsation of the flow ¨. With this method, several frequencies are computed as roots of the indicator curve and are used as initial values for an iterative Newton algorithm. The issued solution corresponds to a bifurcation point. The results exposed in this paper show the efficiency and the accuracy of the studied algorithm. Compared with the previous work, the range for the stationary solution for which Hopf bifurcation points are found is considerably increased. Thus, for the lid-driven cavity with an aspect ratio equal to 1, the first bifurcations are found from an initial Reynolds number close to 3000. Whereas, with the hybrid method proposed in reference [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of Hopf bifurcation points in fluid mechanics[END_REF], this number is nearly 7600. The advantage of the proposed method compared, for example, with an eigenvalue computation is that the bifurcation point is automatically computed without any user intervention. In this study, initial values are obtained using a perturbation method and transmitted to a Newton's algorithm. It could be interesting to develop a new algorithm coupling a direct method (Newton's algorithm) with initial values, which would be isolated by Arnoldi's procedure. Then, it could be possible to compare the CPU times obtained from this method with those obtained from the author's procedure.

Moreover, the calculation times are still significant, and improvements have to be carried out in that way. As the computations are carried out in two steps, the indicator phase and the Newton method, future works will concern these two points. Thus, the computation of the indicator can be realized by using, for example, a linear solver well-adapted to repeated right-hand side problems [START_REF] Cadou | A solver combining reduced basis and convergence acceleration with applications to nonlinear elasticity[END_REF] or with a reduced order model (POD [START_REF] Gervais | Some experiments with stability analysis of discrete incompressible flows in the lid-driven cavity[END_REF], APR [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Verdon | Reduced-order modelling for solving linear and non-linear equations[END_REF]). For the second point, some high order iterative algorithms [START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF][START_REF] Cadou | Projection techniques to improve high-order iterative correctors[END_REF] can be used instead of the Newton method. These methods, based on homotopy and Padé approximants, have the advantage of converging to the solutions with lower CPU times than a classical iterative method such as Newton's method.

Figure 1 .

 1 Figure 1. Geometry for the two-dimensional lid-driven cavity flow (case A D 1.5). In all the calculations of the Reynolds and Strouhal numbers, the reference length is chosen as L D 1.0.

Figure 2 .Figure 3 .

 23 Figure 2. Evolution of the indicator function with respect to the Strouhal number St-Problem of the liddriven cavity-Case A D 1.0 for Re D 6332. The truncature order of asymptotic equations and the parameter ı are, respectively, equal to 20 and 10 6 .

2 φ/dSt 2 Figure 4 .

 224 Figure 4. Evolution of the second derivative function of the indicator with respect to the Strouhal Number St-Problem of the lid-driven cavity-Case A D 1.0 for Re D 6332. The truncature order of asymptotic equations and the parameter ı are, respectively, equal to 20 and 10 6 .

33 Figure 5 .

 335 Figure 5. Streamlines of the real part of the eigenvectors for several critical Reynolds numbers. The lid-driven cavity A D 1.0.

Figure ( 9 ) 33 Figure 6 .

 9336 Figure[START_REF] Baker | Padé approximants: Volume 59 of encyclopaedia of mathematics and its applications[END_REF], for the lid-driven cavity with an aspect ratio equal to 1 (Mesh 2, 13 122 D.o.F.). The evolution of the eigenvalues versus the Reynolds number is quite similar to Figure4of reference[START_REF] Tiesinga | Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method[END_REF]. A comparison of the results obtained with the eigenvalue computations with the critical parameters found with the proposed algorithm is presented in Table VI. The relative percentage error is quite

Figure 7 .

 7 Figure 7. Number of computed roots and determined bifurcation points for different Reynolds numbers-Example of the lid-driven cavity-Case A D 1.0 -Mesh 1.

Figure 8 .

 8 Figure 8. Results of the Hopf bifurcation number for different initial Reynolds numbers (related to the Table V)-Example of the lid-driven cavity-Case A D 1.0 -Mesh 1.

Figure 9 .

 9 Figure 9. Representation in a complex plane of the eigenvalues of the Jacobian matrix for 19 different values of the Reynolds number. Symbol ¢ stands for the complex eigenvalue of the tangent operator K T . Estimation of five critical pulsations corresponding to Re.¢/ D 0: ¨c1 =2183, ¨c2 D 2387, ¨c3 D 2926, ¨c4 D 3424, ¨c5 D 4226-Example of the lid-driven cavity (A D 1.0).

Table I .

 I Characteristics of the smallest finite element used for the meshes -(R) and (I) stand respectively, for regular and irregular meshes.

			A=0.8		A=1.0		A=1.5
	Mesh name	D.o.F.	x	y	x	y	x	y
	Mesh 1	13122 (R)	0.0250	0.020	0.0250	0.0250	0.0250	0.03750
	Mesh 2	13122 (I)	0.0025	0.002	0.0025	0.0025	0.0025	0.00375

Table II .

 II Critical values for the first Hopf bifurcation in the case of the lid-driven cavity for three different geometries defined by A=1.0, A=1.5 and A=0.8.

	Aspect		Critical Reynolds number,	Critical Strouhal number,	Critical pulsation,
	ratio, A	Authors	Re c	St c	¨c
	A=1.0	This study	7890.0	0.440	2185.0
		Boppana

Table III .

 III Numerical results obtained for the problem of the lid-driven cavity-case A=1.0-Critical Reynolds and Strouhal numbers obtained for initial Reynolds and Strouhal numbers evaluated as minima of the indicator curve. Symbol '-' stands for the non-convergence of the calculus.

		Initial values			Critical values
	Reynolds	Pulsation,	Indicator,	Reynolds number,	Pulsation,	Strouhal number,	Number of
	number, Re i	¨i	¥ i	Re c	¨c	St c	iterations
	6332	591	0.745	-	-	-	-
		757	0.744	-	-	-	-
		1085	0.653	-	-	-	-
		1361	1.020	11286	2385	0.33	9
		1799	2.670	7890	2185	0.44	6
		2066	2.860	8829	2925	0.52	8
		2159	3.660	8829	2925	0.52	8
		2323	5.700	8921	3424	0.61	10
		3239	6.680	11417	5636	0.78	8
		3450	6.940	-	-	-	-
		3552	6.830	12625	6877	0.87	16
		3961	6.560	-	-	-	-
		4515	9.030	-	-	-	-
	7257	677	0.654	-	-	-	-
		925	0.667	-	-	-	-
		1231	0.508	-	-	-	-
		1582	0.699	11286	2385	0.33	6
		2021	1.300	7890	2185	0.44	5
		2418	1.610	8829	2925	0.52	7
		2701	2.750	-	-	-	-
		2776	2.660	8921	3424	0.61	6
		3216	4.850	-	-	-	-
		3701	5.740	-	-	-	-
		3759	4.580	-	-	-	-
		4075	6.110	-	-	-	-
		4159	5.690	12625	6877	0.87	10
		4502	6.220	-	-	-	-
	8089	754	0.640	-	-	-	-
		1085	0.606	-	-	-	-
		1360	0.418	-	-	-	-
		1759	0.486	11286	2385	0.33	6
		2236	0.395	7890	2185	0.44	4
		2691	0.669	8829	2925	0.52	5
		3105	1.230	8921	3424	0.61	5
		3614	3.120	9513	4226	0.70	7
		4100	4.630	11417	5636	0.78	9

Table IV .

 IV Example of the lid-driven cavity, value of the correlation number C . The initial Reynolds number is equal to 8089.

	Initial pulsation, ¨i	Critical pulsation, ¨c	Critical Reynolds number, Re c	Correlation factor, C.V i , V c /
	2236	2186	7890	0.97
	3105	3424	8921	0.91
	3614	4226	9513	0.87
	1759	2386	11287	0.86
	2961	2925	8829	0.85
	4100	5637	11418	0.67

no matter the initial Reynolds number included in the range [3000-14 000], at least one bifurcation point is always computed. For example, for a Reynolds number equal to Re 6200, the algorithm retains 16 initial values Re i and computes six critical Reynolds numbers Re c .

Table V .

 V Example of the lid-driven cavity-case A=1.0-Presentation of the 10 computed critical Reynolds and Strouhal numbers.

	Hopf bifurcation number	Critical Reynolds number, Re c	Critical Strouhal number, St c
	1	7890	0.44
	2	8829	0.52
	3	8921	0.61
	4	9513	0.71
	5	11286	0.33
	6	11417	0.78
	7	12625	0.86
	8	14674	0.92
	9	15501	0.99
	10	16170	0.24

Table VI .

 VI Comparison of the critical pulsations obtained by the present algorithm with those obtained by an eigenvalue calculation (Arnoldi method)-Example of the lid-driven cavity-Case A=1.0 -Mesh 1.

	Critical pulsation, ¨AL c obtained for the present algorithm	Critical pulsation, ¨AR c obtained for the Arnoldi method	Relative error calculated as 100 j¨A L c ¨AR c j / ¨AL c
	2181	2183	0.08
	2340	2387	2.00
	2885	2926	1.43
	3419	3424	0.14
	Table VII. Basic central processing unit times for three methods; The
	considered example is the lid-driven cavity with A=1.0 and 13122 D.o.F.
	(Mesh 2). The initial Reynolds number is equal to 8089.
	Method	CPU times (s)

Table VIII .

 VIII Estimated computational times for four numerical methods to detect a Hopf bifurcation point. For the first three methods, the considered central processing unit times are issued from Table VII. The considered example is the lid-driven cavity with A=1.0 and 13 122 D.o.F. (Mesh 2). The initial Reynolds number is Re=8089.

	Method	CPU time, (s)	Number of critical values found	Critical point
	Hybrid method [1]	4500 + 4 60 = 4740	1	accurate
	Bifurcation indicator [2]	4500 3 = 13500	1	estimated
	ARPACK [24]	1440	1	estimated
	Proposed method	4500 + 36 60 = 6660	6	accurate

Table IX .

 IX Presentation of the critical parameters obtained for the lid-driven cavity-Case A=0.8-Mesh 2.

	Critical Reynolds number, Re c	Critical pulsation, ¨c	Critical Strouhal number, St c
	5297	1152	0.346
	5790	1704	0.468
	5966	2208	0.589
	6951	3143	0.720
	8095	4250	0.836
	8934	5358	0.954
	10541	7101	1.072

Table X .

 X Presentation of the critical parameters obtained for the lid-driven cavity-Case A=1.5-Mesh 2.

	Critical Reynolds number Re c	Critical pulsation, ¨c	Critical Strouhal number, St c
	4722	973	0.328
	6465	1715	0.422
	7032	2150	0.486
	7369	2656	0.573
	7903	3267	0.657
	8331	3928	0.750
	9783	5142	0.836


J M Cadou

M Potier-Ferry

A solver combining reduced basis and convergence acceleration with applications to non-linear elasticity

An iterative solver is proposed to solve the family of linear equations arising from the numerical computation of non-linear problems. This solver relies on two quantities coming from previous steps of the computations: the preconditioning matrix is a matrix that has been factorized at an earlier step and previously computed vectors yield a reduced basis. The principle is to define an increment in two sub-steps. In the first sub-step, only the projection of the unknown on a reduced subspace is incremented and the projection of the equation on the reduced subspace is satisfied exactly. In the second sub-step, the full equation is solved approximately with the help of the preconditioner. Last, the convergence of the sequences is accelerated by a well-known method, the modified minimal polynomial extrapolation. This algorithm assessed by classical benchmarks coming from shell buckling analysis. Finally, its insertion in path following techniques is discussed. This leads to non-linear solvers with few matrix factorizations and few iterations.

INTRODUCTION

The algorithms to solve non-linear problems transform a non-linear system into to a family of linear ones. In general, these linear systems cannot be solved simultaneously because the right-hand side vectors depend on the solutions at previous steps. At least two classes of non-linear algorithms can be distinguished. In the first class, the matrix and the right-hand side are updated at each iteration, as for instance in the Newton-Raphson method. In the second group that includes the modified Newton method (MNM) and the asymptotic numerical method (ANM, [START_REF] Cochelin | Asymptotic-numerical method and Padé approximants for non-linear elastic structures[END_REF]), the same matrix holds good for several linear equations. Iterative solvers are convenient for the first class, but direct solvers are more efficient when a matrix is associated with multiple right-hand sides, because the important computation time due to matrix factorization is shared between several linear systems.

Although an iterative method is not a natural choice for linear problems with multiple right-hand sides, there are many papers dealing with this question, see, for example, Reference [START_REF] Gutknecht | See the web site[END_REF]. When the vectors are known since the beginning, all the systems can be solved simultaneously with the so-called block solvers BL-*. The block conjugate gradient method can then be applied [START_REF] Leary | The block conjugate gradient algorithm and related methods[END_REF] if the matrix is symmetric, as it is generally the case in structural mechanics. For non-symmetric matrices, a consequent bibliography can be found in [START_REF] Toutounian | Global least squares method (GL-LSQR) for solving general linear systems with several right-hand sides[END_REF][START_REF] Bouyouli | On block minimal residual methods[END_REF].

Unfortunately in the case of a non-linear problem, the right-hand side vectors are build up recursively and block solvers are not adapted. Another idea is to account for the results of previous linear systems in the numerical treatment of the current system, for instance by using the previously computed vectors. Such techniques have been introduced under various names: deflated conjugate gradient [START_REF] Saad | A deflated version of the conjugate gradient algorithm[END_REF], augmented conjugate gradient [START_REF] Risler | On the use of Ritz vectors for the solution to nonlinear elasticity problems by decomposition methods[END_REF] or subspace projection extrapolation [START_REF] Clemens | Subspace projection extrapolation scheme for transient field simulations[END_REF]. Comparisons between these methods have been done within a computational electromagnetic framework [START_REF] Clemens | Multiple right hand side techniques for numerical simulation of quasistatic electric and magnetic fields[END_REF]. However, these methods have several drawbacks when they are applied with a cheap preconditioning technique such as ILU or Incomplete Cholesky. Indeed a too large number of stored vectors leads to numerical instabilities and many re-used vectors are not necessarily representative of the problem to be solved.

The aim of this paper is to present and discuss two linear iterative solvers that are convenient in the numerical solutions of non-linear problems. These solvers rely on previously computed quantities. Among linear problems, a few are solved by a direct method and most of them are solved iteratively, the last factorized matrix being the preconditioner. The basic idea is more or less the same as in the MNM, but in the present algorithm, a preconditioning matrix is used to solve linear equations and not non-linear ones as in MNM. Further this procedure is coupled with a reduced basis technique, the basis being made of previously computed vectors, as in [START_REF] Saad | A deflated version of the conjugate gradient algorithm[END_REF][START_REF] Risler | On the use of Ritz vectors for the solution to nonlinear elasticity problems by decomposition methods[END_REF][START_REF] Clemens | Subspace projection extrapolation scheme for transient field simulations[END_REF][START_REF] Clemens | Multiple right hand side techniques for numerical simulation of quasistatic electric and magnetic fields[END_REF].

A non-linear solver has been presented recently [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF]. It is based on a homotopy transformation, a perturbation technique and convergence acceleration by Padé approximants (HPP). It improves the high-order Newton algorithm and the high-order iterative algorithms introduced in [START_REF] Damil | An iterative method based upon Padé approximants[END_REF][START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF]. Roughly, one applies the consistent tangent matrix on the reduced subspace and the chosen preconditioner in a correction phase on the whole space. The linear version of this HPP algorithm is one of the two linear solvers evaluated in this paper, see Section 2.3. But it is a bit intricate, because one has to define the homotopy transformation and the recurrence formulae to compute the series. That is why we propose a second linear solver that is an iterative variant of the previous one: it includes the same reduced basis, the same type of reduced problem and the same preconditioner. Convergence acceleration is a family of very efficient techniques to improve the performance of iterative algorithms [START_REF] Brezinski | Padé approximants[END_REF][START_REF] Baker | Padé approximants. Encyclopedia of Mathematics and its Applications[END_REF]. Two acceleration techniques will be used in this paper: the Padé approximants introduced in [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF] to accelerate Taylor series and mainly the modified minimal polynomial extrapolation (MMPE, [START_REF] Jbilou | Vector extrapolation methods, applications and numerical comparison[END_REF]) that is suitable to improve the convergence of sequences and that has about the same efficiency as Padé approximants [START_REF] Jbilou | Vector extrapolation methods, applications and numerical comparison[END_REF][START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF][START_REF] Cadou | Convergence acceleration of iterative algorithms. Applications to thin shell analysis and Navier-Stokes equations[END_REF], see Section 2.2. The resulting algorithm is relatively simple and combines, first a preconditioning by a consistent tangent operator emanating from a previous step, second a projection on a reduced basis, last an acceleration by MMPE.

The paper is organized as follows. In part 2, we state or recall the considered algorithms: the new iterative solver, the linear version of the algorithm from [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF], the acceleration by MMPE and the Newton-Raphson algorithm with arc-length control. In part 3, two classical benchmarks from nonlinear shell analysis permit us to assess the performance of the proposed methods and to compare their efficiency. They will be compared with the preconditioned conjugated gradient associated with various preconditioners. Last, in part 4, the robustness of the proposed procedure is evaluated by considering two path following problems involving the solution of many linear systems.

ALGORITHMS

One deals with the numerical computation of non-linear systems of partial differential equations that depend on a scalar parameter . After discretization, this system can be written in a generic form:

R(U, ) = 0, U ∈ R N , ∈ R (1)
Most of the algorithms to compute non-linear equations as (1) lead to the computation of a family of linear systems in the form:

K i U i = F i with U i ∈ R N , F i ∈ R N , K i ∈ R N * N (2)
In this paper, two algorithms are discussed to solve such a family of linear systems. These algorithms will be based on a preconditioning matrix that is a full tangent matrix coming from a previous step of the computation. We intend to keep the same preconditioner to solve several linear problems.

Since the linear problems are computed by an iterative method, the consistent tangent matrix can be considered without additional cost. That is why we shall choose the Newton-Raphson method to transform the non-near problem (1) into a family of linear problems such as (2): for completeness this non-linear algorithm is recalled in Section 2.4. In the parts 3 and 4, the proposed algorithms will be applied in thin shell analysis that is known to provide ill-conditioned matrices. The chosen class of preconditioning matrices is able to remain very efficient in such a case.

A linear iterative solver involving a reduced subspace technique

We now describe an algorithm to solve one linear problem. For simplicity, the indices are omitted and ( 2) is rewritten as:

K U = F with U ∈ R N , F ∈ R N , K ∈ R N * N (3)
An increment from a trial solution U should satisfy the following incremental equation:

K U = F -K U (4)
The principle of the proposed linear solver is to define an increment in two sub-steps. In the first sub-step, only the projection of the unknown on a reduced subspace is incremented and the projection of the incremental equation ( 4) on the reduced subspace is satisfied exactly. In the second sub-step, the incremental equation ( 4) is satisfied approximately by replacing the consistent matrix K by a given preconditioning matrix K * . Hence the proposed algorithm relies on the choice of a preconditioning matrix and of a reduced subspace. Let E be a given subspace generated by n vectors E i , where n is much smaller than N . We suppose that this basis (E i ) i=1,n is orthonormalized according to a chosen scalar product. From this subspace, one can define a prolongation operator in such a way that a vector of the whole space, X ∈ R N is associated with any vector q of the reduced subspace R n :

X = Pq = n l=1 q l E l (5)
If U is the starting vector of the iteration, the first sub-increment U = P q satisfies a linear problem on the subspace E

k q = f -CU (6)
where k and C are, respectively, a reduced and a coupling matrix and f is a reduced load. All these quantities are defined with the following expressions:

k = t PK P, k ∈ R n * n C = t PK , C ∈ R n * N f = t PF, f ∈ R n (7)
Finally, by considering the initial solution U 0 = 0, we define the iterative method by the following recurrence formulae:

k q = f -CU i-1 K * V = F -K (P q +U i-1 ) U i = U i-1 +P q + V (8)

Convergence acceleration

During the iterations, a sequence of vectors denoted by s has been constructed:

s = {U 1 , . . . ,U i , . . . ,U I } (9)
The previous sequence can slowly converge, or sometimes can be divergent, especially when the preconditioner is not very close to the consistent matrix K . There are several techniques to improve the convergence of sequences: see, the Reference [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF] where a detailed review of these acceleration techniques is presented. Here, we only consider a vector extrapolation method and more precisely the MMPE [START_REF] Jbilou | Vector extrapolation methods, applications and numerical comparison[END_REF]. Within this method, the iterate U I is replaced by the following 'extrapolated' vector:

T I =U 0 + I i=1 a i (U i -U i-1 ) =U 0 + I i=1 a i U i ( 10
)
where the real coefficients a i are computed by solving a small linear system resulting from an orthogonality condition

[M]{a} = {b} (11)
where the matrix [M] and the right-hand side {b} are defined by:

M i j = ( U j+1 -U j )•Y i and b i = -U 1 •Y i with i, j = 1, . . . , I -1 (12)
Within MMPE, the vectors Y i are arbitrary and linearly independent. In our numerical tests, these vectors are the vectors U i .

In practice this acceleration procedure works well, but only if the sequence s is not too large and it has been advised to apply MMPE to a limited number of vectors such as I = 5, 10 or 15 [START_REF] Cadou | Convergence acceleration of iterative algorithms. Applications to thin shell analysis and Navier-Stokes equations[END_REF]. So in this case, once the extrapolated vector has been built with the help of expression [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF] one defines a new linear problem:

K X = F -K T I ( 13
)
Next the new unknown vector X is computed by the same iterative procedure. In the numerical tests, this procedure will be denoted as 'restart algorithm'. Finally, the first presented iterative method combines the recurrence formulae [START_REF] Clemens | Subspace projection extrapolation scheme for transient field simulations[END_REF] and possibly the convergence acceleration method [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF]. In the application, it will be referred as 'proposed method'.

Other acceleration techniques can be deduced from Padé approximants. In the practice, the results from MMPE and Padé approximants are very close [START_REF] Cadou | Convergence acceleration of iterative algorithms. Applications to thin shell analysis and Navier-Stokes equations[END_REF]. It has been established that, in some specific cases, the Padé approximants introduced in [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF] lead to the same mathematical results as MMPE [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF].

A variant based on homotopy and perturbation

In this section, we propose a second linear solver, which is based on homotopy technique, perturbation method and Padé approximants (HPP). This solver is the linear version of an iterative high-order corrector recently introduced in Reference [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF]. For the linear problem (3), the method of Reference [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF] to define a high-order iterative corrector leads to the following system of equations:

-C pen ( t PU -q) = 0 kq +C V -= f K (Pq + V )+P = F U = Pq + V (14)
where is a Lagrange multiplier, ∈ R n , C pen is a pseudo-penalization matrix. As in Reference [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF], this matrix is defined with the following expression:

C pen = Diag(k) (15)
where is a small parameter and is chosen equal to 10 -20 [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF] in all the numerical tests presented in this paper. The reduced matrix k, the coupling matrix C and the reduced load f are the same as in Section 2.1 and are defined by expressions [START_REF] Risler | On the use of Ritz vectors for the solution to nonlinear elasticity problems by decomposition methods[END_REF]. System ( 14) is then modified by introducing an artificial parameter

-C pen ( t PW -q) = 0 kq + C V -= f (1-)K * V + K (Pq + V )+P = F W = Pq + V ( 16
)
where is a real parameter (0 1). The unknowns ( , q, V, W ) are then sought in the form of a truncated integro-power series of the parameter . For example, the asymptotic expansion for vector W is:

W ( ) = W 0 + W 1 + 2 W 2 +•••+ p W p ( 17
)
This series is introduced in Equation ( 16) and equating like power of a sequence of a linear problem is obtained:

Order 0 in 0 = 0 kq 0 = f K * V 0 = F -K Pq 0 W 0 = Pq 0 + V 0 (18) Order p in p = C pen ( t PW p-1 -q p-1 ) kq p = p -C V p-1 K * V p-1 = K * V p-1 -K (Pq p + V p-1 )-P p W p = Pq p + V p (19)
Finally, when all the terms ( i , W i , q i , V i ) are computed, one replaces the asymptotic expansions ( 17) by equivalent rational fractions, called Padé approximants. The method to deduce the Padé approximants is exactly the same as in many previous works (see, for example, [START_REF] Najah | A critical review of asymptotic numerical methods[END_REF]) and this permits to accelerate the convergence of the polynomial approximation [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF].

In this form, the algorithm is not iterative because the series [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF] and the Padé approximants are explicit functions of . The sought solution could be obtained by replacing = 1. Nevertheless, if the radius of convergence is larger than 1, the obtained value U = W ( = 1) can be not satisfactory. The perturbation process can be applied again and this defines an iterative process, see [START_REF] Damil | An iterative method based upon Padé approximants[END_REF][START_REF] Mallil | An iterative process based on homotopy and perturbation techniques[END_REF]. In this paper, we shall not need iterations in the application of this algorithm.

One remarks that the recurrence formulae [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF] are rather similar to formula (8) and to the iterative process of Section 2.1, except for the Lagrange multiplier and the pseudo-penalization that had been introduced in [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF] to stabilize the result of the perturbation process. We no longer discuss this algorithm that is only the linear version of the one presented in [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF]. It is recalled here for comparison with the main algorithm of this paper that is the one of Section 2.1. In the following application, it will be referred as 'HPP algorithm'.

There is a wide recent literature about the association of homotopy and perturbation techniques to solve differential equations or partial differential equations, generally to get analytical solutions and with few terms in the series [START_REF] Andrianov | Asymptotic solution of the theory of shells boundary value problem[END_REF][START_REF] He | Comparison of homotopy perturbation method and homotopy analysis method[END_REF][START_REF] Liao | Beyond perturbation. Introduction to the Homotopy Analysis Method[END_REF]. Homotopy perturbation methods can also be applied in a numerical framework and with many terms in the series.

Applications to non-linear problems

We recall the classical iterative Newton-Raphson method, associated here with an arc-length control. The generic non-linear Equation ( 1) is solved in several steps, each step including a prediction phase and several corrections. In the prediction phase, the increment ( U pred , pred ) is characterized by the following equations:

*R *U 0 U pred + *R * 0 pred = 0 U pred 2 +( pred ) 2 = s 2 (20)
The second one is the arc-length control and s is the given arc-length. In the first one, *R *U 0 is the Jacobian matrix at the starting point of the step. If (U i , i ) is the result of the iteration i, the correction

( U i , i ) is characterized by: *R *U i U i + *R * i i = -R| i U i • U pred + i pred = 0 (21)
The second one is the so-called normal plane condition. As it is well known, this non-linear process requires the solution of N iter +1 linear problems, all the matrices being different.

In this paper, these linear problems will be solved by the iterative methods of Section 2.1 or Section 2.3. The preconditioning matrix will be the Jacobian matrix of a previous step and the same matrix will be used for many linear problems as [START_REF] Andrianov | Asymptotic solution of the theory of shells boundary value problem[END_REF] or [START_REF] He | Comparison of homotopy perturbation method and homotopy analysis method[END_REF]. Several strategies will be discussed to manage the preconditioning matrix: for instance in the next part, the prediction problem [START_REF] Andrianov | Asymptotic solution of the theory of shells boundary value problem[END_REF] will be solved by a direct method and this factorized matrix will be the preconditioner for all the iterations of the same step. Hence this leads to an algorithm with one matrix factorization per step.

NUMERICAL RESULTS

Two numerical tests

The previous linear solvers are now applied to linear problems coming from a predictor-corrector algorithm, the Newton-Raphson method. Two main numerical tests are considered and described in Figures 1(a) and(b). The non-linear response curves are plotted in Figures 2(a) and(b). The first example is a cylindrical shell with two diametrically opposed rectangular cut-outs [START_REF] Riks | Some computational aspects of the stability analysis of nonlinear structures[END_REF]. The corresponding geometric and material characteristics are given in Figure 1. Owing to the symmetries, only one octant of the structure is considered and the mesh involves about 1600 classical triangular DKT18 elements [START_REF] Batoz | Modélisation des structures par éléments finis[END_REF] (5190 d.o.f.). The second example is an open cylinder pulled out by two diametrically opposite point loads [START_REF] Sze | Popular benchmark problems for geometric nonlinear analysis of shells[END_REF]. The geometric and material characteristics are given in Figure 1(b). The finite element for this second example is an eight nodes quadrilateral element that is presented in Reference [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF]. An additional variable accounts for the strain variations throughout the thickness (EAS concept). Owing to the symmetries, only one octant of the cylinder is modelled by 900 elements (30×30). The corresponding number of d.o.f. is then 17 000.

The non-linear response curves have been obtained by the Newton-Raphson method. They are pictured in Figures 2(a) and (b). The prescribed step length of this incremental/iterative method is chosen in such a way that at each step, 2-3 corrections are necessary to get the desired accuracy. For these first tests, the algorithm has been designed in such a way that only one matrix is triangulated at each step. The prediction step is carried out by using a direct solver. The corrections are achieved by the iterative linear solvers presented in Part 2, the preconditioning matrix K * being the triangulated matrix of the prediction step. The prolongation operator P is built with the already computed solutions U , coming from the previous prediction-correction steps. This operator is then upgraded at each prediction-correction step. Consequently, the reduced basis size 'n' increases with the number of prediction-correction steps. In all the numerical tests presented in this study, we consider that the linear solvers have converged when the following criterion is satisfied:

K U -F F < lin ( 22
)
where • is the Euclidian norm and lin is a small parameter (10 -8 lin 10 -4 ).

Evaluation of the proposed solver

The proposed iterative solver is first evaluated by studying its convergence within the first corrective step of a Newton-Raphson algorithm, the preconditioning matrix and the prolongation operator being defined in Section 3.1. The numerical results will be compared with those obtained with the preconditioned conjugate gradient method, the preconditioning matrix K * being the same. In the next Figures 3456, the convergence is represented by curves relating the logarithm of the residual of the linear problem versus the iteration number. The results obtained with the 'proposed method' and with the 'preconditioned conjugate gradient' method will be referred as 'PM' and 'PCG', respectively.

A first example is presented in Figure 3, the physical model being the cut-out cylinder. In this case, we have also considered the linear solver based on the HPP method (Equation ( 19), denoted by 'linear solver' from Reference [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF] in Figure 3). First, one observes that the number of iterations is very small: one can get a residual of 10 -5 after only 4 iterations. This is due to the choice of a good preconditioner and also the efficiency of the proposed algorithm. Second remark: the proposed method converges more rapidly than the preconditioned conjugate method, with the same preconditioning matrix. Third, the proposed method and the HPP solver of the Section 2.3 lead to the same convergence curve: this is in accordance with the results of Reference [START_REF] Damil | Mathematical and numerical connections between polynomial extrapolations and Padé approximants[END_REF], where one has established the mathematical equivalence between a sequence accelerated by MMPE and the 'homotopy-perturbation-Padé' method. Consequently, the HPP will be no longer considered in the following numerical tests. Another example is presented in Figure 4, where the sequence (9) diverges, despite of the great size of the reduced problem (n = 40). The initial sequence has been replaced by the extrapolated vector [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF], which leads to a sequence that converges more quickly than the PCG. Hence, it is very interesting to use a convergence acceleration technique that has moreover a negligible computational cost [START_REF] Cadou | Convergence acceleration of iterative algorithms. Applications to thin shell analysis and Navier-Stokes equations[END_REF].

The influence of the size of the basis (5 n 30) is discussed in Figure 5. In this case, the accuracy is not improved beyond 20 vectors in the basis. A similar result has also been noted in other tests with a larger number of unknowns. This can be explained by the fact that the gap between the operator K * , the prediction matrix and the consistent correction matrix K is small. In this case, this leads to a good preconditioning matrix and the influence of the reduced problem on the convergence is less important.

This figure also illustrates one of the characteristics of the MMPE. Indeed, the accelerated sequence can become unstable or the convergence is not improved when the number of vectors composing the sequence is high. In this case, a solution is to apply the restart algorithm as defined in Section 2.2. In Figure 6, one can see that without using the restart algorithm the proposed The proposed method (PM) is compared with the conjugated gradient associated with three preconditioners: the same as in PM (PCG) and the incomplete Crout triangulations of level 0 and 1 (IC[0] and IC [START_REF] Cochelin | Asymptotic-numerical method and Padé approximants for non-linear elastic structures[END_REF]). Required accuracy lin = 10 -4 for the linear solver and nl = 10 -3 for the non-linear solver. Each Newton-Raphson step includes the prediction and 2 or 3 iterations. In PM, the size of the reduced basis varies between 15 and 37.

Table II. Average number of iterations for the proposed method and the conjugated gradient associated with the same preconditioner. The required accuracies are the same as in Table I. Each Newton-Raphson step includes the prediction and 2 or 3 iterations. The size of the reduced basis varies between 17 and 50.

method does not reach the desired accuracy (≈ 10 -5 ). Whereas, when the restart algorithm is used after 10 or 15 iterations, the proposed linear solver converges.

Two full path following computations

Let us now apply the proposed iterative solver for the whole computation of non-linear solution branches in order to demonstrate its robustness. We consider the two tests described in Section 3.1: cut-out cylinder and pull-out open cylinder pulled by two diametrically opposite sides. The prediction step is done with a direct method and the iterative solver is used in the correction phase. The accuracy on the linear and non-linear problems are, respectively, chosen equal to 10 -4 and 10 -3 . The prediction step length is chosen in such a way that 2 or 3 iterations per step are needed. This leads to 10 steps to solve the cut-out cylinder problem and 23 for the open cylinder (the sought response curve is in the range 0<P R/D<100).

Owing to the efficiency of the preconditioning matrix, the proposed method (PM) and the preconditioned conjugate gradient (PCG) converge in any case. The average number of iterations for these linear solvers is reported in Tables I andII. With the proposed method, the number of iterations is very small, generally lower than 10. The conjugate gradient associated with the same preconditioning matrix (PM) is more expensive: the number of iterations is twice larger in the first example and four times larger in the second one that has a larger number of degrees of freedom.

For comparison, we have also evaluated the conjugate gradient associated with an incomplete Crout factorizations [START_REF] Joly | Quelques méthodes classiques de résolution de systèmes linéaires[END_REF] of the matrix of the prediction step. In Table I, this preconditioner is referred to as IC [m], where m is the 'level' of the incomplete factorization. As it is well known, this technique requires much less memory and much less computation time per iteration than the previous one. The drawbacks are also known and corroborated by our numerical results. Indeed the number of iterations is much higher than with a complete factorization: in the first test (cut-out cylinder), the number of iterations is about 200 with IC[0] and 100 with IC [START_REF] Cochelin | Asymptotic-numerical method and Padé approximants for non-linear elastic structures[END_REF]. In the second test (pull-out cylinder), the algorithm diverges. Likely this could be removed by introducing an appropriate stabilization parameter, but this confirms the lack of robustness of incomplete factorization techniques with ill-conditioned problems as in shell analysis.

PREDICTION-CORRECTION STRATEGIES

Various prediction-correction strategies can be established from the linear iterative solver presented in Part 2. As recalled in Section 2.4, a Newton-Raphson step leads to N iter +1 linear problems. These linear problems will be solved, either by a direct method, either by the iterative solver. In the two first strategies that are discussed in this part, only one matrix per step is triangulated and the corresponding linear problem is solved by the direct method, while all the other problems are solved iteratively, the preconditioner being the last triangulated matrix. More precisely, the prediction problem is solved by the direct method in the first strategy and the linear problem of the first iteration in the second strategy. They will be referred as 'one matrix algorithms', because they require only one matrix triangulation per step. By comparison, the Newton-Raphson method requires classically N iter +1 matrix triangulations or, at least, two in the so-called 'MNM'. With the third strategy, one tries to establish a matrix-less algorithm. More precisely, the direct solver is no longer used at each step, but only when the number of iterations increases too much. In the following numerical experiments, we limit ourselves to 10 iterations before applying the convergence acceleration process (i.e. I = 10 in formula [START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF]) and to 2 cycles '10 iterations + MMPE'. In other words, a new matrix is factorized as soon as the number of iterations is beyond 20. These three strategies are sketched in Table III.

The efficiency of these strategies is assessed by two examples. The first one is the open cylinder (Figure 1). The second one is another classical benchmark: a hemispherical shell with a hole and loaded by pinching forces (F = 1.) [START_REF] Büchter | Shell theory versus degeneration-a comparison in large rotation finite element analysis[END_REF]. The geometric and material descriptions are given in Figure 7. The response curves of the two loaded points are plotted in Figure 7. The mesh is the same as the one used in [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF] (108 eight-node quadrilateral elements with 2166 d.o.f.). In our computations, 30 prediction-correction steps are needed to obtain the non-linear solution branch (Figure 7). The corresponding total number of linear problems to be solved is equal to 139. In Reference [START_REF] Zahrouni | Computing finite rotations of shells by an asymptotic-numerical method[END_REF], this branch had been obtained by the commercial code Abaqus with a cost of 104 tangent matrix triangulations, which is relatively close to our results. For the pull-out cylinder, the number of prediction-correction steps is 23 and 86 linear problems are to be solved. Because the solver is based on a reduced basis, the first three steps are computed by classical direct method, which permits to compute at least 10 basis vectors. The demanded accuracies are the same as in Section 3.3 : 10 -3 for the non-linear iterations and 10 -4 for the linear ones.

The numerical costs depends on the number of iterations and of matrix factorizations. These quantities are reported in Tables IV and V, for the two tests and for the three computational strategies. First, one observes the strategy 2 is much more efficient than the other 'one matrix algorithm' (strategy 1), the number of iterations being about half. This can be explained by the fact that the matrix of the first Newton iteration is a better preconditioner (for the next iterations) than the prediction matrix. Another consequence: the prediction problem ( 20) is also an 'easy' problem that is generally solved with less than 5 iterations. As for the 'matrix-less algorithm', it permits to reduce hugely the number of matrix factorizations: for instance, in the pull-out cylinder problem, only 4 factorizations are necessary, to be compared with 20 for the 'one matrix algorithms' and 76 for the direct solver. The decrease of the number of matrices does not lead to a consequent increase of iterations. More surprisingly, in some cases (Table IV) there are even less iterations than with the strategy 1. This can also be explained by the fact that in strategy 3, the triangulated matrix comes always from the Newton's correction phase, and not from the prediction step, as in strategy 1.

Hence an algorithm based on sparse matrix factorizations, as strategy 3, seems the most attractive technique, especially for large-scale problems since the cost of the matrices is much larger than the cost of one iteration.

CONCLUSION

In this paper, a new linear iterative solver has been proposed, that is designed to solve the numerous linear equations arising in the solution of a non-linear system. It associates a reduced basis technique, a full preconditioning matrix and an accelerating convergence method, the modified minimal polynomial extrapolation (MMPE). It converges more rapidly than the preconditioned conjugate gradient method and as rapidly as a recently introduced algorithm [10] that relies on homotopy transformation, perturbation technique and Padé approximants. The process is robust and leads to convergence with a small number of iterations.

The reliability and the efficiency of the procedure have been proved by several shell buckling non-linear computations within a Newton-Raphson framework. When it is applied to solve nonlinear problems, one gets algorithms that need one matrix factorization per step or one factorization for several Newton-Raphson steps.

This strategy could also be applied in the asymptotic numerical method (ANM) [START_REF] Cochelin | Asymptotic-numerical method and Padé approximants for non-linear elastic structures[END_REF] to compute several steps of the continuation methods with a single matrix factorization. Note that, within ANM, it is very easy to build up a relevant reduced subspace including a rather large number of vectors. There are also possible applications in the computation of non-linear dynamics problems with an implicit scheme.

In this paper, the number of iterations and the speed of convergence have permitted us to discuss the efficiency of the algorithm. Of course, a direct measure of the computation time would be more convincing, but only for tests with many degrees of freedom on parallel computers and with the help of modern direct solvers as MUltifrontal Massively Parallel sparse direct Solver (MUMPS [START_REF]MUltifrontal Massively Parallel sparse direct Solver (MUMPS)[END_REF]). We are now working in this way.

Figure 1 .Figure 2 .

 12 Figure 1. Geometrically non-linear shell problems. Tests 1 and 2: (a) cut-out cylinder and (b) pull-out of an open cylinder.

Figure 3 .

 3 Figure 3. Convergence curve. Comparison between the proposed solver (Equation (8)), the linear solver coming from Reference[START_REF] Cadou | Projections technique to improve high-order iterative correctors[END_REF] and the preconditioned conjugate gradient (PCG). The size of the basis is n = 20. Cut-out cylinder.

Figure 4 .

 4 Figure 4. Convergence curve. Comparison between the extrapolated vector (expression 10) and the initial sequence (expression 9) with n = 40. Cut-out cylinder.

Figure 5 .

 5 Figure 5. Convergence curves for several sizes of the basis. Cut-out cylinder.

Figure 6 .

 6 Figure 6. Comparison between the proposed method+MMPE with or without restart. The size of the reduced basis is n = 5. Cut-out cylinder.

Figure 7 .

 7 Figure 7. Example of the pinched hemisphere: (a) geometry and material descriptions and (b) displacement of points A and B versus the load parameter. The curve is obtained with 30 steps of the Newton-Raphson method.

Table I .

 I Average number of iterations for various linear solvers.

			Cut-out cylinder
			Average number of iterations
	Newton-Raphson Step	PM	PCG	PCG IC[0]	PCG IC[1]
	4	6	11	170	110
	5	6	11	176	111
	6	7	12	178	113
	7	8	17	184	121
	8	9	20	193	124
	9	10	21	198	130
	10	9	21	198	135

Table III .

 III Description of the three prediction-correction strategies. DS and PM design, respectively, a direct solver and the proposed linear solver.

		Strategy 1	Strategy 2	Strategy 3
	Prediction	DS	PM	DS or PM
	Correction	PM	Iteration 1 : DS	DS or PM
			Iteration >1 : PM

Table IV .

 IV Evaluation of the three strategies to update the preconditioner.

				Pinched hemisphere
		Average number of iterations per step		Number of matrices
	Steps	Strategy 1	Strategy 2	Strategy 3	Strategy 1	Strategy 2	Strategy 3
	4-8	24	8	14	5	5	4
	9-13	10	4	13	5	5	2
	14-30	6	2.5	1 2	1 7	1 7	4

Table V .

 V Comparison of the three strategies to update the preconditioner. Average number of iterations per step.

		Pull-out cylinder
		Strategy 1	Strategy 2	Strategy 3
	Iterations per step	7	2.75	13
	Matrices	20	20	4