Learning continuous models for continuous physics

Dynamical systems that evolve continuously over time are ubiquitous throughout science and engineering. Machine learning (ML) provides data-driven approaches to model and predict the dynamics of such systems. A core issue with this approach is that ML models are typically trained on discrete data, using ML methodologies that are not aware of underlying continuity properties. This results in models that often do not capture any underlying continuous dynamics -- either of the system of interest, or indeed of any related system. To address this challenge, we develop a convergence test based on numerical analysis theory. Our test verifies whether a model has learned a function that accurately approximates an underlying continuous dynamics. Models that fail this test fail to capture relevant dynamics, rendering them of limited utility for many scientific prediction tasks; while models that pass this test enable both better interpolation and better extrapolation in multiple ways. Our results illustrate how principled numerical analysis methods can be coupled with existing ML training/testing methodologies to validate models for science and engineering applications.


Introduction
Dynamical systems-systems whose state varies over time-describe many chemical, physical, and biological processes.Thus, understanding and describing these dynamical systems is important for many scientific and engineering applications.Dynamical systems can often be described by differential equations which evolve continuously in time, meaning that the domain of the solution spans a continuum [1].In such systems, the gap between any two timesteps can be subdivided into an infinite number of infinitely smaller timesteps.In practice, these systems are often identified via a finite set of discrete observational data, and there is a long history within scientific computing for dealing with this discrete-to-continuous gap: experimentally measuring scientific data at sufficiently fine timescales to resolve approximately-continuous dynamics of interest; formulating theory within function spaces of sufficient smoothness to guarantee certain continuity requirements; and developing numerical algorithms that come with appropriate stability and convergence guarantees.
In many scientific and engineering applications, we observe measurements that yield a series of discrete data points {x 0 , x 1 , x 2 , . . ., x N }, where each point is spaced apart by some timestep size ∆t.There are many techniques from ML and statistical data analysis to learn data-driven inputoutput mappings (G : x n → x n+1 ) that can provide an approximation for the next discrete timestep.One popular class of data-driven input-output mappings is given by neural networks (NNs).A NN, denoted as N , can be trained to predict x n+1 from x n by learning model parameters θ: However, when considering continuous dynamical systems, there are challenges with this approach.Most obviously, this approach does not learn a continuous function [28,29,30,31]; it simply learns a function that predicts subsequent discrete time steps.This is to be expected, as this model is optimized to make (discrete) point estimates, i.e., to predict solutions at specific (discrete) points.
For this reason, predicting future states of a dynamical system with this approach can result in compounding errors of the dynamics over time [32,19].
A related approach is to assume that the discrete data points can be modeled and described by a continuous differential equation of the form, where F is a function that describes the vector field.In some cases, there is an underlying true F , while in other cases it is simply a modeling assumption.A challenge is that we cannot derive F from first-principles in many situations.Instead, we can use a data-driven approach for modeling F .For instance, an arbitrary NN architecture N can be used to model the vector field F , dx(t) dt = N (x(t); θ).
It is often assumed or simply taken for granted that ODE-Nets and other ML methods for ODEs automatically capture some continuous dynamics, either of the system that generated the data or of some related system [40,41,42,43,44].
However, due to how ODE-Nets are trained, i.e., to predict solutions at specific (discrete) points, these models can easily fail to learn even the simplest continuous dynamical systems [45,46], even when they accurately fit the temporal discretization (i.e., the discrete training points and testing points).An ODE-Net that incorrectly learns a continuous model will simply provide high-quality discrete time predictions -i.e., it is not a ContinuousNet but is simply a very good DiscreteNet [45].
Such a model will fail to extrapolate to new data points outside the temporal discretization, and it will fail to interpolate the solution at timesteps in between the discrete training data.It can also fail   Learning to predict discrete points versus learning continuous dynamics.(a) A NN model learns from a discrete set of points (red dots) that are generated by a dynamical system with an underlying continuous trajectory (dashed black line).(b) After training, the model should be able to predict future data points that are lying on the same or a different trajectory, including data points that are irregularly sampled.(c) A model that has only learned to predict the discrete data points might accurately fit future points on the same trajectory sampled at rate ∆t (shown in red), but fail to predict points sampled at different rates.The blue and green lines evaluated with h < ∆t fall off of the underlying continuous trajectory.(d) A model that has learned an underlying continuous dynamics can converge to a continuous solution as h → 0. In this case, the blue and green lines evaluated with h < ∆t get closer to the underlying continuous trajectory.
to correctly identify qualitative long-term behavior such as bifurcations [47].As we demonstrate later, this Discrete-versus-Continuous distinction affects non-NN ML methods as well, even when they accurately fit the temporal discretization of the data.
Figure 1 illustrates the difference between a model that has learned to predict discrete data points and a model that has learned an underlying continuous dynamics.After training a model at a given discretization ∆t, the trained model can be used to predict trajectories at arbitrary timestep sizes h during inference.For validation, an error metric Error(h) can be defined over a holdout trajectory that allows for evaluation with discretizations that are different than the data spacing.Learning a discrete-only model (a DiscreteNet) means that only the discrete training points-and potentially testing points at the same discretization (i.e., when h = ∆t)-are learned.When evaluating using Error(h = ∆t), the model will appear to perform well.The model may perform well on testing points with a similar discretization, but it will perform poorly for points sampled with other discretizations; that is, Error(h ̸ = ∆t) can be much worse than a discrete-only testing methodology would determine.This will even occur when the discretization h → 0, counter to expectations.In contrast, learning a meaningfully continuous model (a ContinuousNet) means that the model can converge to a smooth solution as the discretization h → 0, or at least that its error will decrease gradually and level off as h → 0. (This will be true regardless of whether the learned continuous model corresponds to the true underlying model, even assuming that such a true model exists and/or is well-conditioned.)In this case, the model will perform well for a broader range of temporal discretizations and thus have a better approximation of the continuous dynamical system.
In this work, we adapt methods from numerical analysis theory to develop a methodology to verify whether an ML model has learned a meaningfully continuous function that describes a dynamical system of interest.Specifically, we introduce a modified convergence test to verify and validate whether a model has learned continuous dynamics for a physical system.Our method allows us to verify that a model approximates a continuous differential operator, rather than only learning discrete points at a given temporal discretization, in the same sense that discrete algorithms from numerical analysis can be said to approximate continuous functions.We also introduce the notion of a ContinuousNet to refer to an ODE-Net model that exhibits the convergence properties that are expected for a continuous time system: trained with a numerical integration scheme is a ContinuousNet if it is convergent to a similar error as the error obtained by using the original training time step, This ContinuousNet convergence criteria is very similar to that of numerical analysis, whereby convergence is judged as h → 0, and can be evaluated with a similar methodology adapted for the ML setting.The criteria of Eq. 4 represents a heuristic that takes into account the error from the learning process that can be observed in the error on the discrete-only validation task, Error(∆t).
(Note that we are not claiming that this method will guarantee that we have learned the true solution-that would require additional assumptions, well-known in scientific computing-simply that we have learned some underlying continuous model of the data.) To illustrate the utility of our approach, we demonstrate how meaningfully continuous models that pass our convergence test enable both better interpolation and better extrapolation in multiple ways.We show that such models can resolve fine-scale features of the solution, despite being trained only on coarse data, including data that are irregularly spaced with non-uniform time intervals; can learn higher resolution solutions through learning continuous temporal dynamics from flow field snapshots; and can correctly predict trajectories starting at different initial conditions on which the model was not trained.We also demonstrate that our convergence test method is generally applicable to ML models.In addition, we derive theoretical error bounds for simple linear ODEs.Our results show promise in bridging between ML methodologies and scientific computing methodologies, by respecting both the fundamentals of ML and the fundamentals of science.

