2.1. The linear dispersion relation as nonlinear correlation

One of the fundamental reasons that Galerkin models of advection-dominated flows tend to be fragile is that they introduce additional variables that do not correspond to physical degrees of freedom. This is perhaps illustrated most clearly by the linear advection equation on a periodic domain

(2.1)\begin{equation} u_t + c u_x = 0, \quad x \in (0, L)\end{equation}

where u is a scalar amplitude, x is the spatial variable, and c is the constant advection speed on a finite domain of length L. For any initial condition $u(x, 0) = u_0(x)$, this equation has the simple travelling wave solution $u(x, t) = u_0(x - ct)$. Given the initial condition, the only effective degree of freedom is the phase $ct/L$. However, the problem could also be solved by means of a Fourier expansion

(2.2)\begin{equation} u(x, t) = \sum_{n={-}\infty}^{\infty} a_n(t) \textrm{e}^{\textrm{i}k_nx} \end{equation}

with $k_n = 2{\rm \pi} n/L$ with time-varying modal coefficients $a_n(t)$. The Galerkin system in this orthogonal basis (see § 5) is

(2.3a,b)\begin{equation} \dot{a}_n(t) ={-}\textrm{i}\omega_n a_n(t), \quad \omega_n = k_n c n = 0, \pm 1, \pm 2, \dots \end{equation}

The relationship between frequency $\omega$ and wavenumber $k_n$ is the dispersion relation; in this case it implies that all scales are carried at the same speed $c$.

With this analytic dispersion relation, (2.3a,b) is equivalent to the travelling wave solution, since $a_n(t) = \textrm {e}^{-\textrm {i}\omega _n t} a_n(0)$

(2.4)\begin{equation} u(x, t) = \sum_{n={-}\infty}^{\infty} a_n(0) \exp({\textrm{i}k_n(x-ct)}) = u_0(x - ct). \end{equation}

However, the Galerkin model has introduced many degrees of freedom in the harmonics $a_n$ by artificially separating space and time. If the projection is approximated numerically or with empirical basis modes, the estimated system may include some error, so that $\omega _n = k_n c + \epsilon _n$. In this case the Galerkin system will be dispersive, i.e. each wavenumber will propagate with a slightly different speed. The travelling wave solution will tend to lose coherence on a time scale $1/\epsilon$, as shown in figure 2.

Figure 2. Linear advection equation with errors $\epsilon _n \sim \mathcal {N}(0, \epsilon ^2)$ in the dispersion relation $\omega _n = c k_n$. The Galerkin model (grey) loses coherence with the exact solution (black) over a time scale $1/\epsilon$. If the polynomial correlations implied by the dispersion relation are enforced explicitly, the model is robust to such errors. Nonlinear correlation in the true system, given by (2.5), appears in the Lissajous-type phase portraits of the Fourier coefficients (b–d). Similar behaviour manifests in Galerkin models of nonlinear advection-dominated flows.

An alternative perspective on the dispersion relation is that it specifies nonlinear correlations between the temporal coefficients $a_n$, removing the spurious degrees of freedom introduced by Galerkin projection. The linear dispersion relation $\omega _n = n k_1 c$ implies the nonlinear relationship for harmonics

(2.5)\begin{equation} a_n(t) = \exp({-\textrm{i}n k_1 c t}) a_n(0) \propto a_1^n, \end{equation}

with the proportionality determined by the initial condition. Then the only degree of freedom is $a_1$, and the travelling wave solution is recovered by the Galerkin model projected onto this mode. In dynamical systems terminology, the solution is restricted to a one-dimensional manifold: a circle representing the phase of the leading Fourier coefficient. In this case the decoherence does not lead to instability because the system is purely linear with purely imaginary eigenvalues, but in nonlinear systems with non-zero linear growth rates the departure from the solution manifold can be catastrophic.

2.2. Triadic interactions and the energy cascade

For more general linear systems the preceding analysis is complicated by non-normality and physical dispersion, and the concept of a dispersion relation is not well defined for nonlinear dynamical systems. Nevertheless, analogous concepts are similarly important in models of nonlinear PDEs. For example, Burgers’ model is a paradigmatic scalar conservation equation illustrating many features of gas dynamics and nonlinear flows more broadly. Burgers’ equation with viscosity $\epsilon$ is

(2.6)\begin{equation} \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \epsilon \frac{\partial^2u}{\partial x^2}, \quad x \in (0, 2{\rm \pi}). \end{equation}

On a periodic domain, we can apply the same Fourier expansion (2.2) with $L=2{\rm \pi}$, leading to the Galerkin ordinary differential equation (ODE) system

(2.7)\begin{equation} \dot{a}_k ={-}\epsilon k^2 a_k -\textrm{i} k \sum_{\ell={-}\infty}^\infty a_\ell a_{k-\ell}, \quad k = 0, \pm 1, \pm 2, \dots. \end{equation}

The two right-hand side terms in (2.7), originating from the viscous and nonlinear PDE terms respectively, capture several key features of the full Navier–Stokes equations. First, the convolution-type sum over wavenumbers $\ell$ includes only pairs that sum to $k$; these are the so-called ‘triadic’-scale interactions. Second, it can be shown that the nonlinear term is energy preserving in the sense that when the energy $a_k a_k^*$ is summed over all wavenumbers the nonlinear term does not contribute to a net change in energy of the system. (A similar result holds for inhomogeneous flows (Schlegel & Noack Reference Schlegel and Noack2015).) This suggests that the only role of nonlinearity is to transfer energy between scales. Meanwhile, the dissipation rate of each mode scales quadratically with wavenumber so that the bulk of dissipation occurs at the smallest scales.

The overall picture of the dynamics in the spectral domain is therefore that the nonlinear term transfers energy from the more energetic large scales to the dissipative small scales. Since (2.7) is very similar to the spectral form of the momentum equations for isotropic turbulence (Tennekes & Lumley Reference Tennekes and Lumley1972), this ‘energy cascade’ is an important feature of real viscous flows as well. The energy cascade points to another often-discussed issue with Galerkin models: if the system is truncated at a wavenumber $r$ which is not sufficiently large to capture the net dissipation rate, the energy cascade is interrupted and the system of ODEs will overestimate the energy, potentially even becoming unstable.

