Sum-of-squares of polynomials approach to nonlinear stability of fluid flows: an example of application

With the goal of providing the first example of application of a recently proposed method, thus demonstrating its ability to give results in principle, global stability of a version of the rotating Couette flow is examined. The flow depends on the Reynolds number and a parameter characterizing the magnitude of the Coriolis force. By converting the original Navier–Stokes equations to a finite-dimensional uncertain dynamical system using a partial Galerkin expansion, high-degree polynomial Lyapunov functionals were found by sum-of-squares of polynomials optimization. It is demonstrated that the proposed method allows obtaining the exact global stability limit for this flow in a range of values of the parameter characterizing the Coriolis force. Outside this range a lower bound for the global stability limit was obtained, which is still better than the energy stability limit. In the course of the study, several results meaningful in the context of the method used were also obtained. Overall, the results obtained demonstrate the applicability of the recently proposed approach to global stability of the fluid flows. To the best of our knowledge, it is the first case in which global stability of a fluid flow has been proved by a generic method for the value of a Reynolds number greater than that which could be achieved with the energy stability approach.


Introduction
Hydrodynamic stability is the field that investigates the transient effects of an initial perturbation of a known steady flow. The area has attracted the attention of many researchers and is closely related to the study of transition to turbulence [1,2].
Using Lyapunov stability theory, a steady flow can be proved to be stable with respect to perturbations of arbitrary amplitude by constructing a Lyapunov functional V [u], which is a positive-definite functional of the velocity perturbation u that decays monotonically on any non-zero solution u(t, x) of the Navier-Stokes equations [3].
Usually, the Lyapunov functional is chosen to be the perturbation energy E = u 2 /2, where the norm is defined as an integral of |u| 2 over the flow domain. Then the problem of proving that E is a Lyapunov functional reduces to a tractable linear eigenvalue problem [4,5]. However, the resulting estimates of the global stability range, usually expressed by the energy stability limit Reynolds number Re E , could be very conservative in the sense that Re E is generally far below the maximum Re for which the flow is globally stable.
Recently, a method was proposed by Goulart & Chernyshenko [6] for exploiting the sumof-squares (SOS) decomposition [7,8] to construct polynomial Lyapunov functionals differing from E, thus extending the range of Re in which the flow can be proved to be globally stable. In this approach, the Navier-Stokes equations are first reduced to a finite-dimensional uncertain dynamical system, that is a system of ordinary differential equations (ODEs) with right-hand side containing terms for which only bounds, but not exact expressions, are available. For incompressible flows both the right-hand side of the ODEs and the bounds are polynomial. The corresponding Lyapunov stability condition can then be reduced to a condition of positive definiteness of other certain polynomials [6]. Noting that the condition of a polynomial being positive-definite can be replaced by a stronger, but more tractable numerically, condition of the polynomial being a SOS of other polynomials, an admissible Lyapunov functionals can be found using the polynomial SOS optimization approach [7,9].
The SOS technique has been applied in stability analysis of constrained ODEs [10], hybrid systems [8] or time-delay systems [11], but there are few results for partial differential equations. The relevant publications are [12] and [13], where SOS-based algorithmic methodologies are presented for the analysis of systems described by certain types of parabolic partial differential equations. It was shown how certain Lyapunov structures could be constructed to prove stability using transformations defined through integration by parts. It is worth noting that in both [12,13], the partial differential equations are considered directly, which is different from the construction process of a Lyapunov functional in [6].
Application of SOS of polynomials to fluid flows goes beyond nonlinear stability. After an overview of nonlinear stability applications including a short announcement of a part of the results of this work, applications for deriving rigorous bounds on time-averaged characteristics of turbulent flows and applications to flow control are discussed in [14].
While [6] provides a full theoretical description of the new approach, it was applied there only to a truncated Galerkin approximation rather than to the full Navier-Stokes equations, thus leaving open the question of whether there is at least one fluid flow for which this method will work. Such an example is given in this paper. To the best of our knowledge, this study is the first case in which global stability of a fluid flow has been proved by a generic method for the value of a Reynolds number greater than that which could be achieved with the energy stability approach.

