2.1. Trajectory data

We assume that (1.1) is observed for $M_2$ time steps, starting at $M_1$ initial conditions. It is helpful to view the trajectory data as an $M_1\times M_2$ matrix

(2.1)\begin{equation} B_{{data}} = \begin{bmatrix} \pmb{x}_0^{(1)} & \cdots & \pmb{x}_{M_2-1}^{(1)}\\ \vdots & \ddots & \vdots \\ \pmb{x}_0^{(M_1)} & \cdots & \pmb{x}_{M_2-1}^{(M_1)}\\ \end{bmatrix}.\end{equation}

Each row of $B_{{data}}$ corresponds to an observation of a single trajectory of the dynamical system that is witnessed for $M_2$ time steps. In particular, $\pmb {x}_{i+1}^{(j)} = F(\pmb {x}_i^{(j)})$ for $0\leq i\leq M_2-2$ and $1\leq j\leq M_1$. For example, these snapshots could be measurements of unsteady velocity field across $M_1$ initial state configurations. One could therefore use just consider $M_2=2$. One could also consider a single trajectory, $M_1=1$, and large $M_2$. The exact structure depends on the data acquisition method. In general, (2.1) could represent measurements of the velocity along $M_1$ trajectories for $M_2$ time steps, a total of $M_1M_2$ measurements. One can also use our algorithms with data consisting of trajectories of different lengths (i.e. general snapshots).

Letting $\{\pmb {x}^{(m)}\}_{m=1}^M$ and $\{\pmb {y}^{(m)}\}_{m=1}^M$ be enumerations of the first $M_2-1$ columns and the second to final columns of $B_{{data}}$, respectively, with $M=M_1(M_2-1)$, we can conveniently represent the data as a finite set of $M$ pairs of measurements of the state

(2.2)\begin{equation} \{\pmb{x}^{(m)},\pmb{y}^{(m)}\}_{m=1}^M \text{ such that } \pmb{y}^{(m)}=F(\pmb{x}^{(m)}),\quad m=1,\ldots,M. \end{equation}

One could also consider measurements of certain observables $g, \{g(\pmb {x}^{(m)}),g(\pmb {y}^{(m)})\}_{m=1}^M$, and interpret the following algorithms as such. EDMD provides a way of using these measurements to approximate the operator $\mathcal {K}$ via a matrix $\mathbb {K}\in \mathbb {C}^{N\times N}$.

2.2. Extended DMD

We work in the Hilbert space $L^2(\varOmega,\omega )$ of square-integrable observables with respect to a positive measure $\omega$ on $\varOmega$. We do not assume that this measure is invariant. A common choice of $\omega$ is the standard Lebesgue measure. This choice is natural for Hamiltonian systems, for which the Koopman operator is unitary on $L^2(\varOmega,\omega )$ (Arnold Reference Arnold1989). Another common choice is for $\omega$ to correspond to an ergodic measure on an attractor (Mezić Reference Mezić2005).

Given a dictionary (a set of basis functions) $\{\psi _1,\ldots,\psi _{N}\}\subset \mathcal {D}(\mathcal {K})\subseteq L^2(\varOmega,\omega )$ of observables, EDMD constructs a matrix $\mathbb {K}\in \mathbb {C}^{N\times N}$ from the snapshot data (2.2) that approximates $\mathcal {K}$ on the finite-dimensional subspace ${V}_{{N}}=\mathrm {span}\{\psi _1,\ldots,\psi _{N}\}$. The choice of the dictionary is up to the user, with some common hand-crafted choices given in Williams et al. (Reference Williams, Kevrekidis and Rowley2015a, table 1). When the state-space dimension $d$ is large, as in this paper, it is beneficial to use a data-driven dictionary. We discuss this in § 3.3, where we present DMD (Kutz et al. Reference Kutz, Brunton, Brunton and Proctor2016a) and kEDMD (Williams et al. Reference Williams, Rowley and Kevrekidis2015b) in a unified framework.

Given a dictionary of observables, we define the vector-valued observable or ‘quasimatrix’

(2.3)\begin{equation} \varPsi(\pmb{x})=[\begin{matrix}\psi_1(\pmb{x}) & \cdots & \psi_{{N}}(\pmb{x}) \end{matrix}]\in\mathbb{C}^{1\times {N}}. \end{equation}

Any new observable $g\in V_{N}$ can then be written as $g(\pmb {x})=\sum _{j=1}^{N}\psi _j(\pmb {x})g_j=\varPsi (\pmb {x})\,\pmb {g}$ for some vector of constant coefficients $\pmb {g}=[g_1\cdots g_N]^\top \in \mathbb {C}^{N}$. It follows from (1.2) that

(2.4)\begin{equation} [\mathcal{K}g](\pmb{x})=\varPsi(F(\pmb{x})) \pmb{g}=\varPsi(\pmb{x})(\mathbb{K} \pmb{g})+R(\pmb{g},\pmb{x}), \end{equation}

where

(2.5)\begin{equation} R(\pmb{g},\pmb{x})=\sum_{j=1}^{N}\psi_j(F(\pmb{x}))g_j- \varPsi(\pmb{x})(\mathbb{K}\,\pmb{g}). \end{equation}

Typically, the subspace $V_{N}$ generated by the dictionary is not an invariant subspace of $\mathcal {K}$, and hence there is no choice of $\mathbb {K}$ that makes the error $R(\pmb {g},\pmb {x})$ zero for all choices of $g\in V_N$ and $\pmb {x}\in \varOmega$. Instead, it is natural to select $\mathbb {K}$ as to minimise

(2.6)\begin{equation} \int_\varOmega \max_{\|\pmb{g}\|=1}|R(\pmb{g},\pmb{x})|^2\,{\rm d} \omega(\pmb{x})=\int_\varOmega \|\varPsi(F(\pmb{x})) - \varPsi(\pmb{x})\mathbb{K}\|^2\,{\rm d} \omega(\pmb{x}),\end{equation}

where $\|\pmb {g}\|$ denotes the standard Euclidean norm of a vector $\pmb {g}$. Given a finite amount of snapshot data, we cannot directly evaluate the integral in (2.6). Instead, we approximate it via a quadrature rule by treating the data points $\{\pmb {x}^{(m)}\}_{m=1}^{M}$ as quadrature nodes with weights $\{w_m\}_{m=1}^{M}$. The discretised version of (2.6) is the following weighted least-squares problem:

(2.7)\begin{equation} \mathbb{K}\in \mathop{\mathrm{argmin}}\limits_{B\in\mathbb{C}^{N\times N}} \sum_{m=1}^{M} w_m\|\varPsi(\pmb{y}^{(m)})-\varPsi(\pmb{x}^{(m)})B\|^2. \end{equation}

We define the following two matrices:

(2.8a,b)\begin{equation} \varPsi_X=\begin{pmatrix} \varPsi(\pmb{x}^{(1)})\\ \vdots\\ \varPsi(\pmb{x}^{(M)}) \end{pmatrix}\in\mathbb{C}^{M\times N},\quad \varPsi_Y=\begin{pmatrix} \varPsi(\pmb{y}^{(1)})\\ \vdots \\ \varPsi(\pmb{y}^{(M)}) \end{pmatrix}\in\mathbb{C}^{M\times N},\end{equation}

