2.1. The Mori–Zwanzig formalism

Different approaches leverage this past information. The Mori–Zwanzig (MZ) formalism, rooted in statistical physics, provides a framework to derive an evolution equation for a set of variables, such as macroscopic observables, associated with a high-dimensional dynamical system (Mori, Reference Mori1965; Zwanzig, Reference Zwanzig1966; Zwanzig et al., Reference Zwanzig, Nordholm and Mitchell1972). This framework is instrumental, for instance, in situations where the full state $ \mathbf{u}(t) $ is unavailable, and where one can only access low-dimensional observations. Similarly, the MZ formalism is also relevant for addressing dimension reduction problems (Zhu, Reference Zhu2019) where the reduced state is a set of variables whose dynamics are affected by the projected-out and nonobserved variables or state components.

In the Mori–Zwanzig framework, the span of $ \mathcal{G} $ of the measurements is considered as resulting from a projection through a projector $ \mathcal{P} $ of the space of $ {L}^2\left(\mathcal{S},\unicode{x211D}\right) $ functions. Similarly, the unobserved functions span a subspace resulting from a projection $ \mathcal{Q} $ , with $ \mathcal{Q}=\mathcal{I}-\mathcal{P} $ . The dynamics of $ \mathbf{y} $ can then be expressed as

(2.3) $$ \frac{\partial }{\partial t}\mathbf{y}\left(\mathbf{u}(t)\right)=\frac{\partial }{\partial t}{e}^{\left(t-{t}_0\right)\mathrm{\mathcal{L}}}\mathbf{y}\left({\mathbf{u}}_0\right)={e}^{\left(t-{t}_0\right)\mathrm{\mathcal{L}}}\left(\mathcal{P}+\mathcal{Q}\right)\mathrm{\mathcal{L}}\hskip0.1em \mathbf{y}\left({\mathbf{u}}_0\right). $$

The influence of past unresolved dynamics on the current observations can be explicitly expressed using the Dyson identity, which states:

(2.4) $$ {e}^{\left(t-{t}_0\right)\mathrm{\mathcal{L}}}={\int}_0^{t-{t}_0}{e}^{\left(t-{t}_0-s\right)\mathrm{\mathcal{L}}}\mathcal{P}\mathrm{\mathcal{L}}\hskip0.1em {e}^{s\mathcal{Q}\mathrm{\mathcal{L}}}\mathrm{d}s+{e}^{\left(t-{t}_0\right)\mathcal{Q}\mathrm{\mathcal{L}}}, $$

so that one can write the dynamics of the current observables $ \mathbf{y}\left(\mathbf{u}(t)\right) $ as:

(2.5) $$ \frac{\partial }{\partial t}\mathbf{y}\left(\mathbf{u}(t)\right)=\mathcal{P}\mathrm{\mathcal{L}}\hskip0.1em \mathbf{y}\left(\mathbf{u}(t)\right)+{\int}_0^{t-{t}_0}\mathcal{P}\mathrm{\mathcal{L}}\hskip0.1em {e}^{s\mathcal{Q}\mathrm{\mathcal{L}}}\mathcal{Q}\mathrm{\mathcal{L}}\hskip0.1em \mathbf{y}\left(\mathbf{u}\left(t-s\right)\right)\mathrm{d}s+{e}^{\left(t-{t}_0\right)\mathcal{Q}\mathrm{\mathcal{L}}}\mathcal{Q}\mathrm{\mathcal{L}}\hskip0.1em \mathbf{y}\left({\mathbf{u}}_0\right). $$

Rewriting Eq. (2.5) in a more compact form, and denoting $ \mathbf{y}\left(\mathbf{u}(t)\right)=\mathbf{y}(t) $ , the dynamics of the vector of observables $ \mathbf{y}=(\hskip0.15em {y}_1,\dots, {y}_p) $ are shown to follow an integro-differential equation (IDE):

(2.6) $$ {\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{\partial}\mathbf{y}(t)}{\mathrm{\partial}t}}=M\hskip0.1em \mathbf{y}(t)+{\int}_{\hskip-5pt 0}^{t-{t}_0}K(\mathbf{y}(t-s),s)\hskip0.1em \mathrm{d}s+F(t,{\mathbf{u}}_0),\end{array}} $$

with $ M $ a Markovian operator and $ K $ an operator applied to past observables and integrated over the whole time span since the initial condition. The so-called “noise” term $ F $ accounts for the impact of the unobserved variables at the initial condition. It can be eliminated from the equations predicting future observations by applying the projector $ \mathcal{P} $ and restricting their dynamics to the linear span of $ \mathcal{P}\mathbf{y} $ . Learning the dynamics of partially observed or known systems hence boils down to estimating each term of the differential equation above, sometimes referred to as the generalized Langevin equation (GLE).

Accounting for the past essentially allows us to isolate the dynamics of the observables. The Mori–Zwanzig framework is general and applies widely, allowing to describe the dynamics of the observables with a non-Markovian model. It similarly provides a principled closure for coarse models, hence missing information, which can be effectively complemented with a history-based term.

2.2. Approximations of integro-differential equations (IDE)

The statement in Eq. (2.6) outlines the structure of the dynamics. However, the operators involved in this dynamical equation are usually poorly known, if any, in particular the integral term. We now discuss several approaches from the literature to estimate these terms.

An approach to approximate the integral consists of studying particular asymptotic regimes. Among the many different models proposed in the literature (Stinis, Reference Stinis2004; Chorin and Stinis, Reference Chorin and Stinis2007; Stinis, Reference Stinis2007), one of the most popular consists in approximating the integral under assumptions of very short memory or very long memory regimes. For example, the t-model (Chorin et al., Reference Chorin, Hald and Kupferman2002), also commonly called slowly decaying memory approximation, leads to Markovian equations with time-dependent coefficients, which can subsequently be modeled using a NODE. However, these remain asymptotic approximations and cannot be extended to intermediate-range memory, in general. Another approach consists of performing Monte Carlo integration (Robert and Casella, Reference Robert and Casella1999). This approach has been extended to a neural network-based formulation (Neural IDE) (Zappala et al., Reference Zappala, de O. Fonseca, Moberly, Higley, Abdallah, Cardin and van Dijk2023), where the memory integrand is decomposed as a product $ {K}_1\left(t,s\right)\hskip0.1em {K}_2\left(\mathbf{y}(s)\right) $ . However, the number of function evaluations required to accurately integrate Equation (2.6) can be large, making the process computationally intensive (cf. Appendix A). In a similar spirit, under the assumption of short memory, one can restrict the integral to a short past and discretize it in time, leading to an equation of the form (Gallage, Reference Gallage2017):

