# Stabilized Neural Ordinary Differential Equations for Long-Time Forecasting
of Dynamical Systems

Alec J. Linot<EMAIL_ADDRESS>Department of Chemical and Biological
Engineering, University of Wisconsin-Madison, Madison WI 53706, USA  Joshua W.
Burby  Los Alamos National Laboratory, Los Alamos, NM 87545, United States of
America  Qi Tang  Los Alamos National Laboratory, Los Alamos, NM 87545, United
States of America  Prasanna Balaprakash  Mathematics and Computer Science
Division, Argonne National Laboratory, Lemont, IL 60439, USA  Michael D.
Graham  Department of Chemical and Biological Engineering, University of
Wisconsin-Madison, Madison WI 53706, USA  Romit Maulik  Mathematics and
Computer Science Division, Argonne National Laboratory, Lemont, IL 60439, USA

###### Abstract

In data-driven modeling of spatiotemporal phenomena careful consideration
often needs to be made in capturing the dynamics of the high wavenumbers. This
problem becomes especially challenging when the system of interest exhibits
shocks or chaotic dynamics. We present a data-driven modeling method that
accurately captures shocks and chaotic dynamics by proposing a novel
architecture,  _stabilized neural ordinary differential equation_ (ODE). In
our proposed architecture, we learn the right-hand-side (RHS) of an ODE by
adding the outputs of two NN together where one learns a linear term and the
other a nonlinear term. Specifically, we implement this by training a sparse
linear convolutional NN to learn the linear term and a dense fully-connected
nonlinear NN to learn the nonlinear term. This is in contrast with the
standard neural ODE which involves training only a single NN for learning the
RHS. We apply this setup to the viscous Burgers equation, which exhibits
shocked behavior, and show better short-time tracking and prediction of the
energy spectrum at high wavenumbers than a standard neural ODE. We also find
that the stabilized neural ODE models are much more robust to noisy initial
conditions than the standard neural ODE approach. We also apply this method to
chaotic trajectories of the Kuramoto-Sivashinsky equation. In this case,
stabilized neural ODEs keep long-time trajectories on the attractor, and are
highly robust to noisy initial conditions, while standard neural ODEs fail at
achieving either of these results. We conclude by demonstrating how
stabilizing neural ODEs provide a natural extension for use in reduced-order
modeling by projecting the dynamics onto the eigenvectors of the learned
linear term.

##  I Introduction

With the ever increasing amount of data available from experiments or high-
resolution simulations it has become common to build models directly and
strictly from data. This tactic is powerful in that no a priori knowledge of
the system is needed and it forces the data to guide the solution.
Unfortunately, models built without any knowledge of the system often lack
robustness to perturbations in the data and do a poor job of generalizing
outside of the training data. In particular, we are interested in the problem
of modeling dynamical systems from time-series data. For example, this could
be predicting the future location of a pendulum given angle and velocity data,
predicting the weather, or predicting chemical concentrations in a reactor.
Time-series predictions are especially susceptible to issues of robustness
because states are repeatedly put through the model and error that pushes
states away from the training data become worse as the state moves further and
further away from the training data.

Here we seek to overcome the issue of robustness, in many data-driven models,
by using the fact that often in dissipative partial differential equations
(PDEs), stability comes from a linear operator that dissipates energy [6]. In
particular, we are concerned with the _nonlinear_ stability of attractors. By
nonlinear stability we mean finite perturbations to states on the attractor
are pushed back to the attractor. This type of stability can be show using the
energy method [32], with Lyapunov functions [18], or for dynamical systems
which possess an inertial manifold [6]. The energy method in particular
highlights the importance of the linear term in stabilizing many dissipative
PDEs. For example, when computing the energy balance for the Viscous Burgers
Equation (VBE) and the Navier-Stokes Equation dissipation comes from the
diffusive term [27] and for the Kuramoto-Sivashinky equation (KSE)
hyperdiffusivity dissipates energy [8]. Below we make this notion of stability
more precise.

In what follows, we briefly review some data-driven modeling approaches and
then outline our stabilizing neural ODE approach. In this approach we train a
linear NN and a nonlinear NN simultaneously to approximate the right hand side
of an ODE, which we find stabilizes the system to perturbations. None of the
existing approaches described explicitly learn a linear term in an ODE from
snapshots of data.

When  modeling time-series data of autonomous systems (i.e. the right hand
side does not explicitly depend on time) two approaches are to either
represent the dynamics as discrete time flows

$u(t+\tau)=F(u(t)),$ (1)

or as an ODE

$\dfrac{du}{dt}=f(u),$ (2)

where $u$ is the state of the system at a given time. Often, snapshots of
data, in the form $\\{u(t_{1}),\ldots,u(t_{N})\\}$, will be available for
training a model. This makes discrete time flows a natural first approach
because that is the form in which the data is available (i.e. the data is in
the precise form that a discrete time flows is capable of predicting).

When the dynamics are simple, a common approach to estimating $F$ is through a
linear map (i.e. $u(t+\tau)=Gu(t)$) using dynamic mode decomposition (DMD)
[19]. DMD works as a predictive tool for systems that decay to a fixed point
or exhibit quasiperiodic dynamics with discrete frequencies. When applying
DMD, or extended DMD (DMD with nonlinear observables), to nonlinear systems
challenges arise in the form of closure of models and robustness to noise
[37]. Lusch et al. [22] showed one way to extend these ideas for a few simple
nonlinear systems is to use an autoencoder to learn a change of basis under
which the dynamics are linear. Moreover, even with more complicated nonlinear
systems there always exists a linear _Koopman operator_ [3]. The Koopman
operator describes the evolution of arbitrary observables, but is infinite-
dimensional. Due to these issues, we do not seek a strictly linear time
evolution method. We instead consider methods that include nonlinearity.

