Artifacts of Numerical Integration in Learning Dynamical Systems

Bing-Ze Lu; Richard Tsai

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.

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

.

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

).

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}

.

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.

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

.

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

.

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

.

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

.

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.

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.

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.

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.

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.

Figure 17: Profiles of

|1+\lambda(x)h|

and

|1+\hat{\lambda}(x)h|

along the upper branch of the limit cycle.

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 .

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.

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

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.

Artifacts of Numerical Integration in Learning Dynamical Systems

Abstract

Figures & tables

Numerical Integration

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

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