2.1.1. The effect of noise

Two common scenarios where equation inference would be of particular value are when the data are generated experimentally (Reinbold et al. Reference Reinbold, Kageorge, Schatz and Grigoriev2021; Joshi et al. Reference Joshi, Ray, Lemma, Varghese, Sharp, Dogic, Baskaran and Hagan2022; Golden et al. Reference Golden, Grigoriev, Nambisan and Fernandez-Nieves2023) or when the data represent coarse graining of the results of direct numerical simulation. An example of the latter is fully kinetic simulations of plasma used to obtain a hydrodynamic description (Alves & Fiuza Reference Alves and Fiuza2022). In both instances, the data will inevitably be noisy, e.g. due to measurement inaccuracies in experiment or fluctuations of the computed coarse-grained fields. To investigate the effects of noise, in addition to the original simulation data downloaded from the turbulence database, we also used synthetic data with varying levels of additive uniform noise. Specifically, we define the noisy data $f_\sigma = f + \sigma \xi _f s_f,$ where $f\in \{\,p,u_x,u_y,u_z\}$ are the hydrodynamic fields, $\sigma$ is the noise level and $\xi _f$ is noise independently sampled from the uniform distribution over $[-1,1]$ at each space–time point.

Parsimonious scalar and vector relations describing velocity and pressure data can be identified by performing sparse regression using the libraries $\mathcal {L}_0$ and $\mathcal {L}_1$, respectively. In the strong form, the terms involving higher-order derivatives, such as $\nabla ^2p$ and $\nabla ^2{\boldsymbol u}$, will be extremely sensitive to noise (Rudy et al. Reference Rudy, Brunton, Proctor and Kutz2017; Reinbold & Grigoriev Reference Reinbold and Grigoriev2019). To make the regression more robust, we use the weak form of both PDEs following the approach introduced in our earlier work (Gurevich et al. Reference Gurevich, Reinbold and Grigoriev2019). Specifically, we multiply each equation by a smooth weight function $w_j({\boldsymbol x},t)$ and then integrate it over a rectangular spatio-temporal domain $\varOmega _i$ of size $H_x\times H_y\times H_z\times H_t$. All side lengths of $\varOmega _i$ are fixed to $H_i=32$ grid points of the numerical grid; this size roughly corresponds to the characteristic length and time scales of the flow field (cf. figure 1).

The derivatives are shifted from the data (${\boldsymbol u}$ and $p$) onto the weight functions $w$ whenever possible via integration by parts, after which the integrals are evaluated numerically using trapezoidal quadratures. (In the few cases where it is not possible to fully integrate a term by parts, remaining derivatives are evaluated using finite differences.) For the scalar library $\mathcal {L}_0$, we use scalar weight functions of the form

(2.6)\begin{equation} \left.\begin{array}{c@{}} w({\boldsymbol x},t) = \tilde{w}(\bar{x})\tilde{w}(\kern0.09em \bar{y})\tilde{w}(\bar{z})\tilde{w}(\bar{t}),\\ \tilde{w}(s) = (1-s^2)^\beta, \end{array}\right\} \end{equation}

where the bar denotes non-dimensionalization using the order-preserving affine map $[x_{min},x_{max}] \to [-1,1]$. In the case of the vector library $\mathcal {L}_1$, we instead use vector weight functions $w({\boldsymbol x}, t){\boldsymbol e}_k$ aligned along each of the coordinate axes. Note that, for the term libraries considered in this problem which involve a Laplacian of $p$ or ${\boldsymbol u}$, we should have $\beta \ge 2$, as this allows us to discard the boundary terms generated during integration by parts. We set $\beta = 8$ in our analysis since this choice (i) ensures that all the boundary terms vanish and (ii) maximizes the accuracy of numerical quadrature along the uniformly gridded dimensions (Gurevich et al. Reference Gurevich, Reinbold and Grigoriev2019). The non-uniform grid in the $y$-direction will control the error of the quadrature and increasing $\beta$ further has no benefit. This is illustrated in figure 2, which shows how the quadrature error scales with $\beta$ for a test function $\cos (2{\rm \pi} y)$. The error decreases quickly with increasing $\beta$ for the uniform grid, plateauing at $\beta \ge 8$. In contrast, the error is found to be almost independent of $\beta$ for the non-uniform grid.

Figure 2.