Two classes of machine learning approaches that have seen success in modeling
nonlinear systems include reservoir networks [21, 26] and recurrent neural
networks (RNN) [13, 35, 34]. Both of these methods work by finding the optimal
parameters to some function that maps $u(t)$ forward to some
$\tilde{u}(t+\tau)$ by using a history of passed states. The parameters are
then trained by minimizing the difference between the prediction and the true
state $\left<||u(t+\tau)-\tilde{u}(t+\tau)||^{2}\right>$. Both methods can
predict multiple Lyapunov times for the chaotic dynamics of the KSE [34]. The
Lyapunov time is the relevant time scale for chaotic systems and will be
discussed in more detail below. Some drawbacks of these methods are they
typically result in an increase in the dimension of the problem, they
typically require evenly spaced data, and the start up typically needs
multiple known states because these methods are non-Markovian (e.g.
predictions of the next state require multiple past states).

Instead of learning a discrete-time representation it is also possible to
learn an ODE. In general, this is less straightforward because the data
available tends to be the state at different points in time. However, if the
time derivative is available Gonzalez Garcia et al. [12] showed that a NN with
states and derivative inputs could learn the ODE. A similar, but more
interpretable approach is “Sparse Identification of Nonlinear Dynamics”
(SINDy)[2]. This approach involves selecting a dictionary of candidate
functions for Eq. 2, and using sparse regression to identify the dominant
terms which most accurately reconstruct the time derivative. Alternatively,
the time derivative can be approximated with a multistep time-integration
scheme and a NN can be trained off of this approximation [29]. A drawback of
all these approaches is the data needs to be close in time to get good
approximations of the derivative, if the time derivative is not known.

The data used can be spaced further in time by using the “neural ODE” approach
developed by Chen et al. [5]. In this approach a trajectory is evolved forward
using a NN for $f$ with some numerical time integrator. Then the NN is trained
to minimize the difference between this prediction and the true state. In
Section II we outline this approach in more detail. The advantages this
approach has over the other methods include data can be spaced further in
time, trajectories can be evolved to arbitrary points in time, and adaptive
time stepping can be used.

Neural ODEs have been successfully used for short-time prediction of Burgers
equation [24], for the evolution of dissipation in decaying isotropic
turbulence [28], and for flow around a cylinder [30]. Additionally, in [20] we
showed that the long-time dynamics of the KSE could be captured using neural
ODEs. However, this only works when the dimension of the problem is reduced.
Without any dimension reduction high wavenumbers are amplified resulting in
trajectories leaving the attractor. In this work we seek to overcome this
problem of high-wavenumbers amplifying not by dimension reduction, but by
explicitly adding a linear term into the neural ODE approach. The properties
of this linear term are described below.

A similar idea to learn a linear term of an ODE is “Linear and Nonlinear
Disambiguation Optimization” (LANDO) [1]. In this approach Eq. 2 is
approximated with a dictionary-based kernel model. Then, the linear term on
the RHS of the ODE is approximated by calculating the Jacobian of $f$ around
the zero base state. There are two clear practical differences between that
approach and the one we present. First, LANDO requires either time derivative
data, or data closely spaced in time, whereas neural ODEs can use data widely
spaced in time. Second, the Jacobian around the base state is only accurate if
the function approximation is accurate around the base state. In our case, the
base state is the zero solution, and the data we use for training never
approaches the base state, so the Jacobian around this state may be poorly
estimated.

In what follows, we show that by adding an explicit linear term into the
neural ODE we can improve long-time predictive capabilities, we can improve
robustness to noise, and we can create a reduced order model from the basis
given by the linear term. Section II describes the framework we take for
constructing and training our models. Then, in Section III we present the
performance of the different models on the VBE (Section III.1) and on the KSE
(Section III.2). The VBE is presented to show how the model can handle shock
behavior and the KSE is presented to show how the model handles prediction of
chaotic dynamics on an attractor. A key result in both situations is for
standard neural ODEs, that lack an explicit linear term, high wavenumbers grow
in time, but by learning a linear term this erroneous behavior is avoidable
and the models become much more robust to noise. Finally, we conclude in
Section IV.

##  II Stabilized Neural ODE Framework

We first introduce our proposed  _stablized neural ODE_ framework, which
builds a dissipative linear mapping into the conventional neural ODE
framework. The approach is motivated by the stabilizing effect the linear term
has in many dynamical systems. We further discuss a natural and interesting
extension of the proposed approach into reduced order modeling (ROM).

For our analysis we consider snapshots of data
$\left\\{u(t_{1}),u(t_{2}),\ldots,u(t_{N})\right\\}$ that come from solutions
to an autonomous ODE with the form

$\dfrac{du}{dt}=f(u)=Au+F(u).$ (3)

Here $u(t_{i})\in\mathbb{R}^{d}$ is the state, $t\in\mathbb{R}$ is time,
$A\in\mathbb{R}^{d\times d}$ is a matrix, and
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and
$F:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ are nonlinear functions. In
particular we investigate systems where $A$ is symmetric and $F$ contains a
quadratic nonlinearity, but the ideas should hold for any system with a
dissipative linear term. For the two equations we consider, the VBE and the
KSE in one spatial dimension, we know, not only that the long-time solutions
lie on attractors, but that the long-time dynamics also collapse onto inertial
manifolds [36, 10]. An inertial manifold is a finite-dimensional invariant
attractive manifold that attracts trajectories exponentially fast [9].

With this data, we seek to train a NN to approximate the RHS of this ODE using
the neural ODE framework. Typically only $f$ is represented by a NN, but here
we show that a better approach is training two NNs _simultaneously_ , a linear
NN – in our case this will be a convolutional NN (CNN) – for $A$ and a
nonlinear NN for $F$. When we approximate one of these terms with a NN, or
make a prediction with these approximations, we will denote the approximation
with a tilde (e.g. $\tilde{f}$).

For our trials we constrain $A$ to be symmetric for the VBE, and apply no
constraints for the KSE. We chose to limit the constraints applied to the
neural ODE because the underlying systems possesses nonlinear stability in the
form of an inertial manifold. We took this approach to show that the model can
learn this property from the data when the neural ODE is provided a linear
term. An alternative future direction is explicitly enforcing the existence of
an inertial manifold. This could be done by enforcing a spectral gap
condition, which guarantees the existence of an inertial manifold [38]. The
spectral gap condition is

