# Data-driven discovery of nonlinear dynamics in complex systems through phase-space informed neural ordinary differential equations

David Garrido González, Nathaniel Saura, Sadruddin Benkadda, Peter Beyer

► To cite this version:

David Garrido González, Nathaniel Saura, Sadruddin Benkadda, Peter Beyer. Data-driven discovery of nonlinear dynamics in complex systems through phase-space informed neural ordinary differential equations. Physics of Plasmas, 2025, 32 (6), 10.1063/5.0270212. hal-05196898

HAL Id: hal-05196898

https://hal.science/hal-05196898v1

Submitted on 5 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Data-driven discovery of nonlinear dynamics in complex systems through phase-space informed neural ordinary differential equations

Special Collection: Papers from the 5th International Conference on Data-Driven Plasma Science

David Garrido González ⚙️; Nathaniel Saura ⚙️; Sadruddin Benkadda 🧑; Peter Beyer ⚙️

<figure><img src="image_2.png" /></figure>

Check for updates

Phys. Plasmas 32, 063904 (2025)

https://doi.org/10.1063/5.0270212

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

## Articles You May Be Interested In

Design of information interactive platform for maintenance application of communication equipment in distribution network based on integrated interface optimization

AIP Advances (October 2025)

Investigating characteristics of the growth rates from QuaLiKiz using an interpretable surrogate model

Phys. Plasmas (May 2025)

Machine learning to analyze images of shocked materials for precise and accurate measurements

J. Appl. Phys. (September 2017)

<figure><img src="image_5.png" /></figure>

**Physics of Plasmas**

Learn more

**Read our Author Testimonials**

*Physics of Plasmas* has a
9.1 author satisfaction rating

<figure><img src="image_6.png" /></figure>

AIP
Publishing
# Data-driven discovery of nonlinear dynamics in complex systems through phase-space informed neural ordinary differential equations

Cite as: Phys. Plasmas 32, 063904 (2025); doi: 10.1063/5.0270212

Submitted: 10 March 2025 · Accepted: 24 May 2025 ·

Published Online: 9 June 2025

<figure><img src="image_7.png" /></figure>

View Online

<figure><img src="image_8.png" /></figure>

Export Citation

<figure><img src="image_9.png" /></figure>

CrossMark

David Garrido González, a) a) Nathaniel Saura, b) Sadruddin Benkadda, b) and Peter Beyer c) c)

## AFFILIATIONS

Aix-Marseille Université, CNRS, PIIM, UMR 7345, 13397 Marseille, France

Note: This paper is part of the Special Topic on Papers from the 5th International Conference on Data-Driven Plasma Science.

a) Electronic mail: david.garrido-gonzalez@univ-amu.fr

b) Author to whom correspondence should be addressed: sadruddin.benkadda@univ-amu.fr

## ABSTRACT

Neural ordinary differential equations (NODEs) offer a data-driven approach to modeling chaotic dynamical systems without prior knowledge of governing equations. We show that enhanced phase-space sampling through diverse initial conditions significantly improves NODEs' ability to capture both local and global system behaviors. Using the Lorenz system, we demonstrate that broader phase-space coverage enhances generalization, reduces overfitting, and improves long-term predictive stability. Evaluations based on Lyapunov exponents, the Kaplan-Yorke dimension, return maps, and probability density functions confirm that increased initial condition diversity leads to more accurate reconstructions of chaotic attractors and system dynamics.

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0270212

## I. INTRODUCTION

Numerical modeling of partial differential equations (PDEs) is an active research area focused on describing physical systems across diverse initial and boundary conditions to analyze their properties and dynamics. However, the computational resources required to fully solve these systems are often substantial, posing a significant challenge particularly for complex, turbulent systems such as those found in fusion plasmas¹,² and neutral fluids.³,⁴ Reduced-order models (ROMs) have been widely investigated to mitigate this challenge, with techniques such as the Galerkin projection frequently employed to reduce PDE dynamics to an ordinary differential equation (ODE) system.⁵,⁶ This approach involves identifying an appropriate basis, often derived from proper orthogonal decomposition (POD), to facilitate the projection.⁷⁻⁹ Despite its effectiveness, this traditional method requires detailed prior knowledge of the physical system and involves complex mathematical operations, such as computing integrals over intricate geometries, which can still be computationally intensive.

Recent advancements in artificial intelligence (AI) have enabled the development of surrogate models that often outperform traditional ROMs in predicting system dynamics, albeit at the cost of interpretability.¹⁰⁻¹⁴ In this context, neural ordinary differential equations

(NODEs) have emerged as a powerful approach for modeling complex dynamical systems directly from available data, even when it is irregularly sampled, without requiring explicit knowledge of the governing ODEs.¹⁵ In classical ordinary differential equations (ODEs), system evolution is defined by its current state and time derivative, forming a vector field in phase space that reveals behaviors like stability points, periodic orbits, and chaotic regions. In NODEs, the right-hand side of the ODE system is parameterized by a neural network (NN) as

$$ \frac{dx}{dt} = f(x, t, \theta), \qquad (1) $$

where $x(t)$ is the system's state and $f$ is the NN-approximated vector field, with $\theta$ as learned parameters. This framework enables the NN to map the system's state into phase space by learning its derivatives, creating a continuous vector field to govern system evolution over time.

In this study, we present a method to improve the predictive capabilities of NODEs in capturing complex dynamical systems, using the Lorenz system as a benchmark. Originally formulated to describe convective fluid dynamics, the Lorenz equations¹⁶ have become a fundamental model for studying chaotic behavior. The governing equations of the Lorenz system are
$$ 
\begin{aligned}
\dot{x} &= \sigma(y - x), \\
\dot{y} &= rx - y - xz, \\
\dot{z} &= xy - bz,
\end{aligned}
\qquad (2) $$

where $\sigma, r$, and $b$ are system parameters. The dynamics of (2) exhibits chaotic behavior, marked by a sensitive dependence on initial conditions and convergence to a strange attractor with a fractal structure in phase space.¹⁷ The trajectories show irregular and unpredictable motion over time, illustrating the “butterfly effect,” where even small changes in initial conditions can lead to markedly different outcomes. Despite this unpredictability, the system remains bounded within a finite region of phase space, with trajectories drawn toward two distinct lobes.