Numerical error $\Delta I$ of the integral $I_\beta = \int _{-1}^1 \cos (2{\rm \pi} y) (1-y^2)^\beta \, {\textrm {d} y}$ evaluated using the trapezoid rule on a 32-point grid. The results are shown for the uniform grid (black squares) and the non-uniform grid (white squares) representing the $y$-coordinates of points taken from the middle of the channel.

By repeating this procedure for a number of integration domains $\varOmega _i$ contained within the full dataset, we construct a feature matrix ${Q} = [{\boldsymbol q}_1 \cdots {\boldsymbol q}_N]$ whose the columns ${\boldsymbol q}_n$ correspond to the various terms in either $\mathcal {L}_0$ or $\mathcal {L}_1$. For instance, for the scalar library $\mathcal {L}_0$,

(2.7a,b)\begin{equation} Q_{ij} = \frac{1}{V_i S_j}\int_{\varOmega_i} w({\boldsymbol x},t)\, f_j({\boldsymbol x},t)\, {\rm d}^3{\boldsymbol x}\, {\rm d}t,\quad {V_i=\int_{\varOmega_i}|w({\boldsymbol x},t)|\, {\rm d}^3{\boldsymbol x}\, {\rm d}t,} \end{equation}

where $S_j = S[\,f_j]$ is the scale of the library term $f_j$ in dimensional units. We estimate the scales of all library terms using the characteristic time scales $T_{\boldsymbol u},T_p$, length scales $L_{\boldsymbol u}, L_p$, mean values $\mu _{\boldsymbol u}, \mu _p$ and standard deviations $\sigma _{\boldsymbol u}, \sigma _p$ of the velocity and pressure fields across the dataset. For instance, for a vector field ${\boldsymbol u}$, the length and time scales can be estimated from the magnitudes of its finite-difference derivatives

(2.8)$$\begin{gather} T_{{\boldsymbol u}} \equiv \frac{\sigma_{\boldsymbol u}}{ \sqrt{ \langle \partial_t {\boldsymbol u} \boldsymbol{\cdot} \partial_t {\boldsymbol u} \rangle} }, \end{gather}$$

(2.9)$$\begin{gather}L_{{\boldsymbol u}} \equiv \frac{\sigma_{\boldsymbol u}}{ \sqrt{ \langle (\boldsymbol{\nabla} {\boldsymbol u}) \boldsymbol{\cdot} (\boldsymbol{\nabla} {\boldsymbol u}) \rangle}}, \end{gather}$$

where the dot products in the denominators are taken over all the indices. We then use the heuristic that the scale of a library term is the product of the scales of its factors: $S[{\boldsymbol u \boldsymbol v}] = S[{\boldsymbol u}] S[{\boldsymbol v}]$, and for $k \geq 1$

(2.10a,b)$$\begin{gather} S[\,p] = \mu_p, \quad S[{\boldsymbol u}] = \mu_{\boldsymbol u}, \end{gather}$$

(2.11a,b)$$\begin{gather}S[\nabla^k p] = L_p^{-k} \sigma_p, \quad S[{\nabla^k {\boldsymbol u}}] = L_{\boldsymbol u}^{-k} \sigma_{\boldsymbol u}, \end{gather}$$

(2.12a,b)$$\begin{gather}S[{\partial_t^k p}] = T_p^{-k}, \sigma_p \quad S[{\partial_t^k {\boldsymbol u}}] = T_{\boldsymbol u}^{-k} \sigma_{\boldsymbol u}. \end{gather}$$

Note that the mass scale does not appear explicitly, as the data are given in units in which the density $\rho =1$. The non-dimensionalization procedure ensures the magnitudes of all columns are comparable, which can dramatically improve the accuracy and robustness of regression. The problem of determining the unknown coefficients ${\boldsymbol c}=[c_1,\ldots,c_N]^\top$ is cast as the solution of an overdetermined linear system of the form

(2.13)\begin{equation} {\boldsymbol{\mathsf{Q}}}{\boldsymbol c} = 0. \end{equation}

In this article, we sample $\varOmega _i$ from a $64^4$-point region of the data lying either in the middle (i.e. the symmetry plane of the flow) or at the edge of the channel (see table 5). Examples of both types of domains are shown in figure 1. We use 256 randomly sampled integration domains $\varOmega _i$ with $32^4$ gridpoints to construct the system (2.13). This yields $M=256$ linear equations on the coefficients of $\mathcal {L}_0$ and $3M=768$ linear equations on the coefficients of $\mathcal {L}_1$. A discussion of how the results are affected by correlated noise, the number of integration domains used and the resolution of the data can be found in Appendix B.