and let $W=\mathrm {diag}(w_1,\ldots,w_M)$ be the diagonal weight matrix of the quadrature rule. A solution to (2.7) is

(2.9)\begin{equation} \mathbb{K}=(\varPsi_X^*W\varPsi_X)^{{{\dagger}}}(\varPsi_X^*W \varPsi_Y)=(\sqrt{W}\varPsi_X)^{\dagger}\sqrt{W}\varPsi_Y, \end{equation}

where ‘${\dagger}$’ denotes the pseudoinverse. In some applications, the matrix $\varPsi _X^*W\varPsi _X$ may be il -conditioned, so it is common to consider a regularisation such as a truncated singular value decomposition (SVD).

2.3. Quadrature and Galerkin methods

We can view EDMD as a Galerkin method. Note that

(2.10a,b)\begin{equation} \varPsi_X^*W\varPsi_X= \sum_{m=1}^{M} w_m \varPsi(\pmb{x}^{(m)})^* \varPsi(\pmb{x}^{(m)}),\quad \varPsi_X^*W\varPsi_Y= \sum_{m=1}^{M} w_m \varPsi(\pmb{x}^{(m)})^*\varPsi(\pmb{y}^{(m)}). \end{equation}

If the quadrature approximation converges, it follows that

(2.11a,b)\begin{equation} \lim_{M\rightarrow\infty}[\varPsi_X^*W\varPsi_X]_{jk} = \langle \psi_k,\psi_j \rangle\quad \text{and}\quad \lim_{M\rightarrow\infty}[\varPsi_X^*W\varPsi_Y]_{jk} = \langle \mathcal{K}\psi_k,\psi_j \rangle, \end{equation}

where $\langle {\cdot },{\cdot } \rangle$ is the inner product associated with the Hilbert space $L^2(\varOmega,\omega )$. Let $\mathcal {P}_{V_{N}}$ denote the orthogonal projection onto $V_{N}$. As $M\rightarrow \infty$, the above convergence means that $\mathbb {K}$ approaches a matrix representation of $\mathcal {P}_{V_{N}}\mathcal {K}\mathcal {P}_{V_{N}}^*$. Thus, EDMD is a Galerkin method in the large data limit $M\rightarrow \infty$. There are typically three scenarios where (2.11a,b) holds:

1. (i) Random sampling: in the initial definition of EDMD, $\omega$ is a probability measure and $\{\pmb {x}^{(m)}\}_{m=1}^M$ are drawn independently according to $\omega$ with the quadrature weights $w_m=1/M$. The strong law of large numbers shows that (2.11a,b) holds with probability one (Klus, Koltai & Schütte Reference Klus, Koltai and Schütte2016, § 3.4), provided that $\omega$ is not supported on a zero level set that is a linear combination of the dictionary (Korda & Mezić Reference Korda and Mezić2018, § 4). Convergence is typically at a Monte Carlo rate of ${O}(M^{-1/2})$ (Caflisch Reference Caflisch1998). From an experimental point of view, an example of random sampling could be $\{\pmb {x}^{(m)}\}$ observed with a sampling rate that is lower than the characteristic time period of the system of interest.

2. (ii) Highorder quadrature: if the dictionary and $F$ are sufficiently regular and we are free to choose the $\{\pmb {x}^{(m)}\}_{m=1}^{M}$, then it is beneficial to select $\{\pmb {x}^{(m)}\}_{m=1}^{M}$ as an $M$-point quadrature rule with weights $\{w_m\}_{m=1}^{M}$. This can lead to much faster convergence rates in (2.11a,b) (Colbrook & Townsend Reference Colbrook and Townsend2021), but can be difficult if $d$ is large.

3. (iii) Ergodic sampling: for a single fixed initial condition $\pmb {x}_0$ and $\pmb {x}^{(m)}=F^{m-1}(\pmb {x}_0)$ (i.e. data collected along one trajectory with $M_1=1$ in (2.1)), if the dynamical system is ergodic, then one can use Birkhoff's ergodic theorem to show (2.11a,b) (Korda & Mezić Reference Korda and Mezić2018). One chooses $w_m=1/M$ but the convergence rate is problem dependent (Kachurovskii Reference Kachurovskii1996). An example of ergodic sampling could be a time-resolved PIV dataset of a post-transient flow field over a long time period.

In this paper, we use (i) random sampling, and (iii) ergodic sampling, which are typical for experimental data collection since they arise from long-time trajectory measurements.

The convergence in (2.11a,b) implies that the EDMD eigenvalues approach the spectrum of $\mathcal {P}_{V_{N}}\mathcal {K}\mathcal {P}_{V_{N}}^*$ as $M\rightarrow \infty$. Thus, approximating the spectrum of $\mathcal {K}, \sigma (\mathcal {K})$, by the eigenvalues of $\mathbb {K}$ is closely related to the so-called finite section method (Böttcher & Silbermann Reference Böttcher and Silbermann1983). Since the finite section method can suffer from spectral pollution (spurious modes), spectral pollution is also a concern for EDMD (Williams et al. Reference Williams, Kevrekidis and Rowley2015a). It is important to have an independent way to measure the accuracy of the candidate eigenvalue–eigenvector pairs, which is what we propose in our ResDMD method presented in § 3.

2.3.1. Relationship with DMD

Williams et al. (Reference Williams, Kevrekidis and Rowley2015a) observed that if we choose a dictionary of observables of the form $\psi _j(\pmb {x})=e_j^*\pmb {x}$ for $j=1,\ldots,d$, where $\{e_j\}$ denote the canonical basis, then the matrix $\mathbb {K}=(\sqrt {W}\varPsi _X)^{\dagger} \sqrt {W}\varPsi _Y$ with $w_m=1/M$ is the transpose of the usual DMD matrix

(2.12)\begin{equation} \mathbb{K}_{{DMD}}=\varPsi_Y^\top\varPsi_X^{\top{\dagger}}= \varPsi_Y^\top\sqrt{W}(\varPsi_X^\top\sqrt{W})^{\dagger}{=} ((\sqrt{W}\varPsi_X)^{\dagger}\sqrt{W}\varPsi_Y)^\top=\mathbb{K}^\top. \end{equation}

Thus, DMD can be interpreted as producing a Galerkin approximation of the Koopman operator using the set of linear monomials as basis functions. EDMD can be considered an extension allowing nonlinear functions in the dictionary.

2.4. Koopman mode decomposition

Given an observable $h:\varOmega \rightarrow \mathbb {C}$, we approximate it by an element of $V_N$ via

(2.13)\begin{equation} h(\pmb{x})\approx \varPsi(\pmb{x})(\varPsi_X^*W\varPsi_X)^{\dagger}\varPsi_X^* W( h(\pmb{x}^{(1)}) \quad \cdots \quad h(\pmb{x}^{(M)}) )^\top. \end{equation}