The Lorenz system has become emblematic in the study of deterministic chaos due to its capacity to capture nonlinear dynamics within a compact, low-dimensional framework. Its significance transcends meteorology, extending into complex systems like plasma physics, where turbulence, transport phenomena, and wave-wave interactions exhibit similarly chaotic behavior. In plasma contexts, reduced-order models derived from high-dimensional nonlinear systems have been shown to exhibit Lorenz-like dynamics, particularly in the modeling of wave coupling,¹⁸,¹⁹ shear flows,²⁰ and intermittent transport events such as edge-localized mode (ELM) dynamics in tokamaks.²¹ They have also been applied to analyze the interaction between vortex motion and shear flow in nonlinear interchange turbulence at the edge of tokamak plasmas.²² In fusion plasmas, the interplay between ion temperature gradient turbulence (ITG) and shear flow regulation has been successfully modeled using Lorenz-type systems, offering insights into transitions between turbulent and regulated transport regimes.²³ These plasma phenomena exhibit chaotic dynamics akin to those modeled by the Lorenz system, making it a valuable testbed for developing and validating modeling techniques. This analogy underscores the structural universality of Lorenz-type attractors in plasma turbulence and motivates their use as surrogate systems in modern machine learning frameworks. The present work applies NODEs to the Lorenz system, leveraging its representative chaotic dynamics to explore the capability of NODEs in capturing and predicting complex behavior from data.

While previous studies have employed NODEs to simulate chaotic dynamics, they have primarily focused on fitting and extending a limited number of known trajectories.²⁴,²⁵ Although these models can successfully reproduce specific chaotic trajectories,²⁶ their ability to generalize across arbitrary initial conditions with the same intrinsic parameters remains largely unexplored. To bridge this gap, we investigate the performance of NODEs trained on multiple initial conditions, demonstrating that a more diverse training set enhances the model's ability to generate Lorenz-like trajectories, even for unseen initial conditions.

To rigorously evaluate NODEs in capturing the chaotic dynamics of the Lorenz system, we compute key dynamical invariants, including Lyapunov exponents and the Kaplan-Yorke dimension,¹⁷ which characterize sensitivity to initial conditions and attractor structure. Additionally, we analyze return maps and lobe residence times to assess how well the models replicate chaotic behavior. Finally, we compute probability density functions (PDFs) of the predicted trajectories to provide a comprehensive statistical validation of the NODE-based approach.

## II. NODES AS AN OPTIMAL CONTROL PROBLEM

NODEs are a class of models that represent the evolution of a hidden state as a continuous process, governed by a differential equation whose vector field is parameterized by a NN. Instead of applying a fixed sequence of transformations, NODEs learn a global function $f(\mathbf{x}(t), t, \theta)$ that continuously governs how the state evolves over time

$$ \frac{d\mathbf{x}}{dt} = f(\mathbf{x}(t), t, \theta), \quad \mathbf{x}(t_0) = \mathbf{x}_0, \quad t \in [t_0, t_f], \qquad (3) $$

where $\mathbf{x}(t) \in \mathbb{R}^N$ is the system state, $\theta \in \mathbb{R}^{N_p}$ are the learnable network parameters, and $f: \mathbb{R}^N \times \mathbb{R}^{N_p} \times \mathbb{R} \to \mathbb{R}^N$ is a NN representing the system's time-dependent vector field.²⁷ These continuous formulations can be seen as the limit of discrete updates of a hidden state, $\mathbf{h}_{t+1} = \mathbf{h}_t + f(\mathbf{h}_t, \theta_t)$, which correspond to an Euler discretization of the underlying dynamics. The state at any time $t$ is obtained by solving this initial value problem using a numerical ODE solver, which adaptively evaluates the dynamics based on the system's complexity.

This continuous-time framework offers several advantages. It allows for smooth modeling of time-evolving processes, naturally supports irregularly sampled data and variable-length sequences, and dynamically adjusts computational cost by letting the solver choose how many evaluations are needed to meet a desired accuracy.

### A. Forward pass

From a control perspective, training NODEs can be seen as solving an optimization problem where the goal is to minimize a loss functional subject to the dynamic constraint (1). Formally, this reads

$$ \min_{\theta} L(\mathbf{x}(t_f)) = L\left(\mathbf{x}(t_0) + \int_{t_0}^{t_f} f(\mathbf{x}(t), t, \theta) dt\right), \qquad (4) $$

which can also be expressed as

$$ L(\text{ODESolve}(\mathbf{x}(t_0), f, t_0, t_f, \theta)). \qquad (5) $$

Here, $L$ measures the discrepancy between the predicted final state $\mathbf{x}(t_f)$ and the observed target $\hat{\mathbf{x}}_f$, while ODESolve denotes any differentiable numerical integration scheme. The training process aims to find the parameters $\theta$ that minimize this loss, effectively recovering the dynamics $f$ that best describe the observed evolution.

This formulation makes explicit the continuous-time nature of NODEs. Instead of applying a discrete sequence of transformations, NODEs learn a continuous transformation by integrating over the state-space trajectory induced by $f$. This allows for flexible temporal modeling, particularly in systems with irregularly spaced data or chaotic sensitivity.

### B. Backward pass

In a regular NN, the process of training involves updating the parameters of the network based on the gradients of the loss function with respect to those parameters. The gradients indicate the direction and magnitude of the changes needed to minimize the loss. In the context of NODEs, we also want to update the parameters of the NN that approximates the ODE dynamics in order to learn the underlying dynamics of the system. The parameters of the NN capture the behavior of the ODE and affect how the system evolves over time.
To compute the gradients of the loss with respect to the NN parameters, we need to perform backpropagation. Backpropagation computes gradients by propagating the error backward through the NN's computational graph. This process allows us to determine how each parameter contributes to the overall loss. In the case of NODEs, the computational graph includes both the NN and the ODE solver. The ODE solver is responsible for integrating the dynamics obtained from the NN to compute the next state. Therefore, in order to compute the gradients of the loss with respect to the NN parameters, we need to account for the contribution of the ODE solver. However, directly differentiating the ODE solver with respect to its parameters is both challenging and memory-intensive.

### C. The adjoint sensitivity method

This primary technical challenge in training was first addressed by Chen et al.¹⁵ using the adjoint sensitivity method²⁸ to compute the gradients. An overview of this process is shown in Fig. 1. This method facilitates the computation of gradients of the loss function L, Eq. (4), concerning both the initial state, $\mathbf{x}(t_0)$, and the parameters, $\theta$, by introducing an adjoint state, $\mathbf{a}(t)$, defined as the derivative of L with respect to $\mathbf{x}(t)$,

$$ \mathbf{a}(t) = \frac{\partial L}{\partial \mathbf{x}}, \qquad (6) $$

and governed by the following differential equation:¹⁵