Problem formulation
Our goal is to demonstrate that the method proposed in [6] can actually be used to prove the stability of a fluid flow for the values of the Reynolds number above the energy stability limit. This section describes briefly the method and the flow to which it is to be applied for achieving this goal.

(a) The method
An unsteady flow of incompressible viscous fluid in a given domain with time-independent boundary conditions is considered. The perturbation velocity u(t, x), defined as the deviation of the instantaneous flow velocity from the steady solutionū, is governed by the Navier-Stokes equations with an additional linear term: The linear operator A depends on the flow in question. It is convenient to leave it in a compact general form here. This formulation is a little more general than in [6], where A had a particular form, but all the results of Goulart & Chernyshenko [6] apply with obvious minor modifications, which are made without further comments in the summary of the method given in this section. The perturbation velocity u is subject to homogeneous boundary conditions. The steady flow can be proved to be stable with respect to perturbations of arbitrary amplitude by constructing a Lyapunov functional V[u].

(i) The Lyapunov functional
For this purpose, in [6] the perturbation velocity is represented as where the finite Galerkin basis fields u i , i = 1, . . . , k, are an orthonormal set of solenoidal vector fields with the inner product w 1 , w 2 defined as the integral of w 1 · w 2 over the flow domain V, the residual perturbation velocity u s is solenoidal and orthogonal to all the u i , and both u i and u s satisfy the homogeneous boundary conditions. The Lyapunov functional is sought in the form 1 V[u] = V(a, q 2 ), where a = (a 1 , . . . , a k ), and q 2 = u s 2 /2 = u s , u s /2. For V[u] to be a Lyapunov functional, the function V(a, q 2 ) : R k × R → R should be positive-definite, and its value should decrease monotonically towards zero along all possible non-zero solutions of (2.1): V(a, q 2 ) > 0 and dV/dt < 0, for all (a, q 2 ) = 0. With the use of (2.1) and (2.2), the latter condition can be rewritten as [6] dV where the components of the vector f are Re u s , ∇ 2 u s − u s , u s · ∇ū +ū · ∇u s (2.5) and χ (u s , a) = u s , g j a j , g j = 2 Re whereū is the steady flow the stability of which is studied, and the vector-valued functional Θ(u s , a) = Θ a (u s ) + Θ b (u s , a) + Θ c (u s ) has the components where h i = (1/Re)∇ 2 u i +ū · ∇u i − ∇ū · u i and h ij = u j · ∇u i − ∇u j · u i . The notation used can be clarified by the Einstein equivalent of the formula for h i : i and x m are the mth components of the vectors h i , u i and x, respectively. Constructing V satisfying (2.3) and V > 0 for all (a, q 2 ) = 0 would prove the global stability of the flow under consideration. However, (2.3) involves u s while V is a function of a and q 2 only. Hence the next step is required [6].