Results and Discussion
In this section, we use the convergence test to demonstrate and identify discrete-overfitting of dynamics models.We start by showing an example of our convergence test on a simple harmonic oscillator system (Example Convergence Method).We then illustrate our convergence test on a variety of different scientific systems, demonstrating that our method can validate whether a trained  The sharp dip at h = 0.1 demonstrates that the model achieves low error when h = ∆t, i.e., it is a good discrete model (DiscreteNet) when the evaluated step size is the same as the ∆t in the training data.However, when the model is evaluated at an h even slightly larger than or smaller than ∆t, the error increases sharply; and thus the Euler-Net model does not pass the convergence test.This is in contrast to numerical integrators, where errors decrease monotonically as h decreases.(c) We visualize a trajectory where the Euler-Net is evaluated at a h one order of magnitude lower (h = 0.01) than the trained ∆t.The resulting trajectory shows that the Euler-Net is unable to follow the baseline numerical Euler solution.(d) Schematic of an ODE-Net block, in which the RK4 numerical integration scheme is used to obtain the next timestep.As before, the RK4-Net is trained on discrete data points.This time, as (e) shows, the error converges monotonically to a fixed value as h decreases; and thus the RK4-Net model passes the convergence test, i.e., it is a good continuous (ContinuousNet) model.The reason the RK4-Net error flattens and converges to a (non-zero but small) fixed value is due to ML-based error sources [48].The RK4 numerical integration scheme also converges to a (non-zero but very small) fixed value-this time, due to numerical-based round-off errors [49].In (f ), the RK4-Net follows behavior similar to the RK4 numerical integration scheme when evaluated at a low h (h < ∆t).
dynamics model has learned (some) meaningfully continuous dynamics (Four Prototypical Dynamical Systems).Next, we show that models that pass this test can predict fine-scale solutions from coarsely spaced data.This includes: predicting continuous temporal dynamics from flow fields (Interpolation: Predicting Fine-scale Solutions from Coarse Training Data); predicting trajectories starting at initial conditions on which the model was not trained (Extrapolation: Predicting Trajectories for New Initial Conditions); and predicting fine-scale solutions from coarse, irregularly spaced data (Irregularly Sampled Training Data).We then show that overfitting to the temporal discretization affects ML methods more generally than just with ODE-Nets (Sparse Identification of Nonlinear Dynamical Systems).
Here, we only consider systems which are non-chaotic, non-divergent, and that are not extremely stiff, such that they can be handled by simple explicit Runge-Kutta integrators.More generally, learning the underlying "true" dynamics would require a test that involves a more sophisticated coupling of numerical and ML methodologies.This is not necessarily needed for many scientific ML tasks, including any of the improvements for predicting fine-scale solutions from coarsely spaced data that we discuss.

Example Convergence Method
We demonstrate our ContinuousNet convergence test on a toy example.We sample discrete training data points from the linear differential equations describing the harmonic oscillator: We show the results in Fig. 2. Two ODE-Nets are trained on this data.Here, Fig. 2(a)(b)(c) use the forward Euler integration scheme, while Fig. 2(d)(e)(f) use the RK4 integration scheme.
Both the Euler-Net (Fig. 2(a)) and RK4-Net (Fig. 2(d)) use the same linear network architecture N (x; θ) = θx to approximate the ODE.In theory, a linear model can exactly represent the linear ODE; however, we will demonstrate that this does not happen.For the example in Fig. 2, the training data was generated from the analytical solution, spaced apart by the training timestep ∆t = 0.1.To measure the performance as a continuous model at inference time, we integrate the model using a range of inference timesteps h.For the Euler-Net, the error when h = ∆t (the step size is equal to the temporal spacing in the training data) is very low, but it increases when h decreases.This is in contrast to the classical Euler numerical integration scheme, where the error decreases as h decreases.Thus, these results for Euler-Net do not pass the convergence test.In contrast, for the RK4-Net, the error decreases as h → 0, and eventually it approaches and levels off at a fixed value.Notably, the error does not increase dramatically as it does with the Euler-Net.In this case, the RK4-Net has learned the right inductive biases to approximate an underlying continuous dynamics for the system.
We illustrate this further by showing an example trajectory at a specific evaluated h.In this case, both trained ODE-Nets are evaluated at h = 0.01 (a 10× increase in resolution in comparison to the training data) up to a final timestep.In Fig. 2(c), the Euler-Net falls off of the true numerical Euler solution.It has clearly not learned the underlying continuous dynamics.In contrast, in Fig. 2(f), the RK4-Net shows good correspondence with the true numerical RK4 solution.

Four Prototypical Dynamical Systems
We consider canonical scientific dynamical systems: the non-linear pendulum, the Lotka-Volterra equations, the Cartesian pendulum, and the double gyre fluid flow.The first two systems are nonlinear dynamical systems; the Cartesian pendulum is a stiff dynamical system (which is difficult to solve with numerical methods without taking very small timesteps); and the double gyre fluid flow consists of vorticity fields describing a stream function.We provide more details about these dynamical systems in Supplementary Note 1: Details about considered dynamical systems.For each system, we sample data points from either the analytical solution or the numerical solution.
The temporal spacing between the discrete data points is denoted as ∆t, while the step size used to evaluate a trained ODE-Net is denoted as h.
Training setup.We train an ODE-Net with a numerical integration scheme (Euler or RK4) for each system.We use simple feed-forward networks with tanh activation functions.See Supplementary Note 2: Model architecture details for details on the architecture used.In every example, the exact same network architecture is used for both the Euler-Net and RK4-Net, respectively.We also include additional results for ODE-Nets trained on training data spaced apart by different ∆t as well as an ODE-Net trained with the Midpoint numerical integration scheme, in Supplementary Note 3: Additional convergence test examples using ODE-Nets.For the double gyre fluid flow, we use a dynamic autoencoder architecture [28,30] to embed the high-dimensional input of flow field snapshots in some latent space.Specifically, we replace the linear discrete map in the architecture proposed by [28] with a linear ODE-block.This means that the model learns to predict the next timestep by integrating forward in latent space (using an Euler or RK4 numerical integration scheme) with step size h = ∆t.Finally, the decoder translates the latent space vectors back to the flow field.
Results.The results of our method are shown in Fig. 3.In each case, the Euler-Net has low error when h = ∆t (i.e., evaluated at the same time spacing as the training data), but it has high error when evaluated at all other h, in particular smaller values of h.Thus, it does not pass the convergence test, and it has not learned a meaningfully continuous dynamics.It is a good discrete model, appropriate for data drawn from the same temporal discretization, but it has overfit to the temporal discretization.In contrast, the error during inference time of the RK4-Net steadily decreases when it is evaluated at lower h, eventually converging to a fixed basal level determined by the model and the noise properties of the data.It has passed the convergence test, and it can be said to have learned a meaningfully-continuous model.We include additional convergence test results in Supplementary Figures 9-12.