Assuming that the quadrature rule converges, the right-hand side converges to the projection $\mathcal {P}_{V_N}h$ in $L^2(\varOmega,\omega )$ as $M\rightarrow \infty$. Assuming that $\mathbb {K}V=V\varLambda$ for a diagonal matrix $\varLambda$ of eigenvalues and a matrix $V$ of eigenvectors, we obtain

(2.14)\begin{equation} h(\pmb{x})\approx \varPsi(\pmb{x})V [V^{{-}1}(\sqrt{W}\varPsi_X)^{\dagger} \sqrt{W}(\begin{matrix} h(\pmb{x}^{(1)}) & \cdots & h(\pmb{x}^{(M)}) \end{matrix})^\top]. \end{equation}

As a special case, and vectorising, we have the Koopman mode decomposition

(2.15)\begin{equation} \pmb{x}\approx \underbrace{\varPsi(\pmb{x})V}_{\text{Koopman e-functions}} \underbrace{[V^{{-}1}(\sqrt{W}\varPsi_X)^{\dagger} \sqrt{W}(\begin{matrix} \pmb{x}^{(1)} & \cdots & \pmb{x}^{(M)} \end{matrix})^\top]}_{N\times d\text{ matrix of Koopman modes}}. \end{equation}

The $j$th Koopman mode corresponds to the $j$th row of the matrix in square brackets, and $\varPsi V$ is a quasimatrix of approximate Koopman eigenfunctions.

3.1. ResDMD and a new data matrix

Whilst EDMD obtains eigenvalue–eigenvector pairs for an approximation of the Koopman operator, it cannot verify the accuracy of the computed pairs. However, we show below how one can confidently identify true physical turbulent flow structures by rigorously rejecting inaccurate predictions (spurious modes). This rigorous measure of accuracy is the linchpin of our new ResDMD method and is shown pictorially in figure 1.

We approximate residuals to measure the accuracy of a candidate eigenvalue–eigenvector pair $(\lambda,g)$ (which could, for example, be computed from $\mathbb {K}$ or some other method). For $\lambda \in \mathbb {C}$ and $g=\varPsi \,\pmb {g}\in V_{N}$, the squared relative residual is

(3.1)\begin{align} &\frac{\displaystyle \int_{\varOmega}|[\mathcal{K}g](\pmb{x})- \lambda g(\pmb{x})|^2\, {\rm d} \omega(\pmb{x})}{\displaystyle \int_{\varOmega} |g(\pmb{x})|^2\, {\rm d} \omega(\pmb{x})}\nonumber\\ &\quad=\frac{\displaystyle \sum_{j,k=1}^{N}\overline{{g}_j}{g}_k [\langle \mathcal{K}\psi_k,\mathcal{K}\psi_j\rangle - \lambda\langle \psi_k,\mathcal{K}\psi_j\rangle -\bar{\lambda}\langle \mathcal{K}\psi_k,\psi_j\rangle+|\lambda|^2\langle \psi_k,\psi_j\rangle]}{\displaystyle \sum_{j,k=1}^{N}\overline{{g}_j}{g}_k\langle \psi_k,\psi_j\rangle}. \end{align}

We approximate this residual using the same quadrature rule in § 2.2

(3.2)\begin{equation} {res}(\lambda,g)^2 = \frac{\pmb{g}^*[\varPsi_Y^*W\varPsi_Y- \lambda[\varPsi_X^*W\varPsi_Y]^* - \bar{\lambda}\varPsi_X^*W\varPsi_Y + |\lambda|^2\varPsi_X^*W\varPsi_X]\pmb{g}}{\pmb{g}^*[\varPsi_X^* W\varPsi_X]\pmb{g}}.\end{equation}

As well as the matrices $\varPsi _X^*W\varPsi _X$ and $\varPsi _X^*W\varPsi _Y$ found in EDMD, (3.2) has the additional matrix $\varPsi _Y^*W\varPsi _Y$. Note that this additional matrix is no more expensive to compute than either of the EDMD matrices, and uses the same data for its construction as EDMD. Since $\varPsi _Y^*W\varPsi _Y= \sum _{m=1}^{M} w_m \varPsi (\pmb {y}^{(m)})^*\varPsi (\pmb {y}^{(m)})$, if the quadrature rule converges then

(3.3)\begin{equation} \lim_{M\rightarrow\infty}[\varPsi_Y^*W\varPsi_Y]_{jk} = \langle \mathcal{K}\psi_k,\mathcal{K}\psi_j \rangle. \end{equation}

Hence, $\varPsi _Y^*W\varPsi _Y$ formally corresponds to an approximation of $\mathcal {K}^*\mathcal {K}$. We also have

(3.4)\begin{equation} \lim_{M\rightarrow\infty} {res}(\lambda,g)^2=\frac{\displaystyle \int_{\varOmega}|[\mathcal{K}g](\pmb{x})-\lambda g(\pmb{x})|^2\, {\rm d} \omega(\pmb{x})}{\displaystyle \int_{\varOmega} |g(\pmb{x})|^2\, {\rm d} \omega(\pmb{x})}. \end{equation}

In Colbrook & Townsend (Reference Colbrook and Townsend2021), it was shown that the quantity ${res}(\lambda,g)$ can be used to rigorously compute spectra and pseudospectra of $\mathcal {K}$. Next, we summarise some of the algorithms and results of Colbrook & Townsend (Reference Colbrook and Townsend2021).

3.2. ResDMD for computing spectra and pseudospectra

3.2.1. Avoiding spurious eigenvalues

Algorithm 1 uses the residual defined in (3.2) to avoid spectral pollution (spurious modes). As is usually done, we assume that $\mathbb {K}$ is diagonalisable. We first find the eigenvalues and eigenvectors of $\mathbb {K}$, i.e. we solve $(\varPsi _X^*W\varPsi _X)^{\dagger} (\varPsi _X^*W\varPsi _Y){\pmb g} = \lambda {\pmb g}$. One can solve this eigenproblem directly, but it is often numerically more stable to solve the generalised eigenproblem $(\varPsi _X^*W\varPsi _Y){\pmb g} = \lambda (\varPsi _X^*W\varPsi _X){\pmb g}$. Afterwards, to avoid spectral pollution, we discard eigenpairs with a relative residual larger than a specified accuracy goal $\epsilon >0$.

The procedure is a simple modification of EDMD, as the only difference is a clean-up step where spurious eigenpairs are discarded based on their residual. This clean-up step is typically faster to execute than the eigendecomposition in step 2 of algorithm 1. The total computational cost of algorithm 1 is ${O}(N^2M+N^3)$, which is the same as EDMD. The clean up avoids spectral pollution and also removes eigenpairs that are inaccurate due to numerical errors associated with non-normal operators, up to the relative tolerance $\epsilon$. The following result (Colbrook & Townsend Reference Colbrook and Townsend2021, theorem 4.1) makes this precise.

Algorithm 1 :