This issue is fundamentally different from the decoherence discussed in the context of the linear advection equation. For example, the issue of fine-scale dissipation is also present in the heat equation, given by (2.6)–(2.7) without the convective nonlinearity. Whereas the Fourier–Galerkin representation of advection introduces spurious degrees of freedom, this discussion suggests that in the representation of the heat equation all coefficients are dynamically important (self-similarity notwithstanding). The Galerkin system is therefore an ideal representation of the parabolic dynamics of the heat equation, where the fundamental assumption of separation of variables is valid.

The inability of Fourier decomposition and POD to produce efficient representations of travelling wave physics has long been recognized. Fundamentally, these decompositions rely on a space–time separation of variables, which is not a valid assumption for travelling waves. Many extensions to POD have been developed for translationally invariant systems and systems with other symmetries (Rowley & Marsden Reference Rowley and Marsden2000; Reiss et al. Reference Reiss, Schulze, Sesterhenn and Mehrmann2018; Rim, Moe & LeVeque Reference Rim, Moe and LeVeque2018; Mendible et al. Reference Mendible, Brunton, Aravkin, Lowrie and Kutz2020).

For general viscous, nonlinear, advection-driven fluid flows, we might expect advection, triadic interactions and small-scale dissipation to all be relevant as a result of the joint hyperbolic–parabolic structure of the Navier–Stokes equations. The intrinsic dimensionality of the system and, conversely, the inaptitude of the Galerkin model, may not be a priori clear as a result of a complex interplay between these mechanisms.

For example, if the leading degree of freedom $a_1$ tends to oscillate at a frequency $\omega _0$, representing either a standing or travelling wave, then the $a_2$ dynamics includes a term of the form $a_1^2 \sim \textrm {e}^{2\textrm {i}\omega _0 t}$. Similarly, $\dot {a}_3 \propto a_1 a_2 \sim \textrm {e}^{3\textrm {i}\omega _0 t}$. In the energy cascade picture, these higher-order modes act as forced, damped oscillators and will tend to respond at the forcing frequencies. In this manner, triadic interactions in the wavenumber domain can also give rise to nonlinear correlations in time and triadic structure in the frequency domain.

This effect is not necessarily limited to systems with spatial or temporal periodicity; Rubini et al. (Reference Rubini, Lasagna and Da Ronch2020) investigated the application of system identification methods to a chaotic lid-driven cavity flow and showed that sparse nonlinear coupling, analogous to triadic interactions, were critical for resolving energy transfers across scales of the flow. They distinguish between this empirical, a posteriori sparsity appearing in the statistics and data-driven models, and the structural, a priori sparsity of systems such as (2.7), which is generally lost in Galerkin models of inhomogeneous flows. Similarly, Schmidt (Reference Schmidt2020) recently proposed a bispectral modal analysis technique that leverages approximate sparsity in frequency interactions.

In analogy with the dispersion relation, these processes may result in latent structure not immediately obvious in the Galerkin representation. If this structure is ignored, the behaviour of the model may depart significantly from that of the underlying system. For example, Majda & Timofeyev (Reference Majda and Timofeyev2000) showed that a truncated Galerkin model of the inviscid Burgers equation tends towards equipartition of energy rather than a physical solution as a result of a catastrophic decoherence mechanism. Consequentially, in the following sections we argue that nonlinear correlation and manifold restriction plays an important role in the stability and accuracy of reduced-order models of advection-dominated fluid flows.

4.1. Proper orthogonal decomposition

One of the most widely used techniques for dimensionality reduction and modal analysis is POD, which solves the optimization problem

(4.2)\begin{equation} \textrm{minimize}~_{\{ \boldsymbol{\psi}_k \} } \left \langle \left\|{ \boldsymbol{u}' - \sum_{k=1}^r \boldsymbol{\psi}_k \left( \boldsymbol{u}' , \boldsymbol{\psi}_k \right)_{\varOmega}} \right\|^2 \right \rangle \quad \textrm{subject to}~ \left( \boldsymbol{\psi}_j, \boldsymbol{\psi}_k \right)_{\varOmega} = \delta_{jk} \end{equation}

in the norm induced by the energy inner product