$$ \frac{\mathrm{d}\mathbf{a}(t)}{\mathrm{d}t} = -\mathbf{a}(t)^T \frac{\partial f(\mathbf{x}(t), t, \theta)}{\partial \mathbf{x}}. \qquad (7) $$

This adjoint state facilitates the computation of gradients through the formulas:¹⁵

$$ \frac{\partial L}{\partial \theta} = - \int_{t_f}^{t_0} \mathbf{a}(t)^T \frac{\partial f(\mathbf{x}(t), t, \theta)}{\partial \theta} dt. \qquad (8) $$

These gradients are efficiently computed by solving the ODE for $\mathbf{a}(t)$ using the same numerical method as in the forward pass. During the forward pass, the ODE solver calculates the solution of the differential

equation $\mathbf{x}(t)$ using the initial state $\mathbf{x}(t_0)$ and the function $f(\mathbf{x}(t), t, \theta)$. During backpropagation, the adjoint state $\mathbf{a}(t)$ is obtained by integrating the differential equation backward in time from $t_f$. Subsequently, this adjoint state is employed to determine gradients for the loss function in relation to the initial state and parameters of the ODE function, which are then utilized to update the model parameters through gradient descent optimization. This approach optimizes gradient computation, enhancing the efficiency and effectiveness of NODE-based models.

### III. METHODOLOGY

For training the models, we used `torchdiffEq`,¹⁵ a powerful library for solving ODEs efficiently and accurately, along with PyTorch Lightning²⁹ for building and training the NNs. The training process begins with constructing a comprehensive dataset consisting of observations over a specified time interval and a separate test set for validation. The primary objective is to minimize the discrepancy between NODE-predicted and expected trajectories using a loss function that balances mean squared and mean absolute errors.

Four models were trained following this procedure, each using one (NODE1), three (NODE3), five (NODE5), or nine (NODE9) initial conditions, respectively. Each NODE model features a shallow, fully connected NN with two dense layers, both employing ReLU³⁰ (rectified linear unit) activations to introduce nonlinearity. Model parameters are optimized using the Adam optimizer³¹ in conjunction with a cyclic learning rate policy,³² which dynamically adjusts the learning rate throughout training. This approach is particularly effective for training NODE models, which are prone to becoming trapped in local minima.³³ A cyclically varying learning rate enhances convergence by allowing the model to escape local minima, thereby facilitating faster and more effective learning.³²

During the learning rate decay phase, model weights are saved whenever validation loss improves. If no improvement occurs after a cooldown,³⁴ the learning rate is increased by a predefined step count, and training resumes from the best-performing model checkpoint. If the loss function stagnates despite learning rate adjustments, we hybridize the current model parameters with those from previous

<figure><img src="image_10.png" /><figcaption>FIG. 1. A schematic illustrating the iterative process of the NODE for updating network parameters. As the loss function in NODEs involves derivatives, it necessitates the use of the adjoint method or the storage of intermediate steps to enable parameter updates.</figcaption></figure>
checkpoints to facilitate escaping larger local minima and continuing optimization.

To further refine NODE training, we mitigate biases arising from underrepresented regions of the system's dynamics. Specifically, because training sub-trajectories of length *p* are randomly sampled from complete system trajectories, certain behaviors may appear as outliers, hindering the model's ability to generalize. To counteract this, we implement a batch selection procedure inspired by genetic algorithms.³⁵ This procedure randomly selects an *n*-sample subset from the training dataset, organizing it into *m* batches, each containing *p* consecutive time steps from the same trajectory. The model prioritizes challenging regions by sampling 70% from these areas.

This dynamic sampling strategy adapts to the model's conver-
gence state, balancing well-understood and challenging regions to miti-
gate catastrophic forgetting.³⁶,³⁷ These techniques collectively reduce
the NODE's tendency to overfit specific trajectories, allowing it to cap-
ture more complex relationships over longer timescales. The combina-
tion of learning rate decay, hybridization, and smart batching, as
illustrated in Fig. 2, facilitates a more comprehensive exploration of the
dataset.

The NN used to parameterize the dynamics function $f(x, t, 0)$ within the neural ODE framework is a fully connected feedforward network with two hidden layers. The input and output sizes are both 3, corresponding to the state variables $(x, y, z)$ and their time derivatives $(\dot{x}, \dot{y}, \dot{z})$, respectively. Each hidden layer contains 75 neurons, with a ReLU activation function applied between layers. This architecture was selected for its simplicity and sufficient expressive power to capture the nonlinear dynamics of the Lorenz system.

We chose a maximum of nine initial conditions to strike a bal-
ance between model performance and computational cost. Prior work
has shown that insufficient coverage of the phase space can limit a
model's ability to generalize beyond its training distribution, particu-
larly in chaotic systems where dynamics vary significantly across
regions.³⁸ Although adding more trajectories can enhance the statisti-
cal fidelity and robustness of the learned dynamics, performance gains
tend to saturate beyond a certain point, making additional simulations
less cost-effective.³⁹

IV. RESULTS

To produce the trajectories used in our experiments, we numeri-
cally integrated the Lorenz system from a range of initial conditions
using the *odeint* solver from the *torchdiffeq* library, employing

<figure><img src="image_11.png" /><figcaption>FIG. 2. Downward staircase learning rate curve on logarithmic scale.</figcaption></figure>

the Runge–Kutta 4 (RK4) method. The Lorenz system was simulated
with its classical parameters: $\sigma = 10$, $r = 28$, and $b = 8/3$. Integration
was carried out over a fixed time interval $[0, t_{\max}] = [0, 5]$, using a tra-
jectory length of $N_t = 800$ uniformly spaced time steps. For each
epoch, $n = 256$ sub-trajectories of length $p = 150$ were sampled from
the full trajectories (with 1, 3, 5, or 9 initial conditions, depending on
the NODE model) and organized into $m = 32$ batches, following the
methodology described in Sec. III. For each initial condition, the full
state trajectory was computed at each time step, forming a dataset rep-
resentative of the system's evolution in state space.

The specific initial conditions used during training are visual-
ized as colored dots in the lower inset of Fig. 3, with each color
corresponding to a different NODE model (NODE{1, 3, 5, 9}).
These initial values, of the form $x_0 = y_0 = z_0$, were selected as
equally spaced points within the interval $-15 < x_0 = y_0 = z_0 < -5$.
The number of initial conditions corresponds to the number of train-
ing trajectories used by each model (1, 3, 5, or 9). This range was
chosen to ensure that the resulting trajectories explore distinct
regions of the Lorenz attractor, specifically, both lobes for more neg-
ative values (i.e., $x_0 \le -7.5$) and primarily one lobe closer to the
upper bound (i.e., $-7.5 \le x_0 < -5$). The number of training initial
conditions per model matches the NODE variant (e.g., NODE3 uses
three equally spaced values within the interval).

