Artifacts of Numerical Integration in Learning Dynamical Systems

Bing-Ze Lu; Richard Tsai

math.NAJul 19, 2025

Artifacts of Numerical Integration in Learning Dynamical Systems

Authors: Bing-Ze Lu, Richard Tsai

Organizations: National Chung Cheng University, Taiwan · The University of Texas at Austin, USA

Abstract

In many applications, one needs to learn a dynamical system from its solutions sampled at a finite number of time points. The learning problem is often formulated as an optimization problem over a chosen function class. However, in the optimization procedure, prediction data from generic dynamics requires a numerical integrator to assess the mismatch with the observed data. This paper reveals potentially serious effects of a chosen numerical scheme on the learning outcome. Specifically, the analysis demonstrates that a damped oscillatory system may be incorrectly identified as having "anti-damping" and exhibiting a reversed oscillation direction, even though it adequately fits the given data points. This paper shows that the stability region of the selected integrator will distort the nature of the learned dynamics. Crucially, reducing the step size or raising the order of an explicit integrator does not, in general, remedy this artifact, because higher-order explicit methods have stability regions that extend further into the right half complex plane. Furthermore, it is shown that the implicit midpoint method can preserve either conservative or dissipative properties from discrete data, offering a principled integrator choice even when the only prior knowledge is that the system is autonomous.

Figures & tables

Figure 1: The gray curve shows the angle profile of the underlying damped pendulum dynamics. The red dots along the curve represent the sampled data, while the black dots depict the fitted data obtained using the fourth-order Runge–Kutta method with a time step size of

h

. The green line shows the high-resolution solution for the learned dynamics’ trajectory. RK4 is chosen here deliberately: its stability region extends into the right half complex plane, which is precisely what permits the learned dynamics to acquire a positive real part and appear expansive even though the sampled data are fitted well. This mechanism is analyzed in Section 2 .

Explicit Runge–Kutta ( $k=1,2,3,4$ )	Implicit Euler	Implicit Midpoint
$p(\xi)=\displaystyle\sum_{j=0}^{k}\frac{\xi^{j}}{j!}$	$p(\xi)=\dfrac{1}{1-\xi}$	$p(\xi)=\dfrac{1+\tfrac{1}{2}\xi}{1-\tfrac{1}{2}\xi}$

Table 1: Stability functions

p(\xi)

for representative one-step schemes: a polynomial for the explicit Runge–Kutta family and a rational function for the implicit schemes shown here. The common symbol

p

is retained throughout for consistency with later sections.

Panel of Figure 2 — Figure 2: Profiles of the learned quantity $\hat{\lambda}h$ using the Explicit Euler, Implicit Euler, RK4, and implicit midpoint rules as $h$ varying from 0 to 0.2. The true $\lambda=-1+4\pi i$ is set. Red curves show how $\hat{\lambda}h$ moves within each method’s stability region; arrows point toward the limit as $h\to 0$ .

Panel of Figure 3 — Figure 3: Convergence of the learned dynamic $\hat{\lambda}(h)$ to the ground truth $\lambda=-1+4\pi i$ (solid triangle) as the step size $h$ decreases from $0.2$ to $0$ for four time-integration schemes: (a) Explicit Euler, (b) RK4, (c) Implicit Euler, and (d) implicit midpoint. The red curve traces the trajectory of $\hat{\lambda}(h)$ in the complex plane, and the arrows indicate the direction of decreasing $h$ .

$y^{\prime}=\lambda y$	$\operatorname{Re}(\lambda)=0$	$\operatorname{Re}(\lambda)<0$
Explicit Euler	$\operatorname{Re}(\hat{\lambda})<0$	$\operatorname{Re}(\hat{\lambda})<0$
RK2	$\operatorname{Re}(\hat{\lambda})<0$	$\operatorname{Re}(\hat{\lambda})<0$
RK3	$\operatorname{Re}(\hat{\lambda})>0$ (for $0<h\ll 1$ )	Cond. Dissipative
Implicit Euler	$\operatorname{Re}(\hat{\lambda})>0$	Cond. Dissipative
Implicit Midpoint	$\operatorname{Re}(\hat{\lambda})=0$	$\operatorname{Re}(\hat{\lambda})<0$

Table 2: The difference between dynamics

\hat{\lambda}

learned from different numerical methods, while the true dynamics

\lambda