ResDMD for computing eigenpairs without spectral pollution. The corresponding Koopman mode decomposition is given in (3.6).

Theorem 3.1 Let $\mathcal {K}$ be the associated Koopman operator of the dynamical system (1.1) from which snapshot data are collected. Let $\varLambda _{M}$ denote the eigenvalues in the output of algorithm 1. Then, assuming convergence of the quadrature rule in § 2.2

(3.5)\begin{equation} \limsup_{M\rightarrow\infty} \max_{\lambda\in \varLambda_{M}}\|(\mathcal{K}-\lambda)^{{-}1}\|^{{-}1}\leq \epsilon. \end{equation}

We can also use algorithm 1 to clean up the Koopman mode decomposition in (2.15). To do this, we simply let $V^{(\epsilon )}$ denote the matrix whose columns are the eigenvectors $\pmb {g}_j$ with ${res}(\lambda _j,g_{(j)})\leq \epsilon$ and compute the Koopman mode decomposition with respect to $\varPsi _X^{(\epsilon )}=\varPsi _X V^{(\epsilon )}$ and $\varPsi _Y^{(\epsilon )}=\varPsi _Y V^{(\epsilon )}$. Since $(\sqrt {W}\varPsi _X V^{(\epsilon )})^{\dagger} \sqrt {W}\varPsi _Y V^{(\epsilon )}$ is diagonal, the Koopman mode decomposition now becomes

(3.6)\begin{equation} \pmb{x}\approx \underbrace{\varPsi(\pmb{x})V^{(\epsilon)}}_{\text{Koopman e-functions}} \underbrace{[(\sqrt{W}\varPsi_XV^{(\epsilon)})^{\dagger} \sqrt{W}(\begin{matrix} \pmb{x}^{(1)} & \cdots & \pmb{x}^{(M)} \end{matrix})^\top]}_{N\times d\text{ matrix of Koopman modes}}. \end{equation}

3.2.2. Computing the full spectrum and pseudospectra

Theorem 3.1 tells us that, in the large data limit, algorithm 1 computes eigenvalues inside the $\epsilon$-pseudospectrum of $\mathcal {K}$. Hence, algorithm 1 avoids spectral pollution and returns reasonable eigenvalues. The $\epsilon$-pseudospectrum of $\mathcal {K}$ is (Trefethen & Embree Reference Trefethen and Embree2005)

(3.7)\begin{equation} {\sigma}_{\epsilon}(\mathcal{K}):=\mathrm{cl} (\{\lambda\in\mathbb{C}:\|(\mathcal{K}-\lambda)^{{-}1}\| > 1/\epsilon\})=\mathrm{cl}(\cup_{\|\mathcal{B}\|< \epsilon}{\sigma}(\mathcal{K}+\mathcal{B})), \end{equation}

where ${\rm cl}$ denotes the closure of a set, and $\lim _{\epsilon \downarrow 0}{\sigma }_{\epsilon }(\mathcal {K})= {\sigma }(\mathcal {K})$. Despite theorem 3.1, algorithm 1 may not approximate the whole of ${\sigma }_{\epsilon }(\mathcal {K})$, even as $M\rightarrow \infty$ and $N\rightarrow \infty$. This is because the eigenvalues of $\mathcal {P}_{V_{N}}\mathcal {K}\mathcal {P}_{V_{N}}^*$ may not approximate the whole of ${\sigma }(\mathcal {K})$ as $N\rightarrow \infty$ (Colbrook & Townsend Reference Colbrook and Townsend2021; Mezic Reference Mezic2022), even if $\cup _{N}V_N$ is dense in $L^2(\varOmega,\omega )$.

To overcome this issue, algorithm 2 computes practical approximations of $\epsilon$-pseudospectra with rigorous convergence guarantees. Assuming convergence of the quadrature rule in § 2.2, in the limit $M\rightarrow \infty$, the key quantity

(3.8)\begin{equation} \tau_N(\lambda) = \min_{\pmb{g}\in\mathbb{C}^{N}} {res}(\lambda,\varPsi\pmb{g}), \end{equation}

is an upper bound for $\|(\mathcal {K}-\lambda )^{-1}\|^{-1}$. The output of algorithm 2 is guaranteed to be inside the $\epsilon$-pseudospectrum of $\mathcal {K}$. As $N\rightarrow \infty$ and the grid of points is refined, algorithm 2 converges to the pseudospectrum uniformly on compact subsets of $\mathbb {C}$ (Colbrook & Townsend Reference Colbrook and Townsend2021, theorems B.1 and B.2). In practice the grid $\{z_j\}_{j=1}^k$ is chosen for plotting purposes (e.g. in figure 3). Strictly speaking, we converge to the approximate point pseudospectrum, a more complicated algorithm leads to computation of the full pseudospectrum – see Colbrook & Townsend (Reference Colbrook and Townsend2021, appendix B). For brevity, we have not included a statement of the results. We can then compute the spectrum by taking $\epsilon \downarrow 0$. Algorithm 2 also computes observables $g$ with ${res}(\lambda,g)< \epsilon$, known as $\epsilon$-approximate eigenfunctions. The computational cost of algorithm 2 is ${O}(N^2M+kN^3)$. However, the computation of $\tau _N$ can be done in parallel over the $k$ grid points, and we can refine the grid adaptively near regions of interest.

Algorithm 2 :

ResDMD for estimating ε-pseudospectra.

3.3. Choice of dictionary

When $d$ is large, it can be impractical to store or form the matrix $\mathbb {K}$ since the initial value of $N$ may be huge. In other words, we run into the curse of dimensionality. We consider two common methods to overcome this issue:

1. (i) DMD: in this case, the dictionary consists of linear functions (see § 2.3.1). It is standard to form a low-rank approximation of $\sqrt {W}\varPsi _X$ via a truncated SVD as

   (3.9)\begin{equation} \sqrt{W}\varPsi_X\approx U_r\varSigma_rV_r^*. \end{equation}
   Here, $\varSigma _r\in \mathbb {C}^{r\times r}$ is diagonal with strictly positive diagonal entries, and $V_r\in \mathbb {C}^{N\times r}$ and $U_r\in \mathbb {C}^{M\times r}$ have $V_r^*V_r=U_r^*U_r=I_r$. We then form the matrix
   (3.10)\begin{equation} \tilde{\mathbb{K}}=(\sqrt{W}\varPsi_XV_r)^{\dagger}\sqrt{W} \varPsi_YV_r=\varSigma_r^{{-}1}U_r^*\sqrt{W}\varPsi_YV_r=V_r^* \mathbb{K}V_r\in\mathbb{C}^{r\times r}. \end{equation}
   Note that to fit into our Galerkin framework, this matrix is the transpose of the DMD matrix that is commonly computed in the literature.