Interpolation: Predicting Fine-scale Solutions from Coarse Training Data
Observational, discrete training data are limited in that they are measured at specific timesteps.To obtain a solution for the system in-between these timesteps, one must retake the data measurements again at finer timesteps.However, selecting a model that has learned meaningfully continuous dynamics should guarantee accurate evaluation at smaller timesteps, despite only training on coarse and/or irregularly spaced temporal data (i.e., measurements taken with large timesteps).By learning continuous dynamics, the trained ODE-Net model can be evaluated at any point in temporal space, and still yield a low error solution.In this case, one would not need to recollect training data with smaller ∆t between data points; the learned ODE-Net can be used instead.Here, we demonstrate that fine-scale evaluation is possible by learning continuous temporal dynamics from flow fields for the double gyre flow example.
Results.We consider two models: the Euler-Net which did not pass the convergence test, and the RK4-Net which did pass the convergence test (see Fig. 3).In Fig. 4, we show the flow field  In each case, for the Euler-Net, the error does not monotonically decrease: error is low when h = ∆t, but high at all other h evaluations.Thus, the Euler-Net does not pass the convergence test.In contrast, the RK4-Net error does monotonically decrease and converges to a fixed value, when the evaluated h → 0. The RK4-Net does pass the convergence test.
snapshots that result from both models being evaluated at different timesteps.The Euler-Net is only able to approximate the true solution at the training data timestep (in this case, h = ∆t = 0.5).It cannot match the true solution at the other timesteps, and it gives a poor approximation that does not capture the flow behavior.In contrast, the RK4-Net has good correspondence to the true solution even when it is evaluated at timesteps that were not in the training data.Thus, our convergence test method has allowed us to choose a model that can recover fine-scale solutions of the system, while only having access to coarse-scale measurements during training.

Extrapolation: Predicting Trajectories for New Initial Conditions
For a given system, temporal trajectories start at some initial condition.Measurements are taken for one trajectory at one initial condition, and then must be taken separately for other trajectories with different initial conditions.Selecting a model that has learned a meaningfully continuous dynamics circumvents this: after training a model on data points sampled from one (or more) trajectories, the model should be able to extrapolate and predict accurate solutions for new initial conditions.Training setup.We look at the non-linear pendulum (described in Eq. 27 in Supplementary Note 1: Details about considered dynamical systems).Here, θ is the initial condition representing the position of the pendulum in time.The phase portrait of this system (representing the true solution trajectories), showing dθ dt against θ, is shown in Fig. 5(a).An Euler-Net and an RK4-Net are trained on trajectories, spaced apart by ∆t = 0.1, starting at certain initial conditions (shown by the black lines in Fig. 5(b)).We then pick a test set of a number of different initial conditions that were not in the training data.The Euler-Net and RK4-Net start at these test initial conditions and are both evaluated at a finer h = 0.001, representing a 100× increase in resolution.Note that we saw in Fig. 3 that the Euler-Net did not pass the convergence test (i.e., it had high error when evaluated at h ≪ ∆t), while the RK4-Net did pass the test.
Results.The results of predicting trajectories starting at different test initial conditions are shown in Fig. 5.The Euler-Net is unable to predict these trajectories and quickly falls off of the phase plot lines corresponding to the true solution.In contrast, the RK4-Net is able to predict the trajectories, starting at different test initial conditions with good correspondence to the true solution.Thus, we see that it is critical to find a model that passes the convergence test and is able to learn a continuous dynamics to succeed at this extrapolation task.

Irregularly Sampled Training Data
It is typically the case that scientific data collection includes measurements that are taken with some amount of imprecision.For example, the measurement of interest is not always taken at the exact same ∆t every time, due to issues such as jitter in the measurement device.Measurements may also be skipped: for example, a measurement is only available at t = ∆t and t = 3∆t because the measurement at t = 2∆t was lost or skipped.Thus, reconstructing the correct trajectory when the measurements are non-uniformly spaced is important in numerous science and engineering problems.Here, we look at an example of using the convergence test to correctly select a meaningfully continuous model in the case of the non-linear pendulum with irregularly spaced temporal data with non-uniform temporal intervals.training.)Both ODE-Nets are then evaluated at a very low h (approximately 100× lower than the general distribution of the training data points) to generate a time series plot.Note, that we run the convergence test on both ODE-Nets.
Results.The Euler-Net quickly falls off of the continuous solution (Fig. 6(b)).Conversely, the RK4-Net follows the continuous solution with good accuracy, including at timesteps not in the training data (Fig. 6(c)).Thus, it is clear that the RK4-Net has learned a meaningfully continuous dynamics while the Euler-Net has not.This is confirmed by RK4-Net passing the convergence test, but Euler-Net not passing it (Fig. 6(d)).The dip for Euler-Net appears at the average ∆t in the training data, which is slightly larger than 0.05 due to the measurement noise.