are either conservative

\operatorname{Re}(\lambda)=0

or dissipative

\operatorname{Re}(\lambda)<0

. The implicit midpoint rule effectively preserves the structural trends of the true dynamics across different methods.

Method/Data type	Data with damping ( $a<0$ )	Purely oscillatory data ( $a=0$ )
Explicit Euler	$\operatorname{Im}(\hat{\lambda}h)=e^{a}\sin\theta$ lagging phase	$\operatorname{Im}(\hat{\lambda}h)=\sin\theta$ small $\|\theta\|\Rightarrow$ small phase error
Implicit Euler	$\operatorname{Im}(\hat{\lambda}h)=e^{-a}\sin\theta$ leading phase	$\operatorname{Im}(\hat{\lambda}h)=\sin\theta$ small $\|\theta\|\Rightarrow$ small phase error
Impl. Midpoint	$\operatorname{Im}(\hat{\lambda}h)=\dfrac{4e^{a}\sin\theta}{\|e^{\lambda h}+1\|^{2}}$ lagging phase	$\operatorname{Im}(\hat{\lambda}h)=\dfrac{4\sin\theta}{\|e^{\lambda h}+1\|^{2}}\approx\theta$ small $\|\theta\|\Rightarrow$ small phase error

Table 3: Sign of

\operatorname{Im}(\hat{\lambda}h)

for one-step methods, inferred from data

e^{\lambda h}=e^{a+i\theta}

$M$	Learned eigenvalue $\hat{\lambda}$	$\|p(\hat{\lambda}h)\|$
Reference (no noise)	$\hat{\lambda}=\phantom{-}0.0000+1.0008i$	$\|z\|=1.00000$
$M=1$	$\hat{\lambda}=-0.0009+1.0018i$	$\|z\|=0.99990$
$M=5$	$\hat{\lambda}=\phantom{-}0.0001+1.0010i$	$\|z\|=1.00001$
$M=100$	$\hat{\lambda}=\phantom{-}0.0000+1.0009i$	$\|z\|=1.00000$

Table 4: Effect of

M

on eigenvalue estimation for

\dot{x}=ix