2. (ii) Kernelised EDMD (kEDMD): kEDMD (Williams et al. Reference Williams, Rowley and Kevrekidis2015b) aims to make EDMD practical for large $d$. Supposing that $\varPsi _X$ is of full rank, kEDMD constructs a matrix with an identical formula to (3.10) when $r=M$, for which we have the equivalent form

   (3.11)\begin{equation} \tilde{\mathbb{K}}=(\varSigma_M^{\dagger} U_M^*)(\sqrt{W}\varPsi_Y \varPsi_X^*\sqrt{W})(U_M\varSigma_M^{\dagger}). \end{equation}
   Suitable matrices $U_M$ and $\varSigma _M$ can be recovered from the eigenvalue decomposition $\sqrt {W}\varPsi _X\varPsi _X^*\sqrt {W}=U_M\varSigma _M^2U_M^*$. Moreover, both matrices $\sqrt {W}\varPsi _X\varPsi _X^*\sqrt {W}$ and $\sqrt {W}\varPsi _Y\varPsi _X^*\sqrt {W}$ can be computed using inner products. The kEDMD applies the kernel trick to compute the inner products in an implicitly defined reproducing Hilbert space $\mathcal {H}$ with inner product $\langle {\cdot },{\cdot }\rangle _{\mathcal {H}}$ (Scholkopf Reference Scholkopf2001). A positive–definite kernel function $\mathcal {S}:\varOmega \times \varOmega \rightarrow \mathbb {R}$ induces a feature map $\varphi :\mathbb {R}^d\rightarrow \mathcal {H}$ so that $\langle \varphi (\pmb {x}), \varphi (\pmb {y})\rangle _{\mathcal {H}}=\mathcal {S}(\pmb {x},\pmb {y})$. This leads to an implicit choice of (typically nonlinear) dictionary $\varPsi (\pmb {x})$, that is dependent on the kernel function, so that $\varPsi (\pmb {x})\varPsi (\pmb {y})^*=\langle \varphi (\pmb {x}),\varphi (\pmb {y})\rangle _{\mathcal {H}}=\mathcal {S}(\pmb {x},\pmb {y})$. Often $\mathcal {S}$ can be evaluated in ${O}(d)$ operations, meaning that $\tilde {\mathbb {K}}$ is constructed in ${O}(dM^2)$ operations.

In either of these two cases, the approximation of $\mathcal {K}$ is equivalent to using a new dictionary with feature map $\varPsi (\pmb {x})V_r\in \mathbb {C}^{1\times r}$. In the case of DMD, we found it beneficial to use the mathematically equivalent choice $\varPsi (\pmb {x})V_r\varSigma _r^{-1}$, which is numerically better conditioned. To see why, note that $\sqrt {W}\varPsi _XV_r\varSigma _r^{-1}\approx U_r$ and $U_r$ has orthonormal columns.

3.3.2. Using two subsets of the snapshot data

A simple remedy to avoid the problem in § 3.3.1 is to consider two sets of snapshot data. We consider an initial set $\{\tilde {\pmb {x}}^{(m)},\tilde {\pmb {y}}^{(m)}\}_{m=1}^{M'}$, that we use to form our dictionary. We then apply ResDMD to the computed dictionary with a second set of snapshot data $\{\hat {\pmb {x}}^{(m)},\hat {\pmb {y}}^{(m)}\}_{m=1}^{M''}$, allowing us to prove convergence as $M''\rightarrow \infty$.

How to acquire a second set of snapshot data depends on the problem and method of data collection. Given snapshot data with random and independent $\{\pmb {x}^{(m)}\}$, one can split up the snapshot data into two parts. For initial conditions that are distributed according to a high-order quadrature rule, if one already has access to $M'$ snapshots, then one must typically go back to the original dynamical system and request $M''$ further snapshots. For ergodic sampling along a trajectory, we can let $\{\tilde {\pmb {x}}^{(m)},\tilde {\pmb {y}}^{(m)}\}_{m=1}^{M'}$ correspond to the initial $M'+1$ points of the trajectory ($\tilde {\pmb {x}}^{(m)}=F^{m-1}(\pmb {x}_0)$ for $m=1,\ldots,M'$) and let $\{\hat {\pmb {x}}^{(m)},\hat {\pmb {y}}^{(m)}\}_{m=1}^{M''}$ correspond to the initial $M''+1$ points of the trajectory ($\hat {\pmb {x}}^{(m)}=F^{m-1}(\pmb {x}_0)$ for $m=1,\ldots,M''$).

Algorithm 3 :

ResDMD with DMD selected observables.

Algorithm 4 :

ResDMD with kEDMD selected observables.

In the case of DMD, algorithm 3 summarises the two-stage process. Often a suitable choice of $N$ can be obtained by studying the decay of the singular values of the data matrix or the condition number of the matrix $\varPsi _X^*W\varPsi _X$. When $M'\leq d$, the computational cost of steps 2 and 3 of algorithm 3 are ${O}(d{M'}^2)$ and ${O}(Nd{M''})$, respectively.

In the case of kEDMD, we follow Colbrook & Townsend (Reference Colbrook and Townsend2021), and algorithm 4 summarises the two-stage process. The choice of kernel $\mathcal {S}$ determines the dictionary and the best choice depends on the application. In the following experiments, we use the Laplacian kernel $\mathcal {S}(\pmb {x},\pmb {y})=\exp (-\gamma {\|\pmb {x}-\pmb {y}\|})$, where $\gamma$ is the reciprocal of the average $\ell ^2$-norm of the snapshot data after it is shifted to have mean zero.

We can now apply the theory of § 3.2 in the limit $M''\rightarrow \infty$. It is well known that the eigenvalues computed by DMD and kEDMD may suffer from spectral pollution. However, crucially in our setting, we do not directly use these methods to compute spectral properties of $\mathcal {K}$. Instead, we only use them to select a good dictionary of size $N$, after which our rigorous ResDMD algorithms can be used. Moreover, we use $\{\hat {\pmb {x}}^{(m)},\hat {\pmb {y}}^{(m)}\}_{m=1}^{M''}$ to check the quality of the constructed dictionary. By studying the residuals and using the error control in ResDMD, we can tell a posteriori whether the dictionary is satisfactory and whether $N$ is sufficiently large.

Finally, it is worth pointing out that the above choices of dictionaries are certainly not the only choices. ResDMD can be applied to any suitable choice. For example, one could use other data-driven methods such as diffusion kernels (Giannakis et al. Reference Giannakis, Kolchinskaya, Krasnov and Schumacher2018) or trained neural networks (Li et al. Reference Li, Dietrich, Bollt and Kevrekidis2017; Murata, Fukami & Fukagata Reference Murata, Fukami and Fukagata2020).

4.1. Spectral measures and Koopman mode decompositions

Given an observable $g\in L^2(\varOmega,\omega )$ of interest, the spectral measure of $\mathcal {K}$ with respect to $g$ is a measure $\nu _g$ defined on the periodic interval $[-{\rm \pi},{\rm \pi} ]_{{per}}$. If $g$ is normalised so that $\|g\|=1$, then $\nu _g$ is a probability measure, otherwise $\nu _g$ is a positive measure of total mass $\|g\|^2$. We provide a mathematical description of spectral measures in § A for completeness. The reader should think of these measures as supplying a diagonalisation of $\mathcal {K}$.