Figure 3 shows the trajectory prediction error, measured as mean squared error (MSE), for test initial conditions spanning the broader range $-20 < x_0 = y_0 = z_0 < 0$. The chart reports MSE values for both interpolated and extrapolated initial conditions, those between and beyond the training points, respectively, enabling assessment of the models’ generalization capabilities.

We restrict the domain to this interval because it allows us to
evaluate the model's ability to generalize across trained and untrained
regions of the attractor. Beyond this range, the MSE tends to be rela-
tively flat, reflecting the chaotic nature of the system: small differences
in initial conditions rapidly lead to diverging trajectories, making long-
term prediction accuracy less meaningful.

As the number of training trajectories increases, the error consis-
tently decreases, especially in regions between and slightly beyond the
training points. This indicates improved generalization and more

<figure><img src="image_12.png" /><figcaption>FIG. 3. Mean squared error (MSE) of predicted trajectories for different initial conditions of the form ($x_0 = y_0 = z_0$), evaluated against the ground truth. Each curve corresponds to a NODE model trained with a different number of initial conditions: NODE{1, 3, 5, 9}. The specific initial conditions used during training for each model are marked as colored dots in the lower inset. As the number of training trajectories increases, the model's ability to generalize across the attractor improves, reducing the average. In the chart, "ICs" refers to initial conditions.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 4. Comparison of predicted (with the NODE trained with nine different initial conditions) and expected trajectories for the Lorenz system. Left: time evolution of x, y, and z. Right: 3D trajectories, with real (red) and predicted (blue) paths.</figcaption></figure>

faithful statistical learning of the underlying dynamics when using a diverse set of initial conditions.

Figure 4 illustrates a comparison between the expected trajectory of the Lorenz system and the prediction generated by the NODE9 model, which was trained using nine different initial conditions. The left panel displays the time evolution of the system’s state variables (x, y, z), demonstrating how the model captures the intricate oscillatory behavior of the Lorenz system. The right panel presents a reconstructed three-dimensional trajectory in phase space, with real trajectories in red and NODE9 predictions in blue. The alignment between the predicted and true attractor structures indicates that the NODE9 model successfully learns and generalizes the underlying chaotic dynamics, accurately reconstructing the Lorenz attractor’s geometry.

To characterize the chaotic nature of the system, several diagnostic measures were employed. Table I presents the estimated Lyapunov exponents for different initial conditions in the chaotic Lorenz system. The bold values represent the expected solutions, while the remaining values are predictions from the NODE models. The numbers in parentheses denote the uncertainty affecting the last two digits of each estimate. The Lyapunov exponents provide insights into the system’s sensitivity to initial conditions, quantifying how nearby trajectories diverge or converge over time. A positive Lyapunov exponent indicates exponential divergence of nearby trajectories, typical of chaotic behavior. In contrast, zero or negative Lyapunov exponents suggest nearby

trajectories either converge or remain stable (for negative or zero values, respectively).

The algorithm for computing Lyapunov exponents is based on QR decomposition (referring to the factorization $A = QR$, where Q is orthogonal and R is upper triangular), a standard method used to estimate a system’s sensitivity to initial conditions.43 This procedure, as outlined in Algorithm 1 in the Appendix, was applied to a sample of 8000 trajectories to ensure robustness.

The performance of the models varies significantly depending on the number of initial conditions used for training. The biggest deviations from the reference values are seen in NODE1, which was trained with only one initial condition. These deviations are particularly noticeable in the first and third Lyapunov exponents. The third exponent, $\lambda_3 = -17.19(11)$, is significantly more negative than the expected value, and the largest exponent, $\lambda_1 = 0.541(94)$, is significantly smaller than the expected value. These differences imply that a single initial condition is not sufficient to fully describe the sensitivity of the system to perturbations, which results in a poor description of the chaotic dynamics.

The models exhibit better alignment with the reference values as the number of initial conditions increases. With five initial conditions, NODE3 gives a closer value to the reference. With five initial conditions, NODE5 exhibits superior performance, particularly for the largest exponent $\lambda_1 = 0.73(13)$. However, the second Lyapunov

<table><caption>TABLE I. Estimated Lyapunov exponents and Kaplan–Yorke dimension for various models using the QR decomposition method. The table includes expected Lyapunov exponents, computed using trajectories with identical initial conditions, along with reference values from the literature for comparison. $\lambda_1$, $\lambda_2$, and $\lambda_3$ represent the first three Lyapunov exponents, while $\sum \lambda_i$ indicates the sum of all exponents. The values in brackets indicate the uncertainty in the last two digits. The Kaplan–Yorke dimension is calculated based on these values.</caption><thead><tr><th>Model</th><th>Kaplan–Yorke</th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th><th>λ<sub>3</sub></th><th>Σ λ<sub>i</sub></th></tr></thead><tbody><tr><td>Reference</td><td>2.06<sup>40,41</sup></td><td>0.9056<sup>42</sup></td><td>0<sup>42</sup></td><td>−14.5721<sup>42</sup></td><td>−13.6665<sup>42</sup></td></tr><tr><td>Expected (QR)</td><td>2.0611(42)</td><td>0.892(60)</td><td>0.000(15)</td><td>−14.592(54)</td><td>−13.700(82)</td></tr><tr><td>NODE1</td><td>2.0337(57)</td><td>0.541(94)</td><td>0.039(29)</td><td>−17.19(11)</td><td>−16.61(15)</td></tr><tr><td>NODE3</td><td>2.076(14)</td><td>0.85(14)</td><td>0.020(47)</td><td>−11.32(45)</td><td>−10.46(48)</td></tr><tr><td>NODE5</td><td>2.0400(98)</td><td>0.73(13)</td><td>−0.073(96)</td><td>−16.40(22)</td><td>−15.74(27)</td></tr><tr><td>NODE9</td><td>2.0555(37)</td><td>0.863(49)</td><td>0.031(16)</td><td>−16.12(54)</td><td>−15.22(54)</td></tr></tbody></table>
exponent, $\lambda_2 = -0.073(96)$, is still negative instead of near zero, suggesting that even though the model does a better job of representing the dynamics of the system, it is still unable to fully capture the intricacy of the phase-space behavior.

With nine initial conditions, NODE9 yields the most accurate results. The Lyapunov exponents for this model, $\lambda_1 = 0.863(49)$ and $\lambda_3 = -16.12(54)$, are closer to the reference values. Moreover, the sum of the Lyapunov exponents, $\sum \lambda_i = -15.22(54)$, is the one that comes closest to the expected value.