2.1.2. Selection of parsimonious relations

Note that the linear system (2.13) is homogeneous and treats all terms in the library on equal terms. This is in contrast to SINDy (Brunton et al. Reference Brunton, Proctor and Kutz2016) and its variants (Messenger & Bortz Reference Messenger and Bortz2021) that solve an inhomogeneous linear system, or many such systems in the case of SINDy-PI (Kaheman et al. Reference Kaheman, Kutz and Brunton2020). Its solutions have a degree of freedom corresponding to the normalization of ${\boldsymbol c}$, which can be eliminated by arbitrarily setting one of the coefficients, say $c_1$, to unity, as done in SINDy, or by fixing the norm of ${\boldsymbol c}$ as in this study. The solutions of a constrained least squares problem

(2.14)\begin{equation} {\boldsymbol c}=\mathop{\rm arg\,min}\limits_{\|{\boldsymbol c}\|=1}\|{\boldsymbol{\mathsf{Q}}}{\boldsymbol c}\|,\end{equation}

are given by the right singular vector of $\boldsymbol{\mathsf{Q}}$ corresponding to the smallest singular value. It is worth noting that, when multiple singular values of $\boldsymbol{\mathsf{Q}}$ are small, there may be several ‘good’ independent solutions for ${\boldsymbol c}$ representing different dominant balances. This is a less restrictive approach compared with SINDy and allows a broader class of functional relations to be identified. It is also more computationally efficient than SINDy-PI, which aims to address the same limitation.

In order to obtain a parsimonious physical relation, we must find a sparse coefficient vector ${\boldsymbol c}^*$ such that the residual $\|\boldsymbol{\mathsf{Q}}{\boldsymbol c}^*\|$ is comparable to the residual $\|\boldsymbol{\mathsf{Q}}{\boldsymbol c}\|$ with dense ${\boldsymbol c}$ given by (2.14). The identified relations either contain a single term or several terms. If the matrix $\boldsymbol{\mathsf{Q}}$ has been properly non-dimensionalized, single-term relations will correspond to columns with small norms. We will use the heuristic that $f_j = 0$ is a valid single-term relation if $\|{\boldsymbol q}_j\| \ll \sqrt {M}$, where $M$ is the number of integration domains. In particular, for the data without added noise, the scalar library $\mathcal {L}_0$ is found to contain terms with $\| {\boldsymbol q}_j \| \approx 10^{-6} \sqrt {M}$, which correspond to the incompressibility condition

(2.15)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol u}=0, \end{equation}

and its trivial corollary $p \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol u} = 0$. Both single-term relations are found using data from the middle of the channel as well as data near the boundary. The single-term relation heuristic crucially relies on proper non-dimensionalization such that velocity gradients are $O(1)$. A more general approach is to compare with the characteristic size of the uncontracted tensor, which in this case is the rate of strain $\nabla _i u_j$. In contrast, prior sparsification algorithms such as SINDy (Brunton et al. Reference Brunton, Proctor and Kutz2016), implicit SINDy (Mangan et al. Reference Mangan, Brunton, Proctor and Kutz2016) and SINDy-PI (Kaheman et al. Reference Kaheman, Kutz and Brunton2020) are unable to identify single-term relations, as these are not examined separately. Note that direct identification of single-term relations is both more robust and more computationally efficient than identification through regression.

Once the library has been pruned, multiple-term relations can be identified by an iterative greedy algorithm. At each iteration, we use the singular value decomposition of ${Q}^{(N)} = [{\boldsymbol q}_1 \cdots {\boldsymbol q}_N]$ to find ${\boldsymbol c}^{(N)}$ as described previously. We also compute the residual $r^{(N)} = \|{Q}^{(N)} {\boldsymbol c}^{(N)}\|$. Next, we consider all of the candidate relations formed by dropping one of the terms and eliminating the corresponding column from ${Q}^{(N)}$. We select the candidate relation with $N-1$ terms that achieves the smallest residual and then repeat until only one term remains. This yields a sequence of increasingly sparse relations described by $N$-dimensional coefficient vectors ${\boldsymbol c}^{(N)}$, forming an approximately Pareto-optimal set (Miettinen Reference Miettinen2012). Note that the use of the absolute residual $r$ guarantees that the residual is a monotonic function of the number of terms $N$, which is not the case for the relative residual $\eta =r/\max _n\|c_n{\boldsymbol q}_n\|$ used by Reinbold et al. (Reference Reinbold, Kageorge, Schatz and Grigoriev2021) and Golden et al. (Reference Golden, Grigoriev, Nambisan and Fernandez-Nieves2023).