Sparse Identification of Nonlinear Dynamical Systems
Here, we demonstrate that ovefitting to the temporal discretization affects ML methods (in the context of dynamical systems) more generally.To illustrate this, we consider the SINDy learning approach, which is a class of methods for system identification [5].
The SINDy method uses the following model structure to represent the dynamics, where Ξ is a matrix of learnable parameters and ϕ(x) is a set of non-linear basis functions which correspond to potential terms in the underlying system.A linear optimizer is used to fit the parameters Ξ to the data, with a sparsity constraint.The sparsity constraint identifies the subset of relevant basis elements in ϕ(x) to reveal an interpretable dynamics model.
In most real applications, the time derivatives, dx n /dt, cannot be measured directly and instead need to be approximated from the observations x n .The common approach in SINDy is to use finite differences to approximate dx n /dt from the data [50,51].(This treatment of time derivatives, where SINDy differentiates the data, is in contrast to the ODE-Net method, which integrates the model.) The finite difference operator F D(x n... ) is applied over the dataset as a pre-processing step.This yields the following set of N discrete equations, which is optimized for Ξ over all observations n: where F D(x n... ) is a finite difference approximation using the region of points around time index n.
Using the series of N equations, Ξ is learned using specialized algorithms designed to seek sparsity, such as LASSO regularization or sequential threshold least squares [5].(This is again in contrast with the ODE-Net method, which uses gradient descent nonlinear optimization.)The discretization order of the finite difference pre-processing is a hyperparameter of SINDy.The first order accurate finite difference (here referred to as FD-1) results in a point-wise approximation of, Note that when rearranged, this is equivalent to the forward Euler integrator: Higher-order finite difference stencils can also be used to increase the accuracy of the time derivative approximation.The second-order stencil (FD-2) (analogous, but not equivalent, to Midpoint) can be expressed as and the fourth-order stencil (FD-4) (analogous, but not equivalent, to RK4) can be expressed as,  an ODE-Net model, SINDy is used instead to learn a sparse polynomial basis.The accuracy of the time derivative approximation is altered by using a different order of central finite difference stencils.Here, FD-1 is the first-order two-point stencil, FD-2 is the second-order three-point stencil, and FD-4 is the fourth-order five-point stencil.In this example, FD-1 is analogous to Euler-Net, FD-2 is analogous to Midpoint-Net, and FD-4 is analogous to RK4-Net.In (a), the convergence test method is shown on the harmonic oscillator, where the convergence test is replicated by training via the SINDy method.In (b) and (c), the convergence test is applied when training on the non-linear pendulum via the SINDy method.We observe similar results as when training using an ODE-Net: the lower-order approximation (FD-1) overfits to the ∆t in the training data, failing the convergence test; and the higher-order approximation (FD-4) passes the convergence test.
Training Setup.We look at the harmonic oscillator described in Eq. 5, and the non-linear pendulum (described in Eq. 27 in Supplementary Note 1: Details about considered dynamical systems).We use the PySINDy implementation [50,51] to train models on these trajectories.We train three different SINDy models on the trajectories, altering the finite difference (FD) approximation order of accuracy: FD-1 is a first-order two-point stencil (analogous to Euler-Net), FD-2 is a second-order three-point stencil (analogous to Midpoint-Net), and FD-4 is a fourth-order five-point stencil (analogous to RK4-Net).The SINDy model is plugged into the convergence test as F , such that the same Runge-Kutta integrators are used for trajectory prediction.For the harmonic oscillator, the training data is spaced apart by ∆t = 0.1, while for the non-linear pendulum, we look at examples where the training data is spaced apart by ∆t = 0.05 and ∆t = 0.1.
Results.We run our convergence test on the different SINDy models.The results of our method are shown in Fig. 7.In each case, the FD-1 model has low error when h = ∆t (i.e., evaluated at the same time spacing as the training data), but it has high error when evaluated at all other h, and especially smaller values of h.Thus, it does not pass the convergence test, and it has not learned a meaningfully continuous dynamics.The FD-1 model shows a sharp dip because it is overfit to forward Euler.In this case, the stencil and integrator correspond to the exact same algebraic structure.Similar to when the models where trained via ODE-Nets, we see that FD-1 has overfit to the temporal discretization.In contrast, the error during inference time of the FD-4 model steadily decreases when it is evaluated at lower h, eventually converging to a fixed basal level.It has passed the convergence test, and it has learned a meaningfully-continuous model.

Conclusion
One of the great challenges in scientific ML is to learn continuous dynamics for physical systemseither "the" underlying continuous dynamics, or "a" continuous dynamics that leads to good predictive results for the spatial/temporal regime of interest for the ML model-such that the learned ML model can be trusted to give accurate and reliable results.ML models are trained on discrete points, and typical ML training/testing methodologies are not aware of the continuity properties of the underlying problem from which the data are generated.Here, we have developed a Contin-uousNet methodology, and we showed that convergence (an important criteria used in numerical analysis) can be used for selecting models that have a strong inductive bias towards learning meaningfully continuous dynamics.Standard ODE-Net approaches, as well as common SINDy methods, both popular in recent years within the ML community, often do not pass this convergence test.In contrast, models that pass this convergence test have favorable properties.For instance, models that learned underlying continuous dynamics can be evaluated at lower or higher resolutions.Our results suggest that principled numerical analysis methods can be coupled with existing ML training/testing methodologies to deliver upon the promise of scientific ML more generally.
Many more concrete directions are of course raised by our methodology.One direction has to do with developing analogous tests to be used for less well-posed dynamical systems.Such systems are of interest in scientific ML, and such tests will be of greatest interest when one needs to obtain "the" correct underlying continuous solution (e.g., to identify correctly qualitative long-term behavior [47]), rather than "a" continuous solution, which is often sufficient for ML prediction tasks.Another direction has to do with whether we can develop analogous tests appropriate for adaptive time-stepping methods, symplectic integrators, and other commonly used numerical simulation methods such as those for optimal control problems using NNs and associated Hamilton-Jacobi partial differential equations [52,53].Work subsequent to the posting of the initial technical report version of this paper has addressed the continuous-discrete equivalence question for learning operators [54,55], and likely our ContinuousNet methodology provides a way to operationalize that in practice.A final direction has to do with whether one can obtain strong theoretical results, e.g., ML-style generalization bounds, to guide the use of methods such as these.Recent theoretical and empirical results suggest that this will be challenging [56,57,58,59,60], at least when using traditional approaches to ML-style generalization bounds.Our success in combining principled numerical analysis methods with existing ML methodologies also leads one to wonder whether we can use a posteriori error bound analysis methods to develop practically useful a posteriori generalization bounds for problems such as those we have considered.

Methods
The basic problem of numerical analysis is to solve problems from continuous mathematics using a discrete computer.The area has a rich history for describing the consistency and convergence behavior of numerical methods for approximating continuous functions [61,62].Here, we expand on the methods we used in Results and Discussion.

Criteria of Classical Numerical Analysis
Given an initial value x(0) = x 0 , we can discretize Eq. 2 along the node points t n = n∆t for n = 0, 1, . . ., N by evaluating the following integral equation: where x n = x(t n ), and ∆t is the discrete timestep.Typically, we are not able to compute an analytic solution for the integral, and thus we rely on numerical schemes to approximate tn+∆t tn F (x(s)) ds.
There are many different types of numerical integration schemes to approximate the integral in Eq. 12.These have different trade-offs between computational efficiency and accuracy.One such scheme, the forward Euler discretization, can be written as: This is a first-order one-step method, where the global error (the error over all of the timesteps) is proportional to the step size, i.e., O(∆t), meaning that the error gets smaller as ∆t decreases.There are also higher-order integration schemes.One popular higher-order scheme is the Runge-Kutta 4 (RK4) discretization, which takes the following form: Here, the global error is proportional to the step size to the fourth power, i.e., O(∆t 4 ); and thus as ∆t gets smaller, the error gets smaller much more quickly than with the forward Euler scheme.
In general, the global error can be written as O(∆t p ), where p denotes the order of accuracy.
Classical numerical integration typically starts by assuming that there exists a true underlying continuous-time system, which is then replaced by a discrete-time problem whose solution approximates that of the continuous problem.However, discretizing the problem introduces an error, and concepts such as stability, convergence, and consistency can be used to quantify the error of the discrete solution [63,64,65].
In the following, we describe the error bounds in the traditional scientific computing context where the system dynamics are known exactly, in which case the only approximation error comes from numerical integration in time.We specifically focus on numerical convergence because this will give us a mechanism to analyze ML models.However, note that stability and consistency are also of interest [62].Let x(t n ) denote the true solution of a dynamical system of interest; and let x∆t n denote a numerical solution after n steps with step size ∆t.We use N to denote the maximum number of time steps such that T = t N = N ∆t is the final time.Decreasing ∆t requires increasing N (the number steps taken to arrive at T ), and vice versa.Then, x(T ) = x(t N ) is the true solution at the final time, while x∆t N is the numerical solution at the final time.Convergence quantifies the global error (the cumulative error of all iterations) of a numerical algorithm.
Definition 2 (Convergent Numerical Approximation.).A numerical one-step method for solving dx(t) dt = F (x(t)), with initial condition x(0) = x 0 , is said to be convergent if and only if the error tends to zero as ∆t goes to zero: Of course, in numerical practice, the error does not converge to zero.Instead, it levels off at some base level determined typically by the level of numerical precision used to describe the data, as observed in 2(e).
The specific metric for quantifying the approximation error across the sequence is somewhat arbitrary.Moreover, there is the problem that the numerical method can potentially converge to the wrong solution [66].Thus, to ensure that a numerical method is not only convergent but also consistent, one can use the mean error, or the maximum error across all N points in time, If, as the step size ∆t decreases, the largest absolute error between the numerical solution x∆t n and the exact solution x(t n ) also decreases, then the numerical approximation converges towards the solution of the continuous system.In the limit of ∆t → 0, the numerical solution converges to the exact solution and the error converges to zero, or to some base level determined by machine precision and numerical round-off noise.This convergence criteria is also a test for continuity in the solution: as ∆t → 0, the time interval between adjacent numerical solutions (e.g., at t n , x n , and at t n+1 , x n+1 ) also decreases towards zero.Thus, the numerical solution collapses onto a continuous solution as ∆t → 0.
Remark 1. Validation of a new integration method involves multiple stages.Consistency, convergence, and stability can be theoretically proven for a rather small class of ODEs (typically only linear ODEs).Thus, the method will be evaluated empirically with a real implementation on a problem of interest.In practice, a convergence test is used, where the numerical integration scheme is used to predict trajectories for a range of ∆t and compared to an analytical solution or to an overrefined solution.The errors are verified to approach zero at the correct rate, at least until they flatten out at some base level.The combination of theoretical proof of consistency, stability, and convergence on simple systems such as linear ODEs, combined with the emprical demonstration of convergence on ODEs of interest, is typically viewed as sufficient to vet the method.Emprirical convergence tests are a standard integration test method for scientific programs, e.g., they are regularly run to automatically catch bugs.