For example, the decomposition in (A3) provides important information on the evolution of dynamical systems. Suppose that there is no singular continuous spectrum, then any $g\in L^2(\varOmega,\omega )$ can be written as

(4.1)\begin{equation} g=\sum_{\lambda\in{\sigma}_{{p}}(\mathcal{K})}c_\lambda \varphi_\lambda+\int_{[-{\rm \pi},{\rm \pi}]_{{per}}}\phi_{\theta,g}\, {\rm d} \theta, \end{equation}

where ${\sigma }_{{p}}(\mathcal {K})$ is the set of eigenvalues of $\mathcal {K}$, the $\varphi _\lambda$ are the eigenfunctions of $\mathcal {K}, c_\lambda$ are expansion coefficients and the density of the absolutely continuous part of $\nu _g$ is $\rho _g(\theta )=\langle \phi _{\theta,g},g\rangle$. One should think of $\phi _{\theta,g}$ as a ‘continuously parametrised’ collection of eigenfunctions. Using (4.1), one obtains the Koopman mode decomposition (Mezić Reference Mezić2005)

(4.2)\begin{equation} g(\pmb{x}_n)=[\mathcal{K}^ng](\pmb{x}_0)=\sum_{\lambda \in{\sigma}_{{p}}(\mathcal{K})}c_{\lambda}\lambda^n \varphi_\lambda(\pmb{x}_0)+\int_{[-{\rm \pi},{\rm \pi}]_{{per}}}\, {\rm e}^{{\rm i}n\theta}\phi_{\theta,g}(\pmb{x}_0)\, {\rm d} \theta. \end{equation}

For characterisations of the dynamical system in terms of these decompositions, see Halmos (Reference Halmos2017), Mezić (Reference Mezić2013) and Zaslavsky (Reference Zaslavsky2002). Typically, the eigenvalues correspond to isolated oscillation frequencies in the fluid flow and the growth rates of stable and unstable modes, whilst the continuous spectrum corresponds to chaotic motion. Computing the measures $\nu _g$ provides us with a diagonalisation of the nonlinear dynamical system in (1.1).

4.2. Spectral measures and autocorrelations

The Fourier coefficients of $\nu _g$ are given by

(4.3)\begin{equation} \widehat{\nu_g}(n):=\frac{1}{2{\rm \pi}}\int_{[-{\rm \pi},{\rm \pi}]_{{per}}} \, {\rm e}^{-{\rm i}n\theta}\,{\rm d} \nu_g(\theta),\quad n\in\mathbb{Z}. \end{equation}

These Fourier coefficients can be expressed in terms of correlations $\langle \mathcal {K}^n g,g\rangle$ and $\langle g,\mathcal {K}^n g\rangle$ (Colbrook & Townsend Reference Colbrook and Townsend2021). That is, for $g\in L^2(\varOmega,\omega )$

(4.4a,b)\begin{equation} \widehat{\nu_g}(n)=\frac{1}{2{\rm \pi}}\langle \mathcal{K}^{{-}n}g,g\rangle, \quad n<0, \quad \widehat{\nu_g}(n)=\frac{1}{2{\rm \pi}}\langle g, \mathcal{K}^ng\rangle,\quad n\geq 0. \end{equation}

From (4.4a,b), we see that $\widehat {\nu _{g}}(-n)=\overline {\widehat {\nu _{g}}(n)}$ for $n\in \mathbb {Z}$. In particular, $\nu _g$ is completely determined by the forward-time dynamical autocorrelations $\langle g,\mathcal {K}^ng\rangle$ with $n\geq 0$. Equivalently, the spectral measure of $\mathcal {K}$ with respect to $g\in L^2(\varOmega,\omega )$ is a signature for the forward-time dynamics of (1.1). This connection allows us to interpret spectral measures as an infinite-dimensional version of power spectra in § 4.5.

4.3. Computing autocorrelations

To approximate spectral measures, we make use of the relation (4.4a,b) between the Fourier coefficients of $\nu _g$ and the autocorrelations $\langle g,\mathcal {K}^ng\rangle$. There are typically three ways to compute the autocorrelations, corresponding to the three scenarios discussed in § 2.3:

1. 1. Random sampling: the autocorrelations can be approximated as

   (4.5)\begin{equation} \langle g,\mathcal{K}^ng\rangle\approx\frac{1}{M} \sum_{m=1}^{M}g(\pmb{x}^{(m)})\overline{[\mathcal{K}^ng](\pmb{x}^{(m)})}. \end{equation}

2. 2. High-order quadrature: the autocorrelations can be approximated as

   (4.6)\begin{equation} \langle g,\mathcal{K}^ng\rangle = \int_{\varOmega}g(\pmb{x}) \overline{[\mathcal{K}^ng](\pmb{x})}\,{\rm d} \omega(\pmb{x}) \approx\sum_{m=1}^{M}w_mg(\pmb{x}^{(m)})\overline{[\mathcal{K}^ng] (\pmb{x}^{(m)})}. \end{equation}

3. 3. Ergodic sampling: the autocorrelations can be approximated as

   (4.7)\begin{equation} \langle g,\mathcal{K}^ng\rangle\approx \frac{1}{M-n} \sum_{m=0}^{M-n-1} g(\pmb{x}_{m})\overline{[\mathcal{K}^ng] (\pmb{x}_{m})}= \frac{1}{M-n} \sum_{m=0}^{M-n-1} g(\pmb{x}_{m})\overline{g(\pmb{x}_{m+n})}. \end{equation}

The first two methods require multiple snapshots of the form $\{\pmb {x}_0^{(m)},\ldots, \pmb {x}_{n}^{(m)}\}$, with $[\mathcal {K}^ng](\pmb {x}_0^{(m)})=g(\pmb {x}_{n}^{(m)})$. Ergodic sampling only requires a single trajectory, and the ergodic averages in (4.7) can be rewritten as discrete (non-periodic) convolutions. Thus, by zero padding, the fast Fourier transform computes all averages simultaneously and rapidly.

4.4. Computing spectral measures

We suppose now that one has already computed the autocorrelations $\langle g,\mathcal {K}^ng\rangle$ for $0\leq n\leq N_{{ac}}$. In practice, given a fixed data set of $M$ snapshots, we choose a suitable $N_{{ac}}$ by checking for convergence of the autocorrelations by comparing with smaller values of $M$. Following Colbrook & Townsend (Reference Colbrook and Townsend2021), we define an approximation to $\nu _g$ as

(4.8)\begin{equation} \nu_{g,N_{{ac}}}(\theta)=\sum_{n={-}N_{{ac}}}^{N_{{ac}}} \varphi \left(\frac{n}{N_{{ac}}}\right)\widehat{\nu_g}(n)\,{\rm e}^{{\rm i}n\theta}. \end{equation}