While Lyapunov exponents tell us how quickly trajectories diverge, they do not provide a complete picture of the system's chaotic behavior, including the overall structure and dimensionality of the attractor. These exponents provide insights into how many degrees of freedom actively shape the system's evolution, which is crucial for understanding long-term behavior. The Kaplan-Yorke dimension⁴⁰ helps quantify the fractal nature of the chaotic attractor, beyond its divergence rate. This dimension obeys the expression

$$D_{KY} = j + \frac{\sum_{i=1}^{j} \lambda_i}{|\lambda_{j+1}|} = 2 + \frac{\lambda_1 + \lambda_2}{|\lambda_3|}, \quad (9)$$

where j is the index at which the sum of Lyapunov exponents up to $\lambda_j$ becomes nonpositive (i.e., when $\sum_{i=1}^{j} \lambda_i \ge 0$ and $\sum_{i=1}^{j+1} \lambda_i < 0$). In this case, $j = 2$.

As the number of initial conditions increases, the models converge toward the expected Kaplan-Yorke dimension of ~2.06. NODE1, trained with a single initial condition, yields a dimension of 2.0337(57). With three initial conditions, NODE3 shows an improved dimension of 2.076(14), though the larger uncertainty and deviation in the sum of the Lyapunov exponents indicate an incomplete representation of the system's chaotic dynamics. NODE5, using five initial conditions, improves to 2.0400(98), but still deviates from the expected value. NODE9, trained with nine initial conditions, achieves the closest result, 2.0555(37), with lower uncertainty. As more initial conditions are accounted for, the model's accuracy improves, with both the dimension and Lyapunov exponents converging toward their true values.

The previous analysis of Lyapunov exponents and the Kaplan-Yorke dimension provided valuable insights into how well the models capture the divergence rates and fractal structure of the chaotic attractor. However, to further evaluate the models' performance, it is crucial to analyze how accurately they reproduce the detailed structure of the attractor. Return maps offer deeper insights into the models' capacity to capture the system's complex dynamics by comparing the predicted and actual maxima of the system's trajectories. This allows for a more granular validation of the attractor's structure and behavior, beyond what is revealed by the Lyapunov exponent and Kaplan-Yorke analyses alone.⁴⁴

Figure 5 shows return maps of successive maxima of $z(t)$ based on Lorenz's method.⁴⁴ These maps were generated by first obtaining $z(t)$ over an extended time interval $0 < t < 200 \times t_{\text{Train}}$ for both predicted and actual time series. The local maxima of $z(t)$ for both the predicted and real data are ordered chronologically as $[z_1, z_2, \dots, z_m]$, and the maximum pairs $[z_i, z_{i+1}]$ are plotted as points. The black points in the panels represent the actual Lorenz system data, while the colored points correspond to predictions from the NODE models. The return maps were constructed for initial

conditions both extrapolated beyond or interpolated within the training set range ($x_0 = y_0 = z_0 \in [-15, -5]$) in three distinct cases: top row ($x_0 = y_0 = z_0 = -17.5$), center row ($x_0 = y_0 = z_0 = -10.0$), and bottom row ($x_0 = y_0 = z_0 = -2.5$). NODE3 and NODE5 models yield accurate results for trajectories with the initial conditions $x_0 = y_0 = z_0 = -17.5$ and $x_0 = y_0 = z_0 = -10.0$. However, they struggle with trajectories where $x_0 = y_0 = z_0 = -2.5$, which exhibit less chaotic behavior. Conversely, NODE9 demonstrates a broader applicability, successfully predicting the correct tendencies even beyond the training set range. Notably, this model performs well in the aforementioned challenging region, offering a more generalized solution with a wider scope of applicability.

The Lorenz attractor, like other chaotic systems, is highly sensitive to initial conditions, particularly during transitions between the attractor's lobes. Studying the time spent in the different lobes offers insights into the system's stability and the frequency of transitions between different states. Figure 6 exhibits the comparison between the expected and predicted residence time spent in the left lobe. The green points represent the ratio of the predicted and the expected residence times in the left lobe for trajectories starting from $-20 < x_0 = y_0 = z_0 < -15$ (out of the training set). The orange and blue points represent the same ratios for trajectories starting from $-15 < x_0 = y_0 = z_0 < -5$ (inside of the training set) and $-5 < x_0 = y_0 = z_0 < 0$ (out of the training set), respectively. NODE9 is the only model that shows a clear clustering around the ideal dashed line, suggesting that it captures the system's long-term stability and transition dynamics more accurately. Despite potential inaccuracies, this prediction is valuable for understanding long-term trends and general system behavior.

Long-term behavior is analyzed through probability density functions (PDFs) of trajectories, providing insights into the system's stability and convergence. The PDFs in Fig. 7 show the range of values for each coordinate in a trajectory starting at $x = y = z = -10$ and extending up to $t = 1000 \times t_{\text{Train}} = 5000$. Trajectories predicted by NODE5 and NODE9 exhibit PDFs comparable to expected solutions. However, NODE5 shows an out-of-scale PDF for the trajectory starting at $x = y = z = -2.5$ (top row) as it converges to a limit cycle, a behavior also observed in NODE3. These convergences were not evident in previous figures because they occur at longer times. Figure 7 highlights that NODE models trained on few trajectories develop a biased understanding of the target system, whereas NODE9 better captures the Lorenz system's global and long-term dynamics.

The enhancement of multiscale learning in chaotic dynamical systems through the increased sampling of initial conditions can be understood as an improvement in phase-space exploration, leading to more robust statistical learning. Chaotic systems exhibit extreme sensitivity to initial conditions, where small perturbations lead to significantly different long-term trajectories.¹⁶ By incorporating a larger set of initial conditions, we effectively sample a broader region of the phase space, capturing a more diverse range of system behaviors.⁴⁵ This improved coverage reduces biases in training datasets and enables machine learning models to learn underlying dynamical structures more accurately.⁴⁶ Furthermore, increasing the number of initial conditions improves generalization by exposing the model to a richer variety of state transitions, thus mitigating overfitting to localized dynamics.⁴⁷
<figure><img src="image_14.png" /><figcaption>FIG. 5. Return maps generated for trajectories with initial conditions extrapolated to the left ($x_0 = y_0 = z_0 = -17.5$), interpolated ($x_0 = y_0 = z_0 = -10.0$), and extrapolated to the right ($x_0 = y_0 = z_0 = -2.5$) of the range of the training set. The results are compared to the real return map (in black).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 6. Comparison of expected vs predicted percentages of time spent in the system's left lobe. The dashed black line indicates perfect alignment. Green points represent extrapolated trajectories (initial conditions left of training set), orange points indicate interpolated trajectories (within training set range), and blue points represent extrapolated trajectories (initial conditions right of training set).</figcaption></figure>

