Symbolic Neural ODEs: Learning interpretable models from time-series data

Nibodh Boddupalli; Jeff Moehlis

cs.LGAug 22, 2026

Symbolic Neural ODEs: Learning interpretable models from time-series data

Authors: Nibodh Boddupalli, Jeff Moehlis

Organizations: Dept. of Mechanical Engineering, University of California, Santa Barbara, CA 93106

Abstract

We present a machine learning framework for identifying sparse, interpretable models of dynamical systems directly from time-series data. Our approach parameterizes the underlying vector field using a neural architecture and trains it by minimizing a multi-step prediction loss over a finite horizon. To ensure numerical tractability, we optimize a mean absolute error objective averaged across prediction steps, and progressively increase the horizon during training. A key feature of this formulation is that it enforces consistency under repeated composition of the learned dynamics. As a result, the identified models exhibit significantly improved stability compared with approaches based on one-step regression of the vector field. When combined with sparsity-promoting regularization, this leads to parsimonious models that generalize beyond the training data. We demonstrate accurate recovery of systems exhibiting a wide range of behaviors, including stable and unstable fixed points, periodic orbits, and chaotic attractors. For chaotic systems, while long-term trajectory prediction is inherently limited by sensitivity to initial conditions, we show that multi-step training yields models with accurate short-term dynamics and strong agreement in long-time statistical properties, including mean, variance, and Lyapunov exponents. Moreover, we establish theoretical bounds linking trajectory error to statistical accuracy, providing a step toward a principled explanation for this behavior.

Figures & tables

Panel of Figure 1 — Figure 1: Architecture of the $l^{\text{th}}$ operational layer in the $k^{\text{th}}$ stack. This is a graphical representation of the primitive operations and functions defined by eq. 10a - eq. 10d in the main text.

Panel of Figure 2 — Figure 2: Architecture of our SymNN, showing $K$ stacks with $L$ layers within each stack.

Panel of Figure 3 — Figure 3: Schematic of the prediction of the future states $\bm{x}[h]$ (far right teal circle) from the initial states $\bm{x}[0]$ (far left teal circle) using the SymNN from fig. 2 . The blue arrows denote the forward computation of the predicted states $\tilde{\bm{x}}[h]$ up to $H$ times (purple circle). The closed loop denotes a recurrence relationship forward in time and the gradients of the loss (red arrow) are propagated backwards in time.

Panel of Figure 4 — Figure 4: Takens–Bogdanov normal form. (a) Phase portrait showing training data (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). The data are predominantly influenced by the unstable manifold of the saddle point. (b) Time-series comparison from the same initial condition. The close agreement demonstrates accurate recovery of both the vector field and transient dynamics despite both state and parametric sensitivity.

Panel of Figure 5 — Figure 5: Hopf normal form with $x=r\cos\theta$ , $y=r\sin\theta$ . (a) Phase portrait illustrating convergence to the limit cycle; training data (teal) are drawn primarily from post-transient dynamics. (b) Time evolution of $r$ and $\theta$ , showing exponential radial convergence and linear phase growth. Despite the non-uniqueness of circular limit-cycle dynamics, the identified model recovers the correct normal form structure.

Panel of Figure 6 — Figure 6: Phase portraits of model approximations for eq. 41a - eq. 41b (top) and eq. 44a - eq. 44b (bottom). Notice how the periodic orbits of the approximate models (dashed red) seem like a slightly perturbed version of the exact model (solid blue). Such models have similar stability properties, dynamic range, and time-periods as explained in section A.2 .

Panel of Figure 7 — Figure 7: Simple pendulum. (a) Phase portrait showing training data (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). The data span multiple energy levels, illustrating nonlinear periodic motion across trajectories of differing amplitudes. (b) Time-series comparison from the same initial condition. The close agreement demonstrates accurate recovery of the underlying Hamiltonian structure and preservation of periodic dynamics over the NODE horizon.

Panel of Figure 8 — Figure 8: Van der Pol oscillator with $\mu=15$ , demonstrating scale separation over length and time. (a) Phase portrait showing training data (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). The trajectories converge to a stable limit cycle, capturing both transient dynamics and nonlinear damping. (b) Time-series comparison from the same initial condition. The agreement illustrates accurate recovery of relaxation oscillations and the stable periodic orbit.

Panel of Figure 9 — Figure 9: Rössler system. (a) Phase portrait showing noise-free training points (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). (b) Time-series comparison, showing exact agreement between the exact and estimated model.

Panel of Figure 10 — Figure 10: Comparison of $m^{*}=2560+H$ data (green dots) from the Rössler attractor with $0-20\%$ NSR. All models are exact as evident from the overlapping true (blue) and estimated (red) trajectories.

Panel of Figure 11 — Figure 11: Phase portraits of models estimated for the Rössler attractor with (a) $25\%$ NSR (corresponding additive $\sigma_{abs}\approx[0.67,0.75,0.7]$ ) and (c) $30\%$ NSR (corresponding additive $\sigma_{abs}\approx[0.75,0.85,0.78]$ ) using horizons of $H=128$ and $256$ time-steps respectively. Corresponding time-series are shown in (b) and (d)