(2.7) $$ {\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{\partial}\mathbf{y}(t)}{\mathrm{\partial}t}}\approx M(\mathbf{y}(t))+\frac{1}{n}\sum \limits_{i=1}^nK(\mathbf{y}(t-{\tau}_i),{\tau}_i),\end{array}} $$

using $ n $ delays $ {\left\{{\tau}_i\right\}}_i $ uniformly spaced over a past horizon instead of sampling them randomly. Such approximations can be improved using high-order discretization schemes, yet, as for Neural IDE, they require an unaffordable number of delays if the integrand varies quickly or if the time interval is too large.

2.3. Exact representation with neural DDE

While Equation (2.7) only provides an approximation of the true dynamics and requires many delays, we show that using a more complex function of a small number of delays, it is actually possible to represent the dynamics exactly. For this, let us consider diffeomorphic dynamical systems, that is, ODEs whose flow is invertible (in the full state space) and smooth in both time directions.

Proposition 2.1 (Exact representation with delays). For any smooth dynamical system ( $ {C}^2 $ is enough), and differentiable observables $ \mathbf{y} $ , there exists almost surely an operator $ M $ of the current observables, a finite number $ n $ of delays $ {\tau}_1,\dots, {\tau}_n>0 $ and a function $ \mathbf{h} $ such that the observables exactly follow the dynamics:

(2.8) $$ {\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{\partial}\mathbf{y}(t)}{\mathrm{\partial}t}}=M(\mathbf{y}(t))+\mathbf{h}(t,\mathbf{y}(t),\mathbf{y}(t-{\tau}_1),\mathbf{y}(t-{\tau}_2),\dots, \mathbf{y}(t-{\tau}_n)).\end{array}} $$

The derivation of the above, deferred to Appendix C, is based on the application of Takens’ theorem (hence the diffeomorphism requirement), which also provides a bound on the required number of delays: at least twice the intrinsic dimension of the manifold $ \mathcal{S} $ in which the full state $ \mathbf{u} $ lives. Note that this evolution equation is exact: approximations may arise from the optimization or the expressivity of the neural networks estimating $ M $ and $ \mathbf{h} $ , but not from the number of delays, provided they reach the Takens’ lower bound. This is in contrast with the discretization of the IDE integral as in Equation (2.7), which becomes asymptotically precise only when the number of delays is large compared to the complexity of the integrand.

3.1. Learning the delays

A constant lag NDDE is part of the larger family of continuous depth models that emerged with NODE (Chen et al., Reference Chen, Rubanova, Bettencourt, Duvenaud, Bengio, Wallach, Larochelle, Grauman, Cesa-Bianchi and Garnett2018). It is defined by:

(3.1) $$ {\displaystyle \frac{\mathrm{\partial}\mathbf{y}(t)}{\mathrm{\partial}t}}={\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),\mathbf{y}(t-{\tau}_1),\dots, \mathbf{y}(t-{\tau}_n)),\hskip2em \mathbf{y}(t\le 0)=\boldsymbol{\psi} (t), $$

with $ \boldsymbol{\psi} $ : $ \unicode{x211D}\to {\unicode{x211D}}^m $ the history function, $ {\tau}_i\in {\unicode{x211D}}^{+} $ a delay constant and $ {\mathbf{h}}_{\boldsymbol{\theta}}:\left[0,T\right]\times {\unicode{x211D}}^m\times \dots \times {\unicode{x211D}}^m\to {\unicode{x211D}}^m $ a neural network.

There are two possible ways of training continuous-depth models: discretize-then-optimize or optimize-then-discretize (Kidger, Reference Kidger2021). In the former, the numerical simulation library’s inherent auto-differentiation capabilities are leveraged. In the latter, the adjoint dynamics are employed to compute the gradient’s loss. The following proposition 3.1 provides the adjoint method for constant lag NDDEs.

Proposition 3.1. Let us consider the Neural DDE model below, where $ \tau $ is a learnable vector parameterized by some components of $ \boldsymbol{\theta} $ :

(3.2) $$ \frac{\partial \mathbf{y}(t)}{\partial t}={\mathbf{h}}_{\boldsymbol{\theta}}\left(t,\mathbf{y}(t),\mathbf{y}\left(t-\tau \right)\right),\hskip2em \mathbf{y}\left(t\le 0\right)=\boldsymbol{\psi} (t),\hskip2em 0\le t\le T, $$

and the following loss function:

$$ L(\mathbf{y})={\int}_{\hskip-5pt 0}^Tl(\mathbf{y}(t))\hskip0.1em \mathrm{d}t. $$

The gradient of the loss with respect to the parameters is given by:

(3.3) $$ {\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{\partial}L}{\mathrm{\partial}\boldsymbol{\theta }}}=-{\int}_{\hskip-5pt 0}^T\boldsymbol{\lambda} (t)\left({\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),\mathbf{y}(t-\tau ))}{\mathrm{\partial}\boldsymbol{\theta }}}-{\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),{\mathbf{y}}_{(\tau )}(t))}{\mathrm{\partial}{\mathbf{y}}_{(\tau )}}\frac{\mathrm{\partial}\mathbf{y}(t-\tau )}{\mathrm{\partial}\boldsymbol{\theta }}}\right)\mathrm{d}t,\end{array}} $$

where $ {\mathbf{y}}_{\left(\tau \right)}(t)\equiv \mathbf{y}\left(t-\tau \right) $ is the third set of variables $ {\mathbf{h}}_{\boldsymbol{\theta}}\left(t,\mathbf{y},{\mathbf{y}}_{\left(\tau \right)}\right) $ depends on. The adjoint dynamics $ \boldsymbol{\lambda} (t) $ are given by another DDE:

(3.4) $$ {\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{d}\boldsymbol{\lambda } (t)}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\partial}l(\mathbf{y}(t))}{\mathrm{\partial}\mathbf{y}}}-\boldsymbol{\lambda} (t){\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),{\mathbf{y}}_{(\tau )}(t))}{\mathrm{\partial}\mathbf{y}}}-\boldsymbol{\lambda} (t+\tau ){\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t+\tau, \mathbf{y}(t+\tau ),{\mathbf{y}}_{(\tau )}(t+\tau ))}{\mathrm{\partial}{\mathbf{y}}_{(\tau )}}},\\ {}\hskip2.5pc \boldsymbol{\lambda} (t\ge T)=\mathbf{0}.\end{array}} $$

