2.1. The original energy method

After projecting the perturbed Navier–Stokes equation on an $L^2$ orthogonal and divergence-free basis, the fluid motion can be described by the following truncated ordinary differential equation system

(2.1) \begin{equation} \frac {\mathrm{d}q_i}{\mathrm{d}t} = A_{ij}\,q_j + Q_{ijk}\, q_j\, q_k, \end{equation}

where $q_i(t)$ represents an $n$ -element vector ( $i=1,\ldots, n$ ) describing the perturbation of the base flow as a function of time $t$ . The coefficients $A_{ij}$ and $Q_{ijk}$ are time-independent arrays characterising the behaviour of the perturbed flow, where $i, j$ and $k$ are running variables ranging from $1$ to $n$ in the Einstein summation notation. The state at the origin is stable if the perturbations ( $q_i$ ) tend to zero as $t \to \infty$ . In cases where the perturbation is assumed to be small ( $q_i\propto \epsilon$ ), neglecting the nonlinear (quadratic) terms in the equation allows for linear stability analysis. This involves examining the eigenvalues of the matrix $A_{ij}$ . However, matrix $A_{ij}$ might be non-normal, that is its eigenvectors might be non-orthogonal. Then perturbations might grow transiently reaching the amplitudes at which nonlinear terms cannot be neglected (Schmid Reference Schmid2007; Kerswell Reference Kerswell2018).

An alternative method of stability analysis involves examining the derivative of the perturbation kinetic energy with respect to time. Assuming the kinetic energy of the perturbations is the inner product of the state vector

(2.2) \begin{equation} e = q_i q_i, \end{equation}

its temporal derivative can be easily obtained from (2.1)

(2.3) \begin{align} \frac {\mathrm{d}e}{\mathrm{d}t} &= 2\, A_{ij}\,q_i\,q_j + 2 \, Q_{ijk}\,q_i\, q_j\,q_k . \end{align}

According to the Reynolds–Orr identity (Orr Reference Orr1907; Schmid & Henningson Reference Schmid and Henningson2001) (utilising Gauss divergence theorem), the nonlinear term does not influence the change in kinetic energy if the perturbations are confined by walls, are periodic or decay to zero in the far field, which are reasonable assumptions in most cases

(2.4) \begin{equation} 2\,Q_{ijk}\,q_i\, q_j\,q_k = 0. \end{equation}

From this point, matrices and vectors are denoted by bold letters to enhance readability. The Einstein summation notation is used when a three-dimensional array appears in an expression or the discussion.

The growth rate of the kinetic energy is

(2.5) \begin{equation} \mu _e = \frac {1}{e}\frac {\mathrm{d}e}{\mathrm{d}t} ,\end{equation}

and using (2.3) and (2.4) the following expression can be derived:

(2.6) \begin{equation} \mu _e = \frac {2\, \boldsymbol{q}^T \textbf {A} \boldsymbol{q}}{\boldsymbol{q}^T \boldsymbol{q}}. \end{equation}

The zero solution of (2.1) is Lyapunov stable if $\mu _e\lt 0$ for any $q_i$ . This statement is equivalent to ensuring that the maximum over any possible state is negative

(2.7) \begin{equation} \mu _{{m},e}= \max\limits _{\boldsymbol{q}} \mu _e(\boldsymbol{q})\lt 0. \end{equation}

The numerator in (2.6) can be written as $2 \,\boldsymbol{q}^T \textbf {A}\,\boldsymbol{q} = \boldsymbol{q}^T (\textbf {A} +\textbf {A}^T)\boldsymbol{q}$ . Moreover, expression (2.6) represents the Rayleigh quotient of $\textbf {A} + \textbf {A}^T$ . As $\textbf {A} + \textbf {A}^T$ is symmetric, the largest Rayleigh quotient corresponds to its largest eigenvalue, which is the maximum possible growth rate of kinetic energy

(2.8) \begin{equation} \max\limits _{\boldsymbol{q}} \mu _e(\boldsymbol{q}) = \lambda _{max }\big (\textbf {A} +\textbf {A}^T\big ). \end{equation}

Therefore, the laminar state is Lyapunov stable if

(2.9) \begin{equation} \lambda _{max }\big (\textbf {A} +\textbf {A}^T\big ) \lt 0. \end{equation}

The critical state, which maximises the growth rate of kinetic energy, is the corresponding eigenvector. Unfortunately, this condition is sufficient but often far from necessary in practice. This analysis is referred to as energy method or nonlinear stability analysis since the results are valid for the nonlinear system due to the nonlinear terms not being assumed zero during the derivation but were eliminated by the Reynolds–Orr identity.

In many fluid dynamics applications, the primary concern is not just whether the flow is stable, but the critical limit at which instability arises. Notably, viscosity or the Reynolds number affects only the linear terms. In certain flow configurations, such as Poiseuille or Couette flows, the non-dimensionalised base flow is independent of viscosity, and the matrix $\textbf {A}$ can be decomposed into components that depend on the Reynolds number and those that do not

(2.10) \begin{equation} \textbf {A}(\textit{Re}) = \textbf {A}_U + \frac {1}{\textit{Re}} \textbf {A}_R. \end{equation}

Considering that the Laplacian term can only dissipate kinetic energy, $\textbf {A}_R$ is a negative definite matrix. The smallest Reynolds number, for which $\mu _{{m},e} = 0$ and $e\neq 0$ , is equal to the smallest Reynolds number, for which $\mu _e$ can be 0, therefore it is equal to the smallest Reynolds number for which $\mu _e=0$ . By substituting (2.10) into (2.6), setting the expression to zero, and subsequently expressing $\textit{Re}$ and calculating its minimum through variation, we arrive at the corresponding Euler–Lagrange equation

(2.11) \begin{equation} {\big(\textbf {A}_R + \textbf {A}^T_R\big )\boldsymbol{q} = \tilde {\textit{Re}}\big (-\textbf {A}_U -\textbf {A}^T_U\big )\boldsymbol{q}.} \end{equation}

This equation represents a general eigenvalue problem where the eigenvalue is the Reynolds number. The smallest eigenvalue, typically denoted as $\textit{Re}_{{E}}$ , is referred to as the global stability limit. If $\textit{Re}\lt \textit{Re}_{{E}}$ , then $\mu _{{m},e}\lt 0$ , signifying unconditional stability.

2.2. The generalised kinetic energy method

The classical energy method often proves to be highly conservative, predicting Reynolds number limits below experimental observations. This issue arises because, at high Reynolds numbers, the $\textbf {A}$ matrix becomes non-normal. In such cases, the eigenvectors are non-orthogonal, and the energy of solutions initialised close to a linearly stable laminar state can grow significantly (Schmid Reference Schmid2007), although they do not necessarily transition to another state. Nerli et al. (Reference Nerli, Camarri and Salvetti2007) defined $\textbf {P}$ in a perturbed form, which is less general than the transformation matrix formulation presented in this work.

A possible way to define a generalised kinetic energy, $h$ , is through the linear transformation of the variables by an invertible $\textbf {S}$ matrix

(2.12) \begin{equation} \boldsymbol{q} =\textbf {S}\,\boldsymbol{r}, \end{equation}

and

(2.13) \begin{equation} h= \boldsymbol{r}^T\boldsymbol{r}. \end{equation}