, identified via the implicit midpoint rule (

h=0.1

\sigma=0.1

Panel of Figure 4 — Figure 4: Learning quality under additive noise ( $\sigma=0.1$ , $h=0.1$ ). The horizontal and vertical axes are the real and imaginary parts of the state; the true solution is the unit circle (dashed). Dots are noise-corrupted observations. The right two panels correspond to $M=100$ initial points used for identification; for visual clarity, 5 trajectories are randomly selected from the pool and displayed in distinct colors. With $M=1$ the predicted trajectory drifts away from the unit circle; with $M=100$ it remains on it.

AB2	AM2	BDF2	Leap-Frog
$\displaystyle\frac{2\zeta^{2}-2\zeta}{3\zeta-1}$	$\displaystyle\frac{12\zeta^{2}-12\zeta}{5\zeta^{2}+8\zeta-1}$	$\displaystyle\frac{3\zeta^{2}-4\zeta+1}{2\zeta^{2}}$	$\displaystyle\frac{\zeta^{2}-1}{2\zeta}$

Table 5: Explicit formula of the learned parameter

\hat{\lambda}h

for two-step methods.

Method	Stable/unstable roots coexist	Repeated roots
	( $\|\zeta_{1}\|\geq 1$ and $\|\zeta_{2}\|\leq 1$ )	( $\zeta_{1}=\zeta_{2}$ )
AB2	Dark Gray region in Figure 5(a)	if $\lambda h=\dfrac{-2\pm 4\sqrt{2}i}{9}$ , $\|\zeta_{1,2}\|<1$
AM2	Dark Gray region in Figure 5(b)	if $\lambda h=\dfrac{-6\pm 4\sqrt{3}i}{7}$ , $\|\zeta_{1,2}\|>1$
BDF2	Dark Gray region in Figure 5(c)	if $\lambda h=-\dfrac{1}{2}$ , $\|\zeta_{1,2}\|<1$
Leap-Frog	always, since $\zeta_{1}\zeta_{2}=-1$	if $\lambda h=\pm i$ , $\|\zeta_{1,2}\|=1$

Table 6: Stable and unstable roots of the characteristic polynomials of two-step methods.

Method	Data with damping ( $a<0$ )	Pure oscillatory data ( $a=0$ )
AB2	Gray region in Figure 6(a)	$\operatorname{Re}(\hat{\lambda}h)<0$ for all $\zeta\in S^{1}$
AM2	Gray region in Figure 6(b)	$\operatorname{Re}(\hat{\lambda}h)<0$ for all $\zeta\in S^{1}$
BDF2	Gray region in Figure 6(c)	$\operatorname{Re}(\hat{\lambda}h)>0$ for all $\zeta\in S^{1}$
Leap-Frog	$\operatorname{Re}(\hat{\lambda}h)>0$ if $\cos(\theta)<0$	$\operatorname{Re}(\hat{\lambda}h)=0$ for all $\zeta\in S^{1}$
	$\operatorname{Im}(\hat{\lambda}h)=\cosh(a)\sin\theta$	small $\|\theta\|\implies$ small phase error

Table 7: Sign of

\operatorname{Re}(\hat{\lambda}h)

for two-step methods, inferred from root

\zeta=e^{a+i\theta}

Panel of Figure 5 — Figure 5: The dark gray region in each panel marks the region $\mathcal{R}_{P}$ . The regions in light gray are the region of absolute stability $\mathcal{R}_{A}$ of each integrator, which is enclosed by the red curve.

Panel of Figure 6 — Figure 6: The dark gray region in each panel marks the set of complex numbers $z=a+i\theta$ where $\operatorname{Re}(\hat{\lambda}h)=\operatorname{Re}({\rho(e^{z})}/{\kappa(e^{z})})>0$ .

Panel of Figure 7 — Figure 7: $\hat{\lambda}h$ versus the region of absolute stability. The gray region in each plot indicates the stability region of the selected method, while the red curve depicts the learning result, $\hat{\lambda}h=\rho(e^{\lambda h})/\kappa(e^{\lambda h})$ for $h\in(0,\pi/4)$ , where $\lambda=-4+2i$ .

Panel of Figure 8 — Figure 8: $\hat{\lambda}\equiv\hat{\lambda}(h)=h^{-1}\rho(e^{\lambda h})/\kappa(e^{\lambda h})$ for $h\in(0,\pi/4)$ . The actual value, $\lambda=-4+2i$ , is indicated by a triangle in each plot. The arrows on selected curves show the trend of $\hat{\lambda}$ as $h\rightarrow 0$ for the respective methods. The blue curves represent the ABk or AMk methods for $k=2$ , the green curves for $k=3$ , and the red curves for $k=4$ . The right plot corresponds to the Leap-Frog scheme.

$N$	$c_{1}$	$c_{2}$	$\|\zeta_{1}\|$	$\|\zeta_{2}\|$
10	$-6.0912\times 10^{-3}$	1.0007	1.1602	0.8619
20	$-1.6658\times 10^{-3}$	0.9962	1.1602	0.8732
50	$-7.8910\times 10^{-6}$	0.9946	1.1602	0.8732

Table 8: Estimated coefficients

c_{1}

and

c_{2}

for different sample sizes

N

Panel of Figure 9 — Figure 9: The $x$ -axis represents the physical time. The blue curves represent the observed data. On the left, the red curves depict the numerical solution of the differential equation $z^{\prime}=\hat{\lambda}z(t)$ obtained using the Leap-Frog scheme with step size $h=0.1$ up to time $t=10$ . On the right, the plots display the evolution of $\zeta_{1}^{n}$ .

Method	$h$	$k=1$		$k=10$
Method	$h$	$\hat{a}$	$\hat{\epsilon}$	$\hat{a}$	$\hat{\epsilon}$
Explicit Euler	$10^{-2}$	$1.9869$	$0.0299$	$1.0281$	$0.0199$
Explicit Euler	$10^{-3}$	$1.9992$	$0.0120$	$1.9181$	$0.0117$
Implicit Euler	$10^{-2}$	$2.0026$	-0.0100	$2.2552$	-0.0133
Implicit Euler	$10^{-3}$	$2.0007$	$0.0080$	$2.0748$	$0.0081$
Implicit Midpoint	$10^{-2}$	$2.0026$	$0.0100$	$2.1765$	$0.0146$

Table 9: Learned coefficients

\hat{a}

and

\hat{\epsilon}

for ( 36 ), true values

a=2

\epsilon=0.01

Panel of Figure 10 — Figure 10: Profiles of the learned convection–diffusion solution at $T=1$ . The blue dashed line denotes the initial condition $\cos(2\pi x)$ , the solid blue line is the true solution, and the red line represents the learned solution.

Panel of Figure 11 — Figure 11: Three Lotka–Volterra orbits are displayed. Filled circles mark samples from the blue orbit. In the time-series panels, the horizontal axis is time; the upper panel plots the $x$ -information and the lower panel plots the $y$ -information of the sampled points.

Panel of Figure 13 — Figure 13: Learning results with timestep $h=0.1$ using the Implicit Euler, implicit trapezoidal, and implicit midpoint methods. The top row presents trajectories generated by the learned dynamics from initial conditions within the given data set. In contrast, the bottom row shows trajectories from initial conditions outside the data set. The purple orbits in (c), (d), (g), and (f) seem to become limiting cycles.

Panel of Figure 15 — Figure 15: Learning results with timestep $h=0.01$ using the Implicit Euler, implicit trapezoidal, and implicit midpoint methods. The top row presents trajectories generated by the learned dynamics from an initial condition within the given data set. In contrast, the bottom row shows trajectories from initial conditions outside the given data set.

Panel of Figure 16 — Figure 16: The left panel shows the trajectory where data is sampled, which also encompasses a limit cycle. The right panel highlights the segment chosen for detailed analysis.

Panel of Figure 17 — Figure 17: Profiles of $|1+\lambda(x)h|$ and $|1+\hat{\lambda}(x)h|$ along the upper branch of the limit cycle.

Panel of Figure 18 — Figure 18: Real and imaginary components of the Jacobian eigenvalues for both the learned and true dynamics along the sampled branch; the x-axis ticks match the labels in Figure 16 .

Panel of Figure 19 — Figure 19: Comparison of the eigenvalues $\lambda$ and $\hat{\lambda}$ of $J_{f}$ and $J_{g_{h}}$ at points $p$ and $q$ . The eigenvalues are scaled by the step size used by the numerical integrator and plotted on the complex plane, with the gray regions indicating the integrator’s stability region.

Panel of Figure 20 — Figure 20: Vector fields using only one trajectory with different sampling timesteps.

Panel of Figure 21 — Figure 21: Amplitude errors of the learned solution component $\vartheta$ (green curves) inferred from applying four different numerical integrators to fit the data (red dots): (a) Explicit Euler, (b) Implicit Euler, (c) Implicit Midpoint, and (d) AB2. For reference, the true solutions are shown by the gray curves.

Symbol	Description
$f:\mathbb{R}^{d}\to\mathbb{R}^{d}$	True (unknown) vector field of the autonomous system $y^{\prime}=f(y)$
$\hat{f}\in\mathcal{F}$	Candidate vector field drawn from function class $\mathcal{F}$
$g$	Global minimizer of the discrete learning problem; in general $g\neq f$
$\varphi_{t}[f]\,y_{0}$	Flow map of $y^{\prime}=f(y)$ : solution at time $t$ starting from $y_{0}$
$y(t;\,y_{0})$	Solution trajectory, i.e. $y(t;\,y_{0})\equiv\varphi_{t}[f]\,y_{0}$
$\lambda\in\mathbb{C}$	Eigenvalue of $f$ ; governs the scalar model $z^{\prime}=\lambda z$

Table 10: Notation: dynamical systems and solutions.

Symbol	Description
$H$	Sampling step size (time gap between consecutive observations)
$t_{n}=nH$	$n$ -th observation time
$\mathcal{D}_{N}=\{(t_{n},y(t_{n}))\}_{n=1}^{N}$	Dataset of $N$ trajectory snapshots
$U_{M}$	Set of $M$ sampled initial conditions
$Z_{n}=e^{\lambda nh}Z_{0}$	Exact solution of the scalar model at step $n$
$z_{n}$	Numerical solution at step $n$ produced by the chosen integrator

Table 11: Notation: data and sampling.

Symbol	Description
$h$	Step size used by the numerical integrator
$m$	Number of integrator steps per observation interval ( $mh=H$ )
$S_{h}[\hat{f}]$	One-step propagator of the integrator applied to $\hat{f}$ with step size $h$
$S_{h}^{mn}[\hat{f}]$	$mn$ -fold composition of $S_{h}[\hat{f}]$ , approximating $\varphi_{mnh}[\hat{f}]$
$p(\xi)$	Stability function of a one-step method (written $R(\xi)$ in Section 4 )
$\xi=\lambda h$	Scaled eigenvalue; the fundamental argument of $p(\xi)$

Table 12: Notation: numerical integrators.

Symbol	Description
$\mathcal{R}_{A}$	Region of absolute stability: all $\xi$ for which every root of $\pi_{\xi}$ lies in the unit disc
$\mathcal{R}_{P}$	Region of partial stability: $\xi$ for which roots lie both inside and outside the unit disc
$\mathcal{F}$	Function class from which the candidate vector field $\hat{f}$ is drawn
$\hat{\xi}=\hat{\lambda}\,h$	Global minimizer of the scalar least-squares problem
$c_{1},\ldots,c_{k}$	Coefficients of the general solution of the linear multistep recurrence

Table 13: Notation: stability regions and optimization.

Numerical Integration

See all 24 papers

May 8, 2026

Learned Lagrangian Models of PDEs via Euler-Lagrange Residual Minimization

We present the first method to directly use a learned continuous Lagrangian to forecast the dynamics of systems governed by partial differential equations, exploiting the inherent conservative structure to achieve stable long-range predictions. We develop an optimization-based integrator that minimizes the squared Euler--Lagrange residual via a mesh-free near-symplectic construction on local space-time patches. Different from integrators for analytical models, integrators for learned models should decouple model error (phase error) from integration error (conservation error). By relying on optimization rather than time-stepping, we bypass the global coupling inherent to fixed discretizations, which slows time- and space-stepping and complicates learning. Our method scales linearly with domain size via Jacobi iteration, and places no structural requirements on the learned network, allowing it to be coupled with existing physics-guided machine learning (ML) methods. We validate our approach on a learned representation of a double pendulum, a one-dimensional wave equation, and a two-dimensional wave equation. Our method achieves error comparable to classical symplectic methods while generalizing to spatially varying dynamics and arbitrary boundary conditions without retraining.

Lyra Zhornyak, Eric Forgoston, M. Ani Hsieh

University of Pennsylvania · Montclair State University

Neural PDE Solvers PDE Surrogate Modeling Physics-Informed ML Lagrangian Neural Networks Numerical Integration

Ask PDF

Aug 22, 2026

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

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.

Nibodh Boddupalli, Jeff Moehlis

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

Nonlinear System Identification System Identification Neural ODEs Sparse Identification of Nonlinear Dynamics Ordinary Differential Equations Chaotic Time-Series Prediction

Ask PDF

Jul 8, 2024

Adversarial dynamical systems characterize when data-driven learning succeeds or fails

Many systems resist analytical modeling, making data-driven inference of dynamics important. Yet data-driven methods can fail to converge or generalize, leaving open a central question: When can system behavior be learned reliably from data, and when is such learning impossible? We answer this question using adversarial dynamical systems to identify the boundary between accessible and inaccessible regimes. In Koopman operator learning, a leading framework for representing nonlinear dynamics through linear spectral objects, we design optimal data-driven spectral algorithms with convergence and certification guarantees under conditions arising broadly in physical systems. This yields a convergence theory for Koopman-operator approximations and resolves a longstanding open problem in Koopman spectral analysis. Conversely, by constructing adversarial systems, we prove matching impossibility results: without these conditions, no single-sequence limiting procedure can guarantee learning, regardless of data quality. These results sharply characterize when data-driven spectral learning can succeed and when it must fail. We validate the framework on oscillators, chaotic fluid flows and Arctic sea ice concentration forecasting. In the latter, we uncover hidden modes of Arctic sea ice decline, deliver long-range forecasts with geographic error bounds, and outperform state-of-the-art dynamical and deep learning models at substantially lower computational cost, enabling real-time deployment on standard CPUs.

Matthew J. Colbrook, Igor Mezić, Alexei Stepanenko

DAMTP, University of Cambridge, Cambridge, CB3 0WA, UK. · University of California, Santa Barbara, CA 93106, USA.

Koopman Operator Learning Spectral Methods Statistical Learning Theory Chaotic Time-Series Prediction

Ask PDF

Abstract

Numerical Integration

Explore similar work

Learned Lagrangian Models of PDEs via Euler-Lagrange Residual Minimization

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

Adversarial dynamical systems characterize when data-driven learning succeeds or fails