V. CONCLUSION

The combination of figure interpretations, metrics, and diagnostics indicates that the NODE9 model can statistically reproduce the Lorenz system's behavior at different levels of detail. Initial investigations into dynamical invariants suggest that the NODE9 model

exhibits a similar sensitivity to initial conditions, divergence rate, and attractor structure as the targeted chaotic system.

The close approximations of the Lyapunov exponents with the expected signs confirm that the NODE9 model accurately captures the system's sensitivity to initial conditions. Furthermore, analysis of
<figure><img src="image_16.png" /><figcaption>FIG. 7. PDF comparison between the models and the true solution for trajectories with initial conditions extrapolated to the left ($x_0 = y_0 = z_0 = -17.5$), interpolated ($x_0 = y_0 = z_0 = -10.0$), and extrapolated to the right ($x_0 = y_0 = z_0 = -2.5$) of the range of the training set. Here, the distributions of values taken by each coordinate in three different trajectories, extending up to $t = 5000$ (a significantly long run), are shown.</figcaption></figure>

return maps and residence times suggests that the NODE9 model exhibits a certain degree of awareness of the statistical properties of the Lorenz system, although interpretability remains a challenge given the black-box nature of NNs. Finally, comparisons of PDFs validated the long-term behavior of the AI-based model for both seen and unseen trajectories. Compared to similar studies in the literature, we focused on modeling the global system dynamics rather than individual trajectories. In our comparisons, we considered the results of the NODE1 model, which is the closest case to those discussed in previous studies. We showed that training a NODE on only one initial condition leads to a model that not only struggles to extrapolate its learning to other trajectories but also fails to reproduce the full spectrum of chaotic behavior.

The statistical representativity gained through this approach also enhances long-term prediction stability, which is a critical challenge in modeling chaotic systems.⁴⁸ Empirical studies in reservoir computing and deep learning have demonstrated that broader phase-space sampling leads to more reliable reconstructions of attractors and better short- and long-term forecasts.⁴⁹,⁵⁰ This aligns with the principle that learning in high-dimensional nonlinear systems benefits from an extensive training distribution to approximate invariant manifolds effectively.⁵¹ Consequently, the strategic selection of diverse initial

conditions is not merely a methodological consideration but an essential step toward enhancing the predictive power of data-driven models in complex dynamical systems.

This study demonstrates that broadening the diversity of initial conditions during NODE training enhances the model's ability to closely approximate the intrinsic dynamics of the Lorenz system. This approach not only improves the temporal extrapolation of known trajectories but also enables the generation of Lorenz-like trajectories for unseen initial conditions, capturing potentially different behaviors. NODE9's consistent performance across key metrics (Lyapunov exponents, Kaplan-Yorke dimension, return maps, and PDFs) demonstrates that a wider range of initial conditions allows the model to capture both the chaotic attractor's local structures and global dynamics, providing accurate predictions across both familiar and novel trajectories.

Increasing the diversity of initial conditions significantly improves the model's ability to capture complex dynamical patterns. However, while this diversification is a powerful technique, it must be carefully balanced to avoid saturation, where additional trajectories no longer provide meaningful improvements, increase computational cost, or even harm the learning process. This is particularly relevant in chaotic systems, where the objective is not to reproduce individual
trajectories, which diverge rapidly, but to capture the system's invariant statistics, such as the attractor structure or Lyapunov exponents, that reflect the underlying physical laws. From this perspective, our results suggest that training with nine well-chosen initial conditions provides sufficient coverage to achieve statistical fidelity without incurring the penalties of oversaturation.

To address training challenges, a custom procedure was proposed, incorporating cyclic learning rate adjustments and a batch resampling strategy (described in Sec. III). These methods help the model converge faster, particularly in cases where traditional optimization approaches struggle with local minima. While these techniques mitigate a common limitation in NODE models, further research is needed to refine them and enhance the efficiency of NODE training.

The ability of NODEs to learn and generalize such complex dynamics from phase-space data thus reinforces the relevance of Lorenz systems as foundational testbeds for developing reduced-order, interpretable surrogates in plasma modeling. This bridging of machine learning with chaotic dynamical systems opens new avenues for real-time predictive modeling, particularly in high-dimensional, data-scarce plasma regimes.

## ACKNOWLEDGMENTS