The $h$ function is the most general one among quadratic energy functions. Olivier Dauchot & Paul Manneville (Reference Dauchot and Manneville1997) defined similarly the ‘exotic’ energy in their paper, where the transformation matrix is the eigenvectors of the linear part of the dynamical system. However, it is not optimal in general, a point that will be discussed later. In another form, $h$ can be expressed as

(2.14) \begin{equation} h = \boldsymbol{q}^T \textbf {P}\, \boldsymbol{q} ,\end{equation}

where

(2.15) \begin{equation} \textbf {P} = \textbf {S}^{-T}\textbf {S}^{-1}, \end{equation}

and the matrix $\textbf {P}$ is clearly positive definite.

Among quadratic energy functions, there is no way to improve the global stability results obtained with kinetic energy (Goulart & Chernyshenko Reference Goulart and Chernyshenko2012), therefore $h$ is constructed to prove the conditional stability of the laminar state. A conditional Lyapunov function must meet the following requirements:

(2.16) \begin{equation} h(\boldsymbol{q}) \gt 0 \; \mathrm{if} \; \boldsymbol{q}\neq \boldsymbol{0} \;\mathrm{and} \; h(\boldsymbol{q}) = 0 \;\mathrm{if} \; \boldsymbol{q} = \boldsymbol{0} ,\end{equation}

and

(2.17) \begin{equation} {\frac {\partial h}{\partial q_i}\frac {\mathrm{d}q_i}{\mathrm{d}t} \lt 0} \; \mathrm{if} \; h(\boldsymbol{q}) \lt \gamma _a^2 .\end{equation}

The region where $h(\boldsymbol{q}) \lt \gamma _a^2$ represents the provable region of attraction (ROA) of the origin, although the true ROA can be larger in general. The definition of $h$ in (2.14) automatically satisfies the first condition (2.16). Our objective is to maximise the size of the ROA satisfying equation (2.17) by constructing an optimal $h$ function. Here, the smallest radius, $\sqrt {e}_{{min}}$ in the original state space is used to describe the region since it makes the results of minimal seed calculations comparable.

There are three essential parts of the investigation. The first is verifying that a given $h$ satisfies equation (2.17). The second involves maximising the ROA and $\gamma _a$ for a given $h$ , where the critical value $\gamma _{{crit}}$ determines $\sqrt {e}_{{min}}$ . The third focuses on finding the optimal $h$ that maximises this region.

2.2.1. Validating the generalised kinetic energy to be a conditional Lyapunov function

For the validation of (2.17), it is convenient to transform the differential equation (2.1) for the new variables

(2.18) \begin{equation} \frac {\mathrm{d} r_i}{\mathrm{d}t} = S_{ij}^{-1} A_{jk}\,S_{kl}\,r_l + S_{ij}^{-1} Q_{jkl}\, S_{km}\, r_m \, S_{lo}\,r_o. \end{equation}

To facilitate this transformation, let us define

(2.19) \begin{align}&\quad \tilde {A}_{ij} = S_{il}^{-1} A_{lk}\,S_{kj}, \end{align}

(2.20) \begin{align}& \tilde {Q}_{ijk} = S_{im}^{-1} Q_{mol}\, S_{oj}\, S_{lk} .\end{align}

These transformations result in a similar ordinary differential to (2.1), if ${A}_{ij}$ and ${Q}_{ijk}$ are replaced by $\tilde {A}_{ij}$ and $\tilde {Q}_{ijk}$ , respectively.

To analyse the system at different perturbation levels, the perturbation magnitude $\gamma$ is introduced by decomposing the transformed state vector into its magnitude, $\gamma = \sqrt {r_i r_i}$ , and a unit vector

(2.21) \begin{equation} r_i = \gamma \tilde {r}_i. \end{equation}

Next, let us define the growth rate of the generalised kinetic energy

(2.22) \begin{equation} \mu _h = \frac {1}{h}\frac {\mathrm{d}h}{\mathrm{d}t}, \end{equation}

which can be calculated as

(2.23) \begin{equation} \mu _h = \frac {2\, \tilde {A}_{ij}\,r_i\,r_j + 2 \, \tilde {Q}_{ijk}\,r_i\, r_j\,r_k}{r_l r_l} = 2\, \tilde {A}_{ij}\,\tilde {r}_i\,\tilde {r}_j + 2\, \gamma \, \tilde {Q}_{ijk}\,\tilde {r}_i\, \tilde {r}_j\,\tilde {r}_k. \end{equation}

The main difference lies in the quadratic term ( $\tilde {Q}_{ijk}$ ) contributing to the growth rate of the generalised kinetic energy ( $h$ ), unlike in the case of the original kinetic energy ( $e$ ). Due to the presence of this term, only conditional but not unconditional stability can be established, and it can be utilised to calculate the threshold amplitude ( $\gamma _{{crit}}$ ).

Let us define the possible maximum growth rate at a given level of perturbation as

(2.24) \begin{equation} \mu _{{max},h} (\textbf {S},\gamma ) = \max\limits _{\tilde {\boldsymbol{r}}} \mu _h(\tilde {\boldsymbol{r}}, \textbf {S},\gamma ). \end{equation}

The generalised kinetic energy, $h$ , should first be investigated at zero perturbation level, $\gamma =0$ . There, similarly to the formula (2.8), the largest growth rate can be determined as

(2.25) \begin{equation} \mu _{h,{lin,max}} =\mu _{{max},h}(\textbf {S},\gamma =0) = \lambda _{{max}}\big (\tilde {\textbf {A}} +\tilde {\textbf {A}}^T\big ), \end{equation}

where $\lambda _{{max}}$ denotes the largest eigenvalue of a matrix. If $\lambda _{{max}} (\tilde {\textbf {A}} +\tilde {\textbf {A}}^T )\geqslant 0$ , then $h$ cannot demonstrate the existence of a ROA, similar to how the original kinetic energy, $e$ , fails to do so if $\textit{Re} \gt \textit{Re}_E$ . However, it is well known that if the equilibrium point is linearly stable, a suitable $h$ function must exist, ensuring the maximal generalised kinetic energy growth rate of the linearised system is negative. For instance, the eigenvectors of $\textbf {A}$ can be used to construct such a function.

If $\lambda _{{max}} (\tilde {\textbf {A}} +\tilde {\textbf {A}}^T )\lt 0$ , $h$ can prove the existence of a ROA. The next step is to increase $\gamma$ and expand the ROA until the point where the maximal growth becomes zero, ( $\gamma = \gamma _{{crit}}$ ). Since the growth rate must not exceed the maximum growth rate, and the maximum remains negative up to this perturbation level, the generalised kinetic energy or h decreases monotonically, if $\gamma \lt \gamma _{{crit}}$ . Therefore, the second requirement (2.17) of the conditional Lyapunov function is fulfilled.

Finding the global maximum of the growth rate at a given perturbation level among the various states ( $\boldsymbol{r}$ ) is essential for proving that $h$ is a conditional Lyapunov function. This will be achieved using robust optimisation techniques, initialised from multiple seed points to increase the likelihood of accurately identifying the global maximum. Unfortunately, the optimisation procedure was not discussed in the paper by Nerli et al. (Reference Nerli, Camarri and Salvetti2007), and thus this aspect of the present study fills an important gap in validating the method.