$\lambda_{N+1}-\lambda_{N}>2L$ (4)

where $\lambda_{i}$ is the $i^{\text{th}}$ eigenvalue of $A$, sorted in
increasing order, and $L$ is the Lipschitz constant of $F$.

To train these NNs, we numerically integrate one of the ODEs forward in time
to estimate the state $\tilde{u}(t_{i}+\tau)$

$\tilde{u}(t_{i}+\tau)=u(t_{i})+\int_{t_{i}}^{t_{i}+\tau}\tilde{h}(u(t);\theta)dt.$
(5)

In what follows we consider three variations on $\tilde{h}$: 1) a nonlinear
NN, $\tilde{h}(u)=\tilde{f}(u;\theta_{1})$, 2) a nonlinear NN with a fixed
linear term, $\tilde{h}(u)=Au+\tilde{F}(u;\theta_{2})$, and 3) a nonlinear NN
with a linear CNN,
$\tilde{h}(u)=\tilde{A}(\theta_{3})u+\tilde{F}(u;\theta_{2})$. Variations 2
and 3 are what we refer to as stabilized neural ODEs. In variation 2 we use
the true linear term of the underlying system, however the framework is
agnostic to this selection. For example, simply replacing $A$ with a scalar
damping term may be useful for model stability. In all of these cases,
$\theta_{i}$ refers to the parameters that are trained. The integration of
$\tilde{h}$ is then used to calculate the loss, which is the difference
between the predicted state, $\tilde{u}(t_{i}+\tau)$, and the known state,
$u(t_{i}+\tau)$

$J=\left<||u(t_{i}+\tau)-\tilde{u}(t_{i}+\tau)||_{1}\right>.$ (6)

We use the $L_{1}$-norm, but other norms are also applicable.

For learning the parameters of these NNs, the gradient of this loss with
respect to all the parameters, $\theta$, must be calculated. This gradient can
either be calculated by backpropagating through the solver with automatic
differentiation or by solving an adjoint problem backwards in time [5]. When
backpropagating all of the data must be stored at each time-step, which
becomes memory intensive when the prediction horizon becomes large. For our
trials we do not consider training data spaced apart far enough in time where
this becomes a problem, and empirically we found using the adjoint method took
longer, so we use backpropagation for training.

In addition to stabilizing the dynamics, which we show in the following
section, the addition of a linear term leads to an interesting and novel ROM
framework, i.e., the stablized neural ODE framework implies a natural basis
for ROM based on the eigenvectors of the linear term. This “flips” the
traditional ROM approach in that we first discovery dynamics and then use this
dynamical system with classical methods for ROM. In particular, we apply
standard Galerkin and nonlinear Galerkin approaches for ROM. These methods
require that the ODE in Eq. 3 is known, but because we learn this equation
from data we can directly apply these methods without knowing the equation a
priori.

Specifically, we perform an eigendecomposition on the linear term, which we
learn from a CNN, to find $\tilde{A}\mathrm{v}_{i}=\lambda_{i}\mathrm{v}_{i}$
and then we project the model approximation of Eq. 3 onto the leading
eigenvectors, giving the  _resolved_ dynamics

$\frac{\mathrm{d}p}{\mathrm{~{}d}t}=\tilde{A}p+P\tilde{F}(p+q),$ (7)

and the trailing eigenvectors, giving the  _unresolved_ dynamics

$\frac{\mathrm{d}q}{\mathrm{~{}d}t}=\tilde{A}q+Q\tilde{F}(p+q),$ (8)

Here $p=Pu$ is the projection onto the leading eigenvectors (resolved
dynamics) and $q=Qu$ is the projection onto the trailing eigenvectors
(unresolved dynamics). For the systems we consider, $\tilde{A}$ is symmetric,
which results in real eigenvalues and orthogonal eigenvectors. This means the
projections can be written $P=V_{p}V_{p}^{T}$ and $Q=V_{q}V_{q}^{T}$, where
$V_{p}=[\mathrm{v}_{0},\dots,\mathrm{v}_{p}]$ and
$V_{q}=[\mathrm{v}_{p+1},\dots,\mathrm{v}_{d}]$. The classic way to sort the
eigenvectors is in the order of decreasing real part of the eigenvalues [23].
In section III.2 we study the best approach to organizing the eigenvectors.

With Eq. 7 and Eq. 8 ROM can be performed by either assuming $q=0$, which is
the standard Galerkin approach, or by assuming $\mathrm{d}q/\mathrm{d}t=0$,
which is the nonlinear Galerkin approach [23, 15]. With the standard Galerkin
approach Eq. 7 can be solved forward in time without any additional equations.
With the nonlinear Galerkin approach $q$ is approximated iteratively by

$q_{i+1}=\tilde{A}^{-1}Q\tilde{F}(p+q_{i}),$ (9)

with an initial guess of $q_{0}=0$. In our trials we only use one iteration to
calculate $q_{1}$, because additional iterations had little affect on the
results. Nonlinear Galerkin is more accurate than Galerkin, but is
computationally more expensive due to solving for $q$ at every timestep.
Another alternative is postprocessing Galerkin where $p$ is solved forward in
time using standard Galerkin and $q$ is approximated for all $p$ at the end
with Eq. 9 [11].

a)b)c)d)e)f)g)h)i)

(a)

(b) (c)

(d) (e)

(f) (g)

(h) (i)

Figure 1: Snapshots of VBE trajectories of true system and model predictions
for different initial conditions. The rows indicate different initial
conditions and the columns indicate solution fields at $t=1$, $3$, and $5$,
respectively.

##  III Results

In the following sections we demonstrate the importance of constructing neural
ODEs with a linear and a nonlinear term on the VBE and the KSE. We examine the
VBE to highlight the ability of this method to accurately capture the shock
behavior seen at low viscosity and to accurately capture the energy spectrum.
Then, we consider the KSE to show how our approach is needed to keep the long-
time trajectories of a chaotic dynamical system on the attractor. This then
leads to a demonstration of reduced order modeling via Galerkin and nonlinear
Galerkin using the eigenvectors of the linear term we discover. In all cases
results are shown for test data that was not used for training any of the
models. Neural ODEs are trained using PyTorch [25] with code modified from
Chen et al. [5].