Panel of Figure 12 — Figure 12: RMSE versus training epochs for $H=$ 8 (red), 16 (blue), 32 (magenta), 64 (green), 128 (black) time-steps. For the same training dataset size of $m=2560$ in $\bm{X}_{0}$ ( $m^{*}=2560+H$ in $\bm{X}$ eq. 4 ), we can observe how the shorter horizons (8, 16) do not converge at all although they might initially have lower RMSE, while the longer horizons (32, 64, 128) not only converge to the true model but also converge earlier (epochs for convergence using $H=128$ are fewer than for $H=64$ , which in turn are fewer than for $H=32$ ).

Panel of Figure 13 — Figure 13: Chua double-scroll. (a) Phase portrait showing training data (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). (b) Time-series comparison from the same initial condition. There is good agreement for short times, and the subsequent divergence is consistent with chaotic dynamics.

Panel of Figure 14 — Figure 14: Chua system with $5\%$ NSR (corresponding to additive $\sigma_{abs}\approx[0.07,0.01,0.09]$ ). (a) Training data in phase portrait (teal dots) seem a lot more scattered due to noise, but trajectory generated by the identified model (dashed red) and the true trajectory (solid blue) are close. (b) Time-series comparison from the same initial condition. There is good agreement for short times, and the subsequent divergence is consistent with chaotic dynamics.

Panel of Figure 15 — Figure 15: Lorenz system without noise. (a) Phase portrait showing training data (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). (b) Time-series comparison, showing exact agreement between the exact and fit model.

Panel of Figure 16 — Figure 16: Lorenz $5\%$ NSR (corresponding to additive $\sigma_{abs}\approx[0.41,0.45,1.28]$ ). (a) Phase portrait showing training data (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). (b) Time-series comparison from the same initial condition. There is good agreement for short times, and the subsequent divergence is consistent with chaotic dynamics.

Panel of Figure 17 — Figure 17: Lorenz $10\%$ NSR (corresponding to additive $\sigma_{abs}\approx[0.79,0.92,2.46]$ ). (a) Phase portrait showing training data (teal), true trajectory (solid blue), and trajectory generated by the identified model (dashed red). (b) Time-series comparison from the same initial condition. There is good agreement for short times, and the subsequent divergence is consistent with chaotic dynamics.

Explore similar work

Turning Time Series into Algebraic Equations: Symbolic Machine Learning for Interpretable Modeling of Chaotic Time Series
Mar 7, 2026Madhurima Panja, Grace Younes, Tanujit ChakrabortyTime Series Forecasting Chaotic Time-Series Prediction
SAFIR, Sorbonne University Abu Dhabi, UAE · Sorbonne Center for Artificial Intelligence, Sorbonne University, Paris, France
Long-Time Trajectory Approximation via SA-NODEs: Model Predictive and Floquet Strategies
Aug 11, 2026Ziqian Li, Nikolaos M. MatzakosDynamical Systems Neural ODEs
Chair for Dynamics, Control, Machine Learning and Numerics, Department of Mathematics, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany · School of Pedagogical & Technological Education (ASPETE), Marousi Attikis, 151 22 Athens, Greece
Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics
Aug 17, 2026Conrad Ainslie, Pedram Hassanzadeh, Michael W. Mahoney +1Neural Surrogate Modeling Chaotic Time-Series Prediction
Department of Applied Mathematics, University of California, Santa Cruz, Santa Cruz, 95064, CA · Department of Geophysical Sciences, University of Chicago, Chicago, IL 60637 · Committee on Computational and Applied Mathematics, University of Chicago, Chicago, IL 60637 +3

Mar 7, 2026J/K move · Enter open · S save

Turning Time Series into Algebraic Equations: Symbolic Machine Learning for Interpretable Modeling of Chaotic Time Series

Madhurima Panja, Grace Younes, Tanujit Chakraborty

SAFIR, Sorbonne University Abu Dhabi, UAE · Sorbonne Center for Artificial Intelligence, Sorbonne University, Paris, France

Time Series Forecasting Chaotic Time-Series Prediction Explainable Time Series Forecasting Symbolic Regression

Open paper AskPDF arXiv

Chaotic time series are notoriously difficult to forecast. Small uncertainties in initial conditions amplify rapidly, while strong nonlinearities and regime dependent variability constrain predictability. Although modern deep learning often delivers strong short horizon accuracy, its black box nature limits scientific insight and practical trust in settings where understanding the underlying dynamics matters. To address this gap, we propose two complementary symbolic forecasters that learn explicit, interpretable algebraic equations from chaotic time series data. Symbolic Neural Forecaster (SyNF) adapts a neural network based equation learning architecture to the forecasting setting, enabling fully differentiable discovery of compact and interpretable algebraic relations. The Symbolic Tree Forecaster (SyTF) builds on evolutionary symbolic regression to search directly over equation structures under a principled accuracy complexity trade off. We evaluate both approaches in a rolling window nowcasting setting with one step ahead forecasting using several accuracy metrics and compare against a broad suite of baselines spanning classical statistical models, tree ensembles, and modern deep learning architectures. Numerical experiments cover a benchmark of 132 low dimensional chaotic attractors and two real world chaotic time series, namely weekly dengue incidence in San Juan and the Nino 3.4 sea surface temperature index. Across datasets, symbolic forecasters achieve competitive one step ahead accuracy while providing transparent equations that reveal salient aspects of the underlying dynamics.

Abstract

Explore similar work

Turning Time Series into Algebraic Equations: Symbolic Machine Learning for Interpretable Modeling of Chaotic Time Series

Long-Time Trajectory Approximation via SA-NODEs: Model Predictive and Floquet Strategies

Eigenanalysis framework for autoregressive neural emulators of multi-scale chaotic dynamics