(4.3)\begin{equation} \left( \boldsymbol{u}, \boldsymbol{v} \right)_{\varOmega} \equiv \int_\varOmega \boldsymbol{u}(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{v}^*(\boldsymbol{x}) \, {\text d} \varOmega, \end{equation}

where $\langle \cdot \rangle$ is an ensemble average, approximated in practice by a time average, $\delta _{jk}$ is the Kronecker delta and the star indicates a complex conjugate. For data on a non-uniform mesh, the inner product is computed with a weighted Riemann sum, approximating ${\text d} \varOmega$ with the mass matrix of the discretization. Thus, the objective is to minimize the residual energy in a linear subspace of $r$ orthonormal modes, providing an optimal low-rank representation of the flow.

This problem can be solved with the calculus of variations, leading to the result that the modes $\{ \boldsymbol {\psi }_k \}$ are eigenfunctions of the correlation tensor $\mathcal {\boldsymbol C}(\boldsymbol {x}, \boldsymbol {x}')$:

(4.4)\begin{equation} \int_\varOmega \mathcal{\boldsymbol C}(\boldsymbol{x}, \boldsymbol{x}') \boldsymbol{\psi}_k(\boldsymbol{x}') \, {\text d} \varOmega' = \sigma_k^2 \boldsymbol{\psi}_k(\boldsymbol{x}), \end{equation}

where $\mathcal {\boldsymbol C}(\boldsymbol {x}, \boldsymbol {x}') = \langle \boldsymbol {u}(\boldsymbol {x}, t) \boldsymbol {u}^*(\boldsymbol {x}', t) \rangle$ and $\{ \sigma _k \}$ are the POD eigenvalues, representing the average fluctuation kinetic energy captured by each mode. The coefficients $\boldsymbol {a}(t)$ can be extracted with the projection $a_k = \left ( \boldsymbol {u}', \boldsymbol {\psi }_k \right )_{\varOmega }$. In practice the correlation tensor is often not feasible to construct, since it scales with the square of the discretized state dimension. Instead it is approximated numerically with either a singular value decomposition (SVD) or the snapshot method (Sirovich Reference Sirovich1987; Holmes et al. Reference Holmes, Lumley and Berkooz1996; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017). In this work we use the snapshot method as implemented in the modred library (Belson, Tu & Rowley Reference Belson, Tu and Rowley2014), since it does not require storing the entire time series of high-dimensional discretized velocity fields in memory.

The method of snapshots is based on simple linear algebraic manipulations of the discretized form of the eigenvalue problem (4.4). We omit a derivation here as it is given in standard references, e.g. Holmes et al. (Reference Holmes, Lumley and Berkooz1996). Rather than form the spatial correlation tensor $\mathcal {\boldsymbol C}(\boldsymbol {x}, \boldsymbol {x}')$, we compute a temporal correlation matrix $\boldsymbol{\mathsf{R}}$ with entries defined by

(4.5)\begin{equation} R_{jk} = \frac{1}{M} \left( \boldsymbol{u}(t_j), \boldsymbol{u}(t_k) \right)_{\varOmega}, \quad j, k = 1, 2, \dots, M. \end{equation}

The temporal correlation matrix $\boldsymbol{\mathsf{R}}$ has dimensions $M \times M$, and is typically much smaller than the discretized spatial correlation tensor. The eigenvalues of $\boldsymbol{\mathsf{R}}$ approximate those of $\mathcal {\boldsymbol C}$, and the modes that solve the discretized form of (4.4) are also given by

(4.6)\begin{equation} \boldsymbol{\psi}_k = \frac{1}{\sigma_k \sqrt{M}} \sum_{j=1}^M \boldsymbol{u}(t_j) W_{jk}, \end{equation}

where $\boldsymbol{\mathsf{R}}\boldsymbol{\mathsf{W}} = \boldsymbol{\mathsf{W}} \boldsymbol {\varSigma }^2$ is the eigendecomposition of $\boldsymbol{\mathsf{R}}$.

The POD has the following useful properties:

1. (i) The spatial modes form an orthonormal set: $\left ( \boldsymbol {\psi }_j, \boldsymbol {\psi }_k \right )_{\varOmega } = \delta _{jk}$.

2. (ii) The temporal coefficients are linearly uncorrelated: $\langle a_j a_k \rangle = \sigma _k^2 \delta _{jk}$.

3. (iii) The modes can be ranked hierarchically by average energy content $\sigma _k^2$.

Since POD can be viewed as a continuous form of the SVD, these properties are analogous to unitarity of the matrices of left and right singular vectors. The singular values quantify the statistical variance captured by the low-rank SVD approximation. As a consequence of the hierarchical ordering, the POD can be computed without a priori specification of the rank $r$, with the truncation determined later by a threshold based on the residual energy. Finally, since the modes are a linear combination of DNS snapshots, the reconstruction (4.1) automatically satisfies the incompressibility constraint and boundary conditions.

We compute the POD from 4000 fields sampled at $\Delta t = 0.05$, approximately ten times the shear layer frequency, using the method of snapshots (Sirovich Reference Sirovich1987). This is around 13 % of the total number of fields retained from the DNS, but is sufficient for statistical convergence of the leading modes. The singular value spectrum and residual energy are shown in figure 5. The singular values converge relatively quickly; the first pair of modes contain 70 % of the fluctuation kinetic energy, the first six account for ${\sim }90\,\%$, and by $r=64$ approximately $99.97\,\%$ of the energy is recovered. We retain 64 modes for further analysis and note that our modelling results are insensitive to moderate changes in truncation.

Figure 5. Singular value spectrum of the quasiperiodic cavity flow. Black dots represent the normalized squared singular values of the snapshot correlation matrix, indicating the fraction of fluctuation kinetic energy resolved by each mode. Red crosses indicate the fraction of residual energy, or normalized cumulative sum of squared singular values. Dashed lines indicate the number of modes retained ($r=64$).

Still, as we will show in § 6, the intrinsic dimensionality of the system is much smaller than that of the linear subspace required for reconstruction. As with the advection system in § 2, this is partly due to the representation of travelling waves, as shown in figure 6. This is made clearer by a DMD analysis.

Figure 6. Harmonic modes identified from POD and DMD analysis. The spatial fields and phase portraits both indicate that certain mode pairs are harmonics arising from the description of wavelike motion in the shear layer and inner cavity. Because DMD is based on both spatial and temporal correlation, this structure is especially pronounced in the DMD coefficients. The vorticity plots are real parts of the DMD modes, but analogous modes exist in the POD basis.

4.2. Dynamic mode decomposition

Although POD is guaranteed to provide an energy-optimal spatial reconstruction of the flow field, it sacrifices all temporal information in the computation of the correlation tensor. The POD basis is therefore purely kinematic and contains no dynamic information. An alternative approach is to compute the discrete Fourier transform of the fields, which suffers from the opposite issue: frequency information is perfectly resolved, but the result is not necessarily associated with a useful reduced-order linear subspace for kinematic representation. DMD, introduced by Schmid (Reference Schmid2010), is a useful compromise between these extremes.

DMD seeks to approximate a discrete-time linear evolution operator defined by

(4.7)\begin{equation} \textrm{minimize}~_{\mathcal{\boldsymbol A}} \left \langle \| \boldsymbol{u}(\boldsymbol{x}, t_{n+1}) - \mathcal{\boldsymbol A} \boldsymbol{u}(\boldsymbol{x}, t_{n}) \|^2 \right \rangle. \end{equation}

The description of nonlinear dynamics in terms a linear evolution operator acting on observables has a deep connection to Koopman theory (Rowley et al. Reference Rowley, Mezić, Bagheri and Schlatter2009b; Mezić Reference Mezić2013; Brunton et al. Reference Brunton, Budišić, Kaiser and Kutz2022). We will discuss this in § 8, but in terms of modal analysis it is more useful to think of DMD as solving an alternative optimization to (4.2).

Given a series of snapshots, a least-squares solution to (4.7) could be found in terms of the pseudoinverse of the snapshot matrix. However, this calculation is typically computationally prohibitive, ill conditioned and disregards low-dimensional structure in the flow. Instead, DMD seeks to approximate the spectral properties of the operator $\mathcal {\boldsymbol A}$ without explicitly forming it. There are a variety of algorithms to compute DMD in conditions with limited or noisy data, but beginning with high-fidelity DNS snapshots we follow the simple exact DMD algorithm introduced by Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014).

Beginning with a truncated POD basis and associated coefficients $\{\boldsymbol {a}(t_n) \}_{n=1}^M$, the coefficients are arranged into time-shifted matrices $\boldsymbol{\mathsf{X}} = \left [\boldsymbol {a}(t_1) \enspace \boldsymbol {a}(t_2)\enspace \cdots \enspace \boldsymbol {a}(t_M-1)\right ]$ and $\boldsymbol{\mathsf{X}}^\prime = \left [\boldsymbol {a}(t_2) \enspace \boldsymbol {a}(t_3)\enspace \cdots \enspace \boldsymbol {a}(t_M) \right ]$. Here, we assume the $\{t_n\}$ are evenly sampled in time, but it is possible to account for situations where this is not the case. A least-squares solution to (4.7) in the POD subspace is $\tilde {\boldsymbol{\mathsf{A}}} = \boldsymbol{\mathsf{X}}^\prime \boldsymbol{\mathsf{X}}^+$, where $\boldsymbol{\mathsf{X}}^+$ is the pseudoinverse. With some assumptions, the spectrum of $\mathcal {\boldsymbol A}$ is approximated by the spectrum of $\tilde {\boldsymbol{\mathsf{A}}}$, which can now be easily computed via an eigendecomposition

(4.8)\begin{equation} \tilde{\boldsymbol{\mathsf{A}}} = \boldsymbol{\mathsf{V}} \boldsymbol{\varLambda} \boldsymbol{\mathsf{V}}^{{-}1}, \end{equation}

with $\tilde {\boldsymbol{\mathsf{A}}}, \boldsymbol{\mathsf{V}}, \boldsymbol {\varLambda } \in \mathbb {C}^{r \times r}$. For details on theory and algorithms of DMD, see Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014) and Kutz et al. (Reference Kutz, Brunton, Brunton and Proctor2016). Based on this eigendecomposition, complex-valued DMD modes $\{ \boldsymbol {\phi }_k(\boldsymbol {x}) \}$ and associated projection coefficients $\boldsymbol {\alpha }(t)$ are linear combinations of the POD modes and coefficients, given by

(4.9a,b)\begin{equation} \boldsymbol{\phi}_k(\boldsymbol{x}) = \sum_{j=1}^r \boldsymbol{\psi}_j(\boldsymbol{x}) V_{jk} \quad \boldsymbol{\alpha}(t) = \boldsymbol{\mathsf{V}}^{{-}1} \boldsymbol{a}(t). \end{equation}

In principle the approximate time evolution is specified by the DMD eigenvalues $\{\lambda _k\} = \mathrm {diag}(\boldsymbol {\varLambda })$, but in terms of reduced-order modelling the decomposition can also be viewed as an alternative expansion to (4.1)

(4.10)\begin{equation} \boldsymbol{u}(\boldsymbol{x}, t) \simeq \bar{\boldsymbol{u}}(\boldsymbol{x}) + \sum_{k=1}^r \boldsymbol{\phi}_k (\boldsymbol{x}) \alpha_k(t). \end{equation}

The DMD coefficient vector $\boldsymbol {\alpha }(t)$ may then be modelled as a time series, as with the POD coefficients $\boldsymbol {a}(t)$.

This representation is essentially a similarity transformation of the POD basis; the two encode the same information and span the same subspace. However, the time dependence in the optimization problem leads DMD to transform the POD basis to modes that tend to have similar frequency content. In terms of the present analysis, figure 6 illustrates the practical relevance of this. Whereas POD happens to identify modes that are roughly coherent in time by coincidence only, the DMD modes are closer to pure harmonics. This perspective on DMD also explains the approximately discrete peaks in figure 4; each DMD eigenvalue (shown in figure 7) can be identified with some combination of the fundamental frequencies of modes $k=1$ and $k=5$.

Figure 7. DMD frequencies $\omega _k$ and average energy $E_k = \langle |\alpha _k|^2 \rangle$ along with vorticity plots for the real part of the most energetic modes. The second mode pair ($k=3, 4$) is a harmonic of the leading pair ($k=1,2$), while the third pair ($k=5, 6$) represents the low-frequency inner-cavity motion. Other modes (e.g. $k=7, 8$) are either harmonics or indicate nonlinear frequency cross-talk between these leading modes, as in figure 4.

Based on the results in §§ 6 and 7, we hypothesize that DMD also filters frequency content and accentuates nonlinear correlations for modes that are not pure harmonics. This approximate nonlinear algebraic dependence clearly indicates a manifold structure of much lower dimensionality than the linear subspace. Given these implications, the results presented below are based on the DMD expansion (4.10).

5.1. POD–Galerkin modelling

In projection-based modelling, the discretized governing equations are projected onto an appropriate modal basis. For simple geometries, this might be done analytically, as for the periodic problems in § 2 and in Noack & Eckelmann (Reference Noack and Eckelmann1994), for instance. Although general and expressive, this approach becomes challenging on complex domains and does not take advantage of structure in the solutions to the particular PDE. As a result, it is increasingly common to project onto an empirical basis, such as POD modes. Assuming that the flow is statistically stationary and the ensemble is sufficiently resolved, this provides optimal kinematic resolution in an orthonormal basis. The following Galerkin projection procedure then leads to a minimum-residual system of ODEs in this basis.

Let $\mathcal {\boldsymbol N}[\boldsymbol {u}] = 0$ be the Navier–Stokes equations in implicit form. By approximating the flow field with a truncated linear combination of basis functions as in (4.1), we expect some residual error $\boldsymbol {r}(x, t)$ in the approximated dynamics defined by

(5.1)\begin{equation} \boldsymbol{r}(\boldsymbol{x}, t) = \mathcal{\boldsymbol N}\left[ \bar{\boldsymbol{u}}(\boldsymbol{x}) + \sum_{\ell = 1}^r \boldsymbol{\psi}_\ell(\boldsymbol{x}) a_\ell(t) \right]. \end{equation}

In order to minimize the residual in this basis, the Galerkin projection condition is that that the residual be orthogonal to each mode

(5.2)\begin{equation} 0 = \left( \boldsymbol{r}, \boldsymbol{\psi}_k \right)_{\varOmega} = \left( \mathcal{\boldsymbol N}\left[ \bar{\boldsymbol{u}}(\boldsymbol{x}) + \sum_{\ell = 1}^r \boldsymbol{\psi}_\ell(\boldsymbol{x}) a_\ell(t) \right], \boldsymbol{\psi}_k \right)_{\varOmega}. \end{equation}

This leads to the linear-quadratic system of ODEs (Holmes et al. Reference Holmes, Lumley and Berkooz1996; Noack et al. Reference Noack, Morzynski and Tadmor2011)

(5.3)\begin{equation} \dot{a}_k = C_k + \sum_{\ell = 1}^r L_{k \ell} a_\ell + \sum_{\ell = 1}^r \sum_{m = 1}^r Q_{k \ell m} a_\ell a_m, \end{equation}

with constant, linear and quadratic terms given by

(5.4a)\begin{gather} C_k = \left( -\boldsymbol{\nabla} \boldsymbol{\cdot} (\bar{\boldsymbol{u}} \otimes \bar{\boldsymbol{u}} ) + \frac{1}{\textit{Re}} \nabla^2 \bar{\boldsymbol{u}} , \boldsymbol{\psi}_k \right)_{\varOmega} \end{gather}

(5.4b)\begin{gather}L_{k \ell} = \left( -\boldsymbol{\nabla} \boldsymbol{\cdot} (\bar{\boldsymbol{u}} \otimes \boldsymbol{\psi}_\ell + \boldsymbol{\psi}_\ell \otimes \bar{\boldsymbol{u} } ) + \frac{1}{\textit{Re}} \nabla^2 \boldsymbol{\psi}_\ell , \boldsymbol{\psi}_k \right)_{\varOmega} \end{gather}

(5.4c)\begin{gather}Q_{k \ell m} = \left( -\boldsymbol{\nabla}\boldsymbol{\cdot} (\boldsymbol{\psi}_m \otimes \boldsymbol{\psi}_\ell + \boldsymbol{\psi}_\ell \otimes \boldsymbol{\psi}_m ) , \boldsymbol{\psi}_k \right)_{\varOmega}. \end{gather}

Note that the constant term vanishes if the flow is expanded about a steady-state solution of the governing equations. Since the mean flow in this case is not a solution, this term represents important mean-flow forcing and is not negligible. Here, we have also neglected the pressure term, though including it does not significantly change any of the results; for detailed discussion of this point see Noack, Papas & Monkewtiz (Reference Noack, Papas and Monkewtiz2005).

In principle, we might expect that the POD–Galerkin system (5.3) leads to approximate solutions with comparable accuracy to the resolution of the expansion basis. However, for reasons introduced in § 2, the long-time behaviour of the reduced-order model may deviate significantly from that of the underlying physical system. In particular, solutions of the model are not constrained to lie on an invariant manifold of the flow. For instance, coefficients associated with shear layer or inner-cavity harmonics evolve independently from the fundamental modes, eventually leading to an unphysical loss of coherence. This effect cannot necessarily be attributed to any particular mode or interaction term, but is related to the structure of the model itself when the POD–Galerkin method is applied to advection-dominated flows.

Figure 8 shows the evolution of the fluctuation kinetic energy as predicted by the POD–Galerkin system for various levels of truncation $r$. Although the estimate does tend to improve with increasing $r$, none of these models capture the quasiperiodic dynamics of the flow, and most exhibit significant instability. This is true despite (and, we will argue, because of) the fact that these models have many more kinematic degrees of freedom than the true dynamics underlying the post-transient cavity flow.

Figure 8. Evolution of the fluctuation kinetic energy predicted by POD–Galerkin reduced-order models of various dimensions along with DNS values. Though all values of $r$ shown here capture sufficient dissipation to remain at finite energy, none resolves the true quasiperiodic dynamics.

Finally, a model similar to (5.3) may also be derived beginning with the DMD expansion (4.10). In this case, since the DMD basis is not orthonormal, the Galerkin projection must be replaced with the more general oblique projection of Petrov–Galerkin methods (Benner et al. Reference Benner, Gugercin and Willcox2015). We implement this with a coordinate transform of the standard POD–Galerkin model based on the DMD eigenvector matrix $\boldsymbol{\mathsf{V}}$, i.e. $\boldsymbol {a} = \boldsymbol{\mathsf{V}} \boldsymbol {\alpha }$. For instance, a DMD–Petrov–Galerkin model with the same truncation as the original POD model has dynamics

(5.5)\begin{equation} \dot{\alpha}_k = \tilde{C}_k + \sum_{\ell = 1}^r \tilde{L}_{k \ell} a_\ell + \sum_{\ell = 1}^r \sum_{m = 1}^r \tilde{Q}_{k \ell m} a_\ell a_m, \end{equation}

with the modified constant, linear and quadratic terms given in tensor summation notation by

(5.6a)\begin{gather} \tilde{C}_k = V^{{-}1}_{k \mu} C_\mu \end{gather}

(5.6b)\begin{gather}\tilde{L}_{k \ell} = V^{{-}1}_{k \mu} L_{\mu \nu} V_{\nu \ell} \end{gather}

(5.6c)\begin{gather}\tilde{Q}_{k \ell m} = V^{{-}1}_{k \mu} Q_{\mu \nu \rho} V_{\nu \ell} V_{\rho m}. \end{gather}

The DMD–Petrov–Galerkin model can also be truncated differently from the POD–Galerkin model by selecting a subset of the DMD eigenvectors, so that $\boldsymbol {a} = \tilde {\boldsymbol{\mathsf{V}}} \boldsymbol {\alpha }$ with $\boldsymbol {a} \in \mathbb {R}^r$, $\boldsymbol {\alpha } \in \mathbb {C}^s$ and $\tilde {\boldsymbol{\mathsf{V}}} \in \mathbb {C}^{r \times s}$, and with $\boldsymbol{\mathsf{V}}^{-1}$ replaced by the pseudoinverse $\tilde {\boldsymbol{\mathsf{V}}}^+$ in (5.6).

Since this transformation is a rotation in the space of modal coefficients, the dynamics and qualitative behaviour of the DMD–Galerkin models do not change compared with figure 8. However, in the following we exploit nonlinear correlations in the coefficients to restrict the dynamics to the manifold of the flow; this is more convenient in the near-harmonic DMD basis, shown in figure 6. Since DMD analysis seeks to approximate the spectrum associated with a linear evolution operator of the flow, this might be considered analogous to the diagonalization step of a centre manifold or normal form analysis.

5.2. Sparse identification of nonlinear dynamics

As an alternative to projection-based reduced-order modelling, a low-dimensional system can be approximated directly from the data in a procedure typically called system identification. In a continuous-time setting, this is done by estimating parameters $\boldsymbol {\varXi }$ for a function $\boldsymbol {f}(\boldsymbol {a}; \boldsymbol {\varXi }) \approx \dot {\boldsymbol {a}}$ that solve the approximation problem

(5.7)\begin{equation} \textrm{minimize}~_{\boldsymbol{\varXi}} \left \langle \| \dot{\boldsymbol{a}} - \boldsymbol{f}(\boldsymbol{a}; \boldsymbol{\varXi}) \|^2 \right \rangle. \end{equation}

This is similar to the residual minimization in Galerkin projection, except that knowledge of the governing equations is neither required nor assumed. Instead, some intuition about the structure of the dynamics is typically encoded in the parameterization of $\boldsymbol {f}$. Different parameterizations lead to symbolic regression (Schmidt & Lipson Reference Schmidt and Lipson2009), operator inference (Peherstorfer & Willcox Reference Peherstorfer and Willcox2016) or deep learning (Vlachas et al. Reference Vlachas, Byeon, Wan, Sapsis and Komoutsakos2018). In the discrete-time counterpart to (5.7), the nonlinear autoregressive moving average model with exogenous inputs (NARMAX) framework provides a powerful approach that can incorporate time delays, stochastic forcing and exogenous inputs (Billings Reference Billings2013).

In this work we apply the SINDy approach to system identification (Brunton et al. Reference Brunton, Proctor and Kutz2016). Let $\boldsymbol {\varTheta }(\boldsymbol {a})$ denote a library of candidate functions of the time series $\boldsymbol {a}(t)$, e.g.

(5.8)\begin{equation} \boldsymbol{\varTheta}\left(\left[a_1 \quad a_2 \right]^\textrm{T} \right) = \left[a_1 \quad a_2 \quad a_1^2 \quad a_1 a_2 \quad a_2^2 \quad \cdots \right]^\textrm{T}. \end{equation}

We seek a sparse approximation $\dot {\boldsymbol {a}}\approx \boldsymbol {f}(\boldsymbol {a}; \boldsymbol {\varXi }) = \boldsymbol {\varXi } \boldsymbol {\varTheta }(\boldsymbol {a})$ in the range of these candidate functions. We can frame this as a linear algebra problem by forming the $r \times M$ data matrix $\boldsymbol{\mathsf{X}}$ as in § 4.2, where each column is a snapshot of modal coefficients in time. Similarly, we estimate the time derivative $\dot {\boldsymbol{\mathsf{X}}}$, in this case with second-order central differences. The SINDy formulation of the optimization problem (5.7) is then

(5.9)\begin{equation} \textrm{minimize}~_{\boldsymbol{\varXi}} \left\|{\dot{\boldsymbol{\mathsf{X}}} - \boldsymbol{\varXi} \boldsymbol{\varTheta}(\boldsymbol{\mathsf{X}}) }\right\|_2^2 + \gamma \|{\boldsymbol{\varXi}}\|_0, \end{equation}

where $\|{\cdot }\|_p$ indicates the $p$-norm and $\gamma$ is some regularization weight. This optimization problem is non-convex and requires a combinatorial search over function combinations. To avoid this, we follow Loiseau (Reference Loiseau2020) and approximate the solution to (5.9) with the greedy forward regression orthogonal least squares (FROLS) algorithm used in NARMAX analysis (Billings Reference Billings2013). See Brunton et al. (Reference Brunton, Proctor and Kutz2016), Loiseau & Brunton (Reference Loiseau and Brunton2018) and Loiseau (Reference Loiseau2020) for details on SINDy reduced-order modelling of fluid flows.

Although the systems obtained from POD–Galerkin projection will be dense in general, we expect that the dynamics can be closely approximated by a sparse combination of candidate functions. This may be justified intuitively by the sparse structure of the triadic interactions in isotropic flow (see § 2.2), where only $r^2$ out of $r^3$ possible interactions are admissible. Moreover, as observed in § 4.2, DMD approximates a diagonalization of the evolution operator, so that by analogy with normal form theory we may reasonably hope for a minimal representation of the dynamics if we work with DMD coefficients $\boldsymbol {\alpha }$ rather than POD coefficients $\boldsymbol {a}$. More generally, sparsity promotion reflects the inductive bias of Occam's razor or Pareto analysis, where we expect that the most important features of the dynamics will be due to a small subset of terms.

For incompressible flows, it has been repeatedly demonstrated that a library of low-order polynomials provides a good basis of functions. Quadratic terms are clearly necessary to capture the advective nonlinearity of the Navier–Stokes equations, but cubic terms allow the model to resolve Stuart–Landau-type nonlinear stability mechanisms (Loiseau & Brunton Reference Loiseau and Brunton2018). From a dynamical systems perspective, higher-order terms may be necessary to describe phenomena such as subcritical bifurcations, but are not necessary to resolve the nonlinear oscillator behaviour in the present case.

For PDEs with more general nonlinearity, low-order polynomials still present an attractive basis for SINDy. In many interesting regimes the effect of the nonlinearity may be relatively weak, so that quadratic and cubic polynomials can be seen as second- or third-order Taylor series approximations to the underlying nonlinearity. Moreover, even strongly nonlinear systems can be lifted with a change of variables to a coordinate system wherein the dynamics is linear quadratic (Rowley, Colonius & Murray Reference Rowley, Colonius and Murray2004; Qian et al. Reference Qian, Kramer, Peherstorfer and Willcox2020).

6.1. A model quasiperiodic cascade

Consider a model system of ODEs including cascading triadic-type interactions of the form of (2.7) where self-sustaining oscillations drive higher-order degrees of freedom

(6.1a)\begin{gather} a_1 = A \textrm{e}^{\textrm{i} \omega_a t}, \quad b_1 = B \textrm{e}^{\textrm{i} \omega_b t} , \end{gather}

(6.1b)\begin{gather}\dot{a}_k = \sum_{|\ell| < k} a_{\ell}a_{k - \ell}, \quad \dot{b}_k = \sum_{|\ell| < k} b_{\ell}b_{k - \ell}, \quad \dot{c}_k = \sum_{|\ell| < k} a_{\ell}b_{k - \ell}, \quad k > 1. \end{gather}

For $k=2$ the only interactions are at $|\ell | = 1$, generating second harmonics for $a_2$ and $b_2$, and oscillations at $\omega _a + \omega _b$ for $c_2$. At $k=3$, the forcing includes the $k=2$ terms, generating third harmonics for $a_3$ and $b_3$ and cross-talk for $c_3$. Up to some complex scaling the solutions take the form

(6.2a)\begin{gather} a_2 = A^2 \exp({2\textrm{i}\omega_a t}), \quad b_2 = B^2 \exp({2\textrm{i}\omega_b t}), \quad c_2 = A B \exp({\textrm{i}(\omega_a + \omega_b) t}), \end{gather}

(6.2b)\begin{gather} \left.\begin{gathered} a_3 = A^3 \exp({3\textrm{i}\omega_a t}), \quad b_3 = B^3 \exp({3\textrm{i}\omega_b t}), \\ c_3 = AB^2 \exp({\textrm{i}(\omega_a + 2 \omega_b) t}) + A^2B \exp({\textrm{i}(2\omega_a + \omega_b) t}). \end{gathered}\right\} \end{gather}

By $k=3$ the response of the $c_k$ cross-talk variable does not represent pure frequency content. If these represented modal coefficients we would expect the dynamics at (for instance) $\omega _a + 2\omega _b$ and $2\omega _a + \omega _b$ corresponds to different spatial structures so $c_3$ and higher-order terms would be separated into distinct coefficients. This serves to emphasize that the temporal coefficients and low-dimensional ODEs are only convenient representations of the full spatio-temporal dynamics and are not fundamental physical quantities.

A global observable such as the energy analogue $\sum _k (|a_k|^2 + |b_k|^2 + |c_k|^2)$ will have frequency content at all integer combinations of $\omega _a$ and $\omega _b$. As with the linear advection equation in § 2.1, these higher-order coefficients can be expressed as algebraic functions of the fundamental oscillators. For instance, the cross-talk variable $c_3$ can be replaced by $a_1 b_1^2 + a_1^2 b_1$, up to scaling. This example shows that periodic oscillation with cascading triadic interactions can generate quasiperiodic time series, point spectra similar to figure 4 and nonlinear algebraic dependence without direct interactions between the oscillators, provided the dominant oscillation frequencies are incommensurate.

6.2. The randomized dependence coefficient

As discussed in § 4, the POD coefficients are guaranteed to be linearly uncorrelated. The same is not necessarily true of DMD coefficients, but in practice they tend to be minimally correlated. However, as in the previous example, a network of triadic interactions forced by a limited number of driving oscillators can exhibit pure algebraic dependence on the active degrees of freedom. In other words, the higher-order variables can have perfect nonlinear correlation, even when uncorrelated in a linear sense.

This is an intuitive result, but challenging to evaluate in a principled way. In a probabilistic setting, mutual information is the most natural metric for generalized correlation, but it requires estimating integrals over conditional probability distributions. This is expensive and difficult for multidimensional signals, and the concept of mutual information itself is not necessarily well suited for purely deterministic systems. To address this issue, various nonlinear generalizations of the standard (Pearson's) linear correlation coefficient have been proposed in the statistics community. Recently, Lopez-Paz et al. (Reference Lopez-Paz, Hennig and Schölkopf2013) proposed the randomized dependence coefficient (RDC) as an efficient and convenient metric for nonlinear correlation that has the properties defined by Rènyi (Reference Rènyi1959) for generalized measures of dependence between variables.

The RDC combines linear canonical correlations analysis with randomized nonlinear projections to estimate nonlinear dependence; details are presented in Lopez-Paz et al. (Reference Lopez-Paz, Hennig and Schölkopf2013). Figure 9 gives examples of both the standard linear correlation coefficient and the RDC for several pairs of DMD coefficients. Coefficient pairs that are pure harmonics tend to score highly on the RDC, while coefficients with nonlinear cross-talk or multiple frequency components do not, even if there is a simple functional dependence on two active coefficients. This limits the potential of the RDC for unravelling the multivariate structure of the manifold, although it is a convenient diagnostic and visualization tool for identifying independent, dynamically active modes.

Figure 9. Selected phase portraits of DMD coefficients along with measures of linear and nonlinear correlation (Pearson's $\rho$ and the randomized dependence coefficient, respectively). While $\alpha _1$ and $\alpha _3$ are linearly uncorrelated ($\rho = 0$), the clear functional relationship between the two is reflected in the randomized dependence coefficient value; physically, $\alpha _3$ is a pure harmonic of $\alpha _1$. On the other hand, $\alpha _{17}$ corresponds to a nonlinear cross-talk mode that has no clear correlation with $\alpha _1$, either linear or nonlinear. Nevertheless, it can be accurately approximated by a simple polynomial function of the active degrees of freedom, as shown by table 1 and figure 11.

For example, figure 10 shows both the phase portraits of leading DMD coefficients (lower triangular portion) and the RDC values (upper triangular). In some cases (coloured outlines), the modes are clearly pure harmonics of one of the two driving mode pairs. This is reflected in the large values of the RDC for the harmonics, indicating that these coefficients can be directly expressed as algebraic functions of one or the other driving modes. However, based on the previous discussion we expect that coefficients representing frequency cross-talk might be multivariate functions of both driving mode pairs. In this case it is clear based on the energy content and harmonic structure of figure 10 that the pairs $(\alpha _1, \alpha _2)$ and $(\alpha _5, \alpha _6)$ are the driving degrees of freedom, but chaotic or turbulent flows might have more opaque causal structure.

Figure 10. Identification of the active degrees of freedom with the RDC. The lower triangular portion of the figure shows phase portraits of the real (horizontal axis) and imaginary (vertical axis) DMD coefficients, while the upper triangular portion depicts the RDC values scaled linearly in colour and radius. Two approximately independent clusters can be identified: the shear layer dynamics (blue) and inner-cavity oscillations (red). Each of these is associated with a dominant mode pair (solid borders) and pure harmonics (dashed borders) that are strongly nonlinearly correlated with the dominant modes. The other modes also have simple polynomial relationships with the active degrees of freedom but include cross terms that break the one-to-one nonlinear correlation (see table 1 and figures 9 and 11).

6.3. Manifold reduction via sparse regression

Based on the RDC analysis of the previous section, it is clear that certain modes are pure algebraic functions of one or the other driving mode pairs. In particular, harmonic modes such as those illustrated in figure 6 are polynomial functions of the fundamental mode. However, as demonstrated by the model in § 6.1, triadic interactions can also generate multivariate nonlinear correlations for modes representing frequency cross-talk. In addition, from a dynamical systems perspective we expect that a quasiperiodic system with two dominant frequencies should have a post-transient attractor described by four degrees of freedom: two generalized amplitude and phase pairs.

If this is the case, the modes that are not pure harmonics may still be approximately polynomial functions of the shear layer and inner cavity modes. We explore this hypothesis with the same approach as outlined in § 5.2 for the SINDy algorithm. Denoting these ‘active’ degrees of freedom by $\hat { \boldsymbol {\alpha }}(t) = \left [\alpha _1 \enspace \alpha _2 \enspace \alpha _5 \enspace \alpha _6 \right ]^\textrm {T}$, the library $\boldsymbol {\varTheta }(\hat {\boldsymbol {\alpha }})$ is defined as in (5.8). We assume the full coefficient vector can be approximated as

(6.3)\begin{equation} \boldsymbol{\alpha} \approx \boldsymbol{\mathsf{H}} \boldsymbol{\varTheta} (\hat{\boldsymbol{\alpha}}), \end{equation}

where the coefficient matrix $\boldsymbol{\mathsf{H}}$ is relatively sparse, as for ${\boldsymbol {\varXi }}$ in the SINDy optimization problem. For rows corresponding to active degrees of freedom, $\boldsymbol{\mathsf{H}}$ are unit vectors that produce an identity map (e.g. $\alpha _1 = \hat {\alpha }_1$).

In this case one motivation for sparse regression is that the library $\boldsymbol {\varTheta }$ tends to be fairly ill conditioned, so that approximation with a sparse combination of polynomials may help avoid overfitting. In addition, based on the preceding discussions about DMD and triadic interactions in frequency space, it is reasonable to expect that relatively few combinations of the driving frequencies will correspond to each DMD coefficient. Once the coefficient matrix $\boldsymbol{\mathsf{H}}$ is identified via a sparse regression algorithm (we use FROLS, as for the SINDy optimization), the functional relationships give a simple nonlinear dimensionality reduction; in this case the ‘latent variables’ are the active degrees of freedom $\hat {\boldsymbol {\alpha }}$ and we can approximate the full coefficient vector with the function $\boldsymbol {\alpha } \approx \boldsymbol {h} (\hat { \boldsymbol {\alpha }}) = \boldsymbol{\mathsf{H}} \boldsymbol {\varTheta } (\hat {\boldsymbol {\alpha }})$.

The key advantage to this representation of the modal coefficients is that it dramatically restricts the dimensionality of the state space. In the case of the quasiperiodic shear-driven cavity, we reduce from the $r=64$-dimensional subspace spanned by the DMD modes to a four-dimensional space of $\hat {\boldsymbol {\alpha }}$. Moreover, as $\alpha _2 = \alpha _1^*$ and $\alpha _6 = \alpha _5^*$, these mode pairs represent two generalized amplitude–phase pairs, as expected for a dynamics with a toroidal attractor. This eliminates the redundant variables introduced by the linear space–time separation of variables via a nonlinear manifold reduction.

The benefit of this reduced state space is immediately clear for system identification; we can apply SINDy to model the evolution of $\hat {\boldsymbol {\alpha }}$ and reconstruct the full state from this minimal representation. However, the manifold representation can also be used to improve the stability and accuracy of the projection-based Galerkin models. For a general nonlinear embedding of the form $\boldsymbol {a} = \boldsymbol {h} (\hat {\boldsymbol {a}} )$ and dynamical system $\dot {\boldsymbol {a}} = \boldsymbol {f} (\boldsymbol {a})$, consistency requires $\dot {\boldsymbol {a}} = \boldsymbol{\mathsf{J}}(\hat {\boldsymbol {a}})\dot {\hat {\boldsymbol {a}}}$, where $\boldsymbol{\mathsf{J}}(\hat {\boldsymbol {a}})$ is the Jacobian of $\boldsymbol {h}$ evaluated at $\hat {\boldsymbol {a}}$. Equivalently,

(6.4)\begin{equation} \frac{\textrm{d}{\hat{\boldsymbol{a}}}}{\textrm{d}t} = \boldsymbol{\mathsf{J}}^+(\hat{\boldsymbol{a}}) \dot{\boldsymbol{a}} = \boldsymbol{\mathsf{J}}^+(\hat{\boldsymbol{a}}) \boldsymbol{f}( \boldsymbol{h} (\hat{\boldsymbol{a}} ) ), \end{equation}

where $\boldsymbol{\mathsf{J}}^+$ is the pseudoinverse of $\boldsymbol{\mathsf{J}}$. This condition defines the reduced dynamics by constraining the velocity of $\boldsymbol {a}$ to the tangent space of the manifold defined by $\boldsymbol {h}$ (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983; Lee & Carlberg Reference Lee and Carlberg2020).

In the case of the linear-quadratic Galerkin dynamics (5.3) with the sparse polynomial manifold equation (6.3), the Jacobian consists of rows of the identity matrix by virtue of the fact that the reduced states $\hat {\boldsymbol {\alpha }}$ are also contained in the full coefficient vector $\boldsymbol {\alpha }$. Then the manifold Galerkin dynamics is

(6.5)\begin{equation} \dot{\alpha}_k = \tilde{C}_k + \sum_{\ell = 1}^r \tilde{L}_{k \ell} h_\ell( \hat{ \boldsymbol{\alpha}}) + \sum_{\ell = 1}^r \sum_{m = 1}^r \tilde{Q}_{k \ell m} h_\ell( \hat{ \boldsymbol{\alpha}}) h_m( \hat{ \boldsymbol{\alpha}}), \quad k = 1, 2, 5, 6, \end{equation}

where the tilde denotes the original POD–Galerkin operators rotated to the DMD coordinates (see §§ 4 and 5). Note that if $\boldsymbol {h}$ contains polynomials up to order $d$, the quadratic interactions in (6.5) lead to an effective nonlinearity of order $2d$. This can be viewed as a generalization of Stuart–Landau-type mean field models (Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003) or centre manifold expansions (Guckenheimer & Holmes Reference Guckenheimer and Holmes1983; Carini et al. Reference Carini, Auteri and Giannetti2015), although in these cases the manifold equation $\boldsymbol {h}$ is usually represented as a Taylor series truncated at second order. This is sufficient near bifurcations, but the more general form enabled by sparse regression allows improved resolution of the manifold structure.