###  III.1 Viscous Burgers Equation

In this section the data comes from numerical solutions of the 1D viscous
Burgers equation

$\dfrac{\partial v}{\partial t}=-v\dfrac{\partial v}{\partial
x}+\nu\dfrac{\partial^{2}v}{\partial x^{2}},$ (10)

in a domain of unit length ($L=1$) with periodic boundary conditions, with
viscosity $\nu=8\cdot 10^{-4}$. This parameter set results in what is known as
Burgers turbulence or _Burgulence_ [24]. We solve this equation using a
pseudospectral Runge-Kutta Crank-Nicolson scheme as described in [31]. The
state $u$ is then represented in $\mathbb{R}^{512}$ by $v$ sampled at
equidistant points in the domain. We construct 1000 random initial conditions
for training the models and 100 random initial conditions for testing model
performance. We then evolve each initial condition over the interval $[0,5]$
sampling every $0.05$ time units. This results in $10^{4}$ total snapshots of
data used for training. All figures compare test data to model reconstruction
of the test data.

When generating training and test datasets, the space of feasible initial
conditions is infinite, so we must first limit this region. Here, we chose to
limit the random initial conditions such that the ensemble-averaged energy
spectrum

$E(k,t)=\left<\dfrac{1}{2}|\hat{u}(k,t)|^{2}\right>$ (11)

at $t=0$ satisfies $\int E(k,0)\,dk=0.5L/(2\pi)$. In Eq. 11 $k$ is the
wavenumber and the Fourier transform is $\hat{u}=\mathcal{F}(u)$. This energy
condition is satisfied by selecting the initial conditions according to

$\hat{u}(k,0)=\sqrt{2E_{0}(k)}\left(\cos(2\pi\Psi(k))-\sin(2\pi\Psi(k)\right),$
(12)

with $E_{0}(k)=Ak^{4}\exp\left(-\left(k/k_{0}\right)^{2}\right)$, and
$\Psi(k)$ coming from a uniform distribution between 0 and 1 for each
wavenumber ($\Psi(k)=\mathcal{U}(0,1)$). The constants $A$ and $k_{0}$ are the
same as in [24]. Selecting the initial conditions in this manner results in
training and testing datasets with a variety of different initial energies.

Table 1: Architectures of NNs and matrices used in Sections III.1. “Shape”
indicates the dimension of each layer, “Activation” the corresponding
activation functions, “Learning Rate” is the learning rate the NN switches to
evenly over the training period, and “Weight Init” is the distribution used to
initialize the weights of the NNs. “relu” is the rectified linear unit [14].
$\tilde{F}_{1}$ is the NN trained with the true $A$ and $\tilde{F}_{2}$ is the
NN trained with $\tilde{A}$.

Function | Shape | Activation | Learning Rate | Weight Init
---|---|---|---|---
$\tilde{f}$ | 512/200/200/200/512 | relu/relu/relu/linear | $[10^{-3},10^{-4},10^{-5}]$ | $\mathcal{N}(0,10^{-2})$
$\tilde{F}_{1}$ | 512/200/200/200/512 | relu/relu/relu/linear | $[10^{-3},10^{-4},10^{-5}]$ | $\mathcal{N}(0,10^{-2})$
$\tilde{F}_{2}$ | 512/200/200/200/512 | relu/relu/relu/linear | $[10^{-3},10^{-4}]$ | $\mathcal{N}(0,10^{-2})$
$\tilde{A}$ | 512/512 | linear | $[10^{0},10^{-1},10^{-2}]$ | $\mathcal{N}(0,10^{4})$

With this data, we train neural ODEs for the three cases listed above: 1)
$\tilde{h}(u)=\tilde{f}(u;\theta_{1})$, 2)
$\tilde{h}(u)=Au+\tilde{F}(u;\theta_{2})$, and 3)
$\tilde{h}(u)=\tilde{A}(\theta_{3})u+\tilde{F}(u;\theta_{2})$, which we will
denote as nonlinear, fixed linear, and CNN. We train these models for $10^{4}$
epochs at which point the error in the loss (Eq. 6) stops decreasing. The
architectures and training parameters for each of these NNs are presented in
Table LABEL:Table. Hyperparameter tuning was performed manually by training
many NNs and selecting the NNs with the best short-time tracking performance.
The linear term in particular showed sensitivity to hyperparemeter selection,
which should be investigated further in future work.

We construct $\tilde{f}$ and $\tilde{F}$ from fully connected NNs, and
$\tilde{A}$ from a linear CNN with stride one, one filter, and a filter width
of 3. This is equivalent to learning a sparse tridiagonal matrix. In the case
of the VBE, we enforce symmetry in $\tilde{A}$ by training a CNN
$B(\theta_{3})$ and letting $\tilde{A}=B(\theta_{3})+B(\theta_{3})^{T}$. With
this tridiagonal structure, the optimal filter for approximating the
derivative is $[210,-420,210]$ (this is the approximation of the second
derivative with central differencing), and the filter learned after training
is $[45,-90,45]$. Although the scale of the filter is off, we will see the
results with this filter and the true filter match well.

In Fig. 1 we show the performance of the three neural ODE approaches on
different initial conditions of the VBE. For all three initial conditions we
see that stabilized neural ODE approaches (fixed linear and CNN) result in
good agreement between the true solution and the model prediction. When,
instead we take the nonlinear approach high-wavenumber behavior appears in the
solution and persists at long-times resulting in a poor prediction.

Figure 2: Ensemble averaged error of model predictions VBE.

The behavior just described becomes even more pronounced when considering the
relative error averaged over initial conditions, shown in Fig. 2. Here, all
methods match at very short times. The fixed linear approach performs the best
followed by the CNN approach, and then by the nonlinear approach. For the two
stabilized neural ODE methods the errors level off, or slightly go down, at
long-times indicating that trajectories tend to separate somewhat and then
decay in a similar fashion causing them to stay the same relative distance
apart. However, this is not the case for the nonlinear method where the high-
wavenumber behavior at long-times continues to push the trajectory further
away increasing the relative error.