2.2.2. Calculation of the maximal growth of the generalised kinetic energy

Proving that a specific $\tilde {\boldsymbol{r}}$ maximises the expression in (2.23), while also adhering to the constraint that ( $\|\tilde {\boldsymbol{r}}\|=1$ ), presents a formidable challenge.

Analytical form of (2.23) helps to derive the gradient and the Hessian, accelerating the optimisation. Among the various methods explored are sequential quadratic programming (SQP), the active set algorithm and the interior-point algorithm (Nocedal & Wright Reference Nocedal and Wright2006), all of which can be efficiently employed using MATLAB’s fmincon function.

First, various dynamical models with increasing degrees of freedom, $n$ , are generated. For $n=2$ , the TTRD’ model is utilised, for $n=4$ the W95A (defined in § 3.1) model is employed, and for $n = {12, 18, 24, 30, 45, 60, 90, 150}$ , a range of Poiseuille flow models (PFMs) are considered. The specifics of these models will be detailed later.

Figure 1. Computation time for calculating the maximal growth rate (a) and the relative frequency of solutions achieving the actual maximum (b) as functions of the dynamical system’s degrees of freedom. The shaded areas represent the standard deviation from the mean values.

The optimisation process was conducted for each dynamical model using 20 randomly generated transformation matrices, across a spectrum of relatively high perturbation levels $\gamma = { 0.1, 0.2, 0.5, 1, 2, 5, 10}$ . For each model, 100 randomly generated seed points – unitary $\tilde {\boldsymbol{r}}$ state vectors – were created. From each seed point, optimisation was executed using three distinct algorithms. The global maximum for a given transformation matrix at a specified perturbation level was determined as the maximum across all optimisation techniques and seed points. The methods were compared based on the convergence rate, defined as the proportion of solutions that converged to the global maximum from differently initialised optimisations.

The convergence rate as a function of the number of degrees of freedom is depicted in figure 1(b), with the thick curve representing the mean value and the dashed line and shadow indicating the standard deviation of the convergence rate across various perturbation levels and transformation matrices. For systems with few degrees of freedom, the convergence rate was approximately 70 %, decreasing to 20 % as $n$ increases. For small systems ( $n\lt 20$ ), no significant difference was observed among the optimisation algorithms; however, for larger systems ( $n\gt 100$ ), the interior-point algorithm clearly outperformed the others. Its mean convergence rate was higher, and its standard deviation was smaller. Although not visualised here, the convergence rate as a function of perturbation level was evaluated. For SQP and active set methods, the convergence rate decreased with increasing perturbation level, whereas for the interior-point algorithm, it remained nearly constant. Based on the convergence rate, the interior-point algorithm is superior, especially for large dynamical systems.

Computational time comparison, shown in figure 1(a), reveals that, for smaller systems ( $n\lt 12$ ), the SQP method is the fastest. There exists a narrow range ( $12\lt n\lt 20$ ) where the active set method is optimal. For larger systems, the interior-point method is significantly faster than the other two methods, primarily because it utilises the analytical Hessian matrix, unlike the other methods.

In summary, the SQP method is recommended for small systems, while the interior-point method is better suited for larger ones. To determine the maximum growth rate for a given transformation matrix and perturbation level, optimisation is repeated from random initial points until at least five solutions converge to the same maximum value, reducing the risk of converging to a local maximum.

Since finding the global maximum is critical, the interior-point method’s convergence rate was analysed for dynamical models with up to 330 degrees of freedom. At this scale, the convergence rate drops to 5.81%. Assuming a sufficiently large sample, this rate can approximate the probability of locating the global maximum with a random initial guess. With 231 random initial guesses, the probability of missing the global maximum decreases to 1 in $10^6$ , making failures highly unlikely. The applied rule requiring at least five consistent results often necessitates even more initial guesses.

2.2.3. Calculation of the critical perturbation level

Figure 2. The maximal growth rate of the generalised kinetic energy ( $\mu _{{max}, h}$ ) as the function perturbation magnitude $\gamma$ in the case of the optimally transformed TTRD’ model at $\textit{Re} = 5$ . The red curve represents the one tenth of the growth rate of the original kinetic energy ( $\mu _e = 0.6$ ), which is independent of the perturbation level.

As the $\mu _{{max},h}$ function is created by a proper optimisation procedure, the next step is to determine the largest possible critical perturbation level below which this growth is negative

(2.26) \begin{equation} \mu _{{max},h} (\textbf {S},\gamma _{{crit}}) = 0. \end{equation}

The corresponding unitary state vector, defined by

(2.27) \begin{equation} \mu _{{max},h} (\textbf {S},\gamma _{{crit}}) = \mu _h(\tilde {\boldsymbol{r}}_{{crit}}, \textbf {S},\gamma _{{crit}}), \end{equation}

can be utilised to obtain the critical state: $\boldsymbol{r}_{{crit}} = \gamma _{{crit}} \tilde {\boldsymbol{r}}_{{crit}}$ .

For illustration, $\mu _{{max},h}$ is plotted in figure 2 for one particular problem. At low $\gamma$ values, the linear part of the dynamical system dominates, where the maximum growth rate is almost constant and equal to the Rayleigh quotient of the $\tilde {A}_{ij}+\tilde {A}_{ji}$ matrix, see equation (2.25). For higher $\gamma$ values, the nonlinearity of the system influences the maximal growth, which tends towards a straight line. The slope of this line corresponds to the maximum of $\lbrace 2\, \, \tilde {Q}_{ijk}\,\tilde {r}_i\, \tilde {r}_j\,\tilde {r}_k\rbrace$ among possible $\tilde {r}_i$ states. Figure 2 illustrates the maximal growth rate of the original kinetic energy, highlighting that the laminar state is unstable according to the regular energy method. Moreover, it demonstrates that the original kinetic energy remains constant and is unaffected by the perturbation amplitude.

Determining the critical perturbation level is relatively simpler compared with the other parts of the method. It involves finding the roots of equation (2.26). The gradient of this expression can be derived easily as

(2.28) \begin{equation} \frac {\partial \mu _{{max},h}}{\partial \gamma } = 2\, \tilde {Q}_{ijk}\,\tilde {r}_{{max},i}\, \tilde {r}_{{max},j}\,\tilde {r}_{{max},k} \end{equation}

where $\tilde {\boldsymbol{r}}_{{max}}$ is the vector that maximises the growth rate at a given perturbation level. By leveraging the analytical gradient, the Newton method, renowned for its second-order convergence rate, is utilised effectively. Typically, after a few iterations (5–6), the solution’s residual falls below the tolerated error threshold of $10^{-10}$ .

Figure 3. The phase space trajectories of the TTRD’ model at $\textit{Re} = 5$ are shown for the original state variables (a) and for the optimally transformed variables (b). Green trajectories converge towards the origin, while red trajectories tend to another equilibrium point, which is not shown. The black curve represents the provable boundary of the ROA utilising the optimal $h$ function. The dashed, blue curve shows the boundary of the true ROA. The red vector is the critical perturbation (2.26), where the growth rate of the generalised kinetic energy was zero at the critical perturbation level (2.27). The yellow vector illustrates the smallest perturbation (2.30) in the original state space (a) whose length is equal to the critical perturbation in the optimally transformed state space (b).