The derivation of Proposition 3.1 is in Appendix B. Algorithm 1 below outlines the training procedure for a Neural DDE with multiple learnable constant delays, defined for convenience as $ {\mathbf{y}}_{\left({\tau}_i\right)}(t)=\mathbf{y}\left(t-{\tau}_i\right) $ .

The developments are packaged in a user-friendly API, developing a numerically robust DDE solver, and implementing the adjoint method in the torchdde package. These advancements allow a seamless integration of DDEs for future users, enhancing reproducibility. Computational details and benchmarks are reported in Appendix A.2.

Algorithm 1 Training a Neural DDE with learnable delays with the adjoint method.

Require: Dataset of one trajectory $ \mathcal{D}=\left\{\left({t}_0,{\mathbf{y}}_0\right),\dots, \left({t}_T,{\mathbf{y}}_T\right)\right\}. $

Require: Initialized model $ {\mathbf{h}}_{\boldsymbol{\theta}}. $

Require: Initialized positive delays $ {\tau}_1,\dots, {\tau}_n $ , handled as additional entries in the parameters vector $ \boldsymbol{\theta} $ .

1: for $ {N}_{\mathrm{epochs}} $ do

2: Set $ {\tau}_{\mathrm{max}}=\max \left\{{\tau}_1,\dots, {\tau}_n\right\} $

3: Create history function interpolation $ \psi $ with data from $ \mathcal{D} $ such that $ t\le {t}_0+{\tau}_{\mathrm{max}} $ .

4: Solve DDE dynamics:

5: $ \left\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}\mathbf{y}(t)}{\mathrm{d}t}}={\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),{\mathbf{y}}_{({\tau}_1)}(t),\dots, {\mathbf{y}}_{({\tau}_n)}(t))\\ {}\mathbf{y}(t\le {t}_0+{\tau}_{\mathrm{max}})=\boldsymbol{\psi} (t)\end{array}\right.\operatorname{}\hskip0.1em $

6: Compute loss $ L\left(\mathbf{y}\right)={\int}_{t_0+{\tau}_{\mathrm{max}}}^{t_T}l\left(\mathbf{y}(t)\right)\hskip0.1em \mathrm{d}t $

7: Solve Adjoint dynamics:

8: $ \left\{\begin{array}{l}\dot{\boldsymbol{\lambda}}(t)={\displaystyle \frac{\mathrm{\partial}l(\mathbf{y}(t))}{\mathrm{\partial}\mathbf{y}}}-\boldsymbol{\lambda} (t){\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),{\mathbf{y}}_{({\tau}_1)}(t),\dots, {\mathbf{y}}_{({\tau}_n)}(t))}{\mathrm{\partial}\mathbf{y}}}\\ {}\hskip1.9pc -\sum \limits_{i=1}^n\boldsymbol{\lambda} (t+{\tau}_i){\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t+{\tau}_i,\mathbf{y}(t+{\tau}_i),{\mathbf{y}}_{({\tau}_1)}(t+{\tau}_i),\dots, {\mathbf{y}}_{({\tau}_n)}(t+{\tau}_i))}{\mathrm{\partial}{\mathbf{y}}_{({\tau}_i)}}}\\ {}\boldsymbol{\lambda} (t\ge {t}_T)=\mathbf{0}.\end{array}\right.\operatorname{}\hskip0.1em $

9: Compute $ \frac{\partial L}{\partial \boldsymbol{\theta}} $ :

10:

$$ {\displaystyle \begin{array}{l}{\displaystyle \frac{\mathrm{\partial}L}{\mathrm{\partial}\boldsymbol{\theta }}}=-{\int}_{\hskip-5pt {t}_0}^T\boldsymbol{\lambda} (t)\hskip2pt ({\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),{\mathbf{y}}_{({\tau}_1)}(t),\dots, {\mathbf{y}}_{({\tau}_n)}(t))}{\mathrm{\partial}\boldsymbol{\theta }}}-\sum \limits_{i=1}^n{\displaystyle \frac{\mathrm{\partial}{\mathbf{h}}_{\boldsymbol{\theta}}(t,\mathbf{y}(t),{\mathbf{y}}_{({\tau}_1)}(t),\dots, {\mathbf{y}}_{({\tau}_n)}(t))}{\mathrm{\partial}{\mathbf{y}}_{({\tau}_i)}}}\operatorname{}\\ {}{\displaystyle \hskip1.8pc \frac{\mathrm{\partial}\mathbf{y}(t-{\tau}_i)}{\mathrm{\partial}\boldsymbol{\theta }}})\mathrm{d}t\end{array}} $$

11: Update $ \boldsymbol{\theta} $

12: end for

3.2. On the importance of relevant delays

Learning the delays $ {\tau}_i $ within NDDEs is crucial for accurately modeling partially observed dynamics, and therefore, delays need to be adapted during training. To illustrate the impact of the relevant delays, we evaluate the information that past observations share with the current dynamics of the observable, and how this varies with the delay. To this end, let us consider a dynamical system evolving on a compact smooth manifold $ \mathcal{S}\subset {\unicode{x211D}}^m $ described by a simple 2-delay dynamical system, assumed to be an attractor, and let us consider a $ {C}^2 $ observable function $ y:\mathcal{S}\to \unicode{x211D} $

(3.5) $$ {\displaystyle \begin{array}{l}y\left(t+\Delta t\right)=\cos \left(y\left(t-{\tau}_1\right)\right)\;\sin \left(y\left(t-{\tau}_2\right)\right)-\alpha \operatorname {sinc}\left(3\hskip0.1em y\left(t-{\tau}_1\right)\right)+\alpha\;\cos \left(y\left(t-{\tau}_2\right)\right),\\ {}y\left(t<0\right)=\psi (t),\end{array}} $$

with $ \alpha =0.2 $ , $ {\tau}_1={p}_1^{\star}\hskip0.1em \Delta t $ , $ {\tau}_2={p}_2^{\star}\hskip0.1em \Delta t $ , $ {p}_1^{\star }=250 $ , $ {p}_2^{\star }=400 $ and $ \operatorname{sinc}(x):= \sin (x)/x $ for $ x\ne 0 $ , $ \operatorname{sinc}(0) $ := 1.