a)b)c)d)e)f)

(a)

(b) (c)

(d) (e)

(f)

Figure 3: Ensemble-averaged energy spectrum of the true system and the model
predictions at different times , with $k^{-2}$ scaling for reference. (a)-(f)
are between $t=0-5$.

We also consider the ensemble-averaged energy spectrum defined in Eq. 11. Fig.
3 shows the energy spectrum at various times. At early times, when the shock
develops, there is an increase in the tail of the energy spectrum, and we see
the development of a $k^{-2}$ scaling. At $t=1$ all methods capture most of
the spectrum with errors appearing at high-wavenumbers. Then, at $t\gtrsim 2$,
both the stabilized methods match the true energy spectrum over full range of
wavenumbers, while the high-wavenumber behavior of the nonlinear model
deviates drastically as time increases. Capturing the dynamics properly at all
scales requires a linear term to damp out high wavenumbers.

a)b)c)

(a)

(b) (c)

Figure 4: Ensemble averaged error of predictions with noisy initial conditions
of $\epsilon=0.1$, $0.3$, and $0.5$ for (a)-(c).

The explicit linear term in our approach also contributes to the robustness
and generalizability of the neural ODE setup. Figure 4 shows the relative
error in predictions for various levels of noise added to the initial
condition. To create this figure, we add random Gaussian noise of the form
$\mathcal{N}(0,\epsilon^{2})$ to each grid location for 100 different initial
conditions. Then those initial conditions are evolved forward with the true
equations and the learned models for comparison. In all these cases the error
in the nonlinear cases rapidly increases, while the stabilized neural ODE
cases exhibit similar behavior as with no noise. This shows that these models
generalize well to noisy initial conditions that are much different from the
training data. This can also be seen in Fig. 4(a), where we plot the ensemble
averaged energy spectrum at different times for the different levels of noise.
The strictly nonlinear model does a poor job of capturing the spectrum at any
wavenumber, while the stabilized neural ODEs accurately capture the spectrum
like they did for the initial conditions with no noise.

a)b)c)d)e)f)g)h)i)

(a)

(b) (c)

(d) (e)

(f) (g)

(h) (i)

Figure 5: Ensemble-averaged energy spectrum of the true system and the model
predictions at different times and noise on the initial conditions.
(a)-(c),(d)-(f), and (g)-(i) are with noise $\epsilon=0.1$, $0.3$, and $0.5$,
at times $t=1$, $3$, and $5$, respectively.

###  III.2 Kuramoto-Sivashinsky Equation

In the previous section we showed that our setup captures the shock dynamics
and transient decay exhibited by the VBE. In that case, all solutions decay to
zero at long times. Now we examine the KSE, where the dynamics are chaotic and
collapse onto a low-dimensional inertial manifold at long-times [10, 33, 16,
38]. A useful model for the KSE should keep long-time trajectories on the
attractor and provide a means for reduced-order modeling to capitalize on the
low-dimensionality of the manifold on which the data lies.

We numerically solve the 1D KSE,

$\dfrac{\partial v}{\partial t}=-v\dfrac{\partial v}{\partial
x}-\dfrac{\partial^{2}v}{\partial x^{2}}-\dfrac{\partial^{4}v}{\partial
x^{4}},$ (13)

with periodic boundary conditions on a domain of size $L=22$. Solutions are
found by performing a Galerkin projection onto Fourier modes and using
exponential time differencing to evolve the ODE forward in time [17]. After
solving in Fourier space, we then transform the data back to physical space
for training the neural ODEs. Code for solving the KSE is available in
Cvitanović et al. [7]. The data used is a single trajectory that has collapsed
onto the inertial manifold evolved over $t=[0,10^{5}]$ sampled every $0.25$
time units. The first 80% of this data is used for training and the remaining
20% is used as test data.

Table 2: Neural ODE details for Section III.2. Labels are the same as in Table
LABEL:Table. “Sig” is the sigmoid activation [14].

Function | Shape | Activation | Learning Rate | Weight Init
---|---|---|---|---
$\tilde{f}$ | 64/200/200/200/64 | sig/sig/sig/linear | $[10^{-3},10^{-4}]$ | $\mathcal{N}(0,10^{-2})$
$\tilde{F}_{1}$ | 64/200/200/200/64 | sig/sig/sig/linear | $[10^{-3},10^{-4}]$ | $\mathcal{N}(0,10^{-2})$
$\tilde{F}_{2}$ | 64/200/200/200/64 | sig/sig/sig/linear | $[10^{-3},10^{-4}]$ | $\mathcal{N}(0,10^{-2})$
$\tilde{A}$ | 64/64 | linear | $[10^{0},10^{-1},10^{-2}]$ | $\mathcal{U}(-\sqrt{1/3},\sqrt{1/3})$

We train neural ODEs for $4\cdot 10^{4}$ epochs, and consider the same three
cases as before. Table LABEL:Table2 contains details on the architectures and
training parameters used for each of these NNs. Small modifications were made
in the architectures due to the different dynamics of the KSE. As before,
$\tilde{f}$ and $\tilde{F}$ are fully connected NN and $\tilde{A}$ is a linear
CNN. For $\tilde{A}$ we selected a wider filter of width 5 because higher-
order gradients need to be approximated. In general the linear term is
unknown, so the filter width can be viewed as an additional tuning parameter
that is similar to selecting the stencil size for a finite difference
approximation. For our system, and filter width, the optimal filter parameters
are $[-72,278,-413,278,-72]$ (this is central differencing of the diffusive
and hyperdiffusive terms), while the parameters learned for the CNN were
$[-0.7,1.8,-2.7,1.8,-0.7]$. Here we also ran trials where we forced symmetry
in the filter, but we found the training procedure resulted in the symmetry
without the need to enforce it. Despite the magnitude of the parameters being
much smaller, we will show the CNN successfully damps out high wavenumbers,
and gives the same eigenvectors.