The investigated region can be envisioned as a multidimensional hypersphere in the $\boldsymbol{r}$ state space around the origin. The radius of this sphere is $\gamma$ . If the radius is smaller than a critical value $\gamma _{{crit}}$ , then $\mu _h\lt 0$ , indicating that the generalised kinetic energy ( $h$ ) is decreasing, and the trajectories move inward, ultimately converging to the origin. At the critical radius, a trajectory becomes tangential to the sphere and may not reach the origin. The hypersphere with radius $\gamma _{{crit}}$ represents the boundary of the provable ROA. However, states outside this sphere can belong to the basin of attraction of either the origin or another attractor. In the case of the two-dimensional problem, the hyperspherical ROA reduces to a circle and is illustrated in figure 3(b) in the case of an optimal transformation matrix and $h$ function.

It is more interpretable to transform the ROA back to the original state space. The linear transformation (scaling and rotating) of the hypersphere results in a hyperellipsoid in the original state space $\boldsymbol{q}$ . This hyperellipsoid defines the boundary of the provable ROA of the origin. A standard quantity describing ROA is its smallest radius, $\sqrt {e_{{min}}}$ .

The value $e_{{min}}$ can be calculated using equations (2.2) and (2.21) as follows:

(2.29) \begin{equation} e_{{min}}(\textbf {S}) = \gamma _{{crit}}^2 (\textbf {S}) \min _{\tilde {\boldsymbol{r}}} \lbrace \tilde {\boldsymbol{r}}^T \textbf {S}^T \textbf {S} \tilde {\boldsymbol{r}} \rbrace . \end{equation}

The argument of the minimum function is the Rayleigh quotient of $\textbf {S}^T \textbf {S}$ , and the minimum value corresponds to the smallest eigenvalue of $\textbf {S}^T \textbf {S}$ , since $\textbf {S}^T \textbf {S}$ is a symmetric matrix

(2.30) \begin{equation} e_{{min}}(\textbf {S}) = \gamma _{{crit}}^2 (\textbf {S}) \;\lambda _{{min}} \big ( {\textbf {S}^T \textbf {S}}\big ) . \end{equation}

The corresponding unitary eigenvector $\tilde {\boldsymbol{r}}_{{min}}$ can be utilised to get the two locations $\boldsymbol{q}_{{min}} = \pm \gamma _{{crit}}\,\textbf {S}\, \tilde {\boldsymbol{r}}_{{min}} = \pm \textbf {S}\, \boldsymbol{r}_{{min}}$ , where the inner hypersphere touches the hyperellipsoid, as illustrated in figure 3(a). This inner hypersphere is a subset of the provable ROA, where $e\lt e_{{min}}$ .

The value of $e_{{min}}$ depends on the definition of $h$ through the transformation matrix $\textbf {S}$ . The aim of the method is to maximise this minimum energy, $e_{{min}}(\textbf {S})$ over all matrices, $\textbf {S}$ . The possible methods for maximisation will be discussed in § 2.3. It is important to note that the value of $\gamma$ is only meaningful for a specific transformation or generalised kinetic energy. Comparisons of $\gamma$ or $\gamma _{{crit}}$ between different $\textbf {S}$ matrices or $h$ functions are not valid.

It should also be mentioned that although the kinetic energy ( $e$ ) can grow significantly inside the provable ROA region – a phenomenon usually called the non-modal growth – stability is guaranteed there due to the exponential decay of $h$ there.

2.3. Optimising the generalised kinetic energy

The last key question is how to determine the optimal $h$ function, represented by the transformation matrix $\textbf {S}_{{opt}}$ , to maximise the size of the provable ROA, $e_{{min}}$ . First, a method utilising the singular value decomposition (SVD) of the eigenvectors is presented. This method cannot provide the largest $e_{{min}}$ , but it is extremely fast compared with other methods. As a second approach, all the elements of the transformation matrix are treated as decision variables, leading to a significantly larger ROA at the cost of higher computational need.

2.3.1. The GKE-SV method utilising the singular value decomposition of the eigenvector matrix

A plausible starting point is to use the eigenvectors of $\textbf {A}$ , $\boldsymbol{\Psi }_c$ , as the transformation matrix, since it diagonalises the $\textbf {A}$ matrix. This transformation redefines the state variables as coefficients of the eigenmodes, addressing non-normality by making $\tilde {\textbf {A}}$ diagonal and its eigenvectors orthogonal. Consequently, the GKE, expressed as the sum of squared eigenmode coefficients, ensures stability at infinitesimally low perturbation levels, as noted by Olivier Dauchot & Paul Manneville (Reference Dauchot and Manneville1997).

Diagonalising $\textbf {A}$ with $\boldsymbol{\Psi }_c$ results in the maximum linear growth rate of

(2.31) \begin{equation} \mu _{h,\mathit{lin,max}} (\boldsymbol{\Psi }_c) = \lambda _{{max}}\big (\tilde {\textbf {A}} +\tilde {\textbf {A}}^T\big ) = 2 {\textit{Re}} \big (\lambda _{{max}} (\textbf {A})\big ), \end{equation}

which must be negative if the equilibrium point is linearly stable and there exists a certain provable ROA. However, it will be shown that $\boldsymbol{\Psi }_c$ used to define $h$ is not optimal, but systematic modifications can yield an improved transformation matrix and $h$ .

As a first step, the eigenvector matrix is transformed into a real-valued matrix, denoted as $\boldsymbol{\Psi }$ . This involves reinterpreting all eigenvector pairs associated with the same complex-conjugate eigenvalue pairs, alternately focusing on the real and imaginary components, and subsequently normalising these to unit vectors. Although this new transformation matrix does not diagonalise the $\textbf {A}$ matrix, the formulation of generalised energy remains analogous and equation (2.31) still holds.

Next, the SVD of $\boldsymbol{\Psi }$ matrix is performed

(2.32) \begin{equation} \boldsymbol{\Psi } = \textbf {U}_\psi \,\textbf {S}_\psi \,\textbf {V}_\psi . \end{equation}

It can be demonstrated that $\textbf {V}_\psi$ does not affect the definition of $h$ in equation (2.15) and $\textbf {V}_\psi$ can therefore be omitted.

For simplicity, the singular values can be normalised by their maximum without sacrificing generality. Furthermore, let the singular values be denoted by $s_i$ , where $s_1$ is the largest and $s_n$ is the smallest. Due to normalisation $s_1 = 1$ and $0\lt s_i\leqslant 1$ . It is well known that, as the Reynolds number increases, the non-normality of A or the condition number of the eigenvector matrix increases. Therefore, the smallest singular value of $\boldsymbol{\Psi }$ , $s_n$ , tends toward zero, leading to two significant issues. One problem is that a decreasing $s_n$ makes the eigenvector matrix ill conditioned and difficult to invert accurately. Additionally, this choice likely reduces the size of ROA, $e_{{min}}$ , since

(2.33) \begin{equation} e_{{min}}\,\propto \, \lambda _{{min}}\big (\textbf {S}^T\textbf {S}\big ) = s_n^2 ,\end{equation}

