  Forecasting high-dimensional dynamical systems is a fundamental challenge in
various fields, such as the geosciences and engineering. Neural Ordinary
Differential Equations (NODEs), which combine the power of neural networks and
numerical solvers, have emerged as a promising algorithm for forecasting
complex nonlinear dynamical systems. However, classical techniques used for
NODE training are ineffective for learning chaotic dynamical systems. In this
work, we propose a novel NODE-training approach that allows for robust learning
of chaotic dynamical systems. Our method addresses the challenges of
non-convexity and exploding gradients associated with underlying chaotic
dynamics. Training data trajectories from such systems are split into multiple,
non-overlapping time windows. In addition to the deviation from the training
data, the optimization loss term further penalizes the discontinuities of the
predicted trajectory between the time windows. The window size is selected
based on the fastest Lyapunov time scale of the system. Multi-step penalty(MP)
method is first demonstrated on Lorenz equation, to illustrate how it improves
the loss landscape and thereby accelerating the optimization convergence. MP
method can optimize chaotic systems in a manner similar to least-squares
shadowing with significantly lower computational costs. Our proposed algorithm,
denoted the Multistep Penalty NODE(MP-NODE), is applied to chaotic systems such
as the Kuramoto-Sivashinsky equation and the two-dimensional Kolmogorov flow.
It is observed that MP-NODE provide viable performance for such chaotic
systems, not only for short-term trajectory predictions but also for invariant
statistics that are hallmarks of the chaotic nature of these dynamics.