white a)white b)white c)white d)

(a)

(b) (c)

(d)

Figure 6: Predictions of trajectories using different models. (a) is the true
trajectory, (b) is the predicted trajectory with the nonlinear model. (c) is
the predicted trajectory with the fixed linear model. (d) is the predicted
trajectory with the CNN model.

Figure 6 compares a true trajectory to predictions from the three methods. The
left side shows short-times, where model predictions should match the true
trajectory. The right side shows long-times, where the model predictions
should have statistics that agree with the true trajectory. In Fig. 6(b) the
nonlinear method prediction is shown. At early times the predicted trajectory
matches the true trajectory, but at long-times erroneous high-wavenumber
behavior appears and the trajectory has moved far away from the attractor.
This result agrees with what was found in [20]. Predictions from the
stabilized neural ODEs appear in Figs. 6(c) and 6(d) for the fixed linear and
CNN methods, respectively. In the linear cases, there is good agreement at
short times, with the fixed linear method tracking slightly better, and at
long-times the trajectories stay on the attractor. Below we show statistics to
further validate this conclusion.

Figure 7: Ensemble averaged error of model predictions for the KSE.

To better quantify the apparent results in Fig. 6 we consider some ensemble
averaged quantities. In Fig. 7 we show the ensemble averaged error for 100
different initial conditions as a function of Lyapunov time normalized by the
difference between states on the attractor sampled at random times $t_{i}$ and
$t_{j}$ ($D=\left<||u(t_{i})-u(t_{j})||^{2}\right>$). The Lyapunov time is the
inverse of the leading Lyapunov exponent and is the relevant timescale for
chaotic systems. For the KSE, with $L=22$, the Lyapunov time is
$\tau_{L}\approx 22$. The error in all three models nearly match over this
time range, with the fixed linear model performing slightly better. All three
models are in agreement because the high-wavenumber behavior of the strictly
nonlinear model does not dominate over these short time ranges. In all three
cases there is good tracking, with trajectories completely diverging at around
$t\sim 4\tau_{L}$.

white a)white b)white c)white d)

(a)

(b) (c)

(d)

Figure 8: (a) is the true joint PDF at $L=22$. (b)-(d) are the joint PDFs for
the nonlinear, fixed linear, and CNN models, respectively.

The advantage stabilized neural ODEs have over the standard neural ODE, for
the KSE, is the prediction at long times. In Fig. 8 we show the true joint PDF
of the first ($u_{x}$) and second ($u_{xx}$) derivatives of the state plotted
for every spatiotemporal point of a trajectory evolved forward for $10^{4}$
time units. Figure 8(b) therein, which corresponds to the nonlinear ODE
prediction, shows the probability of finding a point in the nonlinear
simulations prediction with the correct first and second derivatives is nearly
zero. Both stabilized neural ODEs, shown in Figs. 8(c) and 8(d), accurately
match the true joint PDF.

Similar to the VBE, the presence of a linear term makes the KSE models highly
robust to noisy initial conditions. Figure 9 shows the short-time evolution of
predictions from the CNN approach with varying levels of noise. Here the
initial condition is perturbed with Gaussian noise of magnitude
$\mathcal{N}(0,\epsilon^{2})$ added to Fourier modes of $20^{\text{th}}$ up to
$31^{\text{st}}$. Even with the extreme levels of noise shown in 9(d) the
solution quickly decays back to the attractor. Despite the magnitude of the
CNN linear term being much smaller than the true linear term, the CNN still
rapidly damps high-wavenumbers.

The other useful aspect of the stabilized neural ODE architecture is it lends
itself naturally to reduced-order modeling (ROM). As mentioned in section II,
a standard approach to ROM involves projecting the true ODE (Eq. 3) onto a
basis given by the leading $d_{p}$ eigenvectors of $A$. This is the resolved
dynamics in Eq. 7. Then, the unresolved dynamics, Eq. 8 can be approximated by
setting $q=0$, which is the standard Galerkin approach, or by using Eq. 9 to
solve for $q$, which is the nonlinear Galerkin approach. Here we show we can
successfully apply both these ROM methods to our CNN neural ODE model that we
find strictly from data, whereas typically Eq. 3 must be known. We first
compare the projection operator of the CNN neural ODE to the true projection
operator, and then we compare the long-time statistical performance of the ROM
when varying dimension $d_{p}$.

For these ROM approaches the eigenvector basis is typically selected from
largest to smallest eigenvalues. However, selecting the eigenvectors in this
fashion does not necessarily result in the best agreement with the underlying
assumption of this ROM – that $\mathrm{d}q/\mathrm{d}t=0$. To investigate this
assumption, we consider the dynamics from projecting the equation onto each
eigenvector $d\mathrm{v}_{i}^{T}u/dt=\dot{\mathrm{v}}_{i}$. The ROM assumption
is reasonable when, for some $p$, $\dot{\mathrm{v}}_{i}\approx 0$ for all
$i=p+1,\dots,d$.

white a)white b)white c)white d)

(a)

(b) (c)

(d)

Figure 9: (a)-(d) are trajectories of the CNN model for initial conditions
with $\epsilon=0$, $0.1$, $1$, and $10$, respectively.

In Fig. 10(a) we show the variance in $\dot{\mathrm{v}}_{i}$ sorted from
largest to smallest eigenvalue $\lambda_{i}$. Each eigenvector can be
considered separately because the linear term is symmetric in both the true
system and from the CNN approach resulting in an orthogonal basis. Figure
10(a) shows sorting by eigenvalue does not lead to a monotonic decrease in
$\dot{\mathrm{v}}_{i}$. In Fig. 10(d) we apply this same sorting with the
linear term learned by the CNN. Here the sorting is worse, with the leading
terms showing little variance. In order for the CNN ROM to perform well the
projection matrix must match the projection matrix for the true linear term.
In Fig. 10(b) we show a projection matrix for the true linear term and in Fig.
10(e) we show a projection matrix for the CNN approach. Both of which were
generated using the eigenvalue sorting and the first 10 eigenvectors. These
projection matrices are substantially different, and the CNN projection
matrix, in particular, is not useful for ROM.