A Convergence Test for ODE-Nets
We now describe a convergence test, based on the discussed convergence criteria, to validate properties of an ODE-Net solution.The fact that ODE-Nets are embedded in a numerical integration scheme enables us to use convergence analysis methods that are well-known for studying classical numerical analysis problems.To start, we know that the numerical integrator itself will be convergent if it is given the true f , but we do not know if the ODE-Net will be convergent when an approximate ML model N is used to approximate f .The convergence test is used to determine whether an ODE-Net has learned a meaningfully continuous model for the underlying problem of interest.
Suppose that we are given an ODE-Net N that is trained with a numerical integration scheme (such as Euler or RK4) from t n to t n+1 with stepsize ∆t: Following Definition 2, we compute the global error of the ODE-Net N as it approaches some fixed value b as the time step h goes to zero: Here too, in the ML setting, the error does not necessarily converge to zero.This is analogous to the classical numerical analysis setting, where the numerical analysis test typically converges to a non-zero value determined by the numerical round-off error (e.g., see the floor in Figure 2e).Unlike in the classical numerical analysis setting, even in the absence of numerical errors the b of an ML model will be greater than zero.For an ODE-Net, the numerical value of b depends on the model architecture, integration method, the optimizer, and the noise properties of the data [48].
The value of b will elucidate the convergence properties of the trained ODE-Net.
Computing an error metric in the ML setting requires additional consideration because we do not necessarily have access to the underlying exact solution at arbitrary points in time.Instead, we are restricted to the information that is provided by a given validation set of discrete data points T = {x 0 , x 1 , x 2 , . . ., x N }, where each point is spaced apart by the ∆t between observations.A naive metric to compute the global error in this setting is simply to consider the 2-norm between the end point of the validation trajectory and the predicted value: However, this metric is susceptible to noise and edge cases.Computing the error over all points T is difficult using (19) because the inferred trajectory has a different number of points than the validation trajectory.To mitigate this issue, we suggest to compute the global error on a subset of points S from the validation/test trajectories, which is a set of indexes into the original dataset T ; e.g., S = {0, 9, 18 . . .} for every 9 points spaced by 9∆t.This allows for inferring the trajectory of h computing the error over the subset as follows: where k is the index into the inferred trajectory corresponding to the index into the validation trajectory n such that xk k is the point that lines up at the same time as x n .Note that we can only use certain timesteps h during inference because the solution points must align perfectly with those in the subset trajectory.
Given this setup, we use the term ContinuousNet to refer to an ODE-Net model that also exhibits these convergence properties, as per Definition 1.To evaluate this property, we can apply the same convergence test procedure used in traditional numerical analysis and scientific computing, but with the modifications necessary for it work on training data.Further, it is necessary to apply a weaker heuristic to judge convergence because there will be residual optimization error at Error(∆t), as per Eq. 4. The procedure is as follows.
Given the learned ODE-Net, we first infer a validation/test trajectory on the original stepsize ∆t and evaluate the global error, Error(∆t).This is the standard procedure for evaluating a discrete model, and this error value informs its accuracy at inferring discrete points in a sequence.Then, to further evaluate whether the model is convergent and continuous, we consider a range of h values which are both smaller and larger than the step size ∆t used during training.For example, if the ODE-Net was trained on ∆t = 0.1, we evaluate the ODE-Net on the validation/test trajectory for h ∈ [10 −3 , .., 101 ].Specifically, for a given set of inference timesteps {h 1 , h 2 , . . ., h p }, our proposed convergence test iterates over the elements h i , and executes the following two steps: 1. Evaluate the pre-trained ODE-Net using Eq.16 on the time interval [t 0 , t T ], using step size h i .

Calculate the error between the ODE-Net solution xh i
n and points in the test data, Error(h i ).Algorithms 1 and 2 in Supplementary Note 4: Algorithm for the Convergence Test summarize this procedure.As with many other numerical convergence tests, our proposed algorithm is subject to a heuristic threshold.In practice, we observe that the difference between a model that passes our test condition and one that does not is pronounced.A non-continuous model results in an error that is orders of magnitude larger, as compared to a continuous model, when h is taken smaller than ∆t of the training data.In other words, our convergence test performs a form of model selection by selecting for models that learn inductive biases towards meaningfully continuous dynamics.In the following, we will demonstrate why convergence is an important consideration for ML model selection. 1  Remark 2. In practice, we can also evaluate our convergence test with different starting points x 0 and the respective final time step, x N , and then average the error across the different runs.We highly recommend this to ensure that the same behavior occurs irrespective of start and end point.