according to equation (2.30). Visually, utilising $\boldsymbol{\Psi }$ as the transformation matrix results in an attraction region shaped like a narrow superellipsoid with its smallest axis shrinking to zero. In this context, the concept of ‘exotic’ energy as defined by Olivier Dauchot & Paul Manneville (Reference Dauchot and Manneville1997), while effective for two-dimensional systems at lower Reynolds numbers, is presumably inefficient for larger systems or those at higher Reynolds numbers, where the ratio of the smallest and largest singular values of $\boldsymbol{\Psi }$ becomes large.

A straightforward solution to the problem is proposed by modifying the transformation matrix as

(2.34) \begin{equation} \textbf {S} = \textbf {U}_\psi \,\tilde {\textbf {S}}_\psi, \end{equation}

where $\textbf {U}_\psi$ comes from the SVD decomposition of $\boldsymbol{\Psi }$ and $\tilde {\textbf {S}}$ is the modified matrix of the singular values ( $\textbf {S}_\psi$ ), for which the singular values smaller than a certain threshold, $s_{{min}}$ , are simply set to $s_{{min}}$

(2.35) \begin{equation} \tilde {s}_i = \begin{cases} s_i \,\,\; &\mathrm{if} \, s_i\gt s_{{min}}, \\ s_{{min}} \, &\mathrm{if} \, s_i\leqslant s_{{min}} .\end{cases} \end{equation}

The threshold parameter, $s_{{min}}$ , can be set between $s_n$ and 1. If $s_{{min}} = s_n$ , the transformation matrix is equivalent to the eigenvector matrix in terms of the $h$ definition. Conversely, if $s_{{min}}$ is set to 1, the transformation has no effect, making the GKE equivalent to $e$ . In this case, if $\textit{Re} \gt \textit{Re}_E$ , no provable ROA exists according to $h$ .

By varying $s_{{min}}$ , between the two extrema, $\mu _{h,{lin, max}}$ increases continuously and changes sign. This behaviour is illustrated in figure 4, which shows the maximal growth rate of the linearised system rate after transformation as a function of $s_{{min}}$ for the $n=210$ PFM at $\textit{Re} = 2000$ . Interestingly, for various systems, $s_{{min}}$ can be increased by orders of magnitude above the smallest original singular values without significantly affecting $\mu _{h,{lin,max}}$ . In the presented case, $s_{{min}}$ is increased by two orders of magnitude without altering $\mu _{h,{lin,max}}$ . This demonstrates that the non-normality of $\textbf {A}$ can be effectively managed without relying on an ill-conditioned transformation matrix. Moreover, increasing $s_{{min}}$ significantly enlarges the provable ROA, as $\sqrt {e_{{min}}}$ is proportional to the smallest singular value of the transformation matrix.

Figure 4. The maximal growth rate of the GKE ( $\mu _{h, {lin,max}}$ ) for the linearised system, plotted against the smallest singular value of the modified eigenvector matrix at $\textit{Re} = 2000$ . The system in question is the 210-dimensional PFM. The red curve depicts the maximal growth rate of the original kinetic energy, while the yellow cross marks the critical $s_{{min}}$ value at which the growth rate becomes zero.

During the optimisation of the transformation matrix, the optimal $s_{{min}}$ is determined using a bisection algorithm. The algorithm increases $s_n$ with a relative accuracy of $10^{-3}$ until $\mu _{h,{lin,max}}$ no longer exceeds a relative change of $10^{-5}$ from its minimum possible value.

The optimisation of $\textbf {S}$ such a way offers both benefits and limitations due to its disregard for nonlinear terms. This exclusion streamlines the optimisation, significantly speeding up the process as it necessitates only the eigenvalue analysis of relatively small matrices. After optimising the transformation matrix, a solitary evaluation is sufficient for computing the ROA. However, this method may not achieve optimality, given that it neglects the impact of nonlinear terms during the optimisation of the transformation matrix. This approach has been named the ‘GKE-SV’ algorithm in the context of this study.

2.4. Comparison of the generalised kinetic energy method and sum-of-squares programs

This section introduces a comparison between the proposed procedure for constructing Lyapunov functions and the current state-of-the-art methodologies, particularly those employing SOS programs. A brief description of the applied SOS methods is provided in Appendices A.2.1 and A.2.2.

The first SOS program, SOS1 (Jarvis-Wloszek Reference Jarvis-Wloszek2003; Tan & Packard Reference Tan and Packard2008; Topcu et al. Reference Topcu, Packard and Seiler2008; Khodadadi et al. Reference Khodadadi, Samadi and Khaloozadeh2014; Kalur et al. Reference Kalur, Seiler and Hemati2021), consists of three SOS constraints, one ensures the positive definiteness of the Lyapunov function, the second maximises the size of the ROA while the third mandates that the Lyapunov function decreases along the solution trajectories. The maximal ROA can be obtained by an iterative maximisation of the parameters of the constraints.

The second SOS program, SOS2 (Liu & Gayme Reference Liu and Gayme2020), is similar to the SOS1 program but the second constraint is replaced by an eigenvalue problem, significantly reducing the size of the problem and the necessary number of iterations. Therefore, SOS2 is usually faster but less accurate compared with the SOS1 method.

There is a convex formulation, SOS3, for inner approximating a basin of attraction introduced by Henrion & Korda (Reference Henrion and Korda2014) which avoids the necessary iteration steps of the SOS1 and SOS2 programs. However, in SOS3, time appears as an explicit variable in the polynomials, unlike in SOS1 and SOS2. Additionally, SOS3 consist of four SOS constraints, resulting in significantly larger optimisation problem. Furthermore, it requires the explicit specification of a time horizon and a target domain. These parameters play a critical role in shaping the results. This formulation is excluded from the comparison.

The two SOS algorithms and the two GKE algorithms were evaluated against each other using simple models, focusing on computational time and the size of the ROA as key metrics. For the conversion of SOS programs into a semi-definite program (SDP), two tools, SOSTOOLS (Papachristodoulou et al. Reference Papachristodoulou, Anderson, Valmorbida, Prajna, Seiler and Parrilo2013) and Yalmip (Löfberg Reference Löfberg2004), were employed. However, no significant difference in computational time was observed between them. Consequently, the results presented here utilise the conversion code provided by SOSTOOLS, with the SDP problems solved using MOSEK.

The ROA was calculated for five dynamical systems: TTRD’, W95A, W95B, BT and PFM $_{n=12}$ . These systems have degrees of freedom ( $n$ ) of 2, 4, 4, 4 and 12 respectively. Here, PFM $_{n=12}$ denotes the truncated Galerkin model of Poiseuille flow, while the other models are low-order representations of subcritical transition. These calculations were performed on a computer equipped with an Intel i7-7700 processor running at 3.6 GHz and 32 GB of memory with a speed of 2400 MHz.

The comparative analysis in figure 5 illustrates the performance differences between the GKE and SOS algorithms in determining the ROA for dynamical systems. The GKE-G algorithm identified a ROA radius identical to that of the SOS1 method for low-dimensional models ( $n \leqslant 4$ ) and slightly smaller for the $n = 12$ model, but required more time to do so. The small difference is attributed to the optimisation of the $h$ function in the GKE-G method. When the optimal $h$ definition obtained by the SOS1 algorithm was used to calculate the provable ROA using the GKE-G method, the radius matched that of the SOS1 method. The SOS2 algorithm yielded a slightly smaller ROA, whereas the GKE-SV algorithm predicted significantly smaller regions but demanded substantially less computational resources.