a)b)c)d)e)f)

(a)

(b) (c)

(d) (e)

(f)

Figure 10: (a)-(c) are the variance of the ODE projected onto each
eigenvector, the unsorted projection operator, and the sorted projection
operator, for the true linear term. (d)-(f) are the same plots for the linear
term learned by the CNN. These projection operators are for 10 eigenvectors.

In general, when only data is available there is no immediate way to compute
$\dot{\mathrm{v}}_{i}$, however, because we learn the ODE directly from data
computing $\dot{\mathrm{v}}_{i}$ after the fact is trivial. This motivates
another reason why our stabilized neural ODE approach naturally lends itself
to ROM. Instead of organizing according to eigenvalues, we can, instead, sort
according to the variance of $\dot{\mathrm{v}}_{i}$, which more accurately
reflects the assumption taken in nonlinear Galerkin. In Fig. 10(c) we show the
projection matrix from the true linear term with this sorting and in Fig.
10(f) we show the projection matrix from the CNN linear term. The two
projections matrices match nearly exactly.

With the projection selected according to this variance criterion, we show the
standard Galerkin and nonlinear Galerkin approaches performance in Fig. 11.
For this comparison we consider both the CNN neural ODE, sorted according to
the variance criterion, and the true ODE, sorted by the standard eigenvalue
criterion. We compare ROM performance with the Kullback-Leibler (KL)
divergence of the $u_{x}$ vs $u_{xx}$ joint PDF

$D_{KL}(\tilde{P}||P)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\tilde{P}(u_{x},u_{xx})\text{ln}\dfrac{\tilde{P}(u_{x},u_{xx})}{P(u_{x},u_{xx})}\,du_{x}du_{xx}.$
(14)

In this equation $P$ is the true joint PDF and $\tilde{P}$ is the ROM joint
PDF. When $P=0$ or $\tilde{P}=0$ we adopt the convention of setting the term
in the integral to 0 as was done in [4].

In Fig. 11 we vary the dimension $d_{p}$ (i.e. the number of modes on which we
project) and compute the KL divergence between the true joint PDF and the ROM
joint PDFs, generated the same way as in Fig. 8. In the case of the Galerkin
ROM, we see that the CNN approach does a better job than even the true
Galerkin approach at reconstructing the joint PDF. In the case of nonlinear
Galerkin, the CNN approach shows more error than the true nonlinear Galerkin
approach, but performs better than the standard Galerkin approaches. For all
ROMs we see excellent agreement around a dimension of $d_{p}\sim$18, as can be
seen in the inset, above which there is little change in the KL divergence and
the ROM agree well with the full system. These results show that by using
classical ROM approaches on our data-driven stabilized neural ODEs we can
achieve excellent model performance at far fewer dimensions than the
underlying system.

Figure 11: KL divergence of the joint PDFs generated with either standard
Galerkin “G” or nonlinear Galerkin “NG” when using the CNN neural ODE “CNN” or
the true equations “True”. The inset is the CNN Galerkin ROM at a dimension
$d_{p}=18$.

##  IV Conclusion

Stabilizing neural ODEs with a linear and a nonlinear term is an effective way
to improve the predictive capabilities of these models. By implementing this
structure we find the models are more robust to noise, perform better in
capturing statistics, and provide a method for selecting a natural basis,
different from conventional approaches, for reduced-order models. We test this
setup first on the VBE, where the addition of the linear term results in a
model that effectively captures the shock dynamics and the energy spectrum,
along with generalizing to noisy initial conditions. We then test it on the
KSE, where we find the linear term is necessary to keep trajectories on the
attractor, and again the linear term provides excellent robustness to noise.
With the KSE we further show our stabilized neural ODE method lends itself
naturally to ROM, and that, with proper sorting of the eigenvectors of the
linear term, we can accurately recreate the joint PDF of $u_{x}$ and $u_{xx}$
with far fewer dimensions than required for the full simulation.

The future directions of this research include apply the method to more
complex problems, where the existence of an inertial manifold may be unknown,
and investigating the importance of adding explicity constraints on the NNs
for stability. In particular, we are immediately interested in applying this
method to canonical fluid flow problems (e.g. Couette flow). This modeling
could either be directly data-driven, as in this paper, or for closure of low
resolution physics-based models. Regardless, with problems of increasing
complexity it may require the use of adding additional constraints for
stability. This can be done by enforcing that the linear term has a specific
eigenvalue spectrum, or by tuning the Lipschitz constant of the nonlinearity.
For example, we could enforce the spectral gap condition, which would
guarantee the existence of an inertial manifold.

###### Acknowledgements.

This work was supported by the U.S. Department of Energy (DOE), Office of
Science, Office of Advanced Scientific Computing Research (ASCR), under
Contract No. DE-AC02–06CH11357, at Argonne National Laboratory, and by the
Office of Fusion Energy Sciences and ASCR under the Scientific Discovery
through Advanced Computing (SciDAC) project of Tokamak Disruption Simulation
at Los Alamos National Laboratory. We acknowledge funding support from ASCR
for DOE-FOA-2493 “Data-intensive scientific machine learning”. This research
used resources of the Argonne Leadership Computing Facility, which is a DOE
Office of Science User Facility supported under Contract No. DE-
AC02–06CH11357.