Error Analysis in an Idealized Learning Setting
We provide a theoretical framework for the convergence test, analyzing the discretization error of one-step numerical integration schemes in an idealized setting.Specifically, we consider the problem of learning the simple scalar linear ODE,  where λ ∈ R denotes a scalar parameter, and x ∈ R denotes the state at a given point in time.It is well known that the function x(t) = e λt x 0 is a solution of this system, given the initial condition x(t) = x 0 ∈ R [67].In the following, we assume a similar setting as before: we are given a set of discrete data points D = {x 0 , x 1 , x 2 , . . ., x N } produced by the linear system and spaced by ∆t.Our aim is to learn a scalar ODE-Net model, parameterized by the learnable weight parameter w.Following the ODE-Net process, we discretize the model with a numerical integration scheme and then optimize the squared error loss, min w x∈D We assume that the data are noise-free and can therefore be represented by its analytical solution, x n+1 = e λ∆t x n .When the loss is optimized, the time-discretization step introduces its own unique source of error into the learning process, one that is independent of noise, numerical error, or optimization error.The error stems from the fact that any one-step consistent numerical integration scheme, when applied to a linear ODE, will result in a truncated Taylor series expansion with p terms, where p is the accuracy order of the scheme [62].Thus, the ML model cannot recover the exact parameters of the underlying ODE.The following lemma makes this issue explicit.
The proof is given in Supplementary Note 6: Optimization of w does not learn λ.For finite ∆t, this equation clearly satisfies w ̸ = λ if p ≪ ∞.There are situations where this equation can be solved for values of w that will set the loss in Eq. (22) to zero for all possible data points x n .In the limit as ∆t → 0, there is always a solution at w = λ for any p.Moreover, for practical settings there is at least one root when p is odd for any ∆t; when p = 2 (RK2) and ∆t ≤ log(2)/|λ|; or when p = 4 (RK4) and ∆t < 1.307/|λ|.This equation can be used to find analytical expressions for w for simple integrators; see Supplementary Note 7: Example evaluations of the analytical equations.
From this result, we can characterize the difference between the ML model and the target ODE.In addition, in practice, it is (almost always) only possible to learn a perturbed version, w, of w due to noise, limited numerical precision, and optimization errors.Let ε denote an additive perturbation away from the the optimum of the minimization problem for an observed model due to these sources of errors, w = w + ε.The following theorem bounds the error between the ML and ODE model, due to the overfitting of Eq. ( 23) in presence of additive error sources.
Theorem 1.If optimal weight parameter can be found by Eq. ( 23), the approximation error introduced by scalar ODE-Net is bounded by |w − λ| ≤ c∆t p , where c is a constant proportional to the Lipschitz continuity constants of λx and wx.In the presence of additive numerical error ε, the bound is, The proof is given in Supplementary Note 8: Approximation errors introduced by scalar ODE-Nets.This bound shows that the user needs to increase the order of accuracy of the training scheme p in order to reduce the error between the learned parameter and the true ODE parameter.
In practice, we can only measure Error(h) using a set of data points.Using the above results, we can analyze the expected behavior of the convergence test by bounding the global error using Eq.(24).Fig. 8 plots the global error bound given concrete values of λ, ∆t, and ε.As can be seen, the theoretically derived global error for the scalar ODE exhibits the same behavior as the empirical applications of the convergence test.We can use the global error to further derive expected bounds on the the key points of the convergence test.
Corollary 1.When a scalar ODE-Net is evaluated with the timestep that was used for training, the leading term of the global error is proportional to the optimization error ε for a k ∝ T |x 0 |: The proof is given in Supplementary Note 9: Observable quantities and the convergence test.The c∆t p error term in Eq. ( 24) between w and λ is cancelled out, and the observed value is smaller than | w − λ|.Therefore, the global error only observes the difference between w and w, resulting from the model optimization error.Note that this value can become very small as the optimization error decreases.By applying the convergence test, we are able to extract an estimate of | w − λ|, as given in the following.
Corollary 2. In the limit of decreasing the timestep size during inference, the global error approaches a constant factor (b) based on the bound in Eq. (24).It approaches at a rate of h q , where q is the accuracy order of the ODE integration scheme used at inference time: The proof is given in Supplementary Note 9: Observable quantities and the convergence test.Given these bounds, we can see how using Eq. ( 4) as a threshold yields (b−Error(∆t)) ∝ c∆t p .Therefore, the comparison between b and Error(∆t) allows for the quantitative estimation for the magnitude of the term c∆t p , which describes the error that is induced by the numerical discretization scheme used for training.
In summary, our analysis shows that there are two types of errors in the process of learning the dynamics of a linear scalar ODE using ODE-Nets.These errors can be measured using the data by evaluating Error(∆t) and lim h→0 Error(h).This illustrates the power of our proposed convergence criterion Eq. ( 4).Moreover, even if ε is small, it is required to increase the order of accuracy of the training scheme in order to further decrease the error between w and λ.
The Cartesian pendulum is an example of a stiff ODE, where numerically solving such a system can easily be numerically unstable without special consideration.We generated training data using the analytical solution to in the pendulum in θ(t) from [68].
Double gyre fluid flow.The double gyre fluid flow is a benchmark model for learning temporal dynamics from fluid flow field snapshots.The system dynamics is described by unsteady Stokes equations, and has a closed form analytical solution using the stream function, ϕ(x, y, t) = A sin(πf (x, t)) sin(πy), The velocity field is computed from the stream function, different ∆t.The Cartesian pendulum is an example of a stiff system, so the baseline Euler numerical integrator only starts to converge when the evaluated h is small.Here, we see a similar pattern: at very low ∆t = 0.001, Euler-Net fails the convergence test while Midpoint-Net and RK4-Net pass the test.Once ∆t = 0.1, Midpoint-Net also fails the convergence test.Note that the ODE-Nets, when they pass the convergence test, all have lower error than the numerical integrator for numerous evaluated h.Since the Cartesian pendulum is a stiff system, the ODE-Nets can fail the convergence test at smaller ∆t than the other systems of study.
For a linear ODE, it is possible in certain configurations to obtain a w that sets this residual to zero for all points x n in the dataset: ∃scheme, ∆t, w s.t.0 = e λ∆t x n − scheme(wx, ∆t)(x n ) ∀x n .
Let p denote the order of accuracy of the scheme used for training.A consistent numerical integration scheme will expand a linear ODE into the truncated Taylor series polynomial with p terms.Then, the residual equation is, For example, Euler, p = 1, yields a linear equation; for RK2, p = 2, yields a quadratic equation; and for RK4, p = 4, the equation is a quartic equation.It is thus apparent that it is possible to obtain a w that solves the equation, given certain p and ∆t, for any value of λ.For odd p, there is always such a root.For any p, this equation has the solution w = λ in the limit ∆t → 0. For p = 2 and p = 4, the analytical solutions of the roots yield the following conditions for there to be such a root.For p = 2, the quadratic formula discriminant condition is 1 − e ∆tλ > 0, yielding ∆t ≤ log(2)/|λ|.For p = 4, the discriminant of the quartic solution yields the condition ∆t < 1.307/|λ| (where we approximated the rational expression for brevity).
■ Equation ( 39) can be solved analytically up to p = 4 (see the following); and, beyond that, it can be solved with numerical root finding.It is thus apparent that w is not equal to λ in general, even with error-free optimization in the absence of data noise.That is, λ does not minimize the training loss in the general case.
Note that the solution in Eq. ( 39) can also be adapted for linear ODEs on other domains.For example, for scalar ODEs with complex valued parameters, there is always a solution for w, as complex polynomials always have complex roots.Additionally, for any dimension of vector-valued x with matrix w and λ, there will also always be a solution for p = 1 (Euler).The following results continue to assume real scalar parameters for simplicity of presentation.
yielding the error between learned parameter and ODE parameter, The global convergence error at inference time for any test trajectory for any h is, When h = ∆t, we see that Error euler (∆t) = |ε|.The value for b is, Taking a series expansion in terms of ∆t, We can summarize that term as, Therefore, training using a forward Euler integrator has a convergence error that is linear with respect to both ∆t and ε.Comparing the difference between Error(∆t) and b allows us to estimate empirically that |w − λ| = O(∆t).
Example: RK2 (Explicit Midpoint).For RK2, p = 2, and applying the recursive rule gives us a quadratic polynomial which can be solved analytically, This equation has real roots when ∆t ≤ log(2)/|λ|.Note that there can be two roots; for higher order numerical integrators, there may be multiple parameters that minimize the ML loss function perfectly.The difference between this value and the underling λ can be expanded into a series as, Plugging back into Error(h), we obtain, The error at the training timestep can be expressed by a series expansion, • If λ < 0, then a 0 > 0.Then, all coefficients are positive, and thus, by Descartes' rule of signs, all roots are negative, and w < 0. We will also require a lower bound on w.For odd p, we can obtain a lower bound by observing that there is a sign change in the polynomial residual of Eq. ( 53) between w = 0 and w = λ.When w = 0, the residual is positive, When w = λ, we can use the Taylor remainder term, The sign of this term is opposite the sign of λ p+1 , thus is always negative when p is odd, and positive when p is even.Therefore, for odd p there is a sign flip of the residual between 0 and λ, and we can therefore bound λ < w odd p < 0. For even p, finding such a loose bound is more difficult, as there exist values of λ and ∆t with no solution for w.However, we are assuming combinations of λ, ∆t, and p with a solution for w.Then, for even p, solving for w with p − 1 and p + 1 both satisfy λ < w.By assuming a solution for w even p , we assume that the error |w − λ| solution can be bounded from above by maximum error of the p − 1 or p + 1 cases.
That is, if there is a solution for an even p, we assume that it is not asymptotically less accurate than a constant bound on w p−1 or w p+1 .We can therefore assume the bound aλ < w even p < 0 for a small constant a > 1.Thus, the bound aλ < w < 0 with a small constant a ≥ 1 encompasses both even and odd p.
We will leverage these bounds to make a tighter bound in the following.
The left hand side can be expanded as so, We can use the factorization w i − λ i = (w − λ)(w i−1 + w i−2 λ + • • • + λ i−1 ) to pull out the difference between w and λ, T p (w∆t) − T p (λ∆t) = := K(w, λ)∆t(w − λ), where we define K(w, λ) to shorthand the double summation in the factorization.On the right hand side, we can substitute in the bound of Eq. ( 52), with δ between 0 and 1.Then, we have the expression for the difference between the w and λ, w − λ = e δλ∆t λ p+1 (p + 1)!K(w, λ) ∆t p .(63) We see the main term in ∆t p , but there is still a dependence on w in the denominator.We can now use the above bounds on w as a polynomial root to bound the minimum possible value of |K(w, λ)| to produce the upper bound on |w − λ|.For the two nontrivial (λ ̸ = 0) cases: 1.In case of λ > 0, the lower bound of w > 0 makes K(w, λ) the smallest it could be because it is monotonic, yielding the upper bound on w − λ: |w − λ| < e δλ∆t λ p+1 (p + 1)!K(0, λ) The expression K(w = 0, λ) with w = 0 is We can series expand the upper bound in terms of ∆t to see that this has a leading term in ∆t p , 2. In case of λ < 0, then w < 0. For p odd, we established that |w| < |λ|.For even p, we assume that there exists a constant a that bounds w, if Eq. ( 53) has a solution for w.We thus assume that there exists an a such that |w| ≤ a|λ| for a small constant a.The smallest-magnitude value of K will then occur at some value w = sλ for 0 < s ≤ a. Then we may bound K(w, λ) from the bottom: We can plug this back into expression back into the denominator of Eq. ( 39) to obtain an upper bound for |w − λ|.By series expanding, we can find that the leading term is the same as a above, Proof (Proof of Corollary 1).When we evaluate the global error using the same numerical integration as at training time, h = ∆t and q = p, the global error is dominated by only the numerical error ε.Firstly, we can substitute in the value of h = ∆t and q = p above, Recall that w was obtained by the ML process to solve the first term above exactly.Then, by canceling out a factor ∆t, we have the error bound, Error(∆t) = k|ε| + O(∆twε + ∆tε 2 ).(84) ■ When there is no numerical error, i.e., ε = 0, the observed error will be zero, Error(∆t) = 0, since this is the context in which the model was optimized.
Proof (Proof of Corollary 2).To obtain the fixed constant b to which the global error converges, as the timestep goes to zero, we can take the limit of Error(h) from Eq. ( 77 We can then invoke Taylor's theorem again to bound the remainder term: there exists a δ between 0 and wh such that R q (wh) = e δ (wh) q+1 /(q + 1)!.It follows that lim h→0 Error(h) = lim h→0 k e λh − e wh h + e δ w q+1 (q + 1)! h q . (89) The term multiplying h q is a bounded constant.The limit of (e λh − e wh )/h = λ − w.The limit of the error thus approaches the constant equal to the difference between the learned parameter and ODE parameter at a rate of h q , lim h→0 Error(h) = k|w − λ| + O(h q ).(90) (The rate h q can be observed and measured in the convergence test.)When there is an additional error ε, the limit is the same but for w appearing in the equation.Then, the error separates into two terms, lim h→0 Error(h) = k| w − λ| + O(h q ) ≤ k|w − λ| + k|ε| + O(h q ).(91) Substituting the error bounds of Theorem 2, the limit of the global error can also be represented by lim h→0 Error(h) ≤ kc∆t p + k|ε| + O(h q ).(92)

■
In standard numerical analysis, we would have set w = λ, and we would have proved that our numerical scheme converges to zero at a rate of q for linear ODEs.However, for a learned model, that is not the case.The approached value of b is the same for any q, meaning that any numerical integrator used when performing the convergence test will yield similar results in the limit of the approach value.

Figure 1 :
Figure 1: Learning to predict discrete points versus learning continuous dynamics.(a) A NN model learns from a discrete set of points (red dots) that are generated by a dynamical system with an underlying continuous trajectory (dashed black line).(b)After training, the model should be able to predict future data points that are lying on the same or a different trajectory, including data points that are irregularly sampled.(c) A model that has only learned to predict the discrete data points might accurately fit future points on the same trajectory sampled at rate ∆t (shown in red), but fail to predict points sampled at different rates.The blue and green lines evaluated with h < ∆t fall off of the underlying continuous trajectory.(d) A model that has learned an underlying continuous dynamics can converge to a continuous solution as h → 0. In this case, the blue and green lines evaluated with h < ∆t get closer to the underlying continuous trajectory.

Figure 2 :
Figure 2: Illustration of our convergence test with different ODE-Nets.(a) Schematic of an ODE-Netblock, where the next timestep is obtained with the Euler method.In (b) and (c), an Euler-Net is trained on discrete data points, spaced apart by some ∆t (in this case, ∆t = 0.1).(b) After training at one specific ∆t, Euler-Net is evaluated at many different step sizes, h, both larger and smaller than ∆t.The sharp dip at h = 0.1 demonstrates that the model achieves low error when h = ∆t, i.e., it is a good discrete model (DiscreteNet) when the evaluated step size is the same as the ∆t in the training data.However, when the model is evaluated at an h even slightly larger than or smaller than ∆t, the error increases sharply; and thus the Euler-Net model does not pass the convergence test.This is in contrast to numerical integrators, where errors decrease monotonically as h decreases.(c) We visualize a trajectory where the Euler-Net is evaluated at a h one order of magnitude lower (h = 0.01) than the trained ∆t.The resulting trajectory shows that the Euler-Net is unable to follow the baseline numerical Euler solution.(d) Schematic of an ODE-Net block, in which the RK4 numerical integration scheme is used to obtain the next timestep.As before, the RK4-Net is trained on discrete data points.This time, as (e) shows, the error converges monotonically to a fixed value as h decreases; and thus the RK4-Net model passes the convergence test, i.e., it is a good continuous (ContinuousNet) model.The reason the RK4-Net error flattens and converges to a (non-zero but small) fixed value is due to ML-based error sources[48].The RK4 numerical integration scheme also converges to a (non-zero but very small) fixed value-this time, due to numerical-based round-off errors[49].In (f ), the RK4-Net follows behavior similar to the RK4 numerical integration scheme when evaluated at a low h (h < ∆t).

Figures 2 (
b) and 2(e) plot the results of the convergence test with the Euler-Net and the RK4-Net.

Figure 3 :
Figure 3: Illustration of our convergence test on prototypical dynamical systems.We demonstrate the convergence test on multiple systems: (a) Non-linear pendulum (b) Lotka-Volterra system (c) Cartesian pendulum (d) Double gyre fluid flow.In each case, for the Euler-Net, the error does not monotonically decrease: error is low when h = ∆t, but high at all other h evaluations.Thus, the Euler-Net does not pass the convergence test.In contrast, the RK4-Net error does monotonically decrease and converges to a fixed value, when the evaluated h → 0. The RK4-Net does pass the convergence test.

Figure 4 :
Figure 4: Double gyre fluid flow: Reconstructing fine-scale flow fields from coarse training data.The training data for this problem consists of vorticity field snapshots of the dynamical system taken at ∆t = 0.5.In the images above, the red region is rotating in one direction, and the blue region is rotating in the opposite direction.After training an Euler-Net and an RK4-Net, both models are evaluated at different h (both when h > ∆t and h < ∆t) at a final timestep.We show the evaluation results for this final timestep, T , at h = 0.25, h = 0.5 = ∆t, and h = 0.8.The Euler-Net (which fails the convergence test) approximates a solution close to the true solution at a timestep in the training data but does poorly at the other timesteps; it does not capture the fluid flow behavior and gives a grainer solution.The RK4-Net (which passes the convergence test) is able to output solutions that have close correspondence to the true solution; it successfully interpolates the fine-scale flow fields that are in-between the training data snapshots, resulting in a much higher frame rate solution.

Figure 5 :
Figure 5: Non-linear pendulum: Extrapolation to predict initial condition trajectories on which the model was not trained.ODE-Net models are trained on randomly chosen initial conditions (different θ values), which have temporal points spaced apart by ∆t = 0.1 (see (b)).Each model is then evaluated on a different set of initial conditions (not in the training data) at much smaller step sizes (h = 0.001, 100× higher resolution).After evaluation, the Euler-Net quickly falls off of the phase plot lines corresponding to the true solution for the different test trajectories.The RK4-Net, which has learned a continuous dynamics, is able to extrapolate to different trajectories (starting at different initial conditions), with good correspondence to the phase plot lines of the true solution.
Training setup.An example distribution of irregularly sampled training data is shown in Fig.6(a).The baseline ∆t is 0.05, subject to jitter and frameskipping errors.An Euler-Net and an RK4-Net are both trained on these temporal data points, where the values of ∆t are input into the integration schemes.(That is, at every given measurement, the timestep jiter was also recorded to use in ∆t 2∆t 3∆t 4∆t 5∆t ∆t distribution Probability (a) ∆t of training data distribution : xn + hN (xn; θ) RK4-Net: xn + RK4[N (xn; θ)] (d) Convergence test

Figure 6 :
Figure 6: Learning continuous dynamics from irregularly spaced discrete points.(a)Training data distribution for a scientific problem where temporal data measurements are taken with some amount of imprecision (e.g., the measurement of interest is not always taken at the exact same ∆t) and/or measurements are skipped (e.g., we only have a measurement at ∆t and 3∆t because the measurement at 2∆t failed or was lost).(b) Different ODE-Nets are trained on the irregularly spaced data (indicated by green dots), as denoted by the green circles.Then, the trained models are evaluated at a very low h (where h ≪ ∆t).The RK4-Net (orange) is able to learn a continuous dynamics and follow the continuous solution over time (indicated by the dashed line), while the Euler-Net (blue) does not.(c) The RK4-Net (orange) is able to reconstruct the finescale, high resolution solution (indicated by dashed line), with good correspondence to the continuous solution, from the coarse, irregularly spaced training data (indicated by green dots).(d) The RK4-Net (orange) passes the convergence test, but the Euler-Net (blue) does not.

Figure 7 :
Figure7: SINDy method: additional convergence tests.Instead of training from dynamical data with an ODE-Net model, SINDy is used instead to learn a sparse polynomial basis.The accuracy of the time derivative approximation is altered by using a different order of central finite difference stencils.Here, FD-1 is the first-order two-point stencil, FD-2 is the second-order three-point stencil, and FD-4 is the fourth-order five-point stencil.In this example, FD-1 is analogous to Euler-Net, FD-2 is analogous to Midpoint-Net, and FD-4 is analogous to RK4-Net.In (a), the convergence test method is shown on the harmonic oscillator, where the convergence test is replicated by training via the SINDy method.In (b) and (c), the convergence test is applied when training on the non-linear pendulum via the SINDy method.We observe similar results as when training using an ODE-Net: the lower-order approximation (FD-1) overfits to the ∆t in the training data, failing the convergence test; and the higher-order approximation (FD-4) passes the convergence test.

Figure 8 :
Figure 8: Illustration of our convergence test for linear ODEs.The analytically derived global error estimate, Error(h), using the ideally learned linear ODE.Our results replicate the empirically performed convergence tests.(a) All methods achieve zero global error at h = ∆t, but converge to finite values as h → 0. (b) We introduce the ε error parameter.idpoint and RK4 achieve a finite error of comparable magnitudes at h = ∆t.Overfitting is not observed for RK4 with ε = 10 5 because the |λ − w| error is less significant than |ε|.
previous proof, we can again use the Taylor polynomial T q (wh) and Taylor remainder term R q (wh) to easily compare values with e λh : − e wh + e wh − T q (wh) + 3w 2 ε + 3wε 2 + ε 3 + ... Thus, we see have the original polynomial in w, with additional terms in ε.We truncate the remainder terms of higher powers of ∆t to simplify to obtain, p + pw p−1 ε + • • • + ε p .(80)We can group the terms of w i to obtain, i − ∆tε + O(∆t 2 wε + ε 2 ) .