The Takens’ theorem (Takens, Reference Takens, Rand and Young1981) rigorously discusses the conditions under which a delayed vector of real-valued observables $ \left(y\left(\mathbf{u}(t)\right),y\left(\mathbf{u}\left(t-\tau \right)\right),\dots, y\left(\mathbf{u}\left(t- n\tau \right)\right)\right) $ , $ n\in \unicode{x2115} $ , defines an embedding, a smooth diffeomorphism onto its image. It guarantees a topological equivalence between the original dynamical system and the one constructed from the memory of the observable. The dynamics of the system can then be reformulated on the set $ \left(y\left(\mathbf{u}(t)\right),y\left(\mathbf{u}\left(t-\tau \right)\right),\dots, y\left(\mathbf{u}\left(t- n\tau \right)\right)\right) $ . However, Takens’ theorem, later extended by Sauer et al. (Reference Sauer, Yorke and Casdagli1991), provides a sufficient condition for reconstruction, but does not specify how to choose the time delay $ \tau $ . From a mathematical viewpoint, the delay could be arbitrary, besides some specific values excluded by the theorem. In practice, its value can strongly condition a successful embedding. If too small, entries of the delay vector data are too similar; if too large, the entries tend to be completely uncorrelated and cannot be numerically linked to a consistent dynamical system (Kantz and Schreiber, Reference Kantz and Schreiber2003).

We here illustrate the impact of suitable delays on the relevance of the information available to inform the future evolution of the observable. The relevance of the delays $ \left\{{\tau}_1,{\tau}_2\right\} $ for informing $ y\left(t+\Delta t\right) $ is assessed in terms of the mutual information between $ \left(y\left(t-{\tau}_1\right),y\left(t-{\tau}_2\right)\right) $ and $ y(t) $ , as suggested in Frase and Swinney (Reference Fraser and Swinney1986). The mutual information $ I\left(X,Y\right) $ between two random variables $ X $ and $ Y $ quantifies the amount of information observing a variable brings about the other one. By considering the problem in Equation (3.5), we estimate the information shared between the current observations $ y(t) $ and their past values $ y\left(t-{\tau}_1\right) $ , $ y\left(t-{\tau}_2\right) $ and is shown in Figure 2 as a 2-D map in terms of $ {p}_1 $ and $ {p}_2 $ . The map is symmetric, as expected since $ I\left(\left(y\left(t-{\tau}_1\right),y\left(t-{\tau}_2\right)\right),y(t)\right)=I\left(\left(y\left(t-{\tau}_2\right),y\left(t-{\tau}_1\right)\right),y(t)\right) $ . It can be seen that the amount of information shared between the current observation and a delay vector of the observable significantly varies with the delays and indeed reaches a maximum for $ {\tau}_1={p}_1^{\star}\hskip0.1em \Delta t $ , $ {\tau}_2={p}_2^{\star}\hskip0.1em \Delta t $ . The ability of the present Neural DDE method to learn the delays, in addition to the model $ {\mathbf{h}}_{\boldsymbol{\theta}} $ , is thus key to its performance and wide applicability.

Figure 2.

$ \left\{{\tau}_1={p}_1\hskip0.1em \Delta t,{\tau}_2={p}_2\hskip0.1em \Delta t\right\} $ -map of delayed mutual information, $ I\left(\left(y\left(t-{\tau}_1\right),y\left(t-{\tau}_2\right)\right),y(t)\right) $ . The maximum is exhibited at $ \left({p}_1,{p}_2\right) $ = $ \left(250,400\right) $ and $ \left(\mathrm{400,250}\right) $ , in agreement with $ {p}_1^{\star }=250 $ , $ {p}_2^{\star }=400 $ .

Figure 2. Long description

The heatmap plots p one on the horizontal axis and p two on the vertical axis, both ranging from zero to four hundred fifty. Color intensity represents the value of delayed mutual information, with the scale bar on the far right ranging from zero point five (dark purple) to above three point zero (yellow). A prominent diagonal of high values runs from the bottom left to the top right, indicating symmetry. Brightest maxima are observed at coordinates two hundred fifty, four hundred and four hundred, two hundred fifty, as well as at regular intervals along both axes, forming a grid-like pattern. The axes are labeled in italic font as p one (horizontal) and p two (vertical). The color bar is labeled with values increasing from bottom to top.

4.1. Dynamical systems

In the following, the dynamical systems considered in this article are presented. All datasets are divided into training, validation, and test sets with proportions of 70%, 10%, and 20%, respectively, and the Dormand–Prince (Dopri5) solver was used for data generation and training (Dormand and Prince, Reference Dormand and Prince1980). Table 1 outlines the number of delays employed in NDDE for each experiment.

Table 1.

Number of delays used in NDDE for each experiment

Table 1. Long description

The table contains one header row and one data row. The header row lists three experiments from left to right: KS, Cavity, and Brusselator. The single data row is labeled NDDE at the far left. The NDDE values for each experiment are as follows: KS is 5, Cavity is 1, and Brusselator is 2. All values are centered within their respective columns.

4.1.1. Population dynamics model

As a first example, we consider the dynamics of a scalar-valued system used to model population dynamics in biology (Arino et al., Reference Arino, Hbid and Dads2009; Banks et al., Reference Banks, Banks, Bommarco, Laubmeier, Myers, Rundlöf and Tillman2017). Such a system is formulated through the following DDE:

(4.1) $$ \frac{\mathrm{d}u(t)}{\mathrm{d}t}=u(t)\left(1-u\left(t-\tau \right)\right),\hskip2em u\left(t\le 0\right)=\psi (t), $$

where we integrate from $ t\in \left[0,10\right] $ , $ \tau =1 $ , $ \psi (t)={u}_0 $ , $ {u}_0 $ is sampled from the uniform distribution $ \mathcal{U}\left(\mathrm{2.0,3.0}\right) $ and $ 256 $ trajectories were generated.

4.1.2. Brusselator

A second experiment showcases how NDDEs can effectively model partially observed systems with past state values. We consider the 2-species Belousov–Zhabotinsky kinetic equation (Belousov, Reference Belousov1959; Zhabotinskii, Reference Zhabotinskii1964) that can be modeled by the so-called Brusselator system:

(4.2) $$ \left\{\hskip-1pc ,\begin{array}{ll}& {\displaystyle \frac{\mathrm{d}{u}_1(t)}{\mathrm{d}t}}=A-{Bu}_1-{u}_1+{u}_1^2{u}_2,\\ {}& {\displaystyle \frac{\mathrm{d}{u}_2(t)}{\mathrm{d}t}}={Bu}_1-{u}_1^2{u}_2,\end{array}\right)\operatorname{} $$