The function $\varphi : [-1,1]\rightarrow \mathbb {R}$ is often called a filter function (Tadmor Reference Tadmor2007; Hesthaven Reference Hesthaven2017) and $\varphi (x)$ is close to $1$ when $x$ is close to $0$, but tapers to $0$ near $x=\pm 1$. Algorithm 5 summarises the approach, and the choice of $\varphi$ affects the convergence.

For $m\in \mathbb {N}$, suppose that $\varphi$ is even, continuous and compactly supported on $[-1,1]$ with

1. (a) $\varphi \in \mathcal {C}^{m-1}([-1,1])$;

2. (b) $\varphi (0)=1$ and $\varphi ^{(n)}(0)=0$ for any integer $1\leq n\leq m-1$;

3. (c) $\varphi ^{(n)}(1)=0$ for any integer $0\leq n\leq m-1$,

4. (d) $\varphi |_{[0,1]}\in \mathcal {C}^{m+1}([0,1])$.

Then, we have the following forms of convergence (Colbrook & Townsend Reference Colbrook and Townsend2021, § 3):

1. (i) Under suitable (local) smoothness assumptions, $|\nu _{g,N_{{ac}}}(\theta )-\rho _g(\theta )|= {O}(N_{{ac}}^{-m} \log (N_{{ac}}))$,

2. (ii) A suitable rescaling of $\nu _{g,N_{{ac}}}\ (\theta_0) $ approximates the point spectrum at $\theta _0$,

3. (iii) For any $\phi \in \mathcal {C}^{n,\alpha }([-{\rm \pi},{\rm \pi} ]_{{per}})$,

   (4.9) \begin{equation} \left|\int_{[-{\rm \pi},{\rm \pi}]_{{per}}} \phi(\theta)\nu_{g,N_{{ac}}} (\theta)\,{\rm d} \theta - \int_{[-{\rm \pi},{\rm \pi}]_{{per}}} \phi(\theta)\,{\rm d} \nu_g(\theta)\right| \lesssim \|\phi\|_{\mathcal{C}^{n,\alpha}} (N_{{ac}}^{-(n+\alpha)}+ N_{{ac}}^{{-}m}\log(N_{{ac}})). \end{equation}

This last property is known as weak convergence. One should think of $\nu _{g,N_{{ac}}}$ as a smooth function that approximates the spectral measure $\nu _g$ to order ${O}(N_{{ac}}^{-m}\log (N_{{ac}}))$, with a frequency smoothing scale of ${O}(N_{{ac}}^{-m})$.

Algorithm 5 :

Approximating spectral measures from autocorrelations.

One can also compute spectral measures without autocorrelations. Colbrook & Townsend (Reference Colbrook and Townsend2021) also develop high-order methods based on rational kernels that approximate spectral measures using the ResDMD matrices. Computing suitable residuals allows an adaptive and rigorous selection of the smoothing parameter used in the convolution. In particular, this set-up allows us to deal with general snapshot data $\{\pmb {x}^{(m)},\pmb {y}^{(m)}=F(\pmb {x}^{(m)})\}_{m=1}^M$ without the need for long trajectories. For brevity, we omit the details.

4.5. Interpretation of Koopman spectral measures as power spectra

We can interpret algorithm 5 as an approximation of the power spectrum of the signal $g(\pmb {x}(t))$, given by Glegg & Devenport (Reference Glegg and Devenport2017, chapter 8)

(4.10)\begin{equation} S_{gg}(f)=\int_{{-}T}^TR_{gg}(t)\, {\rm e}^{2{\rm \pi} {\rm i} ft}\, {\rm d} t, \end{equation}

over a time window $[-T,T]$ and for frequency $f$ (measured in Hertz). Here, $R_{gg}(t)$ is the delay autocorrelation function, defined for $t\geq 0$ as

(4.11)\begin{equation} R_{gg}(t)=\langle g,g\circ F_t\rangle, \end{equation}

where $F_t$ is the forward-time propagator for a time step $t$. In particular, we have