## References

  * [1] P. J. Baddoo, B. Herrmann, B. J. McKeon, and S. L. Brunton.  Kernel learning for robust dynamic mode decomposition: Linear and nonlinear disambiguation optimization (lando), 2021.
  * [2] S. L. Brunton, J. L. Proctor, and J. N. Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.
  * [3] M. Budišić, R. Mohr, and I. Mezić.  Applied koopmanism.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(4):047510, 2012.
  * [4] F. Cazáis and A. Lhéritier.  Beyond two-sample-tests: Localizing data discrepancies in high-dimensional spaces.  In 2015 IEEE International Conference on Data Science and Advanced Analytics (DSAA), pages 1–10, 2015.
  * [5] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud.  Neural ordinary differential equations.  arXiv preprint arXiv:1806.07366, 2019.
  * [6] P. Constantin, C. Foias, B. Nicolaenko, and R. Témam.  Spectral barriers and inertial manifolds for dissipative partial differential equations.  Journal of Dynamics and Differential Equations, 1(1):45–73, 1989\.
  * [7] P. Cvitanović, R. Artuso, R. Mainieri, G. Tanner, and G. Vattay.  Chaos: Classical and Quantum.  Niels Bohr Inst., Copenhagen, 2016.
  * [8] P. Cvitanović, D. Borrero-Echeverry, K. M. Carroll, B. Robbins, and E. Siminos.  Cartography of high-dimensional flows: A visual guide to sections and slices.  Chaos, 22(4), 2012.
  * [9] C. Foias, M. S. Jolly, I. G. Kevrekidis, G. R. Sell, and E. S. Titi.  On the computation of inertial manifolds.  Physics Letters A, 131(7-8):433–436, 1988.
  * [10] C. Foias, B. Nicolaenko, G. R. Sell, and R. Temam.  Inertial manifold for the Kuramoto-Sivashinsky equation and an estimate of their lowest dimension.  J. Math. Pure Appl., 67:197–226, 1988.
  * [11] B. García-Archilla, J. Novo, and E. S. Titi.  Postprocessing the galerkin method: a novel approach to approximate inertial manifolds.  SIAM Journal on Numerical Analysis, 35(3):941–972, 1998.
  * [12] R. Gonzalez-Garcia, R. Rico-Martinez, and I. G. Kevrekidis.  Identification of distributed parameter systems: A neural net based approach.  Computers & Chemical EngineeringComputers & Chemical Engineering, 22:S965–S968, 1998.
  * [13] S. Hochreiter and J. Schmidhuber.  Long Short-Term Memory.  Neural Computation, 9(8):1735–1780, 1997.
  * [14] A. C. Ian Goodfellow, Yoshua Bengio.  The Deep Learning Book, volume 521.  2017\.
  * [15] F. Jauberteau, C. Rosier, and R. Temam.  The nonlinear Galerkin method in computational fluid dynamics.  Applied Numerical Mathematics, 6(5):361–370, 1990.
  * [16] M. S. Jolly, R. Rosa, and R. Temam.  Evaluating the dimension of an inertial manifold for the Kuramoto-Sivashinsky equation.  Advances in Differential Equations, 5(1-3):31–66, 2000.
  * [17] A.-K. Kassam and L. N. Trefethen.  Fourth-order time-stepping for stiff PDEs.  SIAM Journal on Scientific Computing, 26(4):1214–1233, 2005.
  * [18] H. Khalil.  Nonlinear Systems.  Pearson Education. Prentice Hall, 2002.
  * [19] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor.  Dynamic Mode Decomposition.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016\.
  * [20] A. J. Linot and M. D. Graham.  Data-driven reduced-order modeling of spatiotemporal chaos with neural ordinary differential equations.  arXiv preprint arXiv:2109.00060, 2021.
  * [21] M. Lukosevicius and H. Jaeger.  Reservoir Computing Approaches to Recurrent Neural Network Training.  2009\.
  * [22] B. Lusch, J. N. Kutz, and S. L. Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature Communications, (2018), 2018.
  * [23] M. Marion and R. Temam.  Nonlinear Galerkin methods: The finite elements case.  Numerische Mathematik, 57(1):205–226, 1989.
  * [24] R. Maulik, A. Mohan, B. Lusch, S. Madireddy, P. Balaprakash, and D. Livescu.  Time-series learning of latent-space dynamics for reduced-order model closure.  Physica D: Nonlinear Phenomena, 405:132368, 2020.
  * [25] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [26] 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.  Physical Review Letters, 120(2):24102, 2018.
  * [27] S. Pope.  Turbulent Flows, 2000.
  * [28] G. D. Portwood, P. P. Mitra, M. D. Ribeiro, T. M. Nguyen, B. T. Nadiga, J. A. Saenz, M. Chertkov, A. Garg, A. Anandkumar, A. Dengel, R. Baraniuk, and D. P. Schmidt.  Turbulence forecasting via neural ODE.  arXiv preprint arXiv:1911.05180, 2019.
  * [29] M. Raissi, P. Perdikaris, and G. E. Karniadakis.  Multistep Neural Networks for Data-driven Discovery of Nonlinear Dynamical Systems.  arXiv preprint arXiv:1801.01236, pages 1–19, 2018.
  * [30] C. J. G. Rojas, A. Dengel, and M. D. Ribeiro.  Reduced-order model for fluid flows via neural ordinary differential equations.  arXiv preprint arXiv:2102.02248, 2021.
  * [31] O. San and A. E. Staples.  Stationary two-dimensional turbulence statistics using a markovian forcing scheme.  Computers & Fluids, 71:1–18, 2013.
  * [32] B. Straughan.  The Energy Method, Stability, and Nonlinear Convection.  Applied mathematical sciences. Springer-Verlag, 1992.
  * [33] R. Temam and X. Wang.  Estimates on the lowest dimension of inertial manifolds for the Kuramoto-Sivasbinsky equation in the general case.  Differential and Integral Equations, 7(3-4):1095–1108, 1994.
  * [34] P. Vlachas, J. Pathak, B. Hunt, T. Sapsis, M. Girvan, E. Ott, and P. Koumoutsakos.  Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics.  Neural Networks, 126:191–217, 2020.
  * [35] 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.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2213):20170844, 2018-05.
  * [36] J. Vukadinovic.  Global dissipativity and inertial manifolds for diffusive burgers equations with low-wavenumber instability.  Discrete & Continuous Dynamical Systems, 29(1):327–341, 2011.
  * [37] Z. Wu, S. L. Brunton, and S. Revzen.  Challenges in dynamic mode decomposition.  Journal of The Royal Society Interface, 18(185):20210686, 2021.
  * [38] S. Zelik.  Inertial manifolds and finite-dimensional reduction for dissipative PDEs.  Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 144(6):1245–1327, 2013.