There are many reasonable ways to select a final relation from this sequence based on the trade-off between their parsimony (i.e. number of terms $N$) and accuracy, quantified by the residuals $r^{(N)}$. For instance, one might select the simplest relation which achieves a relative residual of less than, say, $1\,\%$ or the relation for which discarding a single term results in the largest relative increase in the residual. In this article, we follow Gurevich et al. (Reference Gurevich, Reinbold and Grigoriev2019): specifically, we choose the relation described by the coefficient vector ${\boldsymbol c}^{(N)}$ where $N = \max \{n : r^{(n)}/r^{(n-1)} > \gamma \}$, where the parameter $\gamma = 1.25$ was selected empirically. Once the functional form of a parsimonious relation is determined, the mean values of the coefficients and their standard deviations are computed by subsampling $\boldsymbol{\mathsf{Q}}\ 128$ times, with each member of the ensemble constructed using one half of the rows of $\boldsymbol{\mathsf{Q}}$ selected at random. Finally, to enhance the interpretability of the result, the equation is rescaled by setting the largest coefficient to unity, which defines the mean values $\bar {c}_n$ and their uncertainties $s_n$.

After a sparse single-term or multi-term relation has been identified, one may search for additional sparse relations contained within the same library. Note that the form of our libraries implies a simple connection between a relation and its direct algebraic implications: if ${\boldsymbol c}\boldsymbol {\cdot } {\boldsymbol f}=0$, then ${\boldsymbol c}\boldsymbol {\cdot } (g\,{\boldsymbol f})=0$ and ${\boldsymbol c}\boldsymbol {\cdot } (\partial _s \,f) = 0$ for any term $g$ in any of the libraries and set of partial derivatives $s$. Moreover, iteratively applying these rules produces all implied relations (with a higher complexity) for a given base relation (of a lower complexity). Each implied relation is an equation which allows one of the terms, say the highest-order one, to be eliminated from the respective library without loss of expressivity. (In principle, this procedure could reduce sparsity of future identified equations if applied to many-term relations, but in such cases we find that the most complex term is unlikely to reappear in another independent equation.)

This can be leveraged to devise a simple and efficient algorithm for finding all relations contained within each library $\mathcal {L}_m$. Consider a nested sequence of sub-libraries: $\mathcal {L}^{(1)}_m$ containing terms up to first order, $\mathcal {L}^{(2)}_m$ containing terms up to second order and so on. On each sub-library, we repeatedly run single-term and multi-term regressions. Whenever a new relation is identified, we construct all implied relations in each library and then eliminate from its appropriate library the highest-order term from the base relation as well as from each implied relation before re-running the regression. For instance, the discovery of the incompressibility condition $\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol u} = 0$ from $\mathcal {L}^{(2)}_0$ would lead to the identification of the implied relations $p \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol u} = 0$, ${\boldsymbol u}(\boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol u}) = 0$, and $\boldsymbol {\nabla }( \boldsymbol {\nabla } \boldsymbol {\cdot } {\boldsymbol u}) = 0$ and the elimination of these terms from $\mathcal {L}^{(3)}_0$ and $\mathcal {L}^{(3)}_1$. If no more relations remain after eliminating terms, we advance to the next sub-library (either the same-order sub-library from the next library, or after all sub-libraries of a given order have been exhausted, the next-order sub-library from the first library), continuing until each library has been completely examined. This procedure can be used to identify, for instance, the pressure-Poisson equation and the energy equation from a scalar library of sufficiently high order, as discussed below.

2.1.3. Identified equations and robustness to noise

Regression performed using the vector library $\mathcal {L}^{(3)}_1$ identifies a single relation representing momentum balance

(2.16)\begin{equation} c_1\partial_t {\boldsymbol u} + c_2({\boldsymbol u}\boldsymbol{\cdot} \boldsymbol{\nabla}){\boldsymbol u} + c_3\boldsymbol{\nabla} p + c_4 \nabla^2 {\boldsymbol u}= 0. \end{equation}

Table 1 lists the surviving terms $f_n$, the corresponding coefficients $c_n$, their uncertainties $s_n$ and the respective magnitudes $\chi _n=\|c_n{\boldsymbol q}_n\|$ of the terms which can be used to identify dominant balances in different regions. For data from the edge of the channel, the Navier–Stokes equation with accurate coefficients, including the small viscosity, is identified for all noise levels. In particular, for noiseless data, the magnitude of the viscosity is identified correctly to several significant digits, even though the viscous term involves a second-order derivative, which is the highest in the equation. The corresponding sequence of residuals $r^{(N)}$ can be found in figure 3(c).