Figure 5. Threshold amplitude (a) across five dynamical systems (TTRD’, W95A, W95B, BT, PFM $_{n=12}$ ) at a Reynolds number of $\textit{Re} = 5\textit{Re}_E$ , as determined by four distinct methods, alongside the computational time (b).

For smaller systems, the SOS1 algorithm is recommended due to its precision. The increased computational demand typically does not pose a significant concern for these cases. However, it has been discovered that the SOS1 method becomes impractical for systems with more than 15 degrees of freedom, as the resultant SDP becomes too large for the optimisation program, MOSEK, to solve. Similar limitations were observed with the SOS2 algorithm when the degrees of freedom exceeded 21.

A core issue identified is that, even when restricting the analysis to quadratic Lyapunov functions, the left-hand side of equations (A18–A19) in SOS1 and (A22) in SOS2 manifest as fourth-order polynomials. To ascertain their positive definiteness, the necessary number of monomials is proportional to $n^2$ , implying that the optimisation matrix size – containing the decision variables – scales roughly with $n^4$ . Consequently, the number of decision variables surges as the degrees of freedom ( $n$ ) increase. The SOS2 algorithm accommodates slightly larger systems due to involving only one fourth-order polynomial, in contrast to the two present in SOS1.

Conversely, GKE methods exhibit scalability to systems an order of magnitude larger, as will be demonstrated. Notably, the GKE-G method provides an approximation of the ROA comparable in accuracy to the widely utilised SOS1 method, offering a viable alternative for larger system analysis. The GKE method exemplifies a less conservative approach with better scalability.

3.1. Trefethen’s and Waleffe’s models

Well-known low-order dynamical systems that mimic the subcritical transition of shear flows are investigated in this study. Specifically, the two-dimensional TTRD’ model of Baggett & Trefethen (Reference Baggett and Trefethen1997) and the four-dimensional W95 model of Waleffe (Reference Waleffe1995) are examined. In these models, the linearised part (A) is non-normal, and the condition number of its eigenvector matrix grows with increasing Reynolds number. Meanwhile, the nonlinear part does not affect the growth rate of kinetic energy, as the corresponding matrix remains asymmetric.

The TTRD’ model is represented by the following equation:

(3.1) \begin{align} \frac {\mathrm{d}\boldsymbol{q}}{\mathrm{d}t}= \begin{bmatrix} -\dfrac {1}{\textit{Re}} & 1 \\[8pt] 0 &-\dfrac {1}{\textit{Re}} \\ \end{bmatrix} \boldsymbol{q}+ \begin{bmatrix} 0 & -q_1 \\[3pt] q_1 &0 \end{bmatrix}\boldsymbol{q}, \end{align}

therefore the non-zero elements of $Q_{ijk}$ are

(3.2) \begin{align}& Q_{1\,1\,2} = -1 , \\[-12pt] \nonumber \end{align}

(3.3) \begin{align}&\,\, Q_{2\,1\,1} = 1. \end{align}

The optimal $h$ functions and the corresponding transformation matrices for TTRD’, that maximise $e_{{min}}$ , are calculated at the Reynolds number between 5 and 100 with the step size 2.5 by the proposed GKE-G method. The result is visualised in figure 3 at $\textit{Re}=5$ , where green trajectories converge to the origin and red trajectories to another state. These trajectories are shown in the original variables (figure 3 a) and the transformed variables (figure 3 b). The calculated ROA appears as an ellipse in the original state space, precisely touching the boundary of the true ROA.

Next, the calculated threshold amplitude as the function of Reynolds number are plotted in figure 6 and compared with the findings of Liu & Gayme (Reference Liu and Gayme2020). The cited authors used the quadratic constraint method, which has been proven to be computationally efficient. They treated the nonlinear term as a forcing with an approximated upper bound.

The calculated threshold amplitude ( $\sqrt {e_{{min}}}$ ) decays as a function of the Reynolds number following a power law. The exponents are nearly identical: −3.005 in this study and −3.07 in the work of Liu & Gayme (Reference Liu and Gayme2020). However, the method presented here predicts a ROA with a radius roughly three times larger, indicating a energy level approximately one magnitude higher. This substantial difference arises from their approximation of the nonlinear term, while GKE-G calculation takes into account the exact terms, providing a more precise representation of the system’s behaviour, similarly to SOS methods.

In the next step, the accuracy of the ROA is investigated by solving the ordinary differential equation initialised close to the outside of the calculated ROA. The solutions are initialised from slightly increased threshold state vectors $\boldsymbol{q}_{0} = c_{u} \boldsymbol{q}_{{min}}$ approximating the minimal seed states and computed using the Matlab ode45 Runge–Kutta method.

The average value of $c_{u}$ for such solutions is found to be 1.02, indicating that the proposed method is highly accurate for this small model. The energy level of these approximated minimal seeds are shown in figure 6 but they are almost on the curve representing the radius of ROA.

Figure 6. The size of ROA as the function of Reynolds number in the case of TTRD’ model. The fitted curve from Liu & Gayme (Reference Liu and Gayme2020) using the QC method $(0.912\,\textit{Re}^{-3.07})$ is shown alongside. The best-fitting curve of GKE-G results is $\sqrt {e_{{min}}}\approx 2.228 \,\textit{Re}^{-3.005}$ . The red crosses represent the energy of the approximated minimal seeds. The vertical red line signifies the unconditional stability limit $\textit{Re}_E = 2$ .

In the next step, the GKE-G method is applied to the low-order model of subcritical transition proposed by Waleffe (Reference Waleffe1995). The dynamical system is represented as follows:

(3.4) \begin{equation} \frac {\mathrm{d}\boldsymbol{q}}{\mathrm{d}t}=\frac {1}{\textit{Re}} \begin{bmatrix} -\lambda _w & \textit{Re} & 0 & 0\\[3pt] 0 &-\mu _w & 0 & 0\\[3pt] 0 & 0 &-\nu _w & 0 \\[3pt] 0 & 0 & 0 &-\sigma _w \\ \end{bmatrix} \boldsymbol{q}+ \begin{bmatrix} -\gamma _w q_3^2 + q_2 q_4 \\[3pt] \delta _w q_3^2 \\[3pt] \gamma _w q_3 q_1-\delta _w q_3 q_2 \\[3pt] -q_4 q_2 \\ \end{bmatrix}. \end{equation}

The parameters $\lambda _w, \mu _w, \nu _w, \sigma _w$ represent the decay rates due to viscosity, while $\gamma _w, \delta _w$ describe the nonlinear interaction between rolls ( $q_2$ ) and streaks ( $q_1$ ).

In this study, three different parameter sets are investigated. The first set is characterised by $\lambda _w = \mu _w= \sigma _w = 10,\, \nu _w = 15,\, \delta _w = 1,\, \gamma _w = 0.1$ , denoted as the W95A model (Waleffe Reference Waleffe1995). The parameters of the second set remain the same except $\gamma _w = 0.5$ , and this configuration is denoted as W95B. In the last case, all parameters are set to 1, $\lambda _w = \mu _w= \nu _w =\sigma _w =\delta _w = 1$ , and this configuration is denoted as the BT model (Baggett & Trefethen Reference Baggett and Trefethen1997). It is important to note that these parameter sets significantly influence the system dynamics (Baggett & Trefethen Reference Baggett and Trefethen1997; Kalur et al. Reference Kalur, Seiler and Hemati2021).