with $ {u}_1(t) $ and $ {u}_2(t) $ the two species concentrations at a given time. We integrate in the time domain $ t\in \left[0,25\right] $ . The initial condition $ {u}_1 $ is sampled from the uniform distribution $ \mathcal{U}\left(\mathrm{0,2.0}\right) $ while $ {u}_2=0.0 $ and $ 1024 $ trajectories were generated. We set ourselves in the partially observable case, where we only have access to the dynamics of $ {u}_1 $ , that is, $ y(t)\equiv {u}_1(t) $ , and wish to reconstruct the whole dynamics.

4.1.3. Kuramoto–Sivashinsky (KS) system

This model was originally derived to describe the complex dynamics of flames in a combustion. The Kuramoto–Sivashinsky system is described as:

(4.3) $$ \frac{\partial u}{\partial t}+\frac{\partial^2u}{\partial {x}^2}+\frac{\partial^4u}{\partial {x}^4}+\frac{1}{2}\frac{\partial {u}^2}{\partial x}=0. $$

The system is integrated over the time domain $ t\in \left[0,30\right] $ and its spatial domain $ {D}_x=\left[0,22\right] $ is discretized into 128 points. A total of $ 2048 $ trajectories were generated from the initial condition samples. To put ourselves in the partially observed setting, we choose to observe the solution in time at $ p $ locations uniformly spread across the spatial domain (we retain $ p=5 $ ).

4.1.4. Incompressible open cavity flow

As an experimental demonstrator, we consider the modeling based on time-series derived from wind tunnel experiments of an open cavity flow sketched in Figure 3. The facility and the details of the experimental setup are described in Tuerke et al. (Reference Tuerke, Lusseyran, Sciamarella, Pastur and Artana2020). Open cavity flows have attracted numerous research efforts in the last few decades due to the interesting dynamics at work. The flow is characterized by an impinging shear layer activating a centrifugal instability in a cavity. This interplay, reminiscent of the feedback acoustic mechanisms described in Rossiter (Reference Rossiter1964)), leads to a self-sustained oscillation. A broad range of dynamics is observed, ranging from limit cycles to toroidal and chaotic dynamics. The system is observed through a sensor in $ P $ measuring the local pressure fluctuations, hence $ y(t)\equiv P\left(t;{\mathbf{x}}_{\mathrm{P}}\right) $ . The data are obtained for a Reynolds number $ \mathit{\operatorname{Re}}=9190 $ based on the length L of the cavity.

Figure 3.

Sketch of open cavity flow taken from (Tuerke et al., Reference Tuerke, Lusseyran, Sciamarella, Pastur and Artana2020). A data-acquiring sensor is placed in P. The cavity has a length L and depth H. The incoming laminar boundary layer flow is characterized by the freestream velocity $ {U}_{\infty } $ and the momentum thickness $ {\Theta}_0 $ .

Figure 3. Long description

At the left, horizontal arrows indicate incoming laminar flow with freestream velocity U sub infinity and initial momentum thickness Theta sub 0. The flow enters a rectangular cavity, with the top edge labeled Theta sub effective. Above the cavity, a dashed line marks the shear layer, and a filled circle labeled P denotes the sensor position. Inside the cavity, two large counter-rotating vortices are shown by spiral arrows, labeled k plus and k minus, forming a carousel pattern. The cavity’s horizontal length is labeled L and vertical depth H. All labels are positioned relative to the cavity geometry, with flow direction indicated by arrows.

4.2. Results

We now assess the performance of the models with their ability to predict future measurements of a partially observed system. In this study, LSTM, NODE, ANODE, Latent ODE, and NDDE were selected for comparison. Table 2 displays the test MSE loss over each experiment. Appendix D goes into more detail about each model’s architecture and the training and testing procedure. Every model incorporates a form of “memory” into its architecture, with the exception of NODE. While LSTM and Latent ODE utilize hidden units and ANODE employs its augmented state $ \mathbf{a}(t) $ , NDDE leverages past observations such as $ \mathbf{y}\left(t-\tau \right) $ . In all subsequent figures, the y-axis $ \mathbf{y}(t) $ represents our observables (introduced for each system in Section 4.1), defined as $ \mathbf{y}(t)=y\left(\mathbf{u}(t)\right) $ .

Table 2.

Model performance (MSE) over the test set in each experiment averaged over five runs

Table 2. Long description

The table has three rows for experiments: Brusselator, K S, and Cavity. Each row lists mean squared error values plus or minus standard deviation for five models. For Brusselator: L S T M is 0.0051 plus or minus 0.0031, N O D E is 0.77 plus or minus 0.00080, A N O D E is 0.0050 plus or minus 0.0050, Latent O D E is 0.014 plus or minus 0.0076, N D D E is 0.016 plus or minus 0.0076. For K S: L S T M is 0.77 plus or minus 0.041, N O D E is 0.71 plus or minus 0.10, A N O D E is 0.55 plus or minus 0.027, Latent O D E is 0.43 plus or minus 0.07, N D D E is 0.28 plus or minus 0.024. For Cavity: L S T M is 0.75 plus or minus 0.46, N O D E is 0.96 plus or minus 0.0011, A N O D E is 0.65 plus or minus 0.0090, Latent O D E is 0.25 plus or minus 0.14, N D D E is 0.13 plus or minus 0.012. Lower values indicate better performance, with N D D E generally showing the lowest errors for K S and Cavity, and A N O D E and L S T M performing best for Brusselator.

4.2.1. Population dynamics model