Table 1.

Coefficients of the momentum equation (2.16) in the presence of varying levels of noise and data from (a) the edge of the channel and (b) the middle of the channel. The rows show the mean values of the coefficients $\bar {c}_n$ (normalized by the magnitude of the largest one), their uncertainties $s_n$ and the magnitudes of the terms $\chi _n$ (normalized by the magnitude of the largest term).

Figure 3.

Dependence of the residual $r$ on the number of terms $N$ retained in a given relation in the noiseless case. Black (white) squares represent data collected near the edge (in the middle) of the channel. The identified relations are (a) the energy equation, (b) the pressure equation and (c) the momentum equation.

For data from the middle of the channel, different special cases of (2.16) are identified for different noise levels. For noise levels up to approximately 15 %, sparse regression still identifies the Navier–Stokes equation. For higher noise levels (up to 100 %), the Euler equation is identified instead, also with fairly accurate coefficients. For extreme levels of noise (up to 300 %), the inviscid Burgers equation is identified. Note that the sequence in which the terms in the momentum equation stop being recovered, as the noise level is increased, corresponds to their magnitudes $\chi _n$ in the middle of the channel.

No multi-term relations are found, for any data set, from the scalar library $\mathcal {L}^{(3)}_0$, suggesting it is incomplete. In order to identify any multi-term relations, this library was expanded by including terms quartic in ${\boldsymbol u}$ and/or $\boldsymbol {\nabla }$

(2.17)\begin{align} \mathcal{L}^{(4)}_0=\mathcal{L}^{(3)}_0\cup\{\nabla^2 (u^2),\boldsymbol{\nabla} \boldsymbol{\cdot} [ ({\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}) {\boldsymbol u} ],(\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol u})^2,\boldsymbol{\nabla} {\boldsymbol u}: \boldsymbol{\nabla} {\boldsymbol u}^\top,\boldsymbol{\nabla} {\boldsymbol u}:\boldsymbol{\nabla} {\boldsymbol u},\boldsymbol{\nabla} \boldsymbol{\cdot} ( u^2{\boldsymbol u}),{\boldsymbol u}^4\}.\end{align}

Note that this expanded library is not exhaustive, i.e. it does not include any fourth-order terms which involve either $p$ or $\partial _t$. Whenever possible, the additional terms were written in conservative form to decrease numerical error associated with evaluation of higher-order derivatives in weak form. Two relations are identified from this expanded scalar library via sparsification

(2.18)$$\begin{gather} c_1 \partial_t E + c_2 \boldsymbol{\nabla} \boldsymbol{\cdot} ({\boldsymbol u}E) + c_3 {\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla} p + c_4\nabla^2E +c_5 \boldsymbol{\nabla} {\boldsymbol u} : \boldsymbol{\nabla} {\boldsymbol u} = 0, \end{gather}$$

(2.19)$$\begin{gather}c_6\nabla^2 p + c_7\boldsymbol{\nabla} \boldsymbol{\cdot} [({\boldsymbol u} \boldsymbol{\cdot} \boldsymbol{\nabla}) {\boldsymbol u}] +c_8 = 0, \end{gather}$$

where $E=u^2/2$ is the energy density. The two relations are discovered robustly for both data sampled from the edge and the middle of the channel; however, the order in which they are found depends on the sampled region. The corresponding sequences of residuals $r^{(N)}$ can be found in figure 3(a,b), and the values of the coefficients $c_n$ for different noise levels are summarized in tables 2 and 3.

Table 2.

Coefficients of the energy equation (2.18) in the presence of varying levels of noise and data from (a) the edge and (b) the middle of the channel. The quantities $\bar {c}_n$, $s_n$, and $\chi _n$ are defined in the caption of table 1.

Table 3.

Coefficients of the pressure equation (2.19) in the presence of varying levels of noise and data from (a) the edge and (b) the middle of the channel. The quantities $\bar {c}_n$, $s_n$ and $\chi _n$ are defined in the caption of table 1.