Figure 7. The size of the ROA as the function of the Reynolds number in the case of the W95A model (a), W95B model (b) and BT model (c), computed using the GKE-G method. The QC, SOS, MS curves represent the results of Kalur et al. (Reference Kalur, Seiler and Hemati2021) and MS refers to minimal seeds obtained via direct-adjoint looping. For the W95B model, the minimal seed energy is taken from Cossu (Reference Cossu2005). The red, vertical dashed line indicates the unconditional stability limit, $\textit{Re}_E$ . The best-fitting curves of $\sqrt {e_{{min}}}$ : $\sqrt {e_{{min}}} \approx 176479 \,\textit{Re}^{-2.65089}$ for W95A; $\sqrt {e_{{min}}}\approx 1964.3 \,\textit{Re}^{-2.08506}$ for W95B; $\sqrt {e_{{min}}}\approx 4.1952 \,\textit{Re}^{-1.98875}$ for BT model.

The optimised transformation matrices are calculated for the W95A, W95B and BT models over different ranges of Reynolds numbers: 25–200 for the W95A model, 25–2000 for the W95B model and 5–100 for the BT model by the GKE-G method. For the W95A and BT models, the step size was set to 2.5, while for the W95B model, a logarithmic spacing was applied over 150 steps. Figure 7 shows the largest inner radius of the ROA for the three models.

The proposed GKE-G method is compared with other stability calculation methods. For the W95A and BT models, it yields nearly identical allowable perturbation levels as the SOS1 method employed by Kalur et al. (Reference Kalur, Seiler and Hemati2021), indicating that these are the largest provable ROA utilising quadratic Lyapunov functions.

Kalur et al. (Reference Kalur, Seiler and Hemati2021) also applied the QC method to the same W95A and BT models, predicting significantly smaller regions due to the approximation of nonlinear terms using bounds, despite its lower computational cost. A comparative analysis of the QC method across different model dimensions (TTRD’, W95A, BT) reveals a decreasing accuracy with increasing model complexity.

For the W95B model, the calculated allowable perturbation levels are compared with the results of Nerli et al. (Reference Nerli, Camarri and Salvetti2007). The presented GKE-G implementation slightly outperforms the GKE method of Nerli et al. (Reference Nerli, Camarri and Salvetti2007) due to a more general energy function form. (It is mentioned that Nerli et al. (Reference Nerli, Camarri and Salvetti2007) defined the kinetic energy with a multiplier of 1/2 which was compensated by a factor of $1/\sqrt {2}$ on the plots here.) Additionally, the calculated allowable threshold energy levels by the GKE-G method are close to the energy levels of minimal seeds calculated by Cossu (Reference Cossu2005). Since the minimal seeds represent states with the lowest kinetic energy converging to another attractor, the close agreement between stability threshold energy and minimal seed energy indicates accurate approximation of the stability region boundary. Similar observations hold for the BT model, where minimal seeds (MS) are calculated by Kalur et al. (Reference Kalur, Seiler and Hemati2021) using direct-adjoint looping.

In contrast, the results for the W95A model show different behaviour. Solutions initialised outside the predicted ROA tend to a laminar state, as also observed by Kalur et al. (Reference Kalur, Seiler and Hemati2021). This suggests that the true ROA is significantly larger than the predicted and might be better captured using higher-order Lyapunov functions.

In summary, the GKE-G method accurately predicts perturbation thresholds for the two-dimensional TTRD’ model and the four-dimensional BT and W95B models. However, for the W95A model, the method proves overly conservative, necessitating higher-order Lyapunov functions to capture a larger ROA.

3.2. Poiseuille flow models

To construct a reasonably accurate yet comparatively low-dimensional model of the two-dimensional Poiseuille flow we will use Galerkin projection on the basis of eigenfunctions of the Stokes operator. The construction method is described in Appendix B.

The chosen dimension of the domain is $L_x=2 \pi$ to ensure that the base wavenumber ( $\alpha _0=2 \pi /L_x=1$ ) aligns closely with the critical value identified through linear stability analysis, which is $\alpha = 1.02$ according to Orszag (Reference Orszag1971).

To determine an appropriate number of modes for the analysis, an initial assessment was conducted on the reduced-order model’s linear stability, followed by an examination of the original energy stability.

Initially, the linear stability limit ( $\textit{Re}_{{L}}$ ) is determined by identifying when the first eigenvalue of the linear part ( $\textbf {A}$ ) turns positive. The influence of the number of modes at a particular wavenumber pairs, $N_y$ within the range of 10–50 is depicted in table 1. It becomes apparent that, for fewer than 30 modes, the outcome of the linear stability analysis is markedly sensitive to the mode count, a phenomenon attributable to the high susceptibility of the non-normal linear operator to numerical inaccuracies, as highlighted by Trefethen & Embree (Reference Trefethen and Embree2005).

Table 1. The critical Reynolds numbers, as determined through linear stability analysis ( $\textit{Re}_L$ ) and energy stability analysis ( $\textit{Re}_E$ ), for the two-dimensional PFM with domain size $L_x=2 \pi$ and $N_z=0$ . Values selected for further analysis in this study are highlighted in bold.

Concurrently, the energy stability limit ( $\textit{Re}_{{E}}$ ) exhibits less sensitivity to the choice of modes $N_y$ and numerical errors. Notably, as the number of streamwise modes $N_x$ is increased from 1 to 2, a significant shift in the limit from 121 to 87 is observed. However, further increments in the number of streamwise modes do not influence the outcome. This indicates that the critical perturbation for classical energy analysis corresponds to a wavenumber of $2\alpha _0 = 2$ , whereas, according to linear stability analysis, the streamwise wavenumber of the critical perturbation is $\alpha _0 = 1$ for this configuration.

These two limits serve as pivotal benchmarks for the model and can be calculated with relative ease. Below the energy stability limit, the flow is deemed unconditionally stable, implying an infinite radius for the ROA. Conversely, beyond the linear stability limit, the flow becomes unconditionally unstable, reducing the radius of the ROA to zero. Within the bounds of these two limits, the proposed method offers a viable approach to ascertain the conditional stability threshold.

For the models under consideration, the threshold amplitudes are computed with $N_y=30$ . This provides a balance between the error margin of the linear stability limit and the number of modes.

The evaluation of the models spans Reynolds numbers from 250 to 2000, in increments of 250. Given the complexity associated with a high number of degrees of freedom, which reaches 90 for the smallest model studied, the SOS methods are deemed impractical for estimating the threshold amplitude. Instead, the GKE-G method is utilised for models up to $N_x = 5$ , where the results indicate convergence. For demonstration purposes, the GKE-SV method is applied to the largest systems up to $N_x = 16$ and $n=990$ .

The figure 8(a) presents the allowable perturbation amplitude for various models at $\textit{Re} = 500$ . It depicts the square root of the ratio of allowable perturbation kinetic energy to the base flow kinetic energy ( $E$ ), essentially showcasing the proportionality between the magnitude of perturbation velocity and the base flow velocity. This metric is henceforth referred to as the threshold amplitude ratio.

In the application of the GKE-G method, notable variations in results persist until the number of degrees of freedom reaches 210 ( $N_x=3$ ), beyond which increases to 270 ( $N_x=4$ ) or 330 ( $N_x=5$ ) result in only slight changes. It is important to note that the optimisation process of the transformation matrix does not utilise results from lower-order dynamical systems, rendering these results independently obtained. The convergence of the threshold amplitude ratio, therefore, suggests that the optimisation either identifies the global maximum of the allowable perturbation level or converges to a similar local maximum across various system sizes.

When the transformation matrix is optimised using the GKE-SV method, the predicted allowable perturbation level is lower, as anticipated. This reduction is due to the optimisation focusing solely on the singular values of the $\boldsymbol{\Psi }$ matrix and considering only the linear component of the dynamical system. However, these substantial simplifications in determining the optimal transformation matrix yield a threshold amplitude that is less than one magnitude smaller than that predicted by the GKE-G method. Conversely, the previously described QC method predicts significantly lower threshold amplitude levels, even for smaller dynamical systems.

Figure 8(b) compares the computational times between the two optimisation approaches. The benefits of employing the GKE-SV method are evident; not only is the computational time orders of magnitude lower, but it also facilitates application to larger dynamical systems, underscoring its efficiency and scalability.

Figure 8. (a) The square root of the ratio of allowable perturbation kinetic energy to the base flow kinetic energy, plotted as a function of degrees of freedom for a two-dimensional Poiseuille flow at $\textit{Re} = 500$ . The periodic domain length is set to $L_x = 2\pi$ . (b) The computational time required for the analysis.

Figure 9. The square root of the ratio of allowable perturbation kinetic energy to the base flow kinetic energy plotted against the Reynolds number for various PFMs. The transformation matrix optimisation is carried out using (a) the GKE-G method and (b) the GKE-SV method.

Next, the threshold amplitude ratio as a function of Reynolds number is depicted in figure 9. When the transformation matrix is optimised using the GKE-G method for models with $n = 210, 270, 330$ , the resulting curves closely align with one another, with only minor deviations observed at higher Reynolds numbers. The figure plots the threshold amplitude ratio on a logarithmic y-axis against a linear x-axis. The nearly linear curves suggest a decay rate faster than a power law, consistent with the flow’s linear instability at a finite Reynolds number, where the threshold amplitude approaches zero.

In figure 9(a), the threshold amplitudes obtained using the GKE-SV method are presented. The amplitude ratios are notably lower than those calculated via the GKE-G method, which can be attributed to the simpler optimisation process previously discussed. Additionally, fluctuations are observed around at $\textit{Re} = 1700$ , highlighting the sub-optimality of the GKE-SV method. These fluctuations result from the restrictions imposed on the transformation matrix, which vary with the Reynolds number but remain consistent across the models since they describe the same shear flow configuration. As the Reynolds number increases, the disparity between the results from the two methods widens, yet the difference remains under two magnitudes even at the highest Reynolds number evaluated.

The findings are juxtaposed with the minimal seed calculations conducted by Zhang & Tao (Reference Zhang and Tao2023), who determined that the laminar flow remains globally stable up to a Reynolds number of 2332.5. Zhang calculated the energy of the minimal seed at $\textit{Re}=2333$ to be approximately $e_{\mathit{MS}} = 0.01$ . From this value, the estimated amplitude ratio of the minimal seed, $\sqrt {e_{\mathit{MS}}/E}$ , is roughly 0.007. This represents the amplitude limit beyond which perturbations persist within the system. In our analysis utilising the 330-dimensional model with the GKE-G method, the threshold amplitude – up to which the laminar flow maintains stability – is approximately $\sqrt {e_{{min}}/E}=0.00045$ at $\textit{Re}=2000$ . This value is only one magnitude smaller than the minimal seed estimation. However, it is important to note that their results were obtained in a significantly larger domain of $L_x=400$ , compared with the domain length of $L_x = 2\pi$ employed in our study.

Figure 10. Perturbation field for the truncated PFM with $n=330$ at $\textit{Re}=1000$ . Panels (a) and (b) show the states with the smallest energy ( $e_{{min}}$ ) on the boundary of provable ROA calculated using the GKE-G method. Panels (c) and (d) depict the same states using the GKE-SV method. Panels (e) and (f) present the perturbations with the highest growth rate according to the original energy method. The streamwise and wall-normal velocity components are shown in the left and right columns, respectively.

In figure 10, panels (a–d) present the perturbations with the smallest kinetic energy ( $e_{{min}}$ ) at the boundary of the provable ROA for both the GKE-G and GKE-SV methods. In the case of the GKE-G method, these perturbations manifest as tilted periodic waves localised near the walls, resembling the minimal seed states described by Zhang & Tao (Reference Zhang and Tao2023). However, due to the significantly longer domain (100 units) used in their study, sometimes localised structures in the streamwise direction were observed, depending on the initial perturbation amplitude and Reynolds number. The shorter domain length ( $2 \pi$ ) in our study likely excludes the observation of such structures, resulting in fully periodic structures instead.

For the GKE-SV method, the perturbations appear as stripes localised between the wall and the domain’s centreline. These perturbations are less structured and exhibit asymmetry, likely due to the sub-optimal nature of the GKE-SV method.

For reference, panels (e–f) in figure 10, show the critical perturbation with the highest growth rate ( $\mu _e = 0.7712$ ) according to the original energy method. These perturbations also appear as tilted periodic waves, as demonstrated by Farrell (Reference Farrell1988), but are less localised to the near-wall region and exhibit higher oscillations in the streamwise direction compared with the critical perturbations according to the GKE-G method.

In the following phase of the analysis, power-law functions are fitted to the threshold amplitude as a function of Reynolds number, $\sqrt {e} \propto \textit{Re}^{\eta }$ . This fitting method has previously shown effectiveness for Couette flow, as documented by Duguet et al. (Reference Duguet, Monokrousos, Brandt and Henningson2013), and has been employed in other referenced studies. For Poiseuille flow, the exponents derived from the fitting range between −1.6 and −4.25, according to research conducted by Lundbladh et al. (Reference Lundbladh, Henningson and Reddy1994), Reddy et al. (Reference Reddy, Schmid, Baggett and Henningson1998), Parente et al. (Reference Parente, Robinet, De Palma and Cherubini2022) and Zhang & Tao (Reference Zhang and Tao2023).

The exponents predicted using the GKE-G method surpass these cited figures, yielding an exponent of −1.28 for the model with $n = 330$ dimensions. This discrepancy could stem from several factors. Notably, the current investigation focuses on two-dimensional flow, and the range of Reynolds numbers considered also differs. Additionally, as previously discussed, the relationship between threshold amplitude and Reynolds number in Poiseuille flow is likely to diverge from a simplistic power-law function, further contributing to the observed differences.

Table 2. Exponents of the power law for the threshold amplitude ( $\sqrt {e_{{min}}} \propto \textit{Re}^{-\eta }$ ), corresponding to half of the exponents for the threshold energy.