(4.12) \begin{equation} R_{gg}(n{\rm \Delta} t)=\left\{\begin{array}{@{}ll} \langle g,\mathcal{K}^ng\rangle & \text{if }n\geq 0\\ \overline{\langle g,\mathcal{K}^{{-}n}g\rangle}= \langle\mathcal{K}^{{-}n} g,g\rangle & \text{otherwise} \end{array}\right. . \end{equation}

We apply windowing and multiply the integrand in (4.10) by the filter function $\varphi (t/T)$. We then discretise the resulting integral using the trapezoidal rule (noting that the endpoint contributions vanish) with step size ${\rm \Delta} t=T/N_{{ac}}$ to obtain the approximation

(4.13)\begin{align} \frac{{S}_{gg}(f)}{2{\rm \pi}{\rm \Delta} t}\approx\sum_{n={-}N_{{ac}}}^{N_{{ac}}} \varphi\left(\frac{n}{N_{{ac}}}\right)\frac{R_{gg}(n{\rm \Delta} t)}{2{\rm \pi}} \, {\rm e}^{{\rm i}n(2{\rm \pi} f{\rm \Delta} t)}=\sum_{n={-}N_{{ac}}}^{N_{{ac}}} \varphi\left(\frac{n}{N_{{ac}}}\right)\widehat{\nu_g}(n)\, {\rm e}^{{\rm i}n(2{\rm \pi} f{\rm \Delta} t)}. \end{align}

It follows that $\nu _{g,N_{{ac}}}(2{\rm \pi} f{\rm \Delta} t)$ can be understood as a discretised version of ${S}_{gg}(f)/(2{\rm \pi} {\rm \Delta} t)$ over the time window $[-T,T]$. Taking the limit $N_{{ac}}\rightarrow \infty$ with $(2{\rm \pi} {\rm \Delta} t) f=\theta$, we see that $\nu _g$ is an appropriate limit of the power spectrum with time resolution ${\rm \Delta} t$. There are two key benefits of using $\nu _g$. First, we do not explicitly periodically extend the signal to compute autocorrelations and thus avoid the problem of broadening. Instead, we window in the frequency domain. Second, we have rigorous convergence theory as $N_{{ac}}\rightarrow \infty$. We compare spectral measures with power spectra in § 6.

5.1. Computational set-up

The flow around a circular cylinder of diameter $D$ is obtained using an incompressible, two-dimensional lattice-Boltzmann computational fluid dynamics (CFD) flow solver. The solver uses the two-dimensional 9-velocity lattice model (commonly referred to as D2Q9) lattice and BGKW (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954) collision models to calculate the velocity field. The temporal resolution (time step) of the flow simulations is such that the Courant–Friedrichs–Lewy number is unity on the uniform numerical grid. However, this results in a very finely resolved flow field and hence a large volume of data. An initial down-sampling study revealed that storing 12 snapshots of flow field data within a period of vortex shedding still enables us to use our analysis tools without affecting the results. We also verified our results against higher grid resolution for the CFD solver. The computational domain size is $18D$ in length and $5D$ in height. The cylinder is positioned $2D$ downstream of the inlet at the mid-height of the domain. The cylinder walls and the side walls are defined as bounce-back no-slip walls, and a parabolic velocity inlet profile is defined at the inlet of the domain such that $Re=100$. The outlet is defined as a non-reflecting outflow. For a detailed description of the solver, we refer the reader to Józsa et al. (Reference Józsa, Szőke, Teschner, Könözsy and Moulitsas2016) and Szőke et al. (Reference Szőke, Jozsa, Koleszár, Moulitsas and Könözsy2017).

6.1. Experimental set-up

We consider two cases of turbulent boundary layer flow, a canonical (baseline) and one with an embedded shear layer. Both cases are generated experimentally and detailed in Szőke, Fiscaletti & Azarpeyvand (Reference Szőke, Fiscaletti and Azarpeyvand2020). For clarity, we briefly recall the key features and illustrate the experimental set-up in figure 9. The baseline case consists of a canonical, zero-pressure-gradient turbulent boundary layer (TBL) passing over a flat plate (no flow control turned on). The friction Reynolds number is $Re_\tau = 1400$. In the experiments, the development of a shear layer is triggered using perpendicular flow injection (active flow control) through a finely perforated sheet, denoted by $u_{AFC}$. The velocity field is sensed at several positions downstream of the flow control area by traversing a hot-wire probe across the entire boundary layer. The data gathered using the hot-wire sensor provides a fine temporal description of the flow. For a more detailed description of the experimental set-up and the flow field, see Szőke et al. (Reference Szőke, Fiscaletti and Azarpeyvand2020). A wide range of downstream positions and flow control severity values were considered in the original study to assess flow control's effects on trailing edge noise. Here, the same data are used, but we focus our attention on one streamwise position (labelled as BL3 in figure 9) and consider a flow injection case where a shear layer develops as a result of perpendicular flow injection. Similar qualitative comparisons (against power spectral densities) were found for different sensor locations and different injection angles.

Figure 9.

Schematic of the experimental set-up for a boundary layer with embedded shear profile, and side view of the boundary layer near the flow control region.

7.1. Experimental set-up

Experiments using TR-PIV are performed at the Wall Jet Wind Tunnel of Virginia Tech as schematically shown in figure 13. For a detailed description of the facility, we refer to Kleinfelter et al. (Reference Kleinfelter, Repasky, Hari, Letica, Vishwanathan, Organski, Schwaner, Alexander and Devenport2019). A two-dimensional two-component TR-PIV system is used to capture the wall-jet flow, and the streamwise origin of the field-of-view (FOV) is $\hat {x} = 1282.7$ mm downstream of the wall-jet nozzle. We use a jet velocity of $U_j=50$ m s$^{-1}$, corresponding to a jet Reynolds number of $Re_{jet} = h_{jet} U_j /\nu = 63.5 \times 10^3$. The length and height of the FOV are approximately $75 {\rm mm} \times 40\ {\rm mm}$, and the spatial resolution of the velocity vector field is 0.24 mm. This resolution corresponds to dimension $d = 102,300$ in (1.1). The high-speed cameras are operated in a double frame mode, with a rate of 12 000 frame pairs per second, resulting in a fine temporal resolution of 0.083 ms. For a more detailed description of the experimental set-up and flow field, see Kleinfelter et al. (Reference Kleinfelter, Repasky, Hari, Letica, Vishwanathan, Organski, Schwaner, Alexander and Devenport2019) and Szőke et al. (Reference Szőke, Nurani Hari, Devenport, Glegg and Teschner2021).

Figure 13.

The schematic of the TR-PIV experiments conducted in the Wall Jet Wind Tunnel of Virginia Tech.

The associated flow has some special properties. It is self-similar, and its main characteristics (boundary layer thickness, edge velocity, skin friction coefficient, etc.) can be accurately calculated a priori through power-law curve fits (Kleinfelter et al. Reference Kleinfelter, Repasky, Hari, Letica, Vishwanathan, Organski, Schwaner, Alexander and Devenport2019). The flow consists of two main regions. Within the region bounded by the wall and the peak in the velocity profile, the flow exhibits the properties of a zero-pressure-gradient TBL. Above this fluid portion, the flow is dominated by a two-dimensional shear layer consisting of rather large, energetic flow structures. While the peak in the velocity profile is $y_m \approx 18$ mm from the wall in our case, the overall thickness of the wall-jet flow is of the order of 200 mm. The PIV experiments must compromise between a good spatial resolution and capturing the entire flow field. In our case, the FOV was not tall enough to capture the entire wall-jet flow field. For this reason, the standard DMD algorithm under-predicts the energies corresponding to the shear-layer portion of the wall-jet flow as the corresponding length scales fall outside of the limits of the FOV. We also verified our results, in particular the presence of the modes in figure 16, by comparing with different spatial and temporal resolutions.

8.1. Experimental set-up

We investigate a near-ideal acoustic monopole source that is generated using the laser-optical set-up illustrated in figure 18. When a high-energy laser beam is focused on a point, the air ionises and plasma is generated due to the extremely high electromagnetic energy density (of the order of $10^{12}$ W m$^{-2}$). As a result of the sudden deposit of energy, the air volume undergoes a sudden expansion that generates a shockwave. The initial propagation characteristics can be modelled using von Neumann's point strong explosion theory (Von Neumann Reference Von Neumann1963), which was originally developed for nuclear explosion modelling. For our ResDMD analysis, we use laser-induced plasma (LIP) sound signature data measured using a 1/8 inch, Bruel & Kjaer (B&K) type 4138 microphone operated using a B&K Nexus module Szőke et al. (Reference Szőke, Devenport, Borgoltz, Alexander, Hari, Glegg, Li, Vallabh and Seyam2022). The data from the microphone are acquired using an NI-6358 module at a sampling rate of $f_s=1.25$ MS s$^{-1}$ (million samples per second). With this apparatus, we can resolve the high-frequency nature of the LIP up to 100 kHz. For a detailed description of the laser-optical set-up, experimental apparatus, and data processing procedures, see Szőke et al. (Reference Szőke, Devenport, Borgoltz, Alexander, Hari, Glegg, Li, Vallabh and Seyam2022) and Szőke & Devenport (Reference Szőke and Devenport2021).

The important acoustic characteristic of the LIP is a short time period of initial supersonic propagation speed, shown as Schlieren images taken over a $15\ \mathrm {\mu }$s window in figure 19. When observed from the far field, this initial supersonic propagation is observed as a nonlinear characteristic even though the wave speed is supersonic only until approximately 3–4 mm from the optical focal point. During the experiments, 65 realisations of LIP are captured using microphones. Each realisation of LIP is then gated in time such that only the direct propagation path of the LIP remains in the signal (Szőke & Devenport Reference Szőke and Devenport2021). We use this gated data for our ResDMD analysis.

Figure 19.

Schlieren images of the initial LIP illustrating the shock wave formation and propagation; (a) $t=5\ \mathrm {\mu }{\rm s}$, (b) $t=10\ \mathrm {\mu }{\rm s}$, (c) $t=15\ \mathrm {\mu }{\rm s}$ and (d) $t=20\ \mathrm {\mu }{\rm s}$.