We acknowledge the financial support from the French Federation for Magnetic Fusion Studies (FR-FCM) and the EUROfusion consortium, funded by the Euratom Research and Training Programme under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**David Garrido González:** Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal). **Nathaniel Saura:** Conceptualization (equal); Formal analysis (equal); Methodology (equal). **Sadruddin Benkadda:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **Peter Beyer:** Conceptualization (equal); Methodology (equal); Supervision (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## APPENDIX A: ALGORITHM FOR THE COMPUTATION OF LYAPUNOV EXPONENTS

This paper employs a QR decomposition-based method⁴³ to compute the Lyapunov exponents, an approach used for assessing the sensitivity of dynamical systems to initial conditions. This

<caption>**Algorithm 1.** Computation of Lyapunov exponents using QR decomposition.</caption>

1: **Input:** Trajectories X of shape (N, T, 3), time step dt
2: **Output:** Lyapunov exponents $\lambda_1, \lambda_2, \lambda_3$
3: Initialize $\lambda = 0_{N \times 3}$ (to store final exponents)
4: **for** i = 1 to N **do**
5: &nbsp;&nbsp;&nbsp;&nbsp; ▷ Loop over each trajectory
6: &nbsp;&nbsp;&nbsp;&nbsp; Initialize deformation matrix M ← $I_3$ (identity matrix)
7: &nbsp;&nbsp;&nbsp;&nbsp; Initialize sum of exponents S ← $0_3$
8: &nbsp;&nbsp;&nbsp;&nbsp; **for** k = 1 to T **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Extract current state $x_0 ← X[i, k, :]$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Compute Jacobian matrix J ← Jacobian($x_0$)
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Update deformation matrix $M_n ← (I + dt \cdot J)M$
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Perform QR decomposition: Q, R ← qr($M_n$)
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Accumulate logarithm of diagonal of $R: S ← S + \log[\text{diag}(R)]$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Update deformation matrix M ← Q
15: &nbsp;&nbsp;&nbsp;&nbsp; **end for**
16: **end for**
&nbsp;&nbsp;&nbsp;&nbsp; **return** $\lambda$

technique iteratively tracks the evolution of perturbation vectors to quantify the system's divergence rates along different axes in phase space. By leveraging QR decomposition, the method ensures numerical stability, preventing the exponential growth of perturbation vectors from affecting the accuracy of the computed exponents.

The procedure begins with the initialization of a deformation matrix M, set as the identity matrix, and a sum of exponents vector S, initialized to zero. At each time step, the system's Jacobian matrix J is computed at the current state, representing the local linearization of the dynamics. The deformation matrix is then updated using a first-order approximation

$$M_n ← (I + dt \cdot J)M.$$

Subsequently, QR decomposition is applied to $M_n$, yielding an orthogonal matrix Q and an upper triangular matrix R. The logarithm of the diagonal elements of R is accumulated into S, and the updated Q matrix is reassigned to M to maintain the orthonormality of perturbation vectors throughout the iterations.

Once all time steps have been processed, the Lyapunov exponents $\lambda_i$ are computed as

$$\lambda_i = \frac{S_i}{T \cdot dt},$$

where T is the total integration time and dt is the time step size.

## REFERENCES

¹X. Litaudon, F. Jenko, D. Borba, D. Borodin, B. Braams, S. Brezinsek, I. Calvo, R. Coelho, A. Donné, O. Embréus et al., “Eurofusion-theory and advanced simulation coordination (E-TASC): Programme and the role of high performance computing,” *Plasma Phys. Controlled Fusion* 64, 034005 (2022).

²J. Van Dijk, G. Kroesen, and A. Bogaerts, “Plasma modelling and numerical simulation,” *J. Phys. D* 42, 190301 (2009).

³L. D. Landau and E. M. Lifshitz, *Fluid Mechanics* (Elsevier, 1987), Vol. 6.
<footnote n="4">⁴S. Hickel, C. P. Egerer, and J. Larsson, "Subgrid-scale modeling for implicit large eddy simulation of compressible flows and shock-turbulence interaction," Phys. Fluids 26, 106101 (2014).</footnote>

<footnote n="5">⁵J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, *Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems* (SIAM, 2016).</footnote>

<footnote n="6">⁶C. W. Rowley and S. T. Dawson, "Model reduction for flow analysis and control," *Annu. Rev. Fluid Mech.* 49, 387–417 (2017).</footnote>

<footnote n="7">⁷P. Beyer, S. Benkadda, and X. Garbet, "Proper orthogonal decomposition and Galerkin projection for a three-dimensional plasma dynamical system," *Phys. Rev. E* 61, 813 (2000).</footnote>

<footnote n="8">⁸P. Benner, S. Gugercin, and K. Willcox, "A survey of projection-based model reduction methods for parametric dynamical systems," *SIAM Rev.* 57, 483–531 (2015).</footnote>

<footnote n="9">⁹A. T. Mohan, K. Nagarajan, and D. Livescu, "Learning stable Galerkin models of turbulence with differentiable programming," *arXiv:2107.07559* (2021).</footnote>

<footnote n="10">¹⁰R. Anirudh, R. Archibald, M. S. Asif, M. M. Becker, S. Benkadda, P.-T. Bremer, R. H. S. Budé, C. S. Chang, L. Chen, R. M. Churchill, J. Citrin, J. A. Gaffney, A. Gainaru, W. Gekelman, T. Gibbs, S. Hamaguchi, C. Hill, K. Humbird, S. Jalas, S. Kawaguchi, G.-H. Kim, M. Kirchen, S. Klasky, J. L. Kline, K. Krushelnick, B. Kustowski, G. Lapenta, W. Li, T. Ma, N. J. Mason, A. Mesbah, C. Michoski, T. Munson, I. Murakami, H. N. Najm, K. E. J. Olofsson, S. Park, J. L. Peterson, M. Probst, D. Pugmire, B. Sammuli, K. Sawlani, A. Scheinker, D. P. Schissel, R. J. Shalloo, J. Shinagawa, J. Seong, B. K. Spears, J. Tennyson, J. Thiagarajan, C. M. Ticoş, J. Trieschmann, J. van Dijk, B. Van Essen, P. Ventzek, H. Wang, J. T. L. Wang, Z. Wang, K. Wende, X. Xu, H. Yamada, T. Yokoyama, and X. Zhang, "2022 review of data-driven plasma science," *IEEE Trans. Plasma Sci.* 51, 1750–1838 (2023).</footnote>

<footnote n="11">¹¹J. Han, A. Jentzen, and E. Weinan, "Solving high-dimensional partial differential equations using deep learning," *Proc. Natl. Acad. Sci. U. S. A.* 115, 8505–8510 (2018).</footnote>

<footnote n="12">¹²Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, "Fourier neural operator for parametric partial differential equations," *arXiv:2010.08895* (2020).</footnote>

<footnote n="13">¹³K. Kashinath, M. Mustafa, A. Albert, J. Wu, C. Jiang, S. Esmaeilzadeh, K. Azizzadenesheli, R. Wang, A. Chattopadhyay, A. Singh et al., "Physics-informed machine learning: Case studies for weather and climate modelling," *Philos. Trans. R. Soc. A* 379, 20200093 (2021).</footnote>

<footnote n="14">¹⁴H. Gao, L. Sun, and J.-X. Wang, "PhyGeoNet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state PDEs on irregular domain," *J. Comput. Phys.* 428, 110079 (2021).</footnote>

<footnote n="15">¹⁵R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, "Neural ordinary differential equations," *Adv. Neural Inf. Process. Syst.* 31, 6571–6583 (2018).</footnote>

<footnote n="16">¹⁶E. N. Lorenz, "Deterministic nonperiodic flow," *J. Atmos. Sci.* 20, 130–141 (1963).</footnote>

<footnote n="17">¹⁷S. Lynch, *Dynamical Systems with Applications Using Python* (Springer, 2018).</footnote>

<footnote n="18">¹⁸C. Letellier, L. Aguirre, J. Maquet, and B. Lefebvre, "Analogy between a 1D model for nonlinear wave–wave interaction in a plasma and the 3D Lorenz dynamics," *Physica D* 179, 33–52 (2003).</footnote>

<footnote n="19">¹⁹L. Feng, J. Li, Q.-Y. Feng, Y. Xingyu, X. Liu, and J. Wang, "Lorenz-type chaos (turbulence) in five-wave mixing in a uniform streaming plasma," *Quantum Electron.* 27, 1012 (1997).</footnote>

<footnote n="20">²⁰M. Yagi, S.-I. Itoh, K. Itoh, and A. Fukuyama, "Theoretical study of inhomogeneous plasma turbulence," *J. Plasma Fus. Res.* 2, 113–116 (2000).</footnote>

<footnote n="21">²¹J. Connor, "A review of models for ELMs," *Plasma Phys. Controlled Fusion* 40, 191 (1998).</footnote>

<footnote n="22">²²A. Takayama, T. Unemura, and M. Wakatani, "Comparison of Lorenz type and plasma fluid models for nonlinear-interchange-mode generation of shear flow," *Plasma Phys. Controlled Fusion* 40, 775 (1998).</footnote>

<footnote n="23">²³K. Takeda, S. Benkadda, S. Hamaguchi, and M. Wakatani, "Onset of intermittent thermal transport by ion-temperature-gradient-driven turbulence based on a low-degree-of-freedom model," *Phys. Plasmas* 11, 3561–3571 (2004).</footnote>

<footnote n="24">²⁴Y. D. Zhong, B. Dey, and A. Chakraborty, "Symplectic ode-net: Learning Hamiltonian dynamics with control," *arXiv:1909.12077* (2019).</footnote>

<footnote n="25">²⁵A. Allauzen, T. P. M. Dardis, and H. Plath, "Experimental study of neural ode training with adaptive solver for dynamical systems modeling," *arXiv:2211.06972* (2022).</footnote>

<footnote n="26">²⁶A. J. Linot and M. D. Graham, "Data-driven reduced-order modeling of spatiotemporal chaos with neural ordinary differential equations," *Chaos* 32, 073110 (2022).</footnote>

<footnote n="27">²⁷P. Kidger, "On neural differential equations," *arXiv:2202.02435* (2022).</footnote>

<footnote n="28">²⁸L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, in *Mathematical Theory of Optimal Processes*, edited by L. W. Neustadt (John Wiley & Sons, Nashville, TN, 1962).</footnote>

<footnote n="29">²⁹W. Falcon and T. P. L. Team, *Pytorch Lightning* (Lightning AI, 2019).</footnote>

<footnote n="30">³⁰A. F. Agarap, "Deep learning using rectified linear units (relu)," *arXiv preprint arXiv:1803.08375* (2018).</footnote>

<footnote n="31">³¹D. P. Kingma, "Adam: A method for stochastic optimization," *arXiv:1412.6980* (2014).</footnote>

<footnote n="32">³²L. N. Smith, "Cyclical learning rates for training neural networks," in *IEEE Winter Conference on Applications of Computer Vision (WACV)* (IEEE, 2017), pp. 464–472.</footnote>

<footnote n="33">³³T. Thummerer, J. Stoljar, and L. Mikelsons, "NeuralFMU: Presenting a workflow for integrating hybrid NeuralODES into real-world applications," *Electronics* 11, 3202 (2022).</footnote>

<footnote n="34">³⁴A. cool-down is a waiting period that allows the model to stabilize without overreacting to minor fluctuations in performance.</footnote>

<footnote n="35">³⁵S. Katoch, S. S. Chauhan, and V. Kumar, "A review on genetic algorithm: Past, present, and future," *Multimedia Tools Appl.* 80, 8091–8126 (2021).</footnote>

<footnote n="36">³⁶R. M. French, "Catastrophic forgetting in connectionist networks," *Trends Cognit. Sci.* 3, 128–135 (1999).</footnote>

<footnote n="37">³⁷T. Tadros, G. P. Krishnan, R. Ramyaa, and M. Bazhenov, "Sleep-like unsupervised replay reduces catastrophic forgetting in artificial neural networks," *Nat. Commun.* 13, 7742 (2022).</footnote>

<footnote n="38">³⁸S. Scher and G. Messori, "Generalization properties of feed-forward neural networks trained on Lorenz systems," *Nonlinear Processes Geophys.* 26, 381–399 (2019).</footnote>

<footnote n="39">³⁹D. Chakraborty, S. W. Chung, T. Arcomano, and R. Maulik, "Divide and conquer: Learning chaotic dynamical systems with multistep penalty neural ordinary differential equations," *Comput. Methods Appl. Mech. Eng.* 432, 117442 (2024).</footnote>

<footnote n="40">⁴⁰J. L. Kaplan and J. A. Yorke, in *Functional Differential Equations and Approximation of Fixed Points: Proceedings. Bonn, July 1978* (Springer, 2006), pp. 204–227.</footnote>

<footnote n="41">⁴¹N. V. Kuznetsov, T. N. Mokaev, O. A. Kuznetsova, and E. V. Kudryashova, "The Lorenz system: Hidden boundary of practical stability and the Lyapunov dimension," *Nonlinear Dyn.* 102, 713–732 (2020).</footnote>

<footnote n="42">⁴²B. J. Geurts, D. D. Holm, and E. Luesink, "Lyapunov exponents of two stochastic Lorenz 63 systems," *J. Stat. Phys.* 179, 1343–1365 (2020).</footnote>

<footnote n="43">⁴³Z.-M. Chen, K. Djidjeli, and W. Price, "Computing Lyapunov exponents based on the solution expression of the variational system," *Appl. Math. Comput.* 174, 982–996 (2006).</footnote>

<footnote n="44">⁴⁴J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, "Using machine learning to replicate chaotic attractors and calculate Lyapunov exponents from data," *Chaos* 27, 121102 (2017).</footnote>

<footnote n="45">⁴⁵E. Ott, *Chaos in Dynamical Systems* (Cambridge University Press, 2002).</footnote>

<footnote n="46">⁴⁶J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, "Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach," *Phys. Rev. Lett.* 120, 024102 (2018).</footnote>

<footnote n="47">⁴⁷S. L. Brunton and J. N. Kutz, *Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control* (Cambridge University Press, 2019).</footnote>

<footnote n="48">⁴⁸P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos, "Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks," *Proc. R. Soc. A* 474, 20170844 (2018).</footnote>

<footnote n="49">⁴⁹H. Jaeger and H. Haas, "Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication," *Science* 304, 78–80 (2004).</footnote>

<footnote n="50">⁵⁰S. Pandey, J. Schumacher, and K. R. Sreenivasan, "A perspective on machine learning in turbulent flows," *J. Turbul.* 21, 567–584 (2020).</footnote>

<footnote n="51">⁵¹C. W. Gear, J. M. Hyman, P. G. Kevrekidis, I. G. Kevrekidis, O. Runborg, and C. Theodoropoulos, *Equation-Free, Coarse-Grained Multiscale Computation: Enabling Microscopic Simulators to Perform System-Level Analysis* (Johns Hopkins University, 2003).</footnote>