1. Introduction
Perturbation analyses investigate the response of a dynamical system to infinitesimally small disturbances, providing insights into its short- and long-term behaviours. Classical linear stability theory examines the eigenvalue spectrum of the governing operator to characterise the asymptotic fate of perturbations, with eigenvalues determining whether disturbances decay or grow exponentially over infinite time horizons (Schmid & Henningson Reference Schmid and Henningson2001).
However, significant limitations of asymptotic stability have been identified in practical applications. Experimental and numerical investigations have demonstrated that eigenvalue analyses fail to adequately describe the dynamics of disturbances. These failures were initially attributed to the nonlinear dynamics but were then found to be due to the transient growth of perturbations over a short time interval, which can be substantial even in asymptotically stable flows (Farrell Reference Farrell1988; Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). Further studies demonstrated that transient growth is essential for understanding the sub-critical transition to turbulence in channel flow (Reddy & Henningson Reference Reddy and Henningson1993), bypass transition to turbulence in boundary layer flow (Andersson, Berggren & Henningson Reference Andersson, Berggren and Henningson1999), etc. Consequently, transient growth analyses, focusing on the finite-time amplification of disturbances, have gained attention and offered a more comprehensive understanding of fluid physics. A range of physical mechanisms was then attributed to the transient dynamics. For instance, the ‘lift-up’ effect describes the transient generation of streamwise velocity streaks via the transfer of momentum across the boundary layer by streamwise vortices (Landahl Reference Landahl1980). Similarly, the Orr mechanism elucidates the process by which upstream-oriented perturbations are reoriented into downstream structures under the influence of base shear, thereby amplifying the disturbance energy (Farrell Reference Farrell1988). The anti-lift-up mechanism, in contrast, explains the transient formation of toroidal structures of azimuthal vorticity induced by azimuthal velocity streaks at the periphery of forced vortices (Antkowiak & Brancher Reference Antkowiak and Brancher2007).
From the mathematical perspective, the pseudo-spectrum analysis proves that transient growth is closely tied to the non-normality of the linearised operator, which enables eigenmodes to interact and produce significant short-term growth (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Tumin & Reshotko Reference Tumin and Reshotko2001; Abdessemed, Sherwin & Theofilis Reference Abdessemed, Sherwin and Theofilis2009). Optimal initial perturbations can therefore be interpreted as the superposition of eigenmodes, as demonstrated in Schmid & Henningson (Reference Schmid and Henningson2001). These theoretical advancements were then extended to canonical flows such as the boundary layer flow (Obrist & Schmid Reference Obrist and Schmid2003) and the Batchelor vortex (Mao & Sherwin Reference Mao and Sherwin2012), which reveal that transient growth is often dominated by the continuous part of the spectrum.
Similar to asymptotic instability, the computation of transient growth can be categorised into two branches: direct and iterative methods (Reddy & Henningson Reference Reddy and Henningson1993; Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008). The former requires the explicitly discretised matrix form of the evolution operator of perturbations and is only numerically feasible for low-dimensional problems, such as those in local studies with base flows that are homogeneous in two directions. In the singular value decomposition (SVD) of the evolution operator, the largest singular value corresponds to the square root of the optimal energy growth, while the corresponding singular vectors provide the optimal initial perturbation and its outcome (Schmid & Henningson Reference Schmid and Henningson2001). In more general scenarios, particularly for global analyses involving base flows that are inhomogeneous in multiple directions, the matrix forms of evolution operators are practically unavailable. For instance, the discretisation of the incompressible Navier–Stokes equations over a two-dimensional domain with moderate resolution (e.g. the backward-facing step flow presented in this work) typically leads to state vectors of dimension
$\textit {O}\! (10^{4} \sim 10^{5} )$
, while three-dimensional configurations can exceed
$\textit {O}\! (10^{6} )$
. As such, the full evolution operator would contain
$\sim \textit {O}\! (10^{12} )$
entries, requiring several terabytes of memory just to store, even before any computations are performed. In such conditions, iterative methods, including eigenvalue and optimisation approaches, which are mathematically equivalent (Mao, Blackburn & Sherwin Reference Mao, Blackburn and Sherwin2013), are employed. The former involves constructing a Krylov sequence by iteratively evolving perturbations forwards using the linearised operator and backwards using its adjoint, followed by extracting the dominant eigenvalue and eigenvector of the joint operator (Barkley et al. Reference Barkley, Blackburn and Sherwin2008). The latter determines the optimal initial perturbation by maximising a Lagrangian functional constrained by governing equations, again involving the iterative integration of the linearised operator and its adjoint (Schmid Reference Schmid2007).
The analytical form of the adjoint equations can be derived using integration by parts (Barkley et al. Reference Barkley, Blackburn and Sherwin2008) and is thus called the continuous adjoint. Although widely used in flow stability studies, this form introduces additional boundary terms that complicate implementation, particularly for compressible flows (De Pando et al. Reference De Pando, Sipp and Schmid2012; Raposo, Mughal & Ashworth Reference Raposo, Mughal and Ashworth2019) or in boundary perturbation analyses (Mao Reference Mao2015). Conversely, the discrete adjoint approach derives the adjoint operator directly as the transconjugate of the discretised forward operator. This approach ensures machine-precision consistency between the forward and adjoint operators, rendering it particularly advantageous for gradient-based iterative methods, but it involves complex numerical implementation of the boundary conditions for the adjoint. Considering the reliance on the adjoint for both the discrete and continuous approaches, we propose a data-driven methodology that captures the transient dynamics using only flow field snapshots. By bypassing adjoint equations, this data-driven approach simplifies the analysis process and provides a practical alternative for transient growth studies.
Data-driven analyses of perturbations have been extensively applied in fluid mechanics. One of the well-established techniques is dynamic mode decomposition (DMD) (Schmid Reference Schmid2010), which approximates the eigenvalues and eigenfunctions of a dynamical system from a set of snapshots (Wei & Rowley Reference Wei and Rowley2009; Kutz et al. Reference Kutz, Brunton, Brunton and Proctor2016a ). Over the past decade, DMD has been used to analyse a wide variety of flows, and numerous variations of its original scheme have been developed. These improvements have enabled DMD to handle non-uniform sampling (Tu Reference Tu2013), account for nonlinear observables to perform Koopman analyses for both deterministic and stochastic datasets (Williams, Kevrekidis & Rowley Reference Williams, Kevrekidis and Rowley2015) and improve its robustness to noise using L (Low-rank component) + S (Sparse component) decomposition (Wynn et al. Reference Wynn, Pearson, Ganapathisubramani and Goulart2013; Schmid Reference Schmid2022) or forward–backward mapping (Dawson et al. Reference Dawson, Hemati, Williams and Rowley2016). Furthermore, the processing of large datasets can be facilitated through parallelisation of QR-decomposition (Sayadi & Schmid Reference Sayadi and Schmid2016). By partitioning data sequences based on frequency, DMD also leads to improved resolution of flow features (Kutz et al. Reference Kutz, Fu and Brunton2016b ). Dynamic mode decomposition has also inspired other data-driven perturbation analysis methods. For example, Herrmann et al. (Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021) introduced an approach to perform resolvent analyses and demonstrated its efficiency on three-dimensional transitional channel flow. However, all of these schemes can only approximate a subset of the spectral features of the linear operator. Currently, there is still no robust data-driven framework to perform SVD of the evolution operator from snapshots, which is essential for transient growth analyses.
From a broader perspective, the concept underlying DMD and its related data-driven perturbation analysis techniques can be understood as a form of Galerkin projection (Williams et al. Reference Williams, Kevrekidis and Rowley2015; Colbrook Reference Colbrook2023). For example, standard DMD projects the linear operator onto proper orthogonal decomposition (POD) modes derived from snapshots, while optimal mode decomposition (Wynn et al. Reference Wynn, Pearson, Ganapathisubramani and Goulart2013) projects the same operator onto an optimal low-order subspace determined by a conjugate gradient algorithm. Similarly, data-driven resolvent analysis projects the resolvent operator onto approximated eigenmodes obtained through DMD. These projection methods are used to construct low-dimensional models, where the choice of the subspace significantly influences the accuracy. This is evident from the study in data-driven resolvent analysis, which shows that the resolvent modes extracted from trajectories initialised with random disturbances, optimal forcing and localised impulses can vary considerably, underscoring the sensitivity of the results to the subspace selection (Herrmann et al. Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021).
In the Galerkin framework, POD is often utilised to construct these basis sets since it retains the modes biased toward the large, high-energy scale (Girfoglio, Quaini & Rozza Reference Girfoglio, Quaini and Rozza2021) and the initialisation of trajectories for extracting POD modes often depends on the specific problem context within POD–Galerkin frameworks (Ballarin et al. Reference Ballarin, Manzoni, Quarteroni and Rozza2015; Ballarin & Rozza Reference Ballarin and Rozza2016; Hijazi et al. Reference Hijazi, Stabile, Mola and Rozza2020). In this study, we utilise both one-dimensional (1-D) and two-dimensional (2-D) Hermite polynomials for the initialisation in POD mode extraction, as they serve as the orthogonal basis for Cauchy problems (Williams et al. Reference Williams, Kevrekidis and Rowley2015). Hermite polynomials, originally introduced by Charles Hermite, play a key role in solving the Hermite differential equation and have broad applications in mathematics and physics (Aboites & Ramírez Reference Aboites and Ramírez2019). The 1-D Hermite polynomials describe the eigenfunctions of the quantum harmonic oscillator (Deniz Reference Deniz2019). When extended to two dimensions, Hermite polynomials are widely adopted in the study of quantum and optical systems to analyse and model complex wave patterns (van Dijk & Martens Reference van Dijk and Martens1996; Escalante-Ramírez & López-Caloca Reference Escalante-Ramírez and López-Caloca2008; Marisa et al. Reference Marisa, Niederhauser, Haeberlin, Wildhaber, Vogel, Jacomet and Goette2015). Their orthogonality and functional characteristics are expected to enhance the efficiency of capturing and representing key modes in complex systems. It should be noted that the choice of Hermite polynomials is not unique. Other families of orthogonal polynomials (e.g. Chebyshev or Legendre) could be used. The key requirement is to provide spatially diverse, linearly independent inputs, as the subsequent POD modes define the effective low-dimensional basis for the analysis.
In this work, we present an algorithm designed to extract transient energy growth and associated optimal initial perturbations and their outcomes directly from snapshots of perturbations. This method leverages POD to project the evolution operator onto a reduced-dimensional space. We demonstrate that by selecting appropriate initial conditions, e.g. in the form of Hermite polynomials, the POD modes of their trajectories are sufficient to incorporate optimal initial perturbations. This enables transient growth analysis through direct SVD of the evolution operator in a lower-dimensional space, eliminating the need for deriving and solving adjoint equations.
The structure of this paper is as follows. In § 2, we review the model-based formulation of transient growth analyses, then give a description of the proposed data-driven method in § 3 and demonstrate its effectiveness through three examples and discuss the limitations and potential applications in § 4. We further compare the proposed method with two potential DMD-based approaches in § 5 and clarify why the proposed framework is essential for accurate and general transient growth analyses. Finally, the conclusions are presented in § 6.
2. Model-based transient growth
Before delineating the data-driven transient growth analysis, we first briefly present the algorithms of model-based transient growth.
2.1. Definition of transient growth
For a general description of transient problems, we consider a linear dynamical system in the spatially discretised form
\begin{align} \partial _t \boldsymbol{u} = \boldsymbol{A} \boldsymbol{u}, \end{align}
where
$\boldsymbol{u} \in \mathbb{C}^N$
denotes the state vector with
$N$
representing the dimensionality of the discretised state vectors, and
$\boldsymbol{A} \in \mathbb{C}^{N \times N}$
is the linear dynamics operator. Specifically, the state
$\boldsymbol{u}$
represents the deviation from the base, which can be the steady state of a laminar flow or the temporal mean of a statistically stationary unsteady flow. The elements of
$\boldsymbol{A}$
are typically derived from the underlying physical model, such as the linearised Navier–Stokes (LNS) equations in fluid dynamics. According to (2.1), there is a linear evolution operator that projects the initial condition
$\boldsymbol{u}_0$
to the solution
$\boldsymbol{u}_{\tau }$
\begin{align} \boldsymbol{u}_{\tau } = \mathcal{M}(\tau ) \boldsymbol{u}_0, \end{align}
where
$\mathcal{M}(\tau )= e^{\boldsymbol{A} \tau } \in \mathbb{C}^{N \times N}$
if
$\boldsymbol{A}$
is time-independent.
To facilitate the presentation, we introduce the following inner product defined on the spatial domain:
\begin{align} (\boldsymbol{a}, \boldsymbol{b}) = \boldsymbol{a}^{{T}} \boldsymbol{Q} \boldsymbol{b}, \end{align}
where
$\boldsymbol{a}$
and
$\boldsymbol{b}$
are arbitrary vectors in the same inner product space. The superscript ‘
$\mathrm{T}$
’ denotes the complex conjugate, and the positive–definite matrix
$\boldsymbol{Q}$
accounts for the grid quadrature (Herrmann et al. Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021). Then the energy of the perturbation at time
$\tau$
can be quantified using the
$\boldsymbol{Q}$
-norm of
$\boldsymbol{u}_{\tau }$
\begin{align} E(\tau )=(\boldsymbol{u}_{\tau }, \boldsymbol{u}_{\tau })=\left \|\boldsymbol{u}_{\tau }\right \|_{\boldsymbol{Q}}^{2}. \end{align}
The
$\boldsymbol{Q}$
-norm can be related to the standard
$L_2$
-norm via
$\left \|\boldsymbol{u}\right \|_{\boldsymbol{Q}}^{2} = \left \|\boldsymbol{Fu}\right \|_{{2}}^{2}$
with
$\boldsymbol{F}$
derived from the Cholesky decomposition
$\boldsymbol{Q} = \boldsymbol{F}^{{T}}\boldsymbol{F}$
(Reddy & Henningson Reference Reddy and Henningson1993).
The aim of transient growth analyses is to quantify the maximum energy amplification from a given set of initial conditions of the state over a specified time horizon
$\tau$
. The optimal energy amplification can be expressed as
\begin{align} G(\tau ) = \max _{\boldsymbol{u}_0} \frac {E(\tau )}{E(0)} = \max _{\boldsymbol{u}_0} \frac {\left \|\mathcal{M}(\tau )\boldsymbol{u}_{0}\right \|_{\boldsymbol{Q}}^{2}}{\left \|\boldsymbol{u}_0\right \|_{\boldsymbol{Q}}^{2}} = \sigma _{\textit{max}}^{2} \big(\boldsymbol{F}\mathcal{M}(\tau )\boldsymbol{F}^{-1}\big), \end{align}
where
$\sigma _{\textit{max}}(\boldsymbol{\cdot })$
denotes the largest singular value of a matrix. Transient growth analysis can be performed using adjoint methods, including eigenvalue- and optimisation-based techniques. Mao et al. (Reference Mao, Blackburn and Sherwin2013) have demonstrated the equivalence of these two approaches for both optimal initial and boundary perturbation problems, and therefore we solely focus on the eigenvalue approach in this work.
2.2. Eigenvalue approach
Considering (2.5), when
$\mathcal{M}(\tau )$
is available in a discretised matrix form, transient growth analyses can be conducted using the SVD of the weighted operator
\begin{align} \boldsymbol{F}\mathcal{M}(\tau )\boldsymbol{F}^{-1} \boldsymbol{ {\varPsi }} = \boldsymbol{ {\varPhi }}\boldsymbol{ {\varSigma }} , \end{align}
where
$\boldsymbol{ {\varSigma }}= \mathrm{diag}({\sigma }_1, {\sigma }_2, \ldots , {\sigma }_{N}) \in \mathbb{R}^{N \times N}$
contains real and non-negative singular values and
${\sigma }_1 \geqslant {\sigma }_2 \geqslant , \ldots ,\geqslant {\sigma }_{N} \geqslant 0$
. The matrices
$\boldsymbol{ {\varPhi }}=[\boldsymbol{ {\phi }}_1, \boldsymbol{ {\phi }}_2, \ldots , \boldsymbol{ {\phi }}_{N}] \in \mathbb{C}^{N \times N}$
and
$\boldsymbol{ {\varPsi }}=[\boldsymbol{ {\psi }}_1, \boldsymbol{ {\psi }}_2, \ldots , \boldsymbol{ {\psi }}_{N}] \in \mathbb{C}^{N \times N}$
contain the left- and right-singular vectors, respectively. Clearly, the largest singular value
${\sigma }_1$
corresponds to the square root of the maximum energy amplification, while the associated right- and left-singular vectors,
$\boldsymbol{ {\psi }}_1$
and
$\boldsymbol{ {\phi }}_1$
, represent the optimal initial perturbation and its resulting optimal outcome, respectively.
In practice, however, the explicit matrix form of
$\mathcal{M}(\tau )$
is unavailable. Barkley et al. (Reference Barkley, Blackburn and Sherwin2008) introduced a matrix-free method for transient growth analyses. As defined by Courant & Hilbert (Reference Courant and Hilbert1962), there exists an adjoint operator
$\boldsymbol{L}^*$
for any linear operator
$\boldsymbol{L}$
such that
\begin{align} (\boldsymbol{b}, \boldsymbol{L} \boldsymbol{a}) = (\boldsymbol{L}^* \boldsymbol{b}, \boldsymbol{a}). \end{align}
Then the energy amplification in (2.5) can be expressed as
\begin{align} G(\tau )=\underset {{\boldsymbol{u}_{0}}}{\mathop {\max }}\,\frac {\left ( {\boldsymbol{u}_{0}},\mathcal{M}^*(\tau )\mathcal{M}(\tau ){\boldsymbol{u}_{0}} \right )}{\left ( {\boldsymbol{u}_{0}},{\boldsymbol{u}_{0}} \right )}, \end{align}
where
$\mathcal{M}^*(\tau ) \in \mathbb{C}^{N \times N}$
is the adjoint of
$\mathcal{M}(\tau )$
based on (2.7). The relationship between
$\mathcal{M}^*(\tau )$
and
$\mathcal{M}^{{T}}(\tau )$
is given by
\begin{align} \mathcal{M}^*(\tau ) =\boldsymbol{Q}^{-1} \mathcal{M}^{{T}}(\tau ) \boldsymbol{Q}. \end{align}
The eigenvalue approach involves projecting the joint operator
$\mathcal{M}^*(\tau )\mathcal{M}(\tau )$
onto a low-dimensional space spanned by the Krylov sequence, enabling direct computation of the leading eigenvalues and eigenvectors using the Arnoldi method
\begin{align} \mathcal{M}^*(\tau )\mathcal{M}(\tau )\boldsymbol{ {\varPsi }}_{\kern-1.5pt J} = \boldsymbol{ {\varSigma }}^2_J \boldsymbol{ {\varPsi }}_{\kern-1.5pt J}, \end{align}
where
$\boldsymbol{ {\varSigma }}^2_J = \mathrm{diag}({\sigma }^2_1, {\sigma }^2_2, \ldots , {\sigma }^2_{J}) \in \mathbb{R}^{J \times J}$
contains
$J$
leading eigenvalues of the joint operator, with
$J$
being the length of the Krylov sequence, and
$\boldsymbol{ {\varPsi }}_{\kern-1.5pt J} = [\boldsymbol{ {\psi }}_1, \boldsymbol{ {\psi }}_2, \ldots , \boldsymbol{ {\psi }}_{J}] \in \mathbb{C}^{N \times J}$
is the matrix of corresponding Ritz eigenvectors. This process requires
$J-1$
integrations of the joint operator, i.e. the linearised operator and its adjoint. Since the matrix
$\mathcal{M}^*(\tau )\mathcal{M}(\tau )$
is symmetric and self-adjoint, its eigenvalues are real, and the maximum energy amplification
$G(\tau )$
corresponds to the largest eigenvalue of this joint operator matrix. Here, the matrix
$\boldsymbol{ {\varPsi }}_{\kern-1.5pt J}$
corresponds to the right-singular vectors defined in (2.6), while the left-singular vectors in
$\boldsymbol{ {\varPhi }}_{\kern-1.5pt J}$
can also be obtained by post-processing computation.
Clearly, the process above relies on the adjoint equations. To circumvent this limitation, we project the forward operator onto a lower-dimensional subspace and subsequently bypass the derivation and calculation of the adjoint equations.
3. Data-driven optimal growth analysis
In this section, we outline the algorithm for data-based transient growth calculation, where the initial dataset and low-dimensional space construction are essential.
3.1. Initial matrix construction
This data-driven method employs the Galerkin projection to approximate the operator in a lower-dimensional space, a process that will be further elaborated in § 3.2. In this Galerkin framework, POD is typically used to construct basis functions, as it efficiently captures a low-dimensional representation of the dominant spatio-temporal structures in dynamical systems (Sirovich Reference Sirovich1987; Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993; Lumey Reference Lumey2012).
Supposing a set of initial conditions has been selected, we first capture
$K$
flow trajectories which are solutions to the linearised dynamical equations. The snapshots of the
$k$
th trajectory initialised by the
$k$
th element of a selected set (such as Hermite polynomials or their variations)
$\boldsymbol{q}^k_{t_0}$
, where
$k \in \{1, 2, \ldots , K\}$
, are assembled into the trajectory matrix
$\boldsymbol{D}^k_{1:S} \in \mathbb{C}^{N \times S}$
\begin{align} \boldsymbol{D}^k_{1:S} = \big[\boldsymbol{q}^k_{t_1}, \boldsymbol{q}^k_{t_2}, \boldsymbol{q}^k_{t_3}, \ldots , \boldsymbol{q}^k_{t_{S}}\big], \end{align}
where the subscripts
$t_1, t_2, \ldots , t_{S}$
denote discrete time points. To ensure satisfaction of the continuity constraint in fluid flow examples,
$\boldsymbol{q}^k_{t_0}$
is typically excluded from
$\boldsymbol{D}^k_{1:S}$
. In other words, the perturbed field is first integrated forward to
$t_1$
using the linear solver, and this short-time integration step guarantees that the subsequent velocity field satisfies the continuity equation to numerical accuracy. It is suggested that
$t_{S} \geqslant \tau$
, ensuring that the data contain sufficient transient information at time
$\tau$
where the transient mechanism is of interest. Then, a rank-
$R$
truncated SVD of
$\boldsymbol{D}^k_{1:S}$
yields the matrix
$[\boldsymbol{u}^{(k,1)}, \boldsymbol{u}^{(k,2)}, \boldsymbol{u}^{(k,3)}, \ldots , \boldsymbol{u}^{(k,R)}]$
, containing the first
$R$
left-singular vectors of
$\boldsymbol{D}^k_{1:S}$
, where
$R \leqslant S \ll N$
. These singular vectors are then assembled as columns of
$\tilde {\boldsymbol{X}}$
in the form of
\begin{align} \tilde {\boldsymbol{X}}= \Big [\boldsymbol{u}^{(1,1)}_0, \boldsymbol{u}^{(1,2)}_0, \ldots , \boldsymbol{u}^{(1,R)}_0 \mid \boldsymbol{u}^{(2,1)}_0, \boldsymbol{u}^{(2,2)}_0, \ldots , \boldsymbol{u}^{(2,R)}_0 \mid \ldots \mid \boldsymbol{u}^{(K,1)}_0, \boldsymbol{u}^{(K,2)}_0, \ldots , \boldsymbol{u}^{(K,R)}_0 \Big ], \end{align}
where
$\tilde {\boldsymbol{X}}$
is of size
$N \times (K \times R)$
and the subscript ‘
$0$
’ is added to highlight that they will be used as a pool of initial conditions in the following process.
In practice, the selection of
$K$
or
$R$
may vary considerably, depending on the number of POD modes required to accurately synthesise the optimal modes for the problem. Specifically, the rank truncation parameter
$R$
is determined using a cumulative energy criterion based on the singular values of
$\boldsymbol{D}^k_{1:S}$
. In this study,
$R$
is chosen such that the retained modes capture at least 99 % of the total energy, which has been found to be sufficient to accurately represent the dominant transient dynamics while filtering out low-energy noise. In case studies, a fixed value of
$R$
is selected based on the singular values of
$\boldsymbol{D}^1_{1:S}$
, and the size of
$\tilde {\boldsymbol{X}}$
is increased by expanding
$K$
to improve convergence while retaining the same number of SVD modes for
$\boldsymbol{D}^2_{1:S}, \boldsymbol{D}^3_{1:S}, \ldots , \boldsymbol{D}^K_{1:S}$
. As a practical rule of thumb, as
$K$
or
$R$
increases, the additional vectors in the reduced subspace
$\mathbb{V}$
, spanned by the SVD modes of
$\boldsymbol{X}$
, tend to become linearly dependent if the dynamical system is low rank. This saturation of the subspace can serve as a practical guideline for selecting appropriate values of
$K$
or
$R$
. Other modal decomposition techniques, e.g. the spectral POD (Nekkanti & Schmidt Reference Nekkanti and Schmidt2021), balanced POD (Flinois, Morgans & Schmid Reference Flinois, Morgans and Schmid2015) or even convolutional neural networks (Murata, Fukami & Fukagata Reference Murata, Fukami and Fukagata2020), in addition to snapshot POD, may be used here to reduce the memory requirements of the proposed data-driven algorithm and potentially extend it to more complex scenarios. However, while these methods can increase computational efficiency in subsequent steps of § 3.2, they introduce additional computational costs during the initial matrix construction process or are not completely adjoint free.
To effectively capture transient information, Hermite polynomials are introduced to enhance the POD mode generation process in the selection of
$\boldsymbol{q}^k_{t_0}$
. The details of both 1-D and 2-D forms of Hermite polynomials are given in Appendix A. The POD modes, derived from trajectories which are initialised with Hermite polynomials or their variants, are demonstrated to construct accurate reduced-order models. Here, using the Hermite polynomials as the initial guess significantly improves the efficiency of the scheme, as will be compared against random noise in § 4.2. Moreover, initialising the data matrix with POD-based modes alleviates the ill conditioning that commonly arises in data-driven reduced-order models, an effect that will be further examined in § 5.2.
3.2. Data-driven algorithm
The proposed data-driven transient growth analysis relies on the Galerkin projection of the governing operator onto a space spanned by orthogonal basis vectors, allowing the system state to be approximated as a combination of these vectors. For the convergence test, we establish an initial matrix
$\boldsymbol{X}$
, composed of the first
$n$
column vectors of
$\tilde {\boldsymbol{X}}$
, such that
$\tilde {\boldsymbol{X}} \supseteq \boldsymbol{X} \in \mathbb{C}^{N \times n}$
, with
$1 \leqslant n \leqslant (K \times R)$
and
$n \ll N$
. To facilitate the process of matrix manipulation and low-dimensional space construction, we perform a full-rank SVD of
$\boldsymbol{X}$
\begin{align} \boldsymbol{X}=\boldsymbol{ {V}} \boldsymbol{ {\varLambda }}\boldsymbol{W}^{{T}}, \end{align}
where
$\boldsymbol{ {V}} = [\boldsymbol{v}_1, \boldsymbol{v}_2, \ldots , \boldsymbol{v}_n] \in \mathbb{C}^{N \times n}$
is a unitary matrix and spans a reduced basis space
$ \mathbb{V} = span{\{\boldsymbol{v}_i\}^n_{i=1}}$
. Here,
$\boldsymbol{ {\varLambda }} \in \mathbb{R}^{n \times n}$
is a diagonal matrix of singular values and
$\boldsymbol{W} \in \mathbb{C}^{N \times n}$
is the matrix of right-singular vectors.
By evolving the column vectors of
$\boldsymbol{X}$
forward in time from
$t=0$
to
$ t=\tau$
using (2.1), where the transient growth analysis is to be performed, aligning with (2.2), we obtain
\begin{align} \boldsymbol{Y_{\tau }} = \mathcal{M}(\tau )\boldsymbol{X}, \end{align}
where
$\boldsymbol{Y_{\tau }} \in \mathbb{C}^{N \times n}$
is called the outcome matrix here. For clarity in the subsequent derivation, we use
$\boldsymbol{u}_{0}$
to denote a column vector from
$\boldsymbol{X}$
, and
$\boldsymbol{u}_{\tau }$
to denote the corresponding column in
$\boldsymbol{Y_{\tau }}$
. It should be noted that we can collect the snapshots at an arbitrary time instance
$t \in [0,\tau ]$
during a single forward integration, enabling the computation of transient growth over multiple time instances. This contrasts with model-based methods, which require renewed forward and backward integration for each time instance.
Assuming that
$\boldsymbol{Y_{\tau }}$
can be expressed as linear combinations of columns of
$\boldsymbol{ {V}}$
, which is reasonable given that
$\boldsymbol{X}$
is constructed by POD modes from transient trajectories as described in § 3.1, then the final and initial states can be approximated in the space
$\mathbb{V}$
as
\begin{align} \boldsymbol{u}_{\tau }=\boldsymbol{V} \hat {\boldsymbol{u}}_{\tau }, \quad \boldsymbol{u}_{0}=\boldsymbol{V} \hat {\boldsymbol{u}}_{0}, \end{align}
with
$\hat {\boldsymbol{u}}_{\tau }$
and
$\hat {\boldsymbol{u}}_0$
being the reduced states of size
$n \times 1$
.
Considering the projection in (3.5), we can then construct a reduced-order forward operator
$\hat {\mathcal{M}}(\tau ) = \boldsymbol{V}^{{T}} {\mathcal{M}}(\tau ) \boldsymbol{V}$
, which satisfies
\begin{align} \hat {\boldsymbol{u}}_{\tau }=\hat {\mathcal{M}}(\tau )\hat {\boldsymbol{u}}_{0}. \end{align}
Similar to the approximation process of the linear operator
$e^{\boldsymbol{A} {\Delta t}}$
in DMD studies (Schmid Reference Schmid2022), with
$\Delta t$
being the data sampling interval,
$\hat {\mathcal{M}}(\tau )$
can be computed by substituting the SVD in (3.3) into (3.4) as follows:
\begin{align} \hat {\mathcal{M}}(\tau ) = \boldsymbol{V}^{{T}} \boldsymbol{Y_{\tau }}\boldsymbol{W}\boldsymbol{ {\varLambda }}^{-1}. \end{align}
Since the energy is defined in the discrete form, in this reduced subspace, we need to adjust the inner product to retain the physical meaning in the
$L_2$
-norm framework
\begin{align} \|\boldsymbol{u} \|_{\boldsymbol{Q}}^{2}=\hat {\boldsymbol{u}}^{{T}} \boldsymbol{V}^{{T}} \boldsymbol{Q}\kern-1pt\boldsymbol{V} \hat {\boldsymbol{u}}= \|\kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \hat {\boldsymbol{u}} \|_{2}^2, \end{align}
where the matrix
$\kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \in \mathbb{C}^{n \times n}$
can be derived from the Cholesky factorisation of
$\boldsymbol{V}^{{T}} \boldsymbol{Q}\kern-1pt\boldsymbol{V} = \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}^{{T}} \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}$
. Thus, the optimal energy growth can be expressed as
\begin{align} G(\tau ) = \max _{\boldsymbol{u}_0} \frac {\big \|\boldsymbol{u}_{\tau }\big \|_{\boldsymbol{Q}}^{2}}{\big \|\boldsymbol{u}_0\big \|_{\boldsymbol{Q}}^{2}} = \max _{\boldsymbol{u}_0} \frac {\big \| \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \hat {\boldsymbol{u}}_{\tau }\big \|_{2}^{2}}{\big \| \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \hat {\boldsymbol{u}}_0\big \|_{2}^{2}} = \max _{\boldsymbol{u}_0} \frac {\big \| \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \hat {\mathcal{M}}(\tau ) \hat {\boldsymbol{u}}_{0} \big \|_{2}^{2}}{\big \| \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \hat {\boldsymbol{u}}_0\big \|_{2}^{2}} = {\sigma }_{\textit{max}}^2 \big( \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \hat {\mathcal{M}}(\tau ) \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}^{-1}\big). \end{align}
Then a standard SVD can be performed to obtain the modes and energy growth in the
$n$
-dimensional vector space
$\mathbb{V}$
\begin{align} \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}} \hat {\mathcal{M}}(\tau ) \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}^{-1} \hat {\boldsymbol{ {\psi }}}_1 = \hat {{ {\sigma }}}_1 \kern1.5pt \hat{\boldsymbol{\kern-1.5pt \phi}}_1, \end{align}
where the largest singular value
$\hat { {\sigma }}_1$
gives the square root of the transient energy growth at time
$\tau$
and the modes in original coordinates as (2.6) are represented by
\begin{align} {{\boldsymbol{ {\psi }}}_1 = \boldsymbol{V} \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}^{-1} \hat {\boldsymbol{ {\psi }}}_1, \quad \boldsymbol{ {\phi }}}_1 = \boldsymbol{V} \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}^{-1} \kern1.5pt \hat{\boldsymbol{\kern-1.5pt \phi}}_1. \end{align}
The
$\boldsymbol{Q}$
-norm error
${e}_{\textit{mode}}$
is defined to quantitatively evaluate the performance of this data-driven method in extracting optimal modes and the relative error
${e}_{{ growth}}$
characterises the error in optimal energy growth calculations, as follows:
\begin{align} {e}_{\textit{mode}} (\tau ) = \frac {\big \|\boldsymbol{\zeta } - {\boldsymbol{\zeta }}_{{d}} \big \|_{\boldsymbol{Q}}^{2}}{\left \|\boldsymbol{\zeta }\right \|_{\boldsymbol{Q}}^{2}}, \quad {e}_{ {growth}} (\tau ) = \frac {|{\varepsilon } - {{\varepsilon }}_{{d}}|} {{\varepsilon }}, \end{align}
where
$\varepsilon$
denotes the energy growth and
$\boldsymbol{\zeta }$
represents the optimal initial perturbation or outcome, both serve as the reference from the model-based method here, while the subscript ‘
$\mathrm{d}$
’ refers to results obtained using the data-driven method. Convergence with respect to parameters, particularly the number of initial conditions
$n$
, can then be assessed by comparing the predicted transient growth and associated modes with reference model-based results. The criterion adopted here is that
${e}_{\textit{growth}}$
and
${e}_{\textit{mode}}$
both fall below
$10^{-2}$
.
4. Results
In this section, we apply the established data-driven transient growth analysis algorithm to three representative examples, i.e. the complex Ginzburg–Landau equation with a single velocity component, the backward-facing step flow with two velocity components featuring the Orr transient mechanism (Marquet & Sipp Reference Marquet and Sipp2012) and the Batchelor vortex with three velocity components demonstrating the anti-lift-up transient mechanism (Antkowiak & Brancher Reference Antkowiak and Brancher2007; Mao & Sherwin Reference Mao and Sherwin2012). These examples, with increasing complexity, are well established in the context of model-based hydrodynamic stability and transient growth studies whose results agree well with those of the proposed data-driven approach as will be presented. Moreover, while the proposed method shares similarities with DMD in terms of operator approximation, they differ fundamentally in the objective, construction process and data requirements. These distinctions will be discussed in detail in § 5.
4.1. Complex Ginzburg–Landau equation
Data-driven transient growth analysis of the linearised complex Ginzburg–Landau equation.
$(a)$
The relative error between the model-based and the data-driven optimal energy growth as a function of the number of initial conditions in the dataset
$\boldsymbol{X}$
.
$(b)$
The same as
$(a)$
, but for the
$\boldsymbol{Q}$
-norm error between the optimal initial perturbations (input) and outcomes (output).
$(c)$
Optimal energy growth computed by the data-driven method with
$n=8$
(circle) and the model-based method (line).
$(d)$
Optimal initial perturbations and their corresponding outcomes obtained by the data-driven method with
$n=8$
at
$\tau =5$
,
$8.5$
,
$15$
. Solid and dashed lines denote the real and imaginary parts of the modes, respectively, while the grey lines indicate the model-based results for comparison.
The complex Ginzburg–Landau equation is a widely used model for characterising both the short-term and long-term evolutions of small perturbations in spatially developing flows (Bagheri et al. Reference Bagheri, Henningson, Hoepffner and Schmid2009) and can be written as
\begin{align} \frac {\partial \tilde {u}}{\partial t} = -c_a \frac {\partial \tilde {u}}{\partial x} + c_d\frac {\partial ^2 \tilde {u}}{\partial x^2} + c_s(x)\tilde {u} - c_n|\tilde {u}|^2 \tilde {u}, \end{align}
where the variable
$\tilde {u}(x, t)$
represents the complex amplitude,
$c_n$
the complex weight of nonlinearity,
$c_s(x)$
the spatial amplification parameter,
$c_a$
the complex advection speed and
$c_d$
the complex diffusion coefficient. By assuming the flow is near the equilibrium state
$\tilde {u} = 0$
, the linearised governing equation is given by
\begin{align} \frac {\partial {u}}{\partial t} = \left (-c_a \frac {\partial }{\partial x}+c_d\frac {\partial ^2}{\partial x^2}+c_s(x)\right ) {u}, \end{align}
where
$u$
denotes small deviations from the equilibrium state, and the corresponding adjoint equation can be represented by
\begin{align} -\frac {\partial {u}^*}{\partial t} = \left ({c_a}^{{T}} \frac {\partial }{\partial x} + {c_d}^{{T}} \frac {\partial ^2}{\partial x^2}+{c_s}^{{T}}(x)\right ){u}^*. \end{align}
The parameters in above equations are selected as
$c_a=2+0.4\mathrm{i}, c_d=1-\mathrm{i}$
and
$ c_s(x)= 0.19-0.005{x}^2$
, as adopted by Herrmann et al. (Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021), ensuring linearly stable dynamics. The computational domain extends over the range
$ -85 \lt x \lt 85$
with homogeneous boundary conditions and is decomposed into
$ N = 220$
uniformly spaced points. In a convergence test, the current resolution results in a relative difference below
$ 10^{-6}$
compared with that at
$ N = 440$
. The exponential time differencing method, as outlined in Cox & Matthews (Reference Cox and Matthews2002), is used to evolve both the forward and adjoint equations in the model-based method for comparison.
Following the process outlined in § 3, to construct
$\tilde {\boldsymbol{X}}$
, we employ the first
$K = 6$
Gauss-weighted 1-D Hermite polynomials to generate trajectories as described in (3.1). A total of
$S=100$
snapshots, sampled at intervals of
$\Delta t=0.5$
time units, are recorded for each trajectory. Subsequently,
$R = 5$
POD modes are retained to form
$\tilde {\boldsymbol{X}}$
, which comprises
$30$
column vectors. We then track their evolution to compute both the transient growth and the optimal modes for the time instance
$\tau \in \{1,2,\ldots ,50\}$
.
To assess the convergence performance of the data-driven method, we evaluate the relative error of the optimal energy growth (figure 1
a) and the
$\boldsymbol{Q}$
-norm error of the modes (figure 1
b) as functions of
$n$
, i.e. the number of columns in the initial matrix
$\boldsymbol{X}$
. The results indicate that increasing the number of initial conditions improves the agreement with model-based results. Notably, we can see from these curves that using only
$n=8$
initial conditions is sufficient to perform the data-driven algorithm in this case. Furthermore, since snapshots are collected at several time instances during the assembly of the outcome matrix
$\boldsymbol{Y}_{50}$
, as described in § 3.2, it becomes computationally efficient to compute the envelope of the optimal energy growth, which is in good agreement with that from the model-based method as shown in figure 1
$(c)$
.
The optimal initial perturbations and outcomes at
$\tau =5, 8.5, 15$
, using the data-driven method with
$n=8$
, are shown in figure 1
$(d)$
, where we can see that the data-driven method supports successful extraction of transient modes. These results demonstrate that the data-driven transient growth analysis offers an accurate approximation of the transient behaviour over short time intervals in a dynamical system.
4.2. Flow over a backward-facing step
The second example is the backward-facing step flow, where the fluid motion is governed by the incompressible Navier–Stokes (NS) equations
\begin{align} \partial _t \tilde {\boldsymbol{u}}=-\tilde {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {\boldsymbol{u}}-\boldsymbol{\nabla }\tilde {p} + Re^{-1} {\nabla} ^2 \tilde {\boldsymbol{u}}, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\tilde {\boldsymbol{u}}=0, \end{align}
with
$\tilde {p}$
representing the modified or kinematic pressure and
$\tilde {\boldsymbol{u}}$
denoting the velocity vector field. The Reynolds number
$ \textit{Re}$
is set to
$500$
, where the bifurcation to a 3-D state at
$ \textit{Re}=748$
is avoided to ensure an asymptotically stable flow (Barkley, Gomes & Henderson Reference Barkley, Gomes and Henderson2002). At the inflow boundary
$\partial \varOmega _i$
, a parabolic velocity profile is prescribed. On the solid walls
$\partial \varOmega _w$
, no-slip conditions are enforced. At the downstream outflow boundary
$\partial \varOmega _o$
, a zero-traction condition is applied, with the pressure set to zero. Collectively, the boundary conditions can be expressed as
\begin{align} \begin{aligned} \tilde {\boldsymbol{u}}(\partial \varOmega _i, t) & =\big (4 y-4 y^2, 0,0\big ), \\ \tilde {\boldsymbol{u}}(\partial \varOmega _w, t) & =(0,0,0), \\ \partial _x \tilde {\boldsymbol{u}}(\partial \varOmega _o, t) & =(0,0,0), \quad \tilde {p}(\partial \varOmega _o, t)=0. \end{aligned} \end{align}
By decomposing the flow field into a base flow and a perturbation and omitting the interaction of the perturbation with itself, we obtain the LNS equations that govern the evolution of perturbations
\begin{align} \partial _t \boldsymbol{u} =-(\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}-\left (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{U}-\boldsymbol{\nabla }p+R e^{-1} {\nabla} ^2 \boldsymbol{u}, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} =0, \end{align}
where
$\boldsymbol{u}$
represents the perturbation of velocity with two components
$(u, v)^{{T}}$
in the streamwise and vertical directions. Here,
$\boldsymbol{U}$
is the base flow velocity satisfying the steady NS equations and
$p$
is the pressure perturbation. The adjoint equations, as outlined in Barkley et al. (Reference Barkley, Blackburn and Sherwin2008), take the form
\begin{align} \partial _t \boldsymbol{u}^{*} =(\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{u}^{*}+\left (\boldsymbol{u}^{*} \boldsymbol{\cdot }\boldsymbol{\nabla }\right ) \boldsymbol{U}+\boldsymbol{\nabla }p^{*}-Re^{-1} {\nabla} ^2 \boldsymbol{u}^{*} \quad \text{ with } \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{*} =0, \end{align}
and are subject to the boundary condition
$\boldsymbol{u}^{*}(\partial \varOmega , t) =(0,0,0)$
on the entire domain boundary
$\partial \varOmega$
. This condition, in combination with a homogeneous Dirichlet boundary condition for perturbations, i.e.
$\boldsymbol{u}(\partial \varOmega , t) =(0,0,0)$
, ensures consistency between the forward and adjoint systems (Barkley et al. Reference Barkley, Blackburn and Sherwin2008). The NS (4.4), LNS (4.6) and adjoint (4.7) equations in the computational domain depicted in figure 2 are discretised using the spectral/hp element method (Karniadakis & Sherwin Reference Karniadakis and Sherwin2005) and a polynomial order of 6 is used to ensure computational accuracy (Mao Reference Mao2015).
Spectral elements in the computational domain: (a) overall domain and (b) domain close to the step.
To implement the data-driven approach, we start by generating snapshots with
$K=240$
trajectories, where the initial conditions are from 2-D Hermite polynomials in Cartesian coordinates. Following the formula in (A6), the first
$12$
polynomials are employed with
$C=1$
, which defines the non-zero regions. The polynomials are centred along the
$x$
-direction, extending from
$(x,y)=(-9.5,0.5)$
to
$(x,y)=(9.5,0.5)$
with a spacing of
$1$
. To study the short-term energy growth within
$\tau \in [0, 100]$
which have been studied in Blackburn, Barkley & Sherwin (Reference Blackburn, Barkley and Sherwin2008), each trajectory evolves over a time span of
$t_{S}=100$
with interval
$\Delta t=0.25$
, ensuring that the recorded snapshots capture the key spatial features of the transient dynamics. For ease of implementation in numerical simulations, only the streamwise velocity component is initially perturbed. The field at time
$t_1=0.25$
, denoted by
$\boldsymbol{q}^k_{t_1}$
, satisfies the continuity equation to numerical accuracy, as discussed in § 3.1.
It is crucial to note that the components of the velocity field at each grid point in the
$x$
and
$y$
directions are concatenated into a single column vector, resulting in
$N \gt 5 \times 10^4$
. To optimise the memory usage, only the first
$R=40$
POD modes per trajectory are retained. We then assemble the obtained POD modes into the initial matrix
$\boldsymbol{X}$
, which is arranged according to the rank of Hermite polynomials. The column vectors in
$\boldsymbol{X}$
are then evolved forward in time to generate the outcome matrices
$\boldsymbol{Y}_{\tau }$
at specific time instances for transient growth analyses according to § 3.2. For validation purposes, transient growth analysis is performed at
$\tau =60$
, corresponding to the evolution time of optimal initial perturbations identified in Blackburn et al. (Reference Blackburn, Barkley and Sherwin2008).
To illustrate the importance of the construction of
$\boldsymbol{X}$
, we compare the results using the data generated by Hermite polynomials and random noise, which are zero mean, spatially uncorrelated and with a fixed standard deviation applied to the streamwise velocity component at each mesh node, and present the relative error of energy growth and the
$\boldsymbol{Q}$
-norm error of the modes in figures 3
$(a)$
and 3
$(b)$
, respectively. These figures reveal consistent convergence of both energy growth and modes as the amount of data increases. However, compared with Hermite polynomials, the convergence speed of growth and modes using random noise is much slower and the errors of transient growth and modes are around
$10 \,\%$
even with over
$n=9000$
initial conditions in
$\boldsymbol{X}$
, which is approximately the upper limit of a computer with 64 GB RAM. In contrast, the errors can be reduced to below 1 % using only
$n=4500$
initial conditions generated by the data from Hermite polynomials. Given the original forward operator dimensionality, this represents over two orders of magnitude reduction. Besides, the transient information at any time instance in
$\tau \in (0,60)$
can be obtained by collecting the corresponding outcome matrix
$\boldsymbol{Y}_{\tau }$
as discussed in § 3.2, which can be more efficient than the model-based method when a large number of
$\tau$
are considered.
Comparison of data-driven transient growth analyses of the backward-facing step flow using the dataset obtained from LNS initialised with random noise and Hermite polynomials.
$(a)$
Relative error between the model-based and data-driven transient growth at
$\tau =60$
as a function of the number of initial conditions in the dataset
$\boldsymbol{X}$
.
$(b)$
Corresponding
$\boldsymbol{Q}$
-norm error for optimal initial perturbations (Input) and outcomes (output).
The data-based optimal initial perturbations at
$\tau = 60$
from Hermite polynomials and random noise using different numbers of initial conditions in the dataset
$\boldsymbol{X}$
.
To further illustrate the benefit of employing the Hermite polynomials, the optimal initial perturbations at
$\tau = 60$
using the two forms of basis sets are shown in figure 4. It can be observed that the data generated from random noise are far less efficient in predicting the optimal initial perturbation, particularly around the vicinity of the step edge where the energy is sharply concentrated.
Unlike the complex Ginzburg–Landau equation case, the convergence speed of the optimal initial perturbations is slightly slower than that of the optimal outcomes, which is consistent with observations from the data-driven resolvent analysis (Herrmann et al. Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021). This can be explained by the fact that the spatial support of the POD modes tends to resemble the response modes rather than the forcing modes, which then require more data to converge.
4.3. Axially decomposed Batchelor vortex
To validate the generalisability of the methodology, the Batchelor vortex with three velocity components, which is widely studied in stability and transient growth analyses, is considered. According to Batchelor (Reference Batchelor1964), the non-dimensional base flow velocity
$\boldsymbol{U}=(U, V, W)^{{T}}$
in cylindrical coordinates
$(r,\theta ,z)$
, where
$z$
denotes the axial flow direction, can be described with
\begin{align} \left \{\begin{aligned} U & =0, \\ V & =q\frac {1-e^{-r^2}}{r}, \\ W & =a + e^{-r^2}, \end{aligned}\right . \end{align}
where the viscous decay term is neglected to keep the base in a steady form when considering instability or non-normality problems. The perturbations, which are governed by the LNS equations (see (4.6)), are decomposed in the axial direction in both model-based and data-driven global analyses. We adopt the parameters consistent with those detailed in Mao & Sherwin (Reference Mao and Sherwin2012), where the external non-dimensional free-stream axial velocity is
$a=0$
and the swirl strength
$q=3$
. The Reynolds number is set to
$ \textit{Re}=1000$
and the streamwise wavenumber
$0$
, which maximises the local transient growth, is adopted. The numerical solver employed is consistent with that in § 4.2, and the computational domain and mesh configuration are depicted in figure 5. For convenience, the mesh and subsequent result contours are presented in the Cartesian coordinate system
$(x,y,z)$
, where
$z$
denotes the streamwise direction.
Computational mesh for transient growth calculation in the axially decomposed global analysis of the Batchelor vortex.
$(a)$
Complete mesh;
$(b)$
mesh around the vortex core.
After discretisation, the dimension of the velocity vector exceeds
$2.3\times 10^5$
, an order of magnitude larger than that of the backward-facing step flow case. Due to this high-dimensionality, significantly more data are required to ensure accurate results and to mitigate the curse of dimensionality (Bellman Reference Bellman1966).
One way to minimise the initial conditions required is to pre-determine the azimuthal wavenumber of the optimal initial perturbation since the base flow, in this case, is homogeneous in the axial direction and the optimal mode is with a single axial wavenumber (Bale & Govindarajan Reference Bale and Govindarajan2010). Initially, we select the first
$10$
Hermite polynomials in polar coordinates as in (A7) with the parameter
$C \in \{1,2,3, \ldots , 20\}$
. We generate
$m=200$
trajectories initialised with these polynomials by integrating LNS equations forwards, recording 1000 snapshots per trajectory at time intervals of
$\Delta t=0.6$
. The initial matrix,
$\boldsymbol{X}$
, is constructed using the first
$R=10$
POD modes, resulting in the initial conditions matrix
$\boldsymbol{X}$
with
$n=2000$
. This set-up is sufficient to determine the azimuthal wavenumber of the optimal initial perturbation. Once the wavenumber is obtained, the same procedure is applied using Hermite polynomials with various
$C$
but a fixed azimuthal wavenumber. While this reduces the number of initial conditions, it remains computationally inefficient because the polynomials are not strictly spatially orthogonal within the domain, leading to overlapping regions and thus computational redundancy.
In an alternative approach, we utilise the first
$150$
SVD modes of
$1200$
polynomials in polar coordinates for the selection of
$\boldsymbol{q}^k_{t_0}$
in (3.1). These polynomials correspond to the first 6 Hermite polynomials, with the radius parameter
$C \in \{0.1, 0.2, 0.3, \ldots , 20\}$
, where the smaller interval is used to enhance the resolution for capturing features near the vortex core.
Following the steps in § 3.1,
$K=150$
trajectories are then generated and each is recorded with
$S=1250$
snapshots with time interval
$\Delta t=0.048$
. It should also be noted that only the radial velocity component
$u$
is perturbed, similarly as discussed in § 4.2. We retain
$R=30$
POD modes per trajectory for convergence analyses and the initial matrix
$\boldsymbol{X}$
construction. To ensure reliable convergence, the matrix
$\tilde {\boldsymbol{X}}$
, from which the initial matrix
$\boldsymbol{X}$
is constructed, is organised by ordering the POD modes. This systematic arrangement eliminates possible convergence issues arising from trajectory variations and results in smoother convergence curves, as shown in figure 6
$(a)$
. Then the same methodology in § 3.2 is applied to conduct data-driven transient growth analyses at
$\tau = 15, 30, 45, 60$
.
Figures 6
$(a)$
, 6
$(b)$
and 6
$(c)$
effectively demonstrate the reliability of the data-driven methodology, even when applied to high-dimensional systems. Consistent with that in the backward-facing step case discussed in § 4.2, the optimal initial perturbations converge slightly more slowly than the optimal outcomes. Interestingly, the number of column vectors required in
$\boldsymbol{X}$
is similar to that in § 4.2, despite the fact that the higher-dimensionality typically demands more data. This suggests that the amount of data needed is not only strictly affected by the system’s dimensionality but also influenced by the quality of the POD modes used to synthesise the optimal initial perturbations. This observation highlights the importance of well-chosen modes in ensuring an efficient and accurate representation of the system’s dynamics.
Data-driven transient growth analysis of the axially decomposed Batchelor vortex at
$ \textit{Re}=1000$
and axial wavenumber 0.
$(a)$
The relative error between the model-based and the data-driven optimal energy growth at
$\tau =15$
,
$30$
,
$45$
,
$60$
as a function of the number of initial conditions in the dataset
$\boldsymbol{X}$
.
$(b)$
The same as
$(a)$
, but for the
$\boldsymbol{Q}$
-norm error between the optimal initial conditions.
$(c)$
The
$\boldsymbol{Q}$
-norm error between the optimal outcomes.
$(d)$
The envelope of optima (hollow circle) with curves of linear energy evolution starting from the four optimal initial conditions, which are calculated through the data-driven method. The solid circles mark the points at which the curves of linear energy growth are tangential to the envelope.
Comparison of optimal modes computed by data-driven and model-based methods.
$(a)$
The optimal initial perturbations and outcomes at
$\tau = 15$
, the velocity in radial
$(u)$
, azimuthal
$(v)$
and axial
$(w)$
directions are shown separately.
$(b)$
The streamwise vorticity of optimal initial perturbations at
$\tau = 30, 45, 60$
.
Optimal initial condition at
$\tau =3$
using the data-driven method with varying data sizes and the model-based method with different iteration counts.
We can also note that the error plots for optimal modes in figures 6
$(b)$
and 6
$(c)$
are non-monotonic. This is reasonable because when additional column vectors are included in
$\boldsymbol{X}$
, the corresponding POD subspace may change significantly, leading to a temporary misalignment between the reduced subspace and the true dominant transient dynamics. As shown in figure 8, the azimuthal wavenumber of the computed optimal initial condition using the proposed method changes from 3 to 2 as the amount of data included in
$\boldsymbol{X}$
increases. As an analogy, in model-based analyses, the computed optimal initial condition may only exhibit the correct wavenumber and structure after sufficient iterations. Thus, the observed oscillatory convergence pattern is an expected artefact when enriching the basis.
From the contours of the optimal initial condition at
$\tau =15$
in figure 7
$(a)$
, we observe that the approximation of the axial velocity is less accurate compared with those in radial and azimuthal directions. This discrepancy arises because the axial velocity is typically over two orders of magnitude smaller than the dominant swirl velocity. Additionally, the use of POD modes, which are prioritised based on their energy content, exacerbates this effect. From the contours of the optimal outcome, we see that the initially small radial and axial velocity components grow to become dominant, a phenomenon attributed to the so-called `anti-lift-up’ mechanism (Antkowiak & Brancher Reference Antkowiak and Brancher2007). Figure 6
$(d)$
illustrates the energy evolution of optimal initial perturbations over time, which meets the optimal growth at the corresponding
$\tau$
. Therefore, the streamwise vorticity, related to the two dominant velocity components
$u$
and
$v$
, is presented in figure 7
$(b)$
and the same results have also been shown in Mao & Sherwin (Reference Mao and Sherwin2012). In contrast to the previous example, we observe that achieving mode convergence at later time instances is more challenging. This is because the out-of-core structure of the optimal initial conditions grows as the time instance increases, requiring a larger dataset to fully resolve these perturbations.
5. Correlation with DMD
The method proposed in this work shares similarities with DMD, especially in low-dimensional approximations of the linear operator, as discussed in (3.7). However, the aims and underlying assumptions differ significantly and this section is devoted to clarifying their similarities and differences by first introducing the conventional DMD framework and then discussing its potential extension to transient growth analysis, which can further demonstrate the necessity of our method.
5.1. Recalling DMD
Dynamic mode decomposition is primarily designed to capture the asymptotic behaviour of dynamical systems by approximating the system’s spectral properties. Following the standard formulation (Schmid Reference Schmid2010; Tu Reference Tu2013; Herrmann et al. Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021), DMD begins by assembling two snapshot matrices
\begin{align} \boldsymbol{X}_{\textit{DMD}}, \boldsymbol{Y}_{\textit{DMD}} \in \mathbb{C}^{N \times n_d}, \end{align}
each composed of
$n_d$
state vectors sampled from the system trajectories. These matrices are related through the discrete-time propagator
$\mathcal{M}(\Delta t)$
\begin{align} \boldsymbol{Y}_{\textit{DMD}} = \mathcal{M}(\Delta t) \boldsymbol{X}_{\textit{DMD}} = e^{\boldsymbol{A} \Delta t} \boldsymbol{X}_{\textit{DMD}}, \end{align}
where the sampling interval
$\Delta t$
is typically chosen sufficiently small to resolve all relevant temporal scales of the flow dynamics (Tu Reference Tu2013). Here,
$\boldsymbol{X}_{\textit{DMD}}$
is analogous to
$\boldsymbol{X}$
, and
$\boldsymbol{Y}_{\textit{DMD}}$
corresponds to
$\boldsymbol{Y}_{\Delta t}$
, although they are constructed from consistent trajectories typically. Similar to (3.7), a lower-dimensional surrogate
$\hat {\boldsymbol{A}}$
with rank
$r$
can then be obtained via
\begin{align} \begin{aligned} \boldsymbol{X}_{\textit{DMD}} &= \boldsymbol{ {V}}_r \boldsymbol{ {\varLambda }}_r \boldsymbol{W}^{{T}}_r, \\ e^{\hat {\boldsymbol{A}} \Delta t} &= \boldsymbol{V}^{{T}}_r \boldsymbol{Y}_{\textit{DMD}} \boldsymbol{ {W}}_r \boldsymbol{ {\varLambda }}^{-1}_r. \end{aligned} \end{align}
The eigen-decomposition of
$e^{\hat {\boldsymbol{A}} \Delta t}$
\begin{align} e^{\hat {\boldsymbol{A}} \Delta t} \hat {\boldsymbol{ {\varGamma }}}_r = \hat {\boldsymbol{ {\varGamma }}}_r \boldsymbol{ {\varXi }}_r, \end{align}
yields the DMD eigenvalues
${\xi }_i$
and corresponding modes
$\hat {\boldsymbol{ {\gamma }}}_i$
, with
\begin{align} \boldsymbol{ {\varXi }}_r = \mathrm{diag}({\xi }_1, {\xi }_2, \ldots , {\xi }_{r}), \hat {\boldsymbol{ {\varGamma }}}_r = [\hat {\boldsymbol{ {\gamma }}}_1, \hat {\boldsymbol{ {\gamma }}}_2, \ldots , \hat {\boldsymbol{ {\gamma }}}_{r}] \in \mathbb{C}^{N \times r} \text{and } i \in \{1, 2, 3, \ldots , r\}. \end{align}
The eigenvalues and eigenvectors of
$\boldsymbol{A}$
can be approximated through
$\mu _i \approx {\ln \xi _i}/{\Delta t}$
and
$\boldsymbol{ {\varGamma }}_r = \boldsymbol{ {V}}_r \hat {\boldsymbol{ {\varGamma }}}_r$
, respectively.
As discussed in Schmid (Reference Schmid2022), the system state in (2.1) can be predicted through DMD as
\begin{align} \boldsymbol{u}_{\tau } = \boldsymbol{ {\varGamma }}_r \boldsymbol{ {\varXi }}_r^{\tfrac {\tau }{\Delta t}} \boldsymbol{ {\varpi }} , \end{align}
where the vector of modal amplitudes
$\boldsymbol{ {\varpi }}$
is obtained from the initial condition using
$\boldsymbol{ {\varpi }} = \boldsymbol{ {\varGamma }}_r \backslash \boldsymbol{u}_{0}$
.
From (5.6), we can see a common but potentially not robust approach for estimating transient energy growth using conventional DMD is to define the gain as the squared maximum singular value of the weighted propagator
\begin{align} G(\tau ) = \max _{\boldsymbol{u}_0} \frac {\left \|\boldsymbol{u}_{\tau }\right \|_{\boldsymbol{Q}}^{2}}{\left \|\boldsymbol{u}_0\right \|_{\boldsymbol{Q}}^{2}} \approx \sigma _{\textit{max}}^2 (\kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}_{\textit{DMD}} \boldsymbol{ {\varXi }}_{r}^{\tfrac {\tau }{\Delta t}} \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}^{-1}_{\textit{DMD}}), \end{align}
where the weight matrix
$\kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}_{\textit{DMD}}$
is determined by Cholesky decomposition,
$\boldsymbol{ {\varGamma }}^{{T}}_{r} \boldsymbol{Q} \boldsymbol{ {\varGamma }}_{r} = \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}_{\textit{DMD}}^{{T}} \kern1.5pt \hat {\boldsymbol{\kern-1.5pt F}}_{\textit{DMD}}$
. Since this DMD process through (5.7) focuses on the eigenvalue matrix
$\boldsymbol{ {\varXi }}$
, we refer to it as the
$\boldsymbol{ {\varXi }}$
-based method, while our formulation is characterised as
$\hat {\mathcal{M}}(\tau )$
-based, emphasising its reliance on time-propagation operators.
Since
$\Delta t$
may be chosen arbitrarily large or small, one distinction between the proposed method and the
$\boldsymbol{ {\varXi }}$
-based approach lies in the additional requirement imposed on the construction of the initial matrix
$\boldsymbol{X}$
in § 3.1. This requirement is not only because POD modes are capable of capturing a larger fraction of the transient dynamics for a given data budget. More importantly, employing these modes ensures that the initial matrix is well conditioned, as the associated singular vectors of
$\boldsymbol{X}$
are of comparable magnitude, thereby reducing the impact of spurious modes that often lead to an overestimation of energy growth. This issue is examined in detail in the following section.
It should be noted that, if
$\boldsymbol{X}_{\textit{DMD}} = \boldsymbol{X}$
,
$\boldsymbol{Y}_{\textit{DMD}} = \boldsymbol{Y}_{\Delta t}$
and all SVD modes are retained in (5.3), then
$e^{\hat {\boldsymbol{A}} \Delta t}$
in (5.3) is identical to
$\hat {\mathcal{M}}(\Delta t)$
in (3.7). In this case, the
$\boldsymbol{ {\varXi }}$
-based method is equivalent to the proposed method for computing
$G(\Delta t)$
in steady flows. In this context, another potential approach is evaluating
\begin{align} G(\tau ) = \sigma _{\textit{max}}^2 (\mathcal{M}(\tau )) \approx \sigma _{\textit{max}}^2 (\hat {\mathcal{M}}^{\frac {\tau }{\Delta t}} (\Delta t)), \end{align}
which we refer to as the
$\hat {\mathcal{M}}(\Delta t)$
-based method. This approach, however, relies on the a priori assumption that
$\mathcal{M}(\tau ) \equiv \hat {\mathcal{M}}^{({ {\tau }/{\Delta t}})}(\Delta t)$
. which implicitly presumes a time-invariant evolution operator. This assumption breaks down when analysing transient energy growth in flows with time-dependent base states (Zhao, Ghidaoui & Kolyshkin Reference Zhao, Ghidaoui and Kolyshkin2007; Biau Reference Biau2016; Xu, Song & Avila Reference Xu, Song and Avila2021). In such situations, a more general and robust strategy, as adopted in the present work, is to directly construct the operator
$\boldsymbol{Y}_{\tau }$
for each time horizon
$\tau$
of interest, thereby avoiding any reliance on repeated applications of a fixed-time propagator.
5.2. Ill-condition problem in DMD
It is well known that data collected from system trajectories can be highly ill conditioned, which may lead to numerical issues (Drmač, Mezic & Mohr Reference Drmač, Mezic and Mohr2018). In DMD-related methodologies, this difficulty is commonly alleviated through a rank-
$r$
truncated SVD, as in (5.7). However, the selection of an appropriate truncation rank is non-trivial: a low-rank approximation can lose low-energy unstable modes while a high-rank approximation can be numerically unstable (Drmač & Mezić Reference Drmač and Mezić2023). Although data scaling has been proposed as a preprocessing strategy to improve conditioning (Drmač et al. Reference Drmač, Mezic and Mohr2018), the following investigation indicates that achieving convergence within the
$\boldsymbol{ {\varXi }}$
-based framework remains challenging.
$(a)$
Transient energy growth of the complex Ginzburg–Landau equation estimated using the
$\boldsymbol{ {\varXi }}$
-based approach for various truncation ranks
$r$
.
$(b)$
Contour plot of the MARE of transient energy growth for the Batchelor vortex as a function of data size
$N$
and SVD truncation rank
$r$
using the
$\boldsymbol{ {\varXi }}$
-based method.
$(c)$
Transient energy growth of the Batchelor vortex computed using the
$\boldsymbol{ {\varXi }}$
-based method with
$N=5000$
for different values of
$r$
; results obtained from the model-based approach are shown by the green curve for comparison.
To illustrate this behaviour, we apply the
$\boldsymbol{ {\varXi }}$
-based method to the complex Ginzburg–Landau equation using the same dataset as in Herrmann et al. (Reference Herrmann, Baddoo, Semaan, Brunton and McKeon2021). As shown in figure 9
$(a)$
, the accuracy of the predicted transient energy growth depends sensitively on the choice of the reduced rank
$r$
. When
$r$
is too small, the reduced model fails to capture the essential transient dynamics, whereas excessively large values of
$r$
introduce spurious modes that artificially inflate the predicted energy amplification.
For high-dimensional and strongly non-normal systems, achieving an appropriate balance between the truncation rank and the number of snapshots is particularly challenging and often entails a substantial increase in computational cost. Moreover, in the presence of energy overestimation, it becomes difficult to assess the true magnitude of transient energy growth. To illustrate this issue, we consider the Batchelor vortex case. A total of 5000 snapshots are collected from the first 25 trajectories,
$\boldsymbol{D}_{1:S}^k$
, using a sampling interval of
$\Delta t = 0.3$
, yielding 200 snapshots per trajectory. We note that employing either a smaller or a larger
$\Delta t$
leads to qualitatively similar results. These snapshots are subsequently used to assemble the scaled data matrix
$\boldsymbol{X}_{\textit{DMD}}$
, upon which the
$\boldsymbol{ {\varXi }}$
-based method is applied. As shown in figure 9(b), the mean absolute relative error (MARE) is evaluated as
\begin{align} \mathrm{MARE}=\frac {1}{N_\tau }\sum _{i=1}^{N_\tau } e_{\textit{growth}}(i \Delta t), \end{align}
where
$N_\tau$
is the number of evaluated time instances. The results of figure 9(b) indicate that ill conditioning can lead to non-convergence: across all tested configurations, the minimum MARE remains as large as
$0.097$
, which is significantly higher than that obtained using the proposed method. The best and second-best performances for
$N = 5000$
are achieved with
$r = 750$
and
$r = 3000$
, respectively, and the corresponding transient energy growth curves are shown in figure 9
$(c)$
. Notably, the
$\boldsymbol{ {\varXi }}$
-based method consistently exhibits overestimation at short time horizons, even when convergence has not yet been achieved at longer times. This behaviour persists despite increasing the amount of data.
We also notice that several recent studies have proposed strategies to improve spectral accuracy (Williams et al. Reference Williams, Kevrekidis and Rowley2015; Colbrook & Townsend Reference Colbrook and Townsend2024). Nevertheless, recovering the transient properties of
$\boldsymbol{A}$
remains challenging. It is difficult to distinguish the spurious modes which cause overestimation of transient energy growth because they usually lie in continuous spectrum and DMD is usually robust in approximating part of discrete spectra. It is also generally unclear which modes obtained from the
$\boldsymbol{\varXi }$
-based approach contribute most to transient growth. Filtering out any mode may lead to an underestimation of energy amplification, as demonstrated in the swept Hiemenz flow by Obrist & Schmid (Reference Obrist and Schmid2003) and the vortex flow by Mao & Sherwin (Reference Mao and Sherwin2012). In our study, the number of spurious modes decreases as the condition number of the initial matrix is reduced. In particular, when the condition number of
$\boldsymbol{X}_{\textit{DMD}}$
approaches unity, the eigenvalues (including discrete and continuous parts) approximated by DMD tend to converge towards those obtained using model-based approaches under the same spatial discretisation. This behaviour can be interpreted within the framework of Krylov subspace theory, where the orthogonal property of
$\boldsymbol{X}_{\textit{DMD}}$
reduces the residuals during low-dimensional projection, and it also suggests the potential for developing more efficient eigenvalue computation strategies. A systematic investigation of the spectral properties of reduced-order models is thus of considerable interest, but lies beyond the scope of the present work and will be pursued in future studies.
In the view of the ill-conditioning issues inherent to DMD-related methods, it is advantageous to enforce the condition number of
$\boldsymbol{X}$
to be close to unity, which effectively requires the column vectors of the data matrix to be as nearly orthogonal as possible. At the same time, the right-singular vector matrix
$\boldsymbol{V}_r$
in (5.3) should retain as much dynamical information as possible, since it spans the subspace onto which the evolution operator is projected. These considerations naturally motivate the use of POD modes as a preferred choice for constructing the initial matrix
$\boldsymbol{X}$
. An additional benefit of this construction is that it largely removes sensitivity to other algorithmic parameters commonly encountered in data-driven methods, such as the data sampling frequency and the truncation rank employed in the SVD. It should be noted that, although a rank-
$R$
truncated SVD of
$\boldsymbol{D}^k_{1:S}$
is adopted in the proposed method for computational efficiency, employing a full-rank SVD would not introduce ill conditioning in subsequent computations, owing to the orthogonality of the column vectors in
$\boldsymbol{X}$
. Consequently, within the proposed framework, the only parameter governing convergence is the data size
$N$
of the matrix
$\boldsymbol{X}$
.
6. Summary
In this work, we propose a data-driven methodology for transient growth analyses of linear dynamical systems and validate it through three case studies: the linearised complex Ginzburg–Landau equation, the flow over a backward-facing step and the Batchelor vortex. The data-driven method utilises Hermite polynomials to gather transient information over short time intervals, followed by POD to extract dominant spatial modes from snapshots. By projecting the forward evolution operator onto a reduced space spanned by these modes, the data-driven method accurately captures the spatial distribution of the optimal initial perturbations, with their energy evolutions aligning with the model-based optimal growth. The proposed algorithm also demonstrates robustness in avoiding the influence of spurious global modes, which are often encountered in DMD-based methods.
Furthermore, this study shows that optimal modes can also be interpreted as a combination of POD modes, which are chosen for the construction of the low-dimensional space in this method. In addition, the data-driven approach converges consistently as long as more data are included in the initial matrix
$\boldsymbol{X}$
and outcome matrices
$\boldsymbol{Y}_{\tau }$
. The convergence behaviour also highlights the potential for more efficient strategies to further reduce the computational costs.
A clear advantage of the proposed data-driven approach is that it eliminates the adjoint equations, relying solely on the forward operator/code to perform the analysis. This adjoint-free framework greatly simplifies the implementation and reduces the computational complexity. Contrasted with model-based methods, which require renewed forward and backward integrations for each time instance, the data-driven approach collects the snapshots at arbitrary time instance
$t \in [0,\tau ]$
during a single forward integration for
$\boldsymbol{Y}_{\tau }$
, facilitating the computation of transient growth over multiple time instances. On the other hand, the well-conditioned dataset ensures consistent convergence in computation without parameter selections. Last but not least, this implementation also ensures applicability to systems with time-dependent base flows.
While the present study focuses on linear systems, accounting for nonlinear effects, such as perturbation amplitudes in systems governed by nonlinear equations, represents an interesting direction (Kerswell Reference Kerswell2018). Moreover, applying this method to experimental data that may contain noise could be a valuable avenue for future research. Further development of more robust strategies to mitigate ill conditioning would also be valuable for obtaining more reliable spectral information. We also notice that the majority of wall time is dominated by I/O (reading/writing snapshots). Scalability for higher-dimensional problems can be improved by using parallelised I/O and GPU-accelerated linear algebra libraries.
Funding
The authors would like to acknowledge the support from the National Natural Science Foundation of China (Grant Nos. 92152109).
Declaration of interests
The authors report no conflicts of interest.
Appendix A. Hermite polynomials
The 1-D Hermite polynomials satisfy the differential equation
\begin{align} \frac {d^2 h_{\rho }(x)}{{\rm d}x^2} - \delta x \frac {{\rm d} h_{\rho }(x)}{{\rm d}x} + \delta n h_{\rho }(x) = 0, \end{align}
where
$h_{\rho }(x)$
represents
$\rho$
th polynomial, which can also be defined by the relations
\begin{align} \begin{aligned} h_1(x) &= 1, \\ h_2(x) &= \delta x, \\ h_{m+2}(x) &= \delta x h_{m+1}(x) - m h_{m}(x) \quad \text{ for } m = 1, 2, 3, \ldots . \end{aligned} \end{align}
They satisfy the orthogonality with respect to the Gaussian weight
$e^{-( {\delta }/{2}) x^{2}}$
, namely
\begin{align} \int _{-\infty }^{\infty } h_m(x) h_n(x) e^{-\delta x^{2}} \mathrm{d} x=0 \quad \text{ for } m \neq n, \end{align}
where
$\delta \gt 0$
controls the decay rate of the Gaussian envelope. Those polynomials are complete in the Hilbert space, making them suitable for representing localised initial perturbations in unbounded domains. For enhanced adaptability within the computational domains, we employ Gauss-weighted Hermite functions with
$\delta = 0.24$
, given by
\begin{align} \varPi _{\rho } = h_{\rho }(x)\,e^{-\tfrac {\delta }{2} x^{2}}. \end{align}
The first
$9$
such functions are depicted in figure 10
$(a)$
.
Hermite polynomials.
$(a)$
The first
$9$
Gauss-weighted Hermite functions.
$(b)$
The first
$5$
ranks of 2-D Hermite polynomials in (left) Cartesian form and (right) polar form (Victor et al. Reference Victor, Mechler, Repucci, Purpura and Sharpee2006).
The 2-D Hermite polynomials can be constructed as products of 1-D Hermite functions in the
$x$
and
$y$
directions (Victor et al. Reference Victor, Mechler, Repucci, Purpura and Sharpee2006), satisfying the differential equation
\begin{align} \left ( \frac {\partial }{\partial x} - \frac {\partial }{\partial y} \right )^{\alpha +\beta } H_{\alpha ,\beta }(x,y) = 0. \end{align}
In Cartesian coordinates, the polynomials of rank
$\rho =\alpha +\beta$
, where
$\rho \geqslant \alpha , \beta \geqslant 0$
, are given by
\begin{align} H_{\alpha , \beta }(x, y)=\frac {\varUpsilon }{C} h_\alpha \left (\frac {x}{C}\right ) h_\beta \left (\frac {y}{C}\right ) \exp \left (-\frac {x^2+y^2}{4 C^2}\right ), \end{align}
where
$h_\alpha (u)$
and
$h_\beta (u)$
are 1-D Hermite polynomials. Here,
$C$
represents the scale of the filter set and
$\varUpsilon$
is a normalisation constant. When expressed in polar form, the 2-D Hermite polynomials of the rank
${\rho } = \alpha + 2\beta$
, where
${\rho } \geqslant \alpha \geqslant 0$
and
$\alpha$
,
$\beta$
are positive integers, are structured into ‘cosine’ and ‘sine’ pairs
\begin{align} \left \{\begin{aligned} H_{\alpha , \beta }^{\cos }(r, \theta )=\frac {\varUpsilon }{C} \cos (\alpha \theta )\left (\frac {r}{C}\right )^{\alpha } F_{\alpha , \beta }\left (\frac {r^2}{C^2}\right ) \exp \left (-\frac {r^2}{4 C^2}\right ), \\ H_{\alpha , \beta }^{\sin }(r, \theta )=\frac {\varUpsilon }{C} \sin (\alpha \theta )\left (\frac {r}{C}\right )^{\alpha } F_{\alpha , \beta }\left (\frac {r^2}{C^2}\right ) \exp \left (-\frac {r^2}{4 C^2}\right ), \end{aligned}\right . \end{align}
with
\begin{align} F_{p_c, p_s}(a)=\sum _{i=0}^{p_s}(-2)^{{p_s}-i} \frac {({p_c}+{p_s})!{p_s}!}{(p_c+i)!i!({p_s}-i)!} a^i. \end{align}
The first
$5$
ranks of 2-D Hermite polynomials in both polar and Cartesian coordinates are shown in figure 10
$(b)$
.