For sufficiently low levels of noise (below approximately 10 %), relation (2.18) corresponds to the well-known energy equation, again with fairly accurate coefficients, no matter which region of the flow the data comes from. The accuracy of the coefficients decreases somewhat as the level of noise is increased, as expected. For data from the middle of the channel, small terms such as $\boldsymbol {\nabla }{\boldsymbol u}:\boldsymbol {\nabla }{\boldsymbol u}$ start to disappear as the noise level is increased, similar to what we found for the momentum equation. Relation (2.19) takes the form of the pressure-Poisson equation for moderately noisy data from both the middle and the edge of the channel, also with fairly accurate coefficients. The surprising observation is that, for noiseless data from the middle of the channel, sparse regression reliably identifies a small correction to the pressure-Poisson equation, a term proportional to unity. We will therefore refer to relation (2.19) simply as the pressure equation.

Note that neither the energy equation nor the pressure-Poisson equation is an independent relation; both can be derived from the Navier–Stokes equation and the incompressibility condition. Further, note that SPIDER is superior to alternative approaches to equation inference in both versatility and accuracy: for instance, neither of the two scalar relations would be identified using SINDy, which assumes either $\partial _t{\boldsymbol u}$ or $\partial _tp$ to be present. The energy equation has the form of an evolution equation but involves a temporal derivative of $u^2$ rather than ${\boldsymbol u}$, while the pressure-Poisson equation is a constraint which involves no temporal derivatives at all. Finally, note that one can also heuristically identify the incompressibility condition, the pressure-Poisson equation and the energy equation from the right singular vectors of $\boldsymbol{\mathsf{Q}}$ corresponding to the three smallest singular values without pruning $\mathcal {L}^4_0$. These singular vectors are dense, but only the coefficients associated with the corresponding dominant balances are $O(1)$.

For both sampled regions, all libraries and all noise levels, the residual $r$ asymptotes to a constant value for large $N$. In the noiseless case shown in figure 3, the asymptotic value of the residual is determined by the discretization of the data, both in the numerical simulations and in evaluating the integrals using quadratures, as shown by Gurevich et al. (Reference Gurevich, Reinbold and Grigoriev2019). In the noisy case, shown in figure 4, the asymptotic value of the residual is instead determined by the level of noise and is higher than in the noiseless case, as expected. It should be emphasized that physically meaningful relations can be identified in the presence of very high levels of noise, illustrating the robustness of SPIDER.

Figure 4.

Dependence of the residual $r$ on the number of terms $N$ retained in the momentum equation for synthetic data with added noise for noise levels of (a) 10 % and (b) 100 %. The Navier–Stokes equation is identified in all the cases except for data from the middle of the channel with 100 % noise. In that case, the Euler equation is found instead. The residual for the pressure equation identified from data with 20 % noise (c) and for the energy equations identified from data with 10 % noise (d). Black (white) squares represent data collected near the edge (in the middle) of the channel.

The power of the weak formulation compared with the strong form is vividly illustrated by other physical relations containing Laplacians as well. For instance, the vorticity equation can only be identified correctly in strong form for noise levels up to approximately 1 %, as shown by Rudy et al. (Reference Rudy, Brunton, Proctor and Kutz2017). The pressure-Poisson equation (2.19), which also contains second-order derivatives, can be identified in weak form for noise levels up to 20 % (500 %) for data from the middle (edge) of the domain. The energy equation (2.18) can be identified in weak form using the data from the edge of the channel with up to 10 % noise and using data from the middle of the channel with up to $3\,\%$ noise. The latter equation is more challenging to identify from noisy data due to the presence of the term $\boldsymbol {\nabla }{\boldsymbol u} : \boldsymbol {\nabla } {\boldsymbol u}$ which cannot be integrated by parts and was computed using second-order accurate finite differencing. However, it should be noted that a more sophisticated procedure for numerically differentiating the data would likely allow this equation to be identified even at higher noise levels.

Note that, although we used the absolute residual $r$ in sparse regression, it is the magnitude of the relative residual $\eta$ that better quantifies the accuracy of the identified relations (Reinbold et al. Reference Reinbold, Kageorge, Schatz and Grigoriev2021). The relative residuals are indeed quite low: $4 \times 10^{-4} (2 \times 10^{-2})$ for the energy equation, $2 \times 10^{-3} (1 \times 10^{-4})$ for the two-term pressure-Poisson equation and $4 \times 10^{-5} (5 \times 10^{-5})$ for the Navier–Stokes equation identified using noiseless data from the middle (edge) of the channel. The three-term pressure equation has a relative residual of $8 \times 10^{-4}$ in the middle of the channel, which is less than a half that of the pressure-Poisson equation.