(ii) The bounds
The terms in (2.3) dependent on u s can be bounded by functions of a and q 2 . Note that χ , Θ ai , and Θ bi are linear functionals of u s , while Γ and Θ ci are quadratic functionals of u s .
For Θ ai defined by (2.7) the Cauchy-Schwarz inequality gives |Θ ai (u s )| ≤ u s h i = √ 1/2|q| h i . Obtaining the tight bound requires projecting h i onto the solenoidal subspace and a small modification to account for the boundary conditions [6]. Other linear functionals can be bounded similarly.
The linear functional χ is a special case. If the basis u j is chosen to consist of the eigenfunctions of the classical energy stability problem [15] forū, then χ ≡ 0 [6]. This can become obvious if one notices the similarity between the energy stability operator [15] and (2.6). We will use such a basis in this study. Accordingly, even though we refer to the following formulae as obtained in [6], in fact they are simplified versions that we derived with an additional assumption χ = 0. In most cases, this can be done by simply omitting some terms. The readers wishing to use the method with the basis u i in which χ = 0 should refer to [6] rather than to the formulae below.
The bounds on the quadratic functionals can be obtained by maximizing the functionals subject to a constraint q 2 = 1. This reduces to a linear eigenvalue problem in a fairly standard way. The resulting bound for Γ has the form [6] Γ (u s ) ≤ κ s q 2 . (2.8) Note the link to the energy stability problem following from the form of (2.5): finding κ s reduces to the energy stability problem with an additional constraint u s , u i = 0 ∀ i. When the linear eigenvalue problem has a discrete set of eigenvalues, only a finite number of them can be positive [15] and, hence, κ s can be made negative by selecting a suitable basis u i . Rather than using the generic approach to bounding Θ ci as proposed in [6], in appendix A we derive an explicit tight bound for it.
Putting together all the above gives the following set of bounds where the particular values of κ s and the coefficients of the quadratic polynomial p(a, q 2 ) depend on the particular flow in question and on the selection of the particular eigenfunctions of the corresponding energy stability problem to be used as the finite basis u i . This allows formulating the stability analysis problem as a SOS optimization problem.
(iii) Reduction to a polynomial sum-of-squares optimization Since the term Γ (u s ) in (2.3) is upper-bounded by a negative-definite function κ s q 2 but is not lower-bounded, one has to impose an additional requirement that the candidate function V satisfy The condition (2.11) is difficult to use because both V and p enter it nonlinearly. This can be circumvented via introduction of additional variables. As shown in [6], (2.11) is equivalent to where H(a, q 2 ) = −s 1 (a, q 2 )p(a, q 2 ) s 2 (a, q 2 ) · p(a, q 2 ) s T 2 (a, q 2 ) · p(a, q 2 ) −s 1 (a, q 2 )I (2.13) and In summary, if one is able to construct a function V simultaneously satisfying the three conditions V > 0, (2.10) and (2.12), for all (a, q 2 = 0) then the flow in question is globally stable.
If V is sought for in the form of a polynomial then checking each of these three conditions amounts to checking the global non-negativity of a polynomial function, which is known to be NP-hard; see [9]. However, a sufficient, and numerically tractable, condition for global non-negativity of a polynomial is that it can be written as a SOS of other polynomials. Accordingly, a sufficient condition for joint satisfaction of our three conditions is that where Σ l represents the set of all SOS polynomials in R l , and the positive-definite polynomial functions i (c) = j ij c 2 j , ij ≥ 0, j ij > 0 are used in place of the vector-value conditions z = 0 and/or (a, q 2 ) = 0.
Existence of a function V satisfying the conditions (2.14) can be checked via the polynomial SOS approach of [7,9], which amounts to solving a convex optimization problem in the form of a large semi-definite program. For a prescribed degree of the candidate polynomial V, the solution to that problem either provides an explicit expression for V or states that such V does not exist. The SOS approach is well known and will not be described here.
(iv) Uncertain system interpretation Substituting the partial Galerkin expansion (2.2) into the Navier-Stokes equations (2.1) and projecting onto u i gives da dt = f(a) + Θ, (2.15) and the easily obtainable equation for the energy q 2 of the residual field u s is where Θ, Γ and χ are defined in §2a(i) as functionals of u s and functions of a. One can, however, allow Θ, Γ and χ in (2.15) and (2.16) to assume any values as far as they satisfy the set of known bounds defined in terms of u s and a, such as (2.9). In this sense, ((2.9), (2.15), (2.16)) is an uncertain dynamical system. The solution of this system is therefore not unique. However, if all the solutions of ((2.9), (2.15), (2.16)) tend to zero as time tends to infinity, then the solution of the Navier-Stokes system also tends to zero. The stability conditions described above are in fact the stability condition for this uncertain system. This turns out to be quite useful in understanding and interpretation of the further results. Note that V(a, q 2 ) is a Lyapunov function for the uncertain system and at the same time a Lyapunov functional for the full Navier-Stokes equations.
If one takes V = a 2 /2 + q 2 then (2.11) reduces to the classic energy stability problem. Hence, this approach [6] is guaranteed to give at least as good results as the energy stability approach. In order to demonstrate that it is capable of providing a better result, it has to be applied to a particular flow.

(b) Double-periodic rotating Couette flow
Given the novelty and complexity, and the relatively demanding computational requirements for solving the semi-definite programming problems stemming from polynomial SOS optimization, for the first application of the stability analysis method of Goulart & Chernyshenko [6], it is desirable to select as simple a flow as possible. The particular flow we select is a version of the famous rotating Couette flow between two co-axial cylinders.