Figures 4 and 5, respectively, depict the model’s robust convergence to accurate dynamics and the evolution of the delay $ \tau $ during training over many seeds, showcasing a consequence of Takens’ theorem (Takens, Reference Takens, Rand and Young1981). Using a delay-coordinate map, one can indeed construct a diffeomorphic shadow manifold $ {M}^{\prime } $ from univariate observations of the original system in the generic sense. We here observe the whole state, so that $ y(t)\equiv u(t) $ , and the delay coordinate map is then in terms of $ \left(y(t),y\left(t-\tau \right)\right) $ , with $ \tau >0 $ . The result of Figure 5 shows that the learned delay can lie between about 0.8 and 1.4. In classical approaches (see Tan et al. (Reference Tan, Algar, Corrêa, Small, Stemler and Walker2023), the selected delay with Takens’ theorem for State Space Reconstruction is typically chosen close to the minimum delayed mutual information of the time series.

Figure 4.

The prediction of the DDE model is seen to accurately match the true evolution.

Figure 5.

Toy dataset delay evolution during training for several initial conditions.

Figure 5. Long description

The x axis is labeled Steps, ranging from 0 to 10,000. The y axis is labeled Delay value, spanning 0.5 to 2.5. Multiple colored lines begin at the left edge, each representing a distinct initial condition. Most lines start with rapid changes, then stabilize horizontally between delay values of 0.5 and 2.0 as steps increase. A few lines show higher initial spikes above 2.0 before settling. The majority of trajectories converge to steady values, with minor fluctuations visible throughout. No legend or specific line labels are present.

4.2.2. Brusselator

In the case of this stiff and periodic dynamics, all models demonstrate satisfactory performance except NODE, which essentially predicts a mean trajectory thus highlighting the importance of incorporating memory terms for an accurate prediction. Remarkably, both LSTM and ANODE perform equally well, with NDDE and Latent ODE slightly trailing by a narrow margin as shown in Figure 6. In addition to evaluating the MSE loss performance, Figure 7 demonstrates the stability of each trained model on the Brusselator system over an extended period. After training within a specific time interval, we lengthened the integration period to five times the original training duration to assess the models’ performance. It is observed that NDDE, along with LSTM, NODE, and Latent ODE, are the only models that remain stable throughout this duration, with NDDE exhibiting the best performance over the extended horizon.

Figure 6.

Examples of prediction of the Brusselator dynamics in terms of $ y(t)\equiv {u}_1(t) $ , for initial conditions sampled at random.

Figure 6. Long description

The left panel shows y of t on the vertical axis from 0 to 4 and t on the horizontal axis from 0 to 20. Six lines are plotted: Truth (solid black), NODE (blue dashed), ANODE (orange dashed), NDDE (red dashed), NDDE_FIXED (green dashed), and LSTM forward slash G R U (purple dashed). The Truth, ANODE, NDDE, NDDE_FIXED, and LSTM forward slash G R U lines overlap, showing sharp periodic peaks at t near 10 and 17. The NODE line deviates, remaining nearly flat at y of t equals 1. The right panel repeats the same axes and legend. All models except NODE closely follow the Truth, with three sharp peaks at similar t values. The NODE line remains flat. The legend below the panels identifies each line style and color.

Figure 7.

Phase-portrait of $ \left(y(t),\dot{y}(t)\right) $ . Long-term behavior of each trained model for the Brusselator system.

Figure 7. Long description

The x axis is labeled y of t, ranging from negative zero point two to zero point five. The y axis is labeled y dot of t, ranging from zero to four. Six trajectories are plotted: TRUTH in solid black, NODE in dashed blue, ANODE in dashed orange, NDDE in solid green, NDDE FIXED in dashed red, and LSTM forward slash G R U in dashed purple. All trajectories form closed loops with similar shapes but show visible deviations, especially near the lower left and upper right regions. The TRUTH trajectory is the reference, with other models approximating it to varying degrees. The legend is located in the lower right quadrant.

4.2.3. KS system

This numerical experiment deals with a partially observed system in a chaotic regime, observed through the solution at $ p $ locations uniformly spread across the spatial domain. Figure 8 showcases random test samples from two different initial conditions, highlighting how NDDEs outperform other models struggling with the dynamics of the selected features. In a chaotic setting, the statistics of the dynamics are often more informative than the trajectory itself. We hence focus on the probability distribution of the prediction across a large time horizon for the different $ p $ components of the observables $ {\left\{{y}_i\right\}}_{i=1}^p $ , cf. Figure 9. It is seen that NDDE again outperforms the other approaches, with predictions that are statistically closer to the ground truth. This is further supported by the evolution of the maximum Lyapunov exponent (MLE) of the resulting models. Table 3 displays the MLE estimates for each model, showing that the Neural DDE with learnable delays closely aligns with the ground truth compared to other models.

Figure 8.

Prediction of the KS system from a test sample for different models.

Figure 8. Long description

Starting from the top-left panel and proceeding rightward, each graph plots y of t versus t for a test sample. The x-axis is labeled t, ranging from 0 to about 35. The y-axis is labeled y of t, with varying ranges per panel. Each panel overlays six colored lines: solid blue for Truth, dashed orange for D D E, dashed green for O D E, dashed red for A N O D E, dashed purple for L S T M, and dashed brown for Latent O D E. The legend in the bottom right specifies these mappings. Across panels, the Truth line shows oscillatory behavior, while model predictions diverge in amplitude and phase, with Latent O D E often deviating most. The panels illustrate differences in model accuracy and dynamics for the KS system.

Figure 9.

Probability density functions of the predictions $ {y}_i(t) $ over the time horizon and across several initial conditions for different models. KS system.

Figure 9. Long description

Starting at the top left and moving right, then to the bottom row, each panel shows density on the y axis and X on the x axis. Six lines are plotted per panel: solid blue for Truth, dashed orange for D D E, dashed green for L S T M, dashed red for O D E, dashed purple for A N O D E, and dashed brown for Latent O D E. The legend is in the bottom right panel. Each panel displays different distributions, with the Truth line generally smoother and the model lines showing varying degrees of deviation and multimodality. Peaks and spread differ by panel, reflecting changes in initial conditions. Some models, such as A N O D E and Latent O D E, show sharper or shifted peaks compared to Truth. The density range and x axis limits vary slightly between panels, highlighting the diversity of prediction behaviors across models and conditions.

Table 3.

Estimation of the maximum Lyapunov exponent $ {\boldsymbol{\lambda}}_{\mathrm{max}} $ for the KS system based on the generated trajectories from the test set for each model

Table 3. Long description

Column headers from left to right are Ground truth, N D D E, N O D E, A N O D E, and Latent O D E. The row label is lambda sub max. Ground truth is 0.129, N D D E is 0.128, N O D E is 0.097, A N O D E is 0.120, Latent O D E is 0.035. N D D E is bolded, indicating it is closest to ground truth. Latent O D E is the lowest value. The table compares the estimation of the maximum Lyapunov exponent for the KS system based on generated trajectories from the test set for each model.

4.2.4. Cavity

In this experimental fluid flow configuration, the NDDE formulation again outperforms other models, as illustrated in Figure 10 where the predictions are plotted for two test samples. The benefit of a time-delay model here comes from the fact that this experimental setup involves a large vortex within the cavity, effectively acting like a feedback loop and hence suitably described by governing equations involving delayed contributions. Latent ODE yields acceptable results compared to NODE, which only generates the system’s average trajectory, while LSTM and ANODE capture oscillations, albeit occasionally in conflicting phases. These experiments demonstrate that NDDE can effectively model trajectories even in the presence of noise in the data, such as in this experimental dataset.

Figure 10.

Prediction of the cavity observables from different models for two test samples.

Figure 10. Long description

The left panel shows y of t versus t from 0 to about 33. The black solid line labeled Truth shows the reference data with sharp oscillations. The blue dashed line labeled N O D E remains near zero. The orange dashed line labeled A N O D E shows a large amplitude oscillation out of phase with Truth. The green dashed line labeled N D D E and the red dashed line labeled N D D E underscore F I X E D both follow the Truth curve more closely, capturing peaks and troughs but with some deviation. The purple dashed line labeled L S T M forward slash G R U shows a lower amplitude oscillation. The right panel repeats this structure for a second test sample, with similar trends: Truth (black) is closely followed by N D D E and N D D E underscore F I X E D, while N O D E, A N O D E, and L S T M forward slash G R U diverge more from the reference. The legend below the panels identifies each model by color and line style.

Since the noise in the measurements results from phenomena not related to the physical system at hand, such as, for instance, the electronic noise of the sensors, it is statistically independent from the true observables’ time history. Learning to relate these past data with current measurements then quickly averages out these spurious additional dimensions introduced by the noise as the amount of observations increases. More formally, the noise $ \boldsymbol{\xi} $ being assumed additive, low amplitude, and zero-mean, the expectation of the map $ {\mathbf{h}}_{\boldsymbol{\theta}} $ writes

$$ {\displaystyle \begin{array}{l}\unicode{x1D53C}\left[{\mathbf{h}}_{\boldsymbol{\theta}}\Big(t,\mathbf{y}(t)+\boldsymbol{\xi} (t),\mathbf{y}\left(t-{\tau}_1\right)+\boldsymbol{\xi} \left(t-{\tau}_1\right),\dots, \mathbf{y}\left(t-{\tau}_n\right)+\boldsymbol{\xi} \left(t-{\tau}_n\right)\Big)\right]\\ {}\hskip4em \approx \unicode{x1D53C}\left[{\mathbf{h}}_{\boldsymbol{\theta}}\Big(t,\mathbf{y}(t),\mathbf{y}\left(t-{\tau}_1\right),\dots, \mathbf{y}\left(t-{\tau}_n\right)\Big)\right]+\unicode{x1D53C}\left[\nabla {\mathbf{h}}_{\boldsymbol{\theta}}\Big(t,\mathbf{y}(t),\mathbf{y}\left(t-{\tau}_1\right),\dots, \mathbf{y}\left(t-{\tau}_n\right)\Big)\cdot \boldsymbol{\xi} \right],\\ {}\hskip4em =\unicode{x1D53C}\left[{\mathbf{h}}_{\boldsymbol{\theta}}\Big(t,\mathbf{y}(t),\mathbf{y}\left(t-{\tau}_1\right),\dots, \mathbf{y}\left(t-{\tau}_n\right)\Big)\right]+\unicode{x1D53C}\left[\nabla {\mathbf{h}}_{\boldsymbol{\theta}}\Big(t,\mathbf{y}(t),\mathbf{y}\left(t-{\tau}_1\right),\dots, \mathbf{y}\left(t-{\tau}_n\right)\Big)\right]\cdot \unicode{x1D53C}\left[\boldsymbol{\xi} \right],\\ {}\hskip4em =\unicode{x1D53C}\left[{\mathbf{h}}_{\boldsymbol{\theta}}\Big(t,\mathbf{y}(t),\mathbf{y}\left(t-{\tau}_1\right),\dots, \mathbf{y}\left(t-{\tau}_n\right)\Big)\right],\end{array}} $$

so that the noise is effectively averaged out, provided the empirical mean approximates the expectation well enough.

To illustrate the benefit of learning the delays, instead of a priori setting them, the performance of NDDEs is here illustrated in both situations. A set of delays is considered, ranging from 0.2 to 1, and the mean squared error (MSE) associated with each resulting model is monitored, with and without learning the delays (see Figure 11). The results show that learned delays consistently outperform fixed delays in leading to models associated with a significantly lower MSE.

Figure 11.

Evolution of the MSE train loss of the NDDE model with constant (solid lines) and learnable delays (dashed lines) for different delay initialization values ranging from $ 0.2 $ to $ 1.0 $ .

Figure 11. Long description

The 3D line graph has three axes: x-axis labeled Steps ranging from 0 to 5000, y-axis labeled Loss from 0.5 to 1.3, and z-axis labeled Delay initialization tau from 0.2 to 1.0. Each colored pair of lines corresponds to a specific tau value, with solid lines representing fixed tau and dashed lines representing learnt tau, as indicated in the legend at the top left. For each tau value, the loss generally increases with steps, with a sharp rise near the maximum step value. The separation between solid and dashed lines shows the effect of fixed versus learnt delays, with learnt tau (dashed) often showing lower loss than fixed tau (solid) for the same initialization. The colored lines are stacked along the tau axis, with each pair showing a similar trend but at different loss levels. The legend box in the upper left corner specifies line styles: solid for Fixed tau and dashed for Learnt tau.

4.3. Closure modeling with the ROMs

We now revisit the Kuramoto–Sivashinsky configuration from a closure modeling viewpoint. Conceptually, one assumes an accurate model of the system to be known, but the approximate time-dependent solution to be evaluated with a reduced-order model (ROM) Galerkin approach. Specifically, the solution $ \mathbf{u}(t) $ is approximated as $ \mathbf{u}(t)\approx \sum \limits_{i=1}^{N_{\mathrm{terms}}}{\varphi}_i\left(\mathbf{x}\right)\hskip0.1em {a}_i(t) $ , with $ {\left\{{\varphi}_i\left(\mathbf{x}\right)\equiv {\boldsymbol{\varphi}}_i\right\}}_{i=1}^{N_{\mathrm{terms}}}=: \Phi $ spanning a spatial basis, so that describing the solution boils down to deriving an evolution equation for the coefficients $ \mathbf{a}(t)=\left({a}_1(t),\dots, {a}_{N_{\mathrm{terms}}}(t)\right) $ . Because the approximation basis $ \Phi $ does not span the whole state space whenever $ {N}_{\mathrm{terms}}<\operatorname{card}\left(\mathbf{u}\right) $ , the $ \Phi $ -projected model describing the dynamics of the coefficients $ {\left\{{a}_i(t)\right\}}_i $ , resulting from a Galerkin formulation, is not $ \Phi $ -invariant. To compensate for the loss of information leaking from the subspace spanned by $ \Phi $ , one can complement the reduced-order model with an additional, data-based, closure term in the governing equations. Our experiment employs a Galerkin-projected ROM approach, with modes $ {\left\{{\boldsymbol{\varphi}}_i\right\}}_i $ derived from a proper orthogonal decomposition (POD), with $ {N}_{\mathrm{terms}} $ = 4, 8, and 10 POD modes (low, medium, and high data regimes), respectively accounting for 58%, 87%, and 94% of the system’s energy (i.e., variance). In practice, the more modes selected in our POD Galerkin ROM, the smaller the correction required by the closure term needs to be.

We now consider different models for representing the closure term and compare their performance in terms of the test MSE loss in Table 4. It is seen that the DDE closure term with learnable delays consistently outperforms the ODE closure term across all data regimes (low, medium, and high). This performance difference becomes especially pronounced in the low data regime of the four-mode POD Galerkin ROM, indicating that the Neural ODE closures are not well-suited for such scenarios. Figure 12, along with their respective representations in the original state space in Figure 13, showcase the prediction from two samples of initial conditions from the test dataset for the four-mode POD Galerkin ROM. As the number of modes increases, the discrepancy between ODE and DDE closures decreases quantitatively, as seen in Table 4. This trend is expected since the ROM then captures most of the system’s information. By studying these specific ROMs, we can identify where ODE closure terms fall short and how incorporating past states with the DDE closure term can address ODE’s deficiencies in low data scenarios. Our approach is compared with the closure model CD-ROM (Menier et al., Reference Menier, Bucci, Yagoubi, Mathelin and Schoenauer2023), which considers a continuous embedding of past information in the form of an exponentially decaying dynamics auxiliary term $ \mathbf{z}(t) $ mimicking a memory:

(4.4) $$ \mathbf{z}(t)={\int}_{t_0}^t{e}^{\left(t-\tau \right)\hskip0.1em \Lambda}\hskip0.1em \mathbf{y}\left(\tau \right)\hskip0.1em \mathrm{d}\tau $$

with $ \Lambda $ a diagonal matrix whose elements are associated with the relevant time scales of the system under consideration; see (Menier et al., Reference Menier, Bucci, Yagoubi, Mathelin and Schoenauer2023) for details. Interestingly, compared to a memoryless Markovian model such as NODE, it only improves the quality of the prediction in the low-data regime, and considering a higher-dimensional memory variable $ \mathbf{z}(t) $ does not improve the learning process in high-data regimes.

Table 4.

Summary of model performance metrics, that is, the test MSE loss. CD-ROM, ODE, and DDE closure models are compared across the 4, 8, and 10 POD mode settings

Table 4. Long description

The table has four columns and four rows. The top row lists column headers from left to right as Model, 4 modes, 8 modes, and 10 modes. The first column lists models from top to bottom as C D dash R O M, O D E dash R O M, and D D E dash R O M. For C D dash R O M, the test M S E loss values are 7.94 plus or minus 0.53 for 4 modes, 0.43 plus or minus 0.24 for 8 modes, and 0.099 plus or minus 0.019 for 10 modes. For O D E dash R O M, the values are 13.61 plus or minus 0.34 for 4 modes, 0.44 plus or minus 0.044 for 8 modes, and 0.084 plus or minus 0.0028 for 10 modes. For D D E dash R O M, the values are 3.39 plus or minus 0.034 for 4 modes, 0.18 plus or minus 0.07 for 8 modes, and 0.067 plus or minus 0.0078 for 10 modes. The lowest test M S E loss for each mode count is bolded and occurs for D D E dash R O M at 8 and 10 modes, and for C D dash R O M at 8 modes.

Figure 12.

Examples of predictions of POD Galerkin ROM (4 modes).

Figure 12. Long description

Top-left panel a sub 0 plots t on the x-axis and shows four lines: blue dashed for Truth, orange for D D E R O M, green for C D R O M, red for O D E R O M. All lines rise then fall, with Truth and D D E R O M peaking near t equals 20. Top-right panel a sub 1 shows all lines decreasing then oscillating, with Truth and D D E R O M closely aligned. Bottom-left panel a sub 2 shows all lines rising, with Truth and D D E R O M peaking near t equals 20, C D R O M and O D E R O M diverging. Bottom-right panel a sub 3 shows oscillatory behavior, with Truth and D D E R O M closely following each other, C D R O M and O D E R O M diverging after t equals 20. Legend at bottom-right: Truth (blue dashed), D D E R O M (orange), C D R O M (green), O D E R O M (red). Vertical red line at t equals 30 in all panels.

Figure 13.

Examples of predictions of POD Galerkin ROM (4 modes) with different models from KS testset reconstructed in terms of the full solution field $ u\left(\mathbf{x},t\right) $ (first three columns). The two rightmost columns show the absolute value of the reconstruction error associated with each model.

Figure 13. Long description

The grid contains three rows and five columns. Each row starts with a panel labeled Ground Truth showing the reference solution field with a color scale from minus two to two. The next panel in each row displays a model prediction: D D E dash R O M in the first row, C D dash R O M in the second, O D E dash R O M in the third, each with a similar color scale. The third panel in each row shows R O M predictions with a color scale from minus two hundred to two hundred. The fourth panel in each row is labeled Diff R O M and Truth, showing the absolute error between R O M and ground truth, with a color scale from zero to three. The fifth panel in each row shows the absolute error for the specific model versus ground truth: Diff D D E dash R O M and Truth in the first row, Diff C D dash R O M and Truth in the second, Diff O D E dash R O M and Truth in the third, all using a color scale from zero to three. The spatial order is left to right within each row, top row first, followed by the second and third rows. All panels use vertical axes labeled zero to twenty or zero to thirty, and horizontal axes labeled zero to twenty. Color bars are present below each panel, matching the respective value ranges.

Several numbers of learnable delays were considered for the NDDE approach, ranging from 1 to 3. For this system, the results were essentially unaffected by the number of delays, showing that a single delay is enough to recover most of the necessary information for a good prediction.