(i) Governing equations
The gap between the cylinders is assumed to be much smaller than the cylinder radius. A local Cartesian coordinate system x = (x, y, z) is oriented such that the axis of rotation is parallel to the z-axis, while the circumferential direction corresponds to the x-axis. Only flows independent of x are considered. The flow velocity is represented as (y + u, v, w), so that u = (u, v, w) is the velocity perturbation andū = (y, 0, 0) is the equilibrium flow. Under these assumptions, the governing equations are [16,17] ∂u ∂t where Ω is a non-dimensional parameter characterizing the Coriolis force, Re is the Reynolds number and p is pressure. More compactly, (2.17) can be written in the vector form (2.1) with In §2a, A was assumed to contain only the terms depending on the base flowū, while here the terms with Ω are present. Conveniently, their presence does not affect any of the formulae in §2a, mostly because these terms correspond to the Coriolis force and thus do not enter the energy equation.
For simplicity, the flow is assumed to be 2π -periodic in y and z, u and v are assumed to be odd in y and even in z, while w is assumed odd in z and even in y:

(ii) Stability properties of the flow
We first apply the well-known energy stability approach. Setting the Lyapunov functional as the perturbation energy E = u 2 /2 leads to a linear eigenvalue problem [15]. For (2.17)-(2.18), the resulting eigenfunctions e n,m (x), as can be verified by direct substitution into that eigenvalue problem, are: Note that for this flow, conveniently, neither the eigenfunctions nor the eigenvectors depend on Ω. The eigenvalue λ 1,−1 is positive for Re > 4 √ 2, with all other eigenvalues less than λ 1,−1 . Hence, the energy stability limit is One can show that the flow becomes linearly unstable for 0 < Ω < 1 and Note that Re L = Re E for Ω = 1 2 . For the classical case of no-slip conditions at the wall, it has been proved [18] that the linear stability and the global stability limits coincide. For the double-periodic flow considered here the same is true, too. This can be proved by the same method as in [18], which amounts to selecting the Lyapunov functional in the form V = (λu 2 + v 2 + w 2 ) dy dz (2.22) and adjusting the constant λ to get as best stability limit as possible. Similar results were obtained in [19] using a more systematic approach involving convex optimization technique. These approaches cannot be applied to an arbitrary flow, unlike the method proposed in [6]. On the other hand, the form of the Lyapunov functional V = V(a, q 2 ) might be too restrictive to obtain the same results as with (2.22). Hence, for the flow we are considering both the energy stability limit and the actual global stability limit are known, giving the framework for considering the performance of the method [6].

Application of the method to the double-periodic rotating Couette flow (a) Selection of the Galerkin basis fields u i
In solving the SOS feasibility problem (2.14), the uncertain parts of the system (2.15)-(2.16) are characterized only by their bounds. Further, those bounds are linked directly to the Galerkin basis fields u i , i = 1, . . . , k, as it is evident from (2.5)-(2.7). The choice of the Galerkin basis fields is therefore crucial. If k is too small, then the dynamic model of the system becomes over-simplified and does not adequately capture the salient features of the flow. Consequently, it might be difficult to achieve a better stability result than Re E . If k is too large, the computational cost in SOS analysis will be prohibitively high. We, therefore, aim to select a limited number of Galerkin modes from amongst the eigenfunctions (2.19) of the system (2.1) in such a way that the best stability bound can be obtained.
It is shown in appendix B that for k = 4 with modes e 1,±1 , e 1,±2 there is no polynomial Lyapunov function for the uncertain system (2.15)-(2.16) at Re > Re E . Note that in this case f is linear, while the proof of the existence of a polynomial Lyapunov functional given in [6] explicitly   requires f to be quadratic, which of course can always be achieved by taking more Galerkin modes. For our particular flow, we therefore consider k ≥ 6. It is difficult to foresee the effect of mode selection on the system stability via the changes in f. It is also difficult to foresee the effect of mode selection on p(a, q 2 ). However, minimizing κ s is clearly beneficial. According to the definition (2.8), κ s can be minimized simply by selecting the k modes of the finite basis to consist of the eigenfunctions with the largest eigenvalues λ n,m . Hence, we select the Galerkin modes by following such a criterion.

(b) Global stability of the flow
Proving that the flow is globally stable at a particular value of Re does not prove that it is globally stable for any other Re. However, for the flow in question (see §2b(ii)) there exists a global stability limit Re G (equal to the linear stability limit Re L for this flow) such that the flow is globally stable for all Re < Re G . Hence, any Re for which the flow is globally stable gives a lower bound for Re G .
In particular, the SOS stability bound for the uncertain system Re SOS ≤ Re G . The largest possible value of Re SOS was obtained by trial end error. For a given Re, we try to find V satisfying the SOS constraints (2.14). If this is successful, we increase Re by δRe and repeat the trial. To do this, we assume some partially fixed structure of the candidate function V, that is the degree and the values of a part of the coefficients, and consider the remaining coefficients as the decision variables. We then tune the decision variables to satisfy (2.14), using the SOS package YALMIP [20], under which the decision variables were found using the semi-definite program (SDP) solver MOSEK [21]. Prior to solving the SDP, the SOS problem was pre-processed by the linear program solver GUROBI [22] in order to reduce and simplify the SOS program [23].

(c) The best bound
The best result was achieved using the Galerkin mode set K 1 . For K 1 , (2.4) gives For the flow in question it turns out that Θ a = 0. Further, following the bound evaluation procedure introduced in §2a(ii), we have and The best stability bound was obtained with the following 4th-degree candidate function candidate: V(a, q 2 ) = P(a, c 4 ) + P(a, c 2 )q 2 + αq 4 , (3.4) where α ≥ 0 and P(a, c i ) is a general ith degree polynomial in a, with c i denoting the coefficient vector. Noting that the constant term and the linear terms are obviously redundant in V(a, q 2 ) owing to the first SOS constraint in (2.14), they are eliminated in advance. The Lyapunov functional with the structure of (3.4) can be adjusted by tuning the decision variables c 2 , c 4 and α. Since [c 2 , c T 4 , α] ∈ R 232 and [a, q, z] ∈ R 14 , there are 232 parametric variables and 14 independent variables for the SOS optimization. The parameters ij required to construct the functions i in (2.14) were fixed to 1j = 1 × 10 −5 and 3j = 1.5 × 10 −8 .
We then performed a trial-and-error procedure with δRe = 0.01. Figure 1 shows the result. In the range Ω ∈ (0.2529, 0.7471), Re SOS coincides with Re L . Hence, in that range the method [6] produces the exact result. Outside this range the method gives only the lower bound for the true global stability limit Re G = Re L . This lower bound is still better than the bound Re E obtained by the standard energy stability analysis.
(3.6) Some terms in (3.6) were grouped manually for brevity. In the original expanded form, (3.6) contains only 48 monomial terms in comparison to 232 possible monomials in (3.2). For different values of Ω the monomials themselves and their coefficients show some variation. Note that the first term is proportional to the square of the total energy of the perturbation u. This was not imposed by us but was obtained as a result of the calculations. The likely form of the Lyapunov functional is discussed in [6].
These results demonstrate the feasibility of the method [6], thus achieving the primary goal of the present study.

(d) Further observations
We remind that the stability of the flow at a given Re does not automatically imply stability at all smaller Re, even though such behaviour is typical. To simplify the exposition, we assume that our system does possess this typical behaviour. In fact, our results were obtained for a discrete, although reasonably dense, set of Re values.
The SOS stability limits obtained for the flow in question always turned out to be of the form min{Re L (Ω), Re E + Re}, so that the constant Re can be used as the measure of the quality of the bound. The stability bound Re SOS,T of the truncated system (that is the system (2.15) with Θ = 0) is also presented in figure 1. It is Re SOS,T = min{Re L , Re E + 5.52}. (3.7) One can see that the uncertain term significantly reduces the SOS stability bound, namely, from Re SOS,T = 5.52 in (3.7) to Re SOS = 0.85 given by (3.5). Another immediate observation is that over a certain range of Ω the SOS bound for the stability margin of the truncated system is less than the true global stability margin of the full system, which coincides with Re L (Ω). These differences can be due to truncation, uncertainty and/or the use of SOS relaxation, in particular, the use of a not high-enough degree of the candidate Lyapunov function. In this section, we summarize the related observations. (i) Stability limits of the full system and the truncated system We observed that a better SOS stability limit of the truncated system does not necessarily imply a better SOS stability limit of the uncertain system.
For instance, we considered another set of six Galerkin modes: K 1 = {e 1,0 , e 1,±1 , e 2,0 , e 2,±1 }, which differs from the Galerkin set K 1 in that the mode e 1,−2 is replaced by e 2,1 . The SOS stability limit for the truncated K 1 system was found to be Re SOS,T = Re L (Ω), while it is Re SOS,T = min{Re L (Ω), Re E + 5.52} for K 1 . On the other hand, the SOS stability limits for the corresponding uncertain systems were found to be Re SOS = min{Re L , Re E + 0.23} for K 1 and Re SOS = min{Re L , Re E + 0.85} for K 1 .
This means that the selection of the Galerkin modes cannot be made on the basis of analysing the truncated system alone.
(ii) More Galerkin modes Intuitively, it seems that increasing the number of Galerkin modes explicitly taken into account should improve the resulting SOS stability limit. The observations show a different picture.
For the Galerkin mode set K 2 , which includes the mode set K 1 and two additional modes e 2,−2 and e 2,1 , the uncertain fluid dynamical system is of the 8th order. In this case, (3.10) The expression for f is not shown here due to its large size. For a comparison with the 6mode case, we consider the same 4th-degree Lyapunov function candidate (3.2). Owing to the increase in the number of the modes, a direct solution of the SOS optimization problem, (2.14) requires substantially greater computational effort than in the 6-mode case. For K 2 , there are 532 parametric variables and 18 independent variables in the SOS optimization. Solving the SOS problem (2.14) for K 2 gives the SOS stability limit that is Re SOS,K 2 = 0.50, which is less than Re SOS,K 1 = 0.85. For 10 Galerkin modes set K 3 no feasible Lyapunov function was found for the SOS optimization problem (2.14), even after we decreased Re to Re E . In other words, Re SOS,K 3 = 0.
To understand why increasing the number of Galerkin modes does not necessarily improve the SOS stability limit Re SOS , we revisit the stability condition for the uncertain system: On the one hand, when more modes are taken into account, the negative κ s becomes larger in magnitude, thus increasing the potential for (3.12) to be satisfied. However, using more Galerkin modes changes unfavourably the bounds of the uncertainties Θ b and Θ c , thus increasing p(a, q 2 ). As a result, the potential for (3.12) to be satisfied is reduced. This means that there might be a finite optimum number of Galerkin modes to be included explicitly into the analysis by the method of Goulart & Chernyshenko [6].

(iii) Conservativeness analysis
Everywhere in this section we consider the K 1 system, which gave the best SOS limits we obtained.
The global stability of the uncertain system implies the global stability of the Navier-Stokes system, but not the vice versa. The double-periodic rotating Couette flow (2.17)-(2.18) is globally stable for Re < Re L (Ω), while the best SOS stability limit we could obtain implies global stability for Re < min{Re L (Ω), Re E + 0.85} only. This SOS stability limit can be conservative for two reasons. First, the limit obtained can be less than the actual stability limit for the uncertain system, for example because SOS approach gives only a sufficient condition for the non-negativity of a polynomial, or because the polynomial degree of the candidate Lyapunov function is not large enough. Second, it can be that the global stability limit for the uncertain system is smaller than the global stability limit for the full Navier-Stokes system. The second possibility turns out to be the case here for those values of Ω when Re SOS < Re L (Ω).
To demonstrate this, note that the obtained limits Re SOS and Re SOS,T are independent of Ω for sufficiently small or large Ω, as shown in figure 1. First, we consider Re SOS,T. It is easy to show analytically that the truncated system has five steady solutions: The equilibria O ± 2 exist and are non-zero if Re > 5 √ 5 ≈ Re E + 5.5235. The equilibria O ± 3 exist and are non-zero if Re > Re L . Hence, the truncated system is not globally stable when Re > min{Re L , 5 √ 5}, implying that Re SOS,T ≤ min{Re L , 5 √ 5}. Combining the analytical stability result and the SOS optimization numerical result, one has min{Re L , Re E + 5.52} ≤ Re SOS,T ≤ min{Re L , 5 √ 5} ≈ min{Re L , Re E + 5.5235}.
It can be seen that the stability limit for the truncated system is attained by the SOS optimization analysis, in the sense that the error is less than δRe.
For the case of the uncertain system we could not obtain an analytic solution. Increasing the degree of the Lyapunov function candidate to six by taking V(a, q 2 ) = P(a, c 6 ) + P(a, c 4 )q 2 + P(a, c 2 )q 4 + αq 6 , where c 2 , c 4 , c 6 and α are decision variables, gave the same SOS stability limit as the 4th-degree Lyapunov function. Taking even higher degree polynomial led to so high computational costs that the calculations had to be abandoned.
Fortunately, we can demonstrate numerically that the SOS stability limit of the uncertain system obtained with the Lyapunov function of 4th degree is very close to the actual global stability limit of the uncertain system. For this reason, increasing the degree of Lyapunov function candidates is not helpful.
The idea is to evaluate the lower bound for such Re that there exist a steady non-zero solution of (2.15)-(2.16). This is similar to what we did for the truncated system, but we use numerical rather than analytic solutions. Naturally, the existence of steady non-zero solutions implies that  the zero solution is not globally stable. Starting from any pre-specified Reynolds number for which there exists a non-zero equilibrium of the uncertain system, we decrease Re gradually. The smallest Re for which the uncertain system still has a non-zero equilibrium will be an upper bound of the global stability limit Re U of the uncertain system. Consider the following nonlinear optimization problem: min Re where the constraints F i are The constraints F 1 and F 2 ensure thatȧ =(q 2 ) = 0 in (2.15)-(2.16), and the constraints F 3 and F 4 are the bound for the uncertain terms in (2.15)- (2.16). In other words, we minimize Re subject to a constraint that there exist a, Θ, Γ , q satisfying a steady version of (2.15)-(2.16) and not coinciding with a = 0, q = 0 solution, the stability of which we are investigating. The variable Γ in the program can be eliminated by combining the constraints F 2 and F 3 as aΘ − κ s q 2 ≤ 0. Since the constraints in (3.13) include equalities and inequalities and are highly nonlinear in Re, instead of solving (3.13) directly we consider the following optimization problem for a given Re: where φ(ω 1 , ω 2 , ) = 1 2 ε + (ω 1 − ω 2 ) 2 + 1 2 (ω 1 + ω 2 ), ω 1 = a T Θ − κ s (Re)q 2 , ω 2 = |Θ| 2 − p(a, q 2 ). In (3.14), the small parameter ε > 0 is introduced to smooth the objective function. Note that φ(ω 1 , ω 2 , ) ≥ φ(ω 1 , ω 2 , 0) = 1 2 |ω 1 − ω 2 | + 1 2 (ω 1 + ω 2 ) = max(ω 1 , ω 2 ).
Hence, negative φ(ω 1 , ω 2 , ) implies ω 1 < 0, ω 2 < 0. As such, all the constraints in (3.13) would be satisfied if φ(ω 1 , ω 2 , ) < 0 in (3.14). Now, the lower bound for Re that leads to the non-global stability of system (2.15)-(2.16) can be obtained by decreasing Re and solving (3.14) repeatedly. This procedure is stopped once the minimum of φ(ω 1 , ω 2 , ) is no longer negative. Let = 0.1 in (3.14). Tables 2 and 3 show the optimization results for Ω = 0 and Ω = 0.1 when Re decreases from Re E + 1.00 to Re E + 0.85. The initial searching point in each trial is always set as the stopping point in the previous trial. The sequential quadratic programming method [24] associated with the function NLPSolve in MAPLE optimization toolbox is used to solve (3.14). From the tables, we can see that all the constraints in (3.13) can be satisfied when Re ≥ Re E + 0.86 in the sense that the residual error |f + Θ| 2 is negligible. This implies that the global stability limit  Overall, the analysis of this section shows that for the flow in question further improvement of the SOS stability limit can be achieved only by improving the uncertainty bounds, since increasing the number of Galerkin modes explicitly taken into account gives no improvement, while the SOS stability limit for the uncertain system is already tight.

Discussion and conclusion
The uncertainty bounds obtained by the methods proposed in [6] and appendix A are tight for each of the components of Θ b and Θ c . However, they are attained at different a and u s . Hence, the bound on |Θ| used to obtain the stability limit is not tight. This provides a potential for improving the results.
For all practical purposes, the question of a global stability of a particular flow can often be answered by physical experiment or numerical calculation, at least up to the existence of unstable stationary points or orbits. Theoretical studies of global stability, however, can provide deeper insight into various aspects of the flow and of the methods used. For example, one could hope to gain such insight from the particular form of the Lyapunov functional. It appears, however, that the method [6] generates Lyapunov functionals of a rather complicated form, which is difficult to interpret. In addition, the obtained Lyapunov functional depends on computational parameters such as the number of the Galerkin modes taken explicitly into account and the degree and the form of the polynomial representing the candidate Lyapunov function.
Reduction of the Navier-Stokes equation to an uncertain system, which forms the basis of the method of [6], turns out to be useful for a number of problems, which might be even of larger interest than the global stability problem. For example, it might prove useful in studies of bounds for long-time averages of the characteristics of turbulent flows and other infinite-dimensional systems with complicated behaviour, and in designing control for such systems [14]. For such studies, the observations of the properties of the method of [6] made in the present work might constitute even a greater interest than the central result of this paper.
It remains to give a summary of the obtained results and observations. A new expression for one of the uncertainty bounds required by the method of [6] was obtained (see appendix A). It allows a considerable reduction of the computational cost as compared to the approach proposed in [6].
A systematic approach to selecting the particular Galerkin modes for the uncertain system was proposed.
It was shown that the selection of Galerkin modes leading to a better stability result for the truncated Galerkin system does not necessarily lead to a better stability result for the uncertain system. It was also shown that increasing the dimension of the uncertain system does not necessarily improve the stability bound obtained for the full system, and that increasing the degree of the candidate polynomial Lyapunov function also does not necessarily improve the bounds. This suggests that further progress in this problem is more likely to be achieved by improving the bounds in the uncertain system than by improving SOS optimization.
For the particular version of the double-periodic rotating Couette flow we considered, we demonstrated that (i) for Ω ∈ (0.2529, 0.7471), the SOS stability limit coincides with the actual global stability limit. (ii) for Ω ∈ (0, 0.2529] ∪ [0.7491, 1), the SOS stability limit Re SOS ≥ Re E + 0.85, where Re E is the energy stability limit. This demonstrates the feasibility of the method proposed in [6].
To the best of our knowledge, this is the first case in which a systematic method applicable in principle to any fluid flow was used successfully to prove global stability of a fluid flow for the value of the Reynolds number greater than that which could be achieved with the energy stability approach. maximum of the spectral radii constrained to x ∈ V, and |Θ ci (u s )| ≤ 2|ρ(D i (x))| ∞ q 2 .
( A 3 ) The bound (A 3) can be shown to be tight by constructing a function u s with compact support similar to a delta function centred around the point x m , where ρ(D i (x)) attains its global maximum. The idea here is similar to the one used in [25], while the difference is that u s is a vector and not a scalar, and that it has to be solenoidal, and orthogonal to all the basis field functions u i , i = 1, . . . , k. Moreover, as |a| → ∞ and q → ∞, the leading term of V cannot be proportional to E t . Indeed, let V = E n t + V 1 , such that V 1 E n t as |a| → ∞. Then, since f(a) = O(|a|), (∂V 1 /∂a)f ∼ V 1 (∂E n t /∂a)f = nE n−1 t (λ 1,1 a 2 1 + λ 1,−1 a 2 2 + λ 1,2 a 2 3 + λ 1,−2 a 2 4 ). Hence, for large enough a 2 and a 1 = a 3 = a 4 = q = 0, (∂V/∂a)f > 0, so that V is not a Lyapunov function.
On the other hand, in view of (B 1), p(a, q 2 ) = O(q 2 (|a| 2 + q 2 )). If the main term of V is not a function of E t , then the order of the right-hand side of (B 2) is greater than the order of its lefthand side, and the inequality (B 2) cannot hold true, so that V cannot be a Lyapunov function in this case either.
Overall, when Re > Re E , there does not exist a polynomial Lyapunov function for the uncertain system (2.15)-(2.16), if only the modes e 1,±1 , e 1,±2 are considered.