FNODE: Flow-Matching for data-driven simulation of constrained multibody systems

TL;DR

FNODE introduces a flow-matching neural ODE framework that learns accelerations directly from data, improving efficiency and accuracy in simulating constrained multibody systems while enforcing kinematic constraints.

Contribution

The paper proposes a novel flow-matching approach for neural ODEs that simplifies training and enhances simulation accuracy for constrained multibody dynamics.

Findings

Outperforms existing models in accuracy and efficiency

Enforces kinematic constraints effectively

Reduces training complexity and runtime

Abstract

Data-driven modeling of constrained multibody dynamics remains challenged by (i) the training cost of Neural ODEs, which typically require backpropagation through an ODE solver, and (ii) error accumulation in rollout predictions. We introduce a Flow-Matching Neural ODE (FNODE) framework that learns the acceleration mapping directly from trajectory data by supervising accelerations rather than integrated states, turning training into a supervised regression problem and eliminating the ODE-adjoint/solver backpropagation bottleneck. Acceleration targets are obtained efficiently via numerical differentiation using a hybrid fast Fourier transform (FFT) and finite-difference (FD) scheme. Kinematic constraints are enforced through coordinate partitioning: FNODE learns accelerations only for the independent generalized coordinates, while the dependent coordinates are recovered by solving the position-level constraint equations. We evaluate FNODE on single and triple mass-spring-damper systems, a double pendulum, a slider crank with and without friction, a vehicle model, and a cart-pole, and compare against MBD-NODE, LSTM, and fully connected baselines. Across these benchmarks, FNODE achieves improved prediction accuracy and training/runtime efficiency, while maintaining constraint satisfaction through the partitioning procedure. Our code and scripts are released as open source to support reproducibility and follow-on research.

Figures

Tables

Table 1. Table 1: FNODE Architecture

Layer	Number of Neurons	Activation Function	Initialization
Input Layer	$2n_z$	tanh	Xavier
Hidden Layer 1	$d_{\text{width}}$	tanh	Xavier
Hidden Layer 2	$d_{\text{width}}$	tanh	Xavier
Output Layer	$n_z$	-	Xavier

Table 2. Table 2: Hyper-parameters for the single mass spring system.

Hyper-parameters	Model
	FNODE	HNN	LNN
No. of hidden layers	2	2	2
No. of nodes per hidden layer	256	256	256
Max. epochs	450	30000	400
Initial learning rate	1e-3	1e-3	1e-3
Learning rate decay	0.98	0.98	0.98
Activation function	tanh	tanh	tanh
Loss function	MSE	MSE	MSE
Optimizer	Adam	Adam	Adam

Table 3. Table 3: Hyper-parameters for the single mass-spring-damper system.

Hyper-parameters	Model
	FNODE	MBD-NODE	LSTM	FCNN
No. of hidden layers	2	2	2	2
No. of nodes per hidden layer	256	256	256	256
Max. epochs	450	300	300	450
Initial learning rate	1e-3	1e-3	1e-3	1e-3
Learning rate decay	0.98	0.98	0.98	0.98
Activation function	tanh	tanh	Sigmoid,tanh	tanh
Loss function	MSE	MSE	MSE	MSE
Optimizer	Adam	Adam	Adam	Adam

Table 4. Table 4: Hyper-parameters for the triple-mass-spring-damper system.

Hyper-parameters	Model
	FNODE	MBD-NODE	LSTM	FCNN
No. of hidden layers	2	2	2	2
No. of nodes per hidden layer	256	256	256	256
Max. epochs	450	300	300	450
Initial learning rate	1e-3	1e-3	1e-3	1e-3
Learning rate decay	0.98	0.98	0.98	0.98
Activation function	tanh	tanh	Sigmoid,tanh	tanh
Loss function	MSE	MSE	MSE	MSE
Optimizer	Adam	Adam	Adam	Adam

Table 5. Table 5: Hyper-parameters for the double pendulum system

Hyper-parameters	Model
	FNODE	MBD-NODE	LSTM	FCNN
No. of hidden layers	2	2	2	2
No. of nodes per hidden layer	256	256	256	256
Max. epochs	450	300	300	450
Initial learning rate	3e-3	1e-3	1e-3	3e-3
Learning rate decay	0.98	0.98	0.98	0.98
Activation function	tanh	tanh	Sigmoid,tanh	tanh
Loss function	MSE	MSE	MSE	MSE
Optimizer	Adam	Adam	Adam	Adam

Table 6. Table 6: Hyper-parameters for the slider-crank mechanism.

Hyper-parameters	Model
	FNODE	MBD-NODE	LSTM	FCNN
No. of hidden layers	3	3	3	3
No. of nodes per hidden layer	256	256	256	256
Max. epochs	450	300	300	450
Initial learning rate	5e-4	1e-3	1e-3	1e-3
Learning rate decay	0.99	0.98	0.98	0.98
Activation function	tanh	tanh	Sigmoid,tanh	tanh
Loss function	MSE	MSE	MSE	MSE
Optimizer	Adam	Adam	Adam	Adam

Table 7. Table 7: MSE of various-layers FNODE acceleration prediction for the slider-crank mechanism with friction coefficient μ = 0.6 \mu=0.6

Number of Hidden Layers	MSE
2	3.70e-02
3	1.72e-02
4	1.47e-02
5	1.44e-02
6	1.41e-02

Table 8. Table 8: Hyper-parameters for the cart-pole system.

Hyper-parameters	Model
	FNODE	MBD-NODE	LSTM	FCNN
No. of hidden layers	2	2	2	2
No. of nodes per hidden layer	256	256	256	256
Max. epochs	450	400	300	450
Initial learning rate	1e-3	1e-3	1e-3	1e-3
Learning rate decay	0.98	0.98	0.98	0.98
Activation function	tanh	tanh	tanh	tanh
Loss function	MSE	MSE	MSE	MSE
Optimizer	Adam	Adam	Adam	Adam

Table 9. Table 9: Parameter ranges for the 4DOF vehicle model

Parameter	Range	Explanation
$l$	[2.7, 2.9]	wheelbase
$\beta$	[0.60, 0.75]	steering ratio
$\gamma$	[0.15, 0.25]	gear ratio
$I_{\mathrm{wheel}}$	[0.9, 1.3]	wheel inertia
$R_{\mathrm{wheel}}$	[0.32, 0.35]	wheel radius
$c_0$	[0.015, 0.025]	motor resistance torque coefficient
$c_1$	[0.025, 0.045]	motor resistance torque coefficient
${\tau }_0$	[250, 350]	stalling torque
${\omega }_0$	[1200, 1500]	no-load speed

Table 10. Table 10: Test parameters for the four test cases.

Test	$l$	$R_{\mathrm{wheel}}$	$I_{\mathrm{wheel}}$	${\tau }_0$	${\omega }_0$	$c_0$	$c_1$	$\beta$	$\gamma$
1	2.7	0.32	0.9	250	1200	0.015	0.025	0.6	0.15
2	2.9	0.35	1.3	350	1500	0.025	0.045	0.75	0.25
3	2.726	0.3281	0.952	277	1239	0.0177	0.0276	0.6405	0.163
4	2.846	0.3461	1.192	337	1419	0.0237	0.0396	0.7305	0.223

Table 11. Table 11: Hyper-parameters for the 4DOF vehicle model.

Hyper-parameters	Model
	FNODE
No. of hidden layers	6
No. of nodes per hidden layer	256
Max. epochs	300
Initial learning rate	1e-4
Learning rate decay	0.98
Activation function	tanh
Loss function	MSE
Optimizer	AdamW

Table 12. Table 12: Training time cost for the models. RK4 denotes the 4th-order Runge-Kutta method. The models are trained on Nvidia GeForce 5070Ti GPU with 16 GB memory.

Test Case	Model	Integrator	Time Cost (s)
Single Mass Spring Damper	FNODE	–	112.46
	MBD-NODE	RK4	201.90
	FCNN	–	103.28
	LSTM	–	171.99
Triple Mass Spring Damper	FNODE	–	108.23
	MBD-NODE	RK4	244.80
	FCNN	–	102.78
	LSTM	–	172.05
Double Pendulum	FNODE	–	89.33
	MBD-NODE	RK4	212.09
	FCNN	–	88.04
	LSTM	–	148.41
Slider Crank	FNODE	–	565.01
	MBD-NODE	RK4	1185.64
	FCNN	–	453.50
	LSTM	–	813.11
Cart Pole	FNODE	–	59.46
	MBD-NODE	RK4	144.51
	FCNN	–	59.36
	LSTM	–	96.08
4dof Vehicle	FNODE	–	3304.32
	MBD-NODE	RK4	–
	FCNN	–	–
	LSTM	–	–

Table 13. Table 13: Summary of the numerical MSE for different models.

Test Case	Error
	FNODE	MBD-NODE	LSTM	FCNN
Single Mass-Spring-Damper	7.6e-5	1.4e-3	1.5e-2	1.6e-1
Triple Mass-Spring-Damper	6.1e-3	6.8e-3	2.5e-1	9.0e-2
Double Pendulum	1.4e-1	2.3e-1	7.6e-1	1.9e0
Slider Crank	1.7e-1	8.7e-1	1.4e0	3.1e0
Cart Pole	6.7e-4	1.1e-3	1.6e-1	3.5e-1

Equations

$$\left[\begin{array}{cc}\mathbf{M}& {\mathbf{\Phi }}_{\mathbf{q}}^T\\ {\mathbf{\Phi }}_{\mathbf{q}}& 0\end{array}\right]\left[\begin{array}{c}\ddot{\mathbf{q}}\\ \lambda \end{array}\right]=\left[\begin{array}{c}{\mathbf{F}}_{\mathbf{e}}\\ {\gamma }_{\mathbf{c}}\end{array}\right],$$

$$\frac{d\mathbf{z}(t)}{dt}=f_{\Theta }(\mathbf{z}(t),t),$$

$$\mathbf{z}(t)=\mathbf{z}(0)+{\int }_0^t{f}_{\Theta }(\mathbf{z}(\tau ),\tau ),d\tau =\Phi (\mathbf{z}(0),f_{\Theta },t),$$

$$\ddot{\mathbf{z}}(t,\mathbf{\mu })=f_{\Theta }(\mathbf{Z}(t,\mathbf{\mu })),$$

$$\mathbf{Z}(0,\mathbf{\mu })=({\mathbf{z}}^T(0,\mathbf{\mu }),{\dot{\mathbf{z}}}^T(0,\mathbf{\mu }){)}^T,$$

$$\mathbf{\mu }=({\mu }_1,{\mu }_2,\dots ,{\mu }_{n_{\mathbf{\mu }}}{)}^T\in {\mathbb{R}}^{n_{\mathbf{\mu }}}.$$

$$\mathbf{z}(t,\mathbf{\mu })=(z^1(t,\mathbf{\mu }),...,z^{n_z}(t,\mathbf{\mu }){)}^T\in {\mathbb{R}}^{n_z},$$

$$\dot{\mathbf{z}}(t,\mathbf{\mu })=({\dot{z}}^1(t,\mathbf{\mu }),\dots ,{\dot{z}}^{n_z}(t,\mathbf{\mu }){)}^T\in {\mathbb{R}}^{n_z}.$$

$$f_{\Theta }(\mathbf{z}(t),\mathbf{\mu }):{\mathbb{R}}^{2n_z}\times {\mathbb{R}}^{n_{\mathbf{\mu }}}\to {\mathbb{R}}^{n_z}.$$

$$\mathbf{Z}(t,\mathbf{\mu })=\mathbf{Z}(0,\mathbf{\mu })+{\int }_0^t{f}_{\Theta }(\mathbf{Z}(\tau ,\mathbf{\mu }))d\tau =\Phi (\mathbf{Z}(0,\mathbf{\mu }),t)\circ f_{\Theta }={\Phi }_{\Theta }(\mathbf{Z}(0,\mathbf{\mu }),t).$$

$$\mathbf{z}(t)=\sum _{k=-\infty }^{\infty }{\mathbf{Z}}_k{e}^{2\pi ikt/L}.$$

$$\frac{d\mathbf{z}}{dt}=\sum _{k=-\infty }^{\infty }\frac{2\pi ik}{L}{\mathbf{Z}}_k{e}^{2\pi ikt/L}.$$

$${\Delta }_{boundary}=z((N-1)\Delta t)-z(0)\ne 0.$$

$$Z(k)\sim \frac{{\Delta }_{boundary}}{2\pi ik/N}.$$

$$Z^{\prime }(k)=\frac{2\pi ik}{N\Delta t}\cdot Z(k)\sim \frac{{\Delta }_{boundary}}{2\Delta t}.$$

$${ϵ}_{\text{Gibbs}}\sim \frac{|{\Delta }_{\text{boundary}}|}{2\Delta t}.$$

$$\ell (t)=a+bt,\left[\begin{array}{cc}a& b\end{array}\right]=\underset{a,b}{\mathrm{argmin}}\sum _i(z_i-(a+b{t}_i){)}^2,$$

$$z_d(t)=z(t)-\ell (t).$$

$${\Delta }_{new}=z_d(N-1)-z_d(0)={\Delta }_{\text{boundary}}-bL.$$

$$z_e(t)$$

$$\tau (t)$$

$$w(t)=\{\begin{array}{ll}\frac{1}{2}[1+\cos (\pi \frac{\alpha L/2-t}{\alpha L/2}\left)\right],& 0\le t<\frac{\alpha L}{2},\\ 1,& \frac{\alpha L}{2}\le t\le L-\frac{\alpha L}{2},\\ \frac{1}{2}[1+\cos (\pi \frac{t-(L-\alpha L/2)}{\alpha L/2}\left)\right],& L-\frac{\alpha L}{2}<t\le L,\end{array}\alpha =0.2$$

$$G(k)=\exp [-\frac{1}{2}(\frac{k}{\sigma }{)}^2],\sigma \approx \frac{N}{20},$$

$${\mathbf{Z}}_e(k)=\sum _{n=0}^{N-1}{\mathbf{z}}_e[n]w[n]e^{-j2\pi kn/N}.$$

$${\mathbf{Z}}_e^{\prime }(k)=j{\omega }_k{\mathbf{Z}}_e(k)G(k),$$

$$\dot{\mathbf{z}}[n]={\mathcal{F}}^{-1}[{\mathbf{Z}}_e^{\prime }(k)]=\frac{1}{N}\sum _{k=-N/2}^{N/2-1}{\mathbf{Z}}_e^{\prime }(k)e^{j2\pi kn/N}.$$

$${\rho }_{\delta }(r)=\{\begin{array}{ll}\frac{1}{2}{r}^2,& |r|\le \delta ,\\ \delta \left(|r|-\frac{1}{2}\delta \right),& |r|>\delta .\end{array}$$

$${\dot{z}}_t^{(5)}=\frac{-z(t+2\Delta t)+8z(t+\Delta t)-8z(t-\Delta t)+z(t-2\Delta t)}{12\Delta t},$$

$$Z_{n+1}=\Phi (Z_n),{\hat{Z}}_{n+1}={\Phi }_{\Theta }({\hat{Z}}_n),$$

$$E_{n+1}={\hat{Z}}_{n+1}-Z_{n+1}={\Phi }_{\Theta }({\hat{Z}}_n)-\Phi (Z_n).$$

$$E_{n+1}={\Phi }_{\Theta }(Z_n+E_n)-\Phi (Z_n).$$

$${\Phi }_{\Theta }(Z_n+E_n)={\Phi }_{\Theta }(Z_n)+\frac{\partial {\Phi }_{\Theta }}{\partial Z}{|}_{Z_n}{E}_n+\mathcal{O}(\Vert E_n{\Vert }^2),$$

Full text

K****eywords Multibody dynamics $\cdot$ Neural ODE $\cdot$ Constrained dynamics $\cdot$ Friction dynamics $\cdot$ Scientific machine learning $\cdot$ Flow-Matching

MSC2020: 68T07, 68T45

1 Introduction

Multibody dynamics is relevant in simulation-based engineering, with applications in numerous fields, e.g., robotics [29], vehicle and train dynamics [2, 47], biomechanics [1], aircraft dynamics [53] and molecular dynamics [41], to name a few. Traditionally, modeling multibody systems has relied on first-principles approaches such as Newton-Euler formulations [35], Lagrangian mechanics [4], and Hamiltonian mechanics [6]. While these classical methods provide accurate and interpretable models, they are challenged in scenarios where the underlying dynamics are partially unknown, highly complex, or subject to unmodeled effects such as friction, contact forces, or material nonlinearities [10].

Machine learning has opened the door to data-driven simulation of dynamical systems, offering the potential to learn complex behaviors directly from observational data [32, 20]. Early approaches in this domain employed classical machine learning architectures such as LSTM [43, 19] and FCNN [34] to learn discrete-time mappings from current states to future states. In addition, recent work [12, 13] explored fixed-time increment DNN methods combined with dimensionality reduction techniques for efficient multibody dynamics simulation. While these methods have shown promise in short-term prediction tasks, they demonstrated limited out-of-distribution generalization. Specifically, these black-box predictors learn fixed-timestep transitions $f:s_{t}\mapsto s_{t+\Delta t}$, which limits their temporal resolution and also causes error accumulation that grows exponentially with the prediction horizon [54]. Furthermore, such discrete-time approaches often fail to enforce fundamental conservation laws when applied in conjunction with conservative systems [15]. Physics-informed approaches have been proposed to address these limitations by incorporating physical constraints and laws into the model, presenting improved reliability through enforcement of physical principles in neural network training [39, 40, 50]. Neural Ordinary Differential Equations (NODEs) [8] marked a major step in continuous-time learning by representing system dynamics as neural vector fields integrated with standard ODE solvers. This formulation handles irregular time sampling and directly connects learning with established numerical integration theory [21, 31].

Meanwhile, hybrid methods combining first-principle models with NODEs have shown promise in specific applications through automatic parameter tuning [55]. However, these approaches typically assume that the functional form of the governing equations (or physics structure) is known a priori, which can limit their ability to generalize when the underlying mechanisms are partially unknown or misspecified. To improve physical consistency, several structure-preserving learning directions have been proposed, including energy-based formulations such as Hamiltonian and Lagrangian Neural Networks [14, 9], symplectic-geometry-based approaches such as Symplectic ODE-Net and Symplectic Machine Learning (Sym-ML) frameworks [52, 38], and methods that target numerical stability and explicit constraint satisfaction, e.g., Stabilized Neural Differential Equations [48]. For flow-matching-based methods, recent work proposes a contrastive feature-alignment remedy for systematic failure modes that arise in the low-noise regime [49]. Separately, another work develops improved training strategies for rectified flows [24]. In a different modeling paradigm, operator-learning methods such as DeepONets [30] and Fourier Neural Operators [26] learn solution operators to generalize across parameterized families of dynamical systems. However, despite their strengths, NODEs suffer from a bottleneck: the adjoint sensitivity method for backpropagation requires solving an additional backward ODE, which increases computational cost and memory use, and limits scalability for stiff or high-dimensional systems [11, 45].

Recent advances in flow-matching algorithms [28, 27] show that vector fields can be learned by directly supervising derivatives at sampled points, avoiding backward ODE solves [22]. The approach of learning state-acceleration mapping relationship has been successfully applied to simulate flexible body dynamics [37]. Applying these ideas to constrained multibody dynamics, we train neural networks to approximate accelerations directly, eliminating costly integration in the learning loop. We propose FNODE, a framework for data-driven modeling of multibody systems. Instead of predicting future states via time integration, FNODE trains a neural network to learn directly the acceleration vector field – the natural quantity in mechanical systems governed by force laws [3]. Kinematic constraints are enforced through a coordinate-partitioning approach [46], in which FNODE learns the acceleration field only for the independent generalized coordinates, while the dependent coordinates are recovered implicitly by enforcing the kinematic constraint equations at the position level. This reformulation removes the need for ODE solvers in backpropagation, improving computational efficiency while preserving accuracy. By learning accelerations directly, the model captures the instantaneous dynamics of the system, leading to better generalization and more stable long-term predictions. Moreover, FNODE naturally accommodates various numerical differentiation schemes for computing acceleration targets from position data, including finite differences and FFT-based spectral methods, allowing practitioners to balance between noise robustness and accuracy based on their data characteristics [25, 7].

Directly learning accelerations requires accurate derivative information, which can be difficult to obtain from noisy trajectory data. To address this, we use a hybrid scheme that combines the accuracy of FFT-based spectral differentiation with the stability of FD near trajectory boundaries. Our numerical experiments show that, with proper differentiation techniques, the advantages of acceleration-based learning outweigh the challenges of derivative estimation.

The following are our main contributions in this paper:

We introduce FNODE, an approach for data-driven simulation of dynamical systems that learns acceleration vector fields directly from trajectory data. This removes the backpropagation bottleneck of conventional NODEs and reduces training cost while preserving accuracy. 2. 2.

We present a numerical differentiation framework for constructing high-quality acceleration targets, combining FFT-based spectral differentiation with finite differences to balance accuracy and stability. 3. 3.

We evaluate FNODE on a diverse set of constrained multibody dynamics benchmarks, showing consistent improvements over baselines—including MBD-NODE, LSTMs, and FCNNs—in long-term prediction accuracy, out-of-distribution generalization, and computational efficiency. 4. 4.

To support reproducibility, we provide an open-source implementation of FNODE together with benchmark problems, implementations of baseline methods, data generation tools, and evaluation protocols.

2 Methodology

2.1 The Multibody Dynamics Problem

We employ an MBD formulation in redundant generalized coordinates in which Lagrange multipliers are used to account for the presence of kinematic constraints in the equations of motion [36].

[TABLE]

Above, the mass matrix $\mathbf{M}$ encodes the system’s inertia, while $\mathbf{\Phi_{q}}$ is the Jacobian of the kinematic constraints. The generalized coordinates are collected in $\mathbf{q}$, with $\ddot{\mathbf{q}}$ denoting accelerations. The Lagrange multipliers $\lambda$ correspond to the constraint forces. The right-hand side includes $\mathbf{F_{e}}$, the vector of external and velocity-dependent forces, and $\mathbf{\gamma_{c}}$, the contribution from the second time derivative of the kinematic constraints. For details, see [36, 17].

2.2 Neural Ordinary Differential Equations for Multibody System Dynamics

2.2.1 Neural Ordinary Differential Equation (NODE)

For a general dynamical system, the NODE framework assumes the existence of an underlying continuous-time process that governs the evolution of the system states. To that end, we assume a hidden state $\mathbf{z}(t)\in\mathbb{R}^{n_{z}}$ with dynamics

[TABLE]

where $f_{\Theta}:\mathbb{R}^{n_{z}}\times\mathbb{R}^{+}\to\mathbb{R}^{n_{z}}$ is a neural network parameterized by $\mathbf{\Theta}$. The future state $\mathbf{z}(t)$ is obtained by solving the initial value problem

[TABLE]

with $\Phi$ denoting the chosen numerical ODE solver. The parameters $\mathbf{\Theta}$ are identified from trajectory data ${\mathbf{z}(t_{i})}_{i=1}^{n}$. Since time steps need not be uniform, adaptive and variable-step integrators can be used, which is particularly advantageous in multibody dynamics where time scales vary widely.

2.2.2 Flow-Matching in NODE

FNODE operates on an augmented state space defined as $\bm{Z}(t)=(\bm{z}^{T}(t),\dot{\bm{z}}^{T}(t))^{T}$ that is used to model the acceleration directly:

[TABLE]

where the initial conditions for positions and velocities are defined as:

[TABLE]

and the parameter vector capturing problem-specific characteristics such as material properties, geometric parameters, or external loading conditions is

[TABLE]

The generalized coordinates and velocity are

[TABLE]

[TABLE]

The neural network function processes this enriched input space that is modulated by the parameter vector $\bm{\mu}$ to predict acceleration values directly

[TABLE]

In the inference phase, the forward integration process maintains the same mathematical structure as traditional NODEs with the definition from Eq. (3) and Eq. (4).

[TABLE]

Figure 1 shows the discretized forward pass of FNODE. The network maps the system state to accelerations using the three-layer architecture in Table 1, with tanh activations and Xavier initialization.

2.3 Calculation of Acceleration

Two primary approaches are employed depending on the characteristics of the available dataset: FFT spectral differentiation and FD methods.

2.3.1 FFT Spectral Differentiation

Continuous Fourier Representation

For dynamic systems which satisfy periodicity assumptions, FFT-based spectral differentiation provides high accuracy in time derivative calculation. We begin by presenting the ideal mathematical framework for a smooth periodic function $z(t)$ with time interval $L=N\Delta t$ and uniform time step $\Delta t$. The exact Fourier series representation is:

[TABLE]

The time derivative computed by the exact spectral representation assumes then the expression:

[TABLE]

Discrete Implementation and Artifacts

Eqs. (10)–(11) describe the ideal continuous-time Fourier representation for a smooth periodic trajectory, and motivate spectral differentiation through a frequency-domain multiplier. In practice, however, we only observe a finite set of samples and compute derivatives using the discrete Fourier transform (DFT)/FFT.

In this discrete setting, sampling (finite $\Delta t$) introduces a Nyquist limit $f_{N}=\tfrac{1}{2\Delta t}$ and thus possible aliasing of unresolved high-frequency content, while finite-length windowing and imperfect periodicity introduce spectral leakage. These artifacts are not properties of the exact Fourier series itself, but become central once differentiation is implemented on sampled data. Moreover, since differentiation multiplies each mode by its (angular) frequency, high-frequency components—including aliased noise—are amplified in the derivative estimates. Accordingly, in the sampled DFT/FFT implementation, aliasing and finite-window effects (including Gibbs-type ringing) may introduce spurious high-frequency oscillations.

One particularly important manifestation of these discrete-sampling artifacts is the Gibbs error, which arises from boundary discontinuity in the periodic expansion:

[TABLE]

When the boundary discontinuity of magnitude $\Delta_{\text{boundary}}$ exists as defined in Eq. (12), the Fourier coefficients exhibit asymptotic behavior

[TABLE]

The derivative operation multiplies each Fourier coefficient in Eq. (13) by $\frac{2\pi ik}{N\Delta t}$, resulting in

[TABLE]

This implies that, in the presence of a boundary discontinuity, the differentiated Fourier modes do not decay with $|k|$ in the high-frequency range; instead, their magnitude remains $O\!\left(\tfrac{|\Delta_{\text{boundary}}|}{\Delta t}\right)$ as $|k|\to\infty$ (cf. Eq. (14)). As a result, the truncated Fourier reconstruction of the derivative can exhibit spurious high-frequency oscillations (Gibbs-type ringing), with a characteristic (order-of-magnitude) error scale

[TABLE]

The error is independent of $N$ but inversely proportional to the time step, which means it would be particularly problematic for fine temporal resolution.

To address the Gibbs phenomenon problem described in Eq. (15) and other boundary discontinuities, least-square detrending is applied:

[TABLE]

where $\ell(t)$ is the linear trend function to approximate the trend of the raw dataset. The detrended dataset is formulated using Eq. 16 as

[TABLE]

The new discontinuity value is formulated from Eq. (17) as

[TABLE]

Even though the detrend narrows the gap of the two ends of dataset as described in Eq. (18), the first-order derivative of the dataset is still not continuous. In order to further reduce Gibbs error while maintain the spectrum information as much as possible, we adopt mirror reflection and cosine taper.

[TABLE]

Above, $\tau(i)$ is the cosine window weight and $M$ is the length of Mirror Reflection. The amplitude of the dataset can be smoothly transited through this process. A Tukey window is then applied, tapering the signal to zero at both ends so that it appears periodic and suitable for FFT-based differentiation:

[TABLE]

Detrending, mirror reflection, and Tukey windowing smooth structural discontinuities, but FFT differentiation still amplifies residual noise and high-frequency artifacts. To suppress these effects, we apply a Gaussian low-pass filter in the frequency domain, defined as

[TABLE]

where $k$ is the frequency index and $\sigma$ controls the filter bandwidth. The choice of $\sigma\approx N/20$ ensures that the filter effectively suppresses high-frequency noise while preserving the essential spectral content of the signal.

Complete FFT-Based Differentiation Pipeline

The complete FFT-based differentiation process is formulated as follows. Let $t_{n}=n\Delta t$ for $n=0,1,\dots,N-1$, and denote the sampled (preprocessed) sequence by $\mathbf{z}_{e}[n]:=\mathbf{z}_{e}(t_{n})$. First, we compute the discrete Fourier transform (DFT) of the preprocessed signal $\mathbf{z}_{e}[n]$ with the applied window function $w[n]$ from Eq. (20):

[TABLE]

This approach combines spectral differentiation with noise suppression. The FFT-based pipeline is most reliable for smooth, near-periodic trajectories after preprocessing; for non-periodic or non-smooth signals (e.g., impacts/contact events), Fourier truncation can lead to Gibbs oscillations and boundary artifacts, and any explicit differentiation operator (FFT or FD) amplifies measurement noise. Therefore, when FNODE is trained using explicit acceleration labels obtained via differentiation, the current differentiation choices primarily limit applicability to smooth regimes. Note that practical error diagnostics for differentiated signals, including a reconstruction-consistency check based on integrating the derivatives and comparing against the original trajectory, are described in Section 2.3.2.

Since acceleration labels are obtained via explicit differentiation, occasional high-frequency artifacts (e.g., due to boundary treatment, FFT/finite-difference effects, or nonsmooth events such as frictional transitions) can manifest as outliers in the target residual. In such cases, replacing the mean-squared error (MSE) with a robust regression loss can reduce sensitivity to outliers without changing the FNODE architecture. A common choice is the Huber loss applied to the component-wise residual $r=\hat{a}-a$:

[TABLE]

For small residuals the loss behaves like MSE, while large residuals are down-weighted (linear growth), making training less sensitive to impulsive artifacts. In this manuscript we focus on noise-free trajectories for clarity and do not perform a systematic study over injected noise types/levels (e.g., Gaussian or impulse noise); however, robust losses such as (23) are a standard mitigation when differentiated targets contain occasional spikes.

2.3.2 Finite Difference Method

For non-periodic datasets, finite-difference (FD) schemes provide a simple and robust alternative to FFT-based differentiation. Let $\dot{z}_{t}$ denote the FD estimate of $\dot{z}(t)$. Using the central-difference formula, $\dot{z}_{t}=\tfrac{z(t+\Delta t)-z(t-\Delta t)}{2\Delta t}$, and thus $\dot{z}(t)=\dot{z}_{t}+\mathcal{O}(\Delta t^{2})$ [18]. At the boundaries, we use one-sided formulas, $\dot{z}_{t}=\tfrac{z(t+\Delta t)-z(t)}{\Delta t}$ and $\dot{z}_{t}=\tfrac{z(t)-z(t-\Delta t)}{\Delta t}$, and hence $\dot{z}(t)=\dot{z}_{t}+\mathcal{O}(\Delta t)$ in either case. Although less accurate than spectral differentiation, FD methods are more stable for short trajectories and at dataset boundaries, making them complementary in our hybrid FFT–FD approach.

Higher-order FD discretizations are also available. For example, let $\dot{z}^{(5)}_{t}$ denote the five-point central FD estimate,

[TABLE]

so that $\dot{z}(t)=\dot{z}_{t}^{(5)}+\mathcal{O}(\Delta t^{4})$, together with matching higher-order one-sided stencils near boundaries. In our setting, trajectories are short and derivative labels are obtained from sampled data; wider stencils increase the number of samples affected by boundary treatment and can amplify high-frequency noise. We therefore adopt the three-point central stencil in the interior and one-sided stencils at the endpoints as a conservative accuracy–robustness tradeoff.

To obtain a practical estimate of differentiation error on a given dataset, one may use two complementary consistency checks. First, compute derivative estimates with two stencil orders and define a relative disagreement indicator $\eta^{\prime}=\tfrac{\|\dot{z}_{t}^{(5)}-\dot{z}_{t}^{(3)}\|_{2}}{\|\dot{z}_{t}^{(5)}\|_{2}}$, where $\dot{z}_{t}^{(3)}$ and $\dot{z}_{t}^{(5)}$ denote the three-point and five-point FD estimates of $\dot{z}(t)$ (and analogously for $\ddot{z}(t)$), and $\|\cdot\|_{2}$ denotes the discrete $\ell_{2}$ norm over the sampled trajectory. Second, define a reconstruction-consistency indicator by integrating the differentiated signals and comparing to the original trajectory, e.g., $\tilde{z}_{t+\Delta t}=\tilde{z}_{t}+\Delta t\,\dot{z}_{t}$ with $\tilde{z}_{0}=z_{0}$, and form the normalized reconstruction residual $\eta_{z}=\tfrac{\|z_{t}-\tilde{z}_{t}\|_{2}}{\|z_{t}\|_{2}}$. These diagnostics provide a practical way to quantify sensitivity to the differentiation scheme and to detect differentiation artifacts when constructing acceleration targets for FNODE training; a systematic evaluation of $\eta^{\prime}$ and $\eta_{z}$ across datasets is deferred to future work and falls outside the scope of the present contribution.

2.4 Qualitative Analysis of Error Propagation in Integration

While FNODE avoids numerical integration during training, its long-term accuracy depends on how errors propagate through the ODE solver. To carry out a basic error analysis following a classical framework [16, 5], let the true system state at step $n$ be $Z_{n}$, advanced by the exact evolution operator $\Phi$, and let $\hat{Z}_{n}$ denote the predicted state advanced by the learned solver $\Phi_{\Theta}$. The global error is $E_{n}=\hat{Z}_{n}-Z_{n}$. The exact and predicted updates are

[TABLE]

so that

[TABLE]

Substituting $\hat{Z}_{n}=Z_{n}+E_{n}$ gives

[TABLE]

A first-order Taylor expansion around $Z_{n}$ yields

[TABLE]

where $\tfrac{\partial\Phi_{\Theta}}{\partial Z}\big|_{Z_{n}}$ denotes the Jacobian of the learned flow map $\Phi_{\Theta}$ with respect to the state, evaluated at $Z_{n}$. Consequently, one obtains

[TABLE]

Qualitatively, each step combines a local error (discrepancy between $\Phi_{\Theta}$ and $\Phi$) and a propagated error (amplification of the previous error by the solver Jacobian), showing how small inaccuracies in the learned acceleration field can accumulate over time, and thus underscoring the need for accurate acceleration targets during training.

To connect integration error with training error and differentiation-induced label error, we decompose the one-step defect. Let the true augmented dynamics be $\ddot{z}=f(Z)$ and let $\hat{f}_{\Theta}(Z)$ denote the learned vector field (FNODE trained on numerically differentiated acceleration labels). Let $\Psi_{\Delta t}(\cdot;f)$ denote one step of a numerical integrator with step size $\Delta t$ applied to the vector field $f$. In this view, $\Phi(Z)=\Psi_{\Delta t}(Z;f)$ and $\Phi_{\Theta}(Z)=\Psi_{\Delta t}(Z;\hat{f}_{\Theta})$ represent the corresponding one-step maps. The true and learned rollouts satisfy $Z_{n+1}=\Psi_{\Delta t}(Z_{n};f)$ and $\widehat{Z}_{n+1}=\Psi_{\Delta t}(\widehat{Z}_{n};\hat{f}_{\Theta})$, and the global error $E_{n}=\widehat{Z}_{n}-Z_{n}$ evolves as $E_{n+1}=\Psi_{\Delta t}(Z_{n}+E_{n};\hat{f}_{\Theta})-\Psi_{\Delta t}(Z_{n};f)$.

Adding and subtracting $\Psi_{\Delta t}(Z_{n};\hat{f}_{\Theta})$ yields $E_{n+1}=\delta_{n}+\bigl(\Psi_{\Delta t}(Z_{n}+E_{n};\hat{f}_{\Theta})-\Psi_{\Delta t}(Z_{n};\hat{f}_{\Theta})\bigr)$, where $\delta_{n}:=\Psi_{\Delta t}(Z_{n};\hat{f}_{\Theta})-\Psi_{\Delta t}(Z_{n};f)$ is the one-step defect at the true state. Under standard smoothness assumptions on the true dynamics and assuming the use of a consistent order-$p$ numerical integrator whose induced numerical flow satisfies a Lipschitz-type stability bound over the rollout region, the one-step defect can be bounded as: $\|\delta_{n}\|\leq C_{\mathrm{int}}\Delta t^{p+1}+C_{\mathrm{fld}}\Delta t\,\|\hat{f}_{\Theta}(Z_{n})-f(Z_{n})\|$, i.e., it combines the local truncation error of the time integrator and the effect of vector-field mismatch over one step.

The learned-field error has two contributions. FNODE is trained on numerically differentiated acceleration labels $\tilde{\ddot{z}}=\ddot{z}+e_{\mathrm{diff}}$, so $\hat{f}_{\Theta}-f=(\hat{f}_{\Theta}-\tilde{f})+(\tilde{f}-f)$, where $\tilde{f}$ is the vector field consistent with the differentiated labels. The first term is controlled by the training objective; in particular, if $L(\Theta)=\tfrac{1}{N}\sum_{i}\|\hat{f}_{\Theta}(Z_{i})-\tilde{f}(Z_{i})\|_{2}^{2}$, then $\varepsilon_{\mathrm{train}}:=\sqrt{L(\Theta)}$ provides an RMS proxy for the learned-field mismatch on the training distribution. The second term reflects differentiation artifacts and can be assessed through the diagnostics introduced in Section 2.3.2, such as the stencil-disagreement indicator $\eta^{\prime}$ and the trajectory reconstruction residual $\eta_{z}$.

Finally, the second term in $E_{n+1}$ captures error propagation through integration. If the integrator is stable and satisfies $\bigl\|\tfrac{\partial}{\partial Z}\Psi_{\Delta t}(Z;\hat{f}_{\Theta})\bigr\|\leq 1+C\Delta t$ over the rollout region, then $\|\Psi_{\Delta t}(Z_{n}+E_{n};\hat{f}_{\Theta})-\Psi_{\Delta t}(Z_{n};\hat{f}_{\Theta})\|\leq(1+C\Delta t)\|E_{n}\|$, which implies the recursion $\|E_{n+1}\|\leq(1+C\Delta t)\|E_{n}\|+\|\delta_{n}\|$. Consequently, $\|E_{n}\|$ grows at most like $(1+C\Delta t)^{n}$, and for $T=n\Delta t$ one obtains $\|E(T)\|\lesssim\exp(CT)\bigl(\|E(0)\|+\sum_{k}\|\delta_{k}\|\bigr)$. This shows explicitly how integration error accumulation couples to (i) the integrator order and step size through $C_{\mathrm{int}}\Delta t^{p+1}$, (ii) the training error through $\varepsilon_{\mathrm{train}}$, and (iii) differentiation-induced label error through the diagnostics $\eta^{\prime}$ and $\eta_{z}$. For the three error sources, (i) could be illustrated by the trajectory error, while (ii) and (iii) represent the vector-field mismatch in acceleration. These three error sources are further analyzed below.

2.4.1 Integrator error

For a $p$th-order accurate time integrator, the global trajectory-level MSE is represented as:

[TABLE]

where $\Phi(Z_{n})$ is the one-step update produced using the ground truth acceleration with step size $\Delta t$. As mentioned above, the single-step truncation error is $\mathcal{O}(\Delta t^{p+1})$, so the error over $t\in[0,T]$ is $\mathcal{O}(\Delta t^{p})$ uniformly, then

[TABLE]

Note that since RK4 has $p=4$, the global discretization error (in the state) over a fixed time horizon typically scales as $\mathcal{O}(\Delta t^{4})$ under standard smoothness assumptions.

2.4.2 Neural Network learned field mismatch

We train FNODE with analytical accelerations $\ddot{z}$ and acceleration calculated by FD $\tilde{\ddot{z}}$ as targets and evaluate the predicted accelerations $\mathring{\ddot{z}}_{n}$ (trained on $\ddot{z}$) and $\hat{\ddot{z}}$ (trained on $\tilde{\ddot{z}}$) on the ground truth $\ddot{z}_{n}$, reporting $\mathring{\mathrm{MSE}}_{\ddot{z}}(\Delta t)$ and $\widehat{\mathrm{MSE}}_{\ddot{z}}(\Delta t)$:

[TABLE]

[TABLE]

2.4.3 Differentiation-induced label error

To interpret FD slopes, we partition indices into an interior set $\mathcal{I}$ and a boundary set $\mathcal{B}=\{0,N-1\}$, where $|\mathcal{B}|=2$ is a constant independent of $N$ (the samples where one-sided stencils are used) and $\mathcal{I}=\{1,\dots,N-2\}$, where $|\mathcal{I}|=N-2$. Define the pointwise FD acceleration error $\tilde{e}_{n}=\tilde{\ddot{z}}_{n}-\ddot{z}_{n}$. Let $q_{\mathcal{I}}$ and $q_{\mathcal{B}}$ denote the order of the finite-difference method used to compute accelerations in the interior and boundary regions, respectively, i.e., $\|\tilde{e}_{n}\|=\mathcal{O}(\Delta t^{q_{\mathcal{I}}})$ for $n\in\mathcal{I}$ and $\|\tilde{e}_{n}\|=\mathcal{O}(\Delta t^{q_{\mathcal{B}}})$ for $n\in\mathcal{B}$. Here $\mathcal{I}$ contains $N-\mathcal{O}(1)$ interior time indices, whereas $\mathcal{B}$ contains only $\mathcal{O}(1)$ boundary indices (e.g., the first/last few samples required by the stencil). Over a fixed time horizon $T$, the number of time samples satisfies $N\sim T/\Delta t$.

Then the interior and boundary contributions to acceleration MSE scale as

[TABLE]

[TABLE]

Since $N-1=\frac{T}{\Delta t}$ and the physical time is not extremely long ($10^{6}$ or above), $\frac{1}{N-1}=\frac{\Delta t}{T}\sim\mathcal{O}(\Delta t)$. Besides, as $N\gg 1$,

[TABLE]

Then we can rewrite Eq.(25) as:

[TABLE]

The global FD acceleration MSE is:

[TABLE]

The overall MSE in Eq. (27) is dominated, for sufficiently small $\Delta t$, by the term with the smaller exponent in $\Delta t$ (i.e., the term that decays more slowly and therefore yields a larger error). Therefore, according to Eq. (24) and Eq. (26), if $2q_{\mathcal{I}}<2q_{\mathcal{B}}+1$, the interior error dominates, and vice versa.

To simplify step-size convergence order comparisons, the root mean squared error (RMSE) can also be used to quantify the error. Although RMSE makes it difficult to directly separate the contributions from $\mathcal{I}$ and $\mathcal{B}$, its effect on the observed convergence order is universal. Therefore, the order of $\widetilde{\mathrm{RMSE}}_{\ddot{z},\mathcal{I}}(\Delta t)$ is $q_{\mathcal{I}}$, whereas the order of $\widetilde{\mathrm{RMSE}}_{\ddot{z},\mathcal{B}}(\Delta t)$ is $q_{\mathcal{B}}+\tfrac{1}{2}$.

For a fixed-order time integrator with a fixed step size, any acceleration error, regardless of its source, enters the same numerical update operator at each step. Therefore, when comparing rollouts generated with the same integrator and $\Delta t$, differences observed at the trajectory level primarily reflect differences in the underlying acceleration (vector-field) errors, rather than differences in the integrator. Consequently, to assess the robustness of our numerical differentiation, it is necessary to compare, under identical integration settings, the learned vector-field mismatch against the numerical-differentiation error.

The quantitative error-propagation analysis that isolates the contribution of each source is provided in Section 3.2.

2.5 Loss Function and Optimization

The FNODE loss directly supervises accelerations. Given predicted accelerations $\hat{\ddot{\mathbf{z}}}_{i}=f(\mathbf{z}_{i},\dot{\mathbf{z}}_{i};\mathbf{\Theta})$ and reference accelerations $\ddot{\mathbf{z}}_{i}$ obtained from trajectory data via FFT- or FD-based differentiation, we minimize the mean squared error (MSE)

[TABLE]

The optimal parameters are obtained as

[TABLE]

We train FNODE using the Adam optimizer [23]. The training algorithm is shown in A. A discussion of the training costs is deferred to Section 4.

3 Numerical Experiments

3.1 Single-Mass-Spring System

An undamped single mass–spring system serves as an energy-conserving benchmark for long-horizon prediction. Its dynamics are governed by

[TABLE]

where $x$ denotes the displacement from equilibrium, $m=10$ kg, and $k=50$ N/m. The system is initialized with $x(0)=1$ m and $\dot{x}(0)=0$, and simulated with step size $\Delta t=0.01$ s for 3000 steps (30 s). The first 300 steps are used for training and the remaining 2700 steps are reserved for extrapolation.

To highlight the system’s energy-conserving structure, we also express the dynamics in Hamiltonian form. Let $q=x$ be the generalized position and $p=m\dot{q}$ be the generalized momentum. The total energy (Hamiltonian) is

[TABLE]

which induces Hamilton’s equations

[TABLE]

We compare FNODE against structure-preserving baselines, including the Hamiltonian Neural Network (HNN) and Lagrangian Neural Network (LNN) [45]. HNN learns a scalar Hamiltonian $H_{\theta}(p,q)$ and obtains time derivatives through Hamilton’s equations, while LNN learns a Lagrangian $L_{\theta}(q,\dot{q})=T(\dot{q})-V(q)$ and predicts accelerations via the Euler–Lagrange equation. We also consider the Runge-Kutta 4th order (RK4) integrator as a baseline, which infers the future steps through the exact dynamics from Eq. (30) without learning.

Figure 2 compares long-term rollout behavior over the 30-second trajectory. In phase space, all methods remain close to the ground truth closed orbit, indicating that they capture the periodic motion. However, the energy trajectories highlight marked differences in stability. Specifically, structure-aware modeling and numerical integration, e.g., use of the Hamiltonian or Lagrangian structure or using a symplectic integrator (leap-frog in our implementation – FNODE_LF), mitigate energy drift, see Fig. 2. HNN and LNN maintain near-constant energy throughout the rollout. In contrast, the RK4 integrator shows tangible energy drift, both when used directly on the equations of motion or in conjunction with the FNODE model (see FNODE_RK4).

3.2 Single-Mass-Spring-Damper System

The equation of motion for the single-mass-spring-damper system is given by

[TABLE]

where $m=10$ kg, $k=50$ N/m, and $d=2$ Ns/m. The system is initialized with $x(0)=1$ m, $\dot{x}(0)=0$, and solved by a Runge-Kutta 4th order (RK4) formula at $\Delta t=0.01$ s. A trajectory of 300 steps was used for training, with an additional 100 steps reserved for extrapolation. The model hyperparameters are listed in Table 3.

Figure 3 compares predicted displacement and velocity. All methods fit well within the training interval, but differences emerge in extrapolation. FNODE maintains close agreement with the ground truth while MBD-NODE degrades more quickly. LSTM fails to capture the dissipative behavior, and FCNN diverges entirely. Phase-space trajectories in Figure 4 confirm FNODE’s stable long-term prediction, in contrast to the drift or stagnation observed in the baselines.

We also investigate the sources of rollout error in FNODE on the single-mass-spring-damper system. All errors are calculated over a fixed physical time horizon $T=2s$ by MSE. Extract 15 samples of time step $\Delta t\in[5\times 10^{-4},1\times 10^{-2}]$ and define the uniform grid $t_{n}=n\Delta t$ for $n=0,\dots,N$, where $N-1=\frac{T}{\Delta t}$. Define $Z_{n}=Z(t_{n})\in\mathbb{R}^{2d}$ as the accurate reference solution, and let $\widetilde{Z}_{n}=\widetilde{Z}(t_{n})\in\mathbb{R}^{2d}$ denote the system state (e.g., position and velocity). Let $f$ denote the ground-truth vector field $\ddot{z}=f(Z_{n})$, and let $f_{\Theta}$ denote the learned vector field $\mathring{\ddot{z}}=f_{\Theta}(Z_{n})$, see Eq. (8).

Figure 5 reports the decomposition of the acceleration MSE into boundary and interior contributions, together with the overall MSE. The log–log slope of each curve is also provided, and the ground truth is given by the analytical acceleration solution of the system. Consistent with the analysis in Section 2.4, the contribution with the lower FD order dominates the overall MSE. In particular, the boundary-point error (1st order FD) exhibits a slope close to $3$, matching the theoretical scaling $\mathcal{O}(\Delta t^{3})$ for the MSE of 1st order FD. The interior-point error has a slope close to $4$, consistent with the second-order central finite difference used for interior-point acceleration calculation. Accordingly, the overall MSE follows the scaling of the dominant term in the regime shown. The same figure also compares the MSE of the learned vector field $f_{\Theta}$ trained on the analytical acceleration solution with that of $f_{\Theta}$ trained on the finite difference acceleration result. In both cases, the learned acceleration error is significantly larger than the FD error, and the overall MSE is dominated by the learned-field mismatch. This suggests that, in FNODE, the acceleration MSE is dominated by the error in the learned vector field.

Figure 6 shows the log–log slope results for trajectory MSE of the RK4 integrator. The measured slope is close to $8$, which is consistent with the theoretical order of accuracy of RK4. This confirms that, when the learned vector field is accurate, the trajectory MSE is dominated by the integrator error.

3.3 Multiscale Triple-Mass-Spring-Damper System

We next consider a three-mass-spring-damper system with strong scale separation (Figure 7). The masses are $m_{1}=100$ kg, $m_{2}=10$ kg, and $m_{3}=1$ kg. All springs have stiffness $k_{i}=50$ N/m and all dampers $d_{i}=2$ Ns/m. The governing equations are

[TABLE]

The system was simulated with RK4 at a fixed step of $\Delta t=0.01$ s. Training trajectories consisted of 300 steps from the initial state $x_{1}=1$, $x_{2}=2$, $x_{3}=3$, $\dot{x}_{1}=\dot{x}_{2}=\dot{x}_{3}=0$ (SI units). Performance was then evaluated by extrapolation for an additional 100 steps. Model hyperparameters are listed in Table 4.

Figure 8 shows the trajectories of $x_{1},x_{2},x_{3}$ and $\dot{x}_{1},\dot{x}_{2},\dot{x}_{3}$ during training and testing. All models fit the training data, though FCNN exhibits minor spurious oscillations. In extrapolation ($t>3$), FNODE achieves the best agreement with ground truth ($\epsilon=7.4\text{e-}4$). MBD-NODE remains reasonable but drifts due to solver error accumulation. LSTM largely repeats training patterns and fails to generalize, while FCNN diverges under multiscale dynamics.

Figure 9 shows phase-space trajectories on the test set. FNODE tracks the dynamics well – the small discrepancies for the lighter masses are consistent with increased sensitivity to local inaccuracies in the learned acceleration field and their subsequent accumulation during rollout, as discussed in Section 2.4. MBD-NODE performs poorly for Mass 1, where low variance in the signal reduces its contribution to the loss. LSTM largely repeats training patterns and fails to generalize, while FCNN is inaccurate both in-distribution and out-of-distribution.

3.4 Double Pendulum

We next evaluate FNODE on a classic example of a chaotic system – the double pendulum. The system is described by the Hamiltonian formulation (see B for details). We used $L_{1}=L_{2}=1$ m, $m_{1}=m_{2}=1$ kg, and $g=9.81$ m/s2. The initial state was $\theta_{1}=3\pi/7$, $\theta_{2}=3\pi/4$, and $\dot{\theta}_{1}=\dot{\theta}_{2}=0$. Trajectories were generated using RK4 with $\Delta t=0.01$ s for 300 training steps and 100 extrapolation steps. The model hyperparameters are listed in Table 5.

Figure 11 shows the double pendulum trajectories. All models fit the training data, but diverge in extrapolation. FNODE, despite error accumulation in the ODE solver, achieves the best accuracy ($\epsilon=1.4\text{e-}1$).

Phase-space plots in Figure 12 confirm this trend: FNODE reproduces the main chaotic patterns, with only minor discrepancies. MBD-NODE loses accuracy in testing, LSTM largely repeats training patterns without generalization, and FCNN fails to capture the dynamics outside the training regime.

For the double pendulum, which exhibits chaotic behavior in the regime considered here, “long-term” pointwise trajectory accuracy must be interpreted relative to the system’s intrinsic predictability horizon. In chaotic dynamics, small initial-condition or modeling errors often grow approximately exponentially (on average) at a rate set by the largest Lyapunov exponent $\lambda_{\max}$, motivating the Lyapunov time scale $T_{L}:=\tfrac{1}{\lambda_{\max}}$. Beyond horizons on the order of a few $T_{L}$, pointwise trajectory comparisons become increasingly stringent because even very small discrepancies are rapidly amplified by the dynamics. The rollout horizons reported in this section should therefore be interpreted as finite-horizon predictions (short to intermediate horizons) rather than claims of sustained pointwise accuracy arbitrarily far into the chaotic regime. Estimating $\lambda_{\max}$ for the specific operating regime (and reporting horizons in normalized units $t/T_{L}$) is an important direction for future work. Finally, even when long-horizon prediction is inherently challenging, FNODE (and related learned dynamics models) can still be useful in downstream tasks that rely on short-horizon predictive accuracy, e.g., when the prediction horizon of interest is short relative to $T_{L}$.

3.5 Slider Crank with Coulomb friction

Next, FNODE’s performance is evaluated on a constrained slider-crank mechanism (Figure 13) under both frictionless and frictional conditions. The mechanism consists of three rigid bodies with a slider connected to a spring and damper, forming a constrained three-body problem. The generalized coordinates are $q=(x_{1},y_{1},\theta_{1},x_{2},y_{2},\theta_{2},x_{3},y_{3},\theta_{3})$; see C for details. The mechanism is defined by the following three bodies:

A crank (body 1) with mass $m_{1}=1$ kg and moment of inertia $I_{1}=0.1\text{ kg}\cdot\text{m}^{2}$ is connected to the ground by a revolute joint. Its length is $l=2$ m. 2. 2.

A connecting rod (body 2) with mass $m_{2}=1$ kg and moment of inertia $I_{2}=0.1\text{ kg}\cdot\text{m}^{2}$ is attached to the crank via a revolute joint. Its length is $r=1$ m. 3. 3.

A slider (body 3) with mass $m_{3}=1$ kg and moment of inertia $I_{3}=0.1\text{ kg}\cdot\text{m}^{2}$ is joined to the connecting rod by a revolute joint. It is also connected to a fixed wall via a spring with a torsional constant of $k=1$ Nm/rad. There is friction between the slider and the ground, which is modeled as Coulomb friction with a coefficient of $\mu$.

FNODE learns a mapping from the full 18-dimensional state space $(q,\dot{q})$ to accelerations $\ddot{q}$. Since the mechanism has only one DOF, the dynamics can be represented compactly with the minimal coordinates $(\theta_{1},\dot{\theta}_{1})$, with all other states dependent. Following the coordinate-partitioning approach [46], once FNODE produces the acceleration of the crank, numerical integration is used to compute the orientation and the angular velocity of the crank, while the kinematic constraint equations are used to recover the dependent generalized coordinates. The model hyperparameters are summarized in Table 6.

3.5.1 Non-friction Case

Figure 14 shows the slider-crank responses. All models fit the training data, but diverge in extrapolation. FNODE maintains accurate long-term predictions. MBD-NODE, while structurally similar to FNODE, is less stable and incurs higher computational cost due to its solver-based training. In contrast, LSTM and FCNN fail to generalize beyond the training window. Figures 15–18 illustrate the remaining state variables, which inherit errors from the minimal coordinates $(\theta_{1},\dot{\theta}_{1})$.

3.5.2 Friction Case

In this experiment, the friction coefficient $\mu$ is added as an input to the FNODE model to test the interpolation in the parameter space. The friction coefficient is sampled from $[0.0,0.6]$ uniformly for training data generation. For testing, we use interpolated friction coefficients to test the generalization ability of FNODE in parameter space.

Figure 19 shows FNODE-predicted phase-space trajectories of the slider-crank mechanism for different friction coefficients. The FNODE model was trained on coefficients sampled uniformly from $[0.0,0.6]$ and evaluated on interpolated values. This test case indicates that FNODE can learn the frictional dynamics and approximate the resulting non-smooth acceleration. To improve fidelity near these non-smooth transitions, one can increase FNODE model complexity (e.g., deeper/wider hidden layers) to better represent piecewise-smooth behavior.

A central challenge for FNODE modeling in the presence of friction is accurately predicting the resulting non-smooth acceleration field. A common way to improve accuracy in such settings is to increase network depth, which enhances the model’s ability to represent sharp, piecewise-smooth transitions. Figure 20 shows FNODE acceleration approximations for the slider-crank mechanism at $\mu=0.6$ using models with different numbers of hidden layers. As depth increases, the approximation consistently improves, as evidenced by the decreasing MSE values reported in Table 7.

3.6 Cart-Pole System

We conclude with the cart-pole system, considered here because it is a canonical control benchmark (Figure 21). The cart can move horizontally on a frictionless track, while the pendulum swings freely about its pivot, with motion restricted to $[-\pi/2,\pi/2]$. The system state is described by four variables: cart position $x$, velocity $\dot{x}$, pendulum angle $\theta$, and angular velocity $\dot{\theta}$. The dynamics are governed by nonlinear ODEs:

[TABLE]

The cart-pole was first tested in free evolution (no external input), initialized at

[TABLE]

Trajectories were generated with a midpoint integrator at $\Delta t=0.01$ s. The first 200 steps (2 s) were used for training and the following 50 steps (0.5 s) for testing. Model hyperparameters are listed in Table 8.

Figures 22-23 compare model performance. FNODE and MBD-NODE both match the ground truth closely, with low MSEs ($6.7\!\times\!10^{-4}$ and $1.1\!\times\!10^{-3}$, respectively). MBD-NODE is slightly more accurate, but FNODE remains competitive while training far more efficiently. LSTM largely reproduces training patterns, which fit the periodic pattern of cart-pole system, so it still provides relatively accurate predictions. FCNN diverges significantly in extrapolation.

We also test FNODE in a model predictive control (MPC) task, where an external force $u$ is applied to stabilize the pole and return the cart to the origin. The MPC strategy, for control purposes, linearizes the dynamics and then solves a quadratic program over a finite horizon [33]. For a linearized system $\dot{z}=Az+Bu$, the optimization problem is formulated as:

[TABLE]

At each step, the system state is $z_{k}=(\theta_{k},x_{k},\dot{\theta}_{k},\dot{x}_{k})^{T}$ with control input $u_{k}$. The MPC problem is solved over a horizon of length $N$, subject to state and input constraints $\mathcal{Z}$ and $\mathcal{U}$. The cost function weights are $Q$ and $R$, both chosen as identity matrices. $A$ and $B$ are the linearized system matrices, obtained from the Jacobian of the dynamics.

To solve the control problem, FNODE learns the mapping from state and control input to accelerations, $f(\theta_{k},x_{k},\dot{x}_{k},\dot{\theta}_{k},u)\mapsto(\ddot{\theta}_{k},\ddot{x}_{k})$. We use automatic differentiation and the inherent differentiable nature of a Neural Network to obtain the linearization matrices $A$ and $B$ required by the MPC framework. The standard (non-NN based) MPC for Cart-Pole is also provided as a baseline for comparison, see D for the expression of the linearization matrices.

Figure 24 shows the MPC results for the cart-pole. The five subplots depict the pole angle $\theta$, cart position $x$, angular velocity $\dot{\theta}$, cart velocity $\dot{x}$, and control input $u$. The FNODE-based controller (red dotted) closely matches the analytical baseline (blue solid), demonstrating that FNODE learns the dynamics with sufficient accuracy for supporting an MPC controller.

3.7 Parameterized Vehicle Model

We finally consider the simplified 4DOF vehicle model proposed in [51, 42]. The model has a simple powertrain component and adopts throttle and steering as control inputs. The state and control vectors are:

[TABLE]

where $(x,y)$ is the vehicle position in the global frame, $\theta$ is the heading (yaw) angle, $v$ is the longitudinal speed, $\alpha$ is the normalized throttle, and $\delta$ is the steering angle.

The state transition model is

[TABLE]

where $l$ is the wheelbase, $\beta$ maps steering command $\delta$ to a physical wheel steering angle, $\gamma$ is the gear ratio, and $I_{\mathrm{wheel}}$ and $R_{\mathrm{wheel}}$ are the wheel inertia and radius. The motor torque is modeled as

[TABLE]

where $\tau_{0}$ and $\omega_{0}$ are the stalling torque and no-load speed of the motor, and $c_{0},c_{1}$ model the motor resistance torque.

This example contains nine scalar parameters

[TABLE]

each sampled on a grid with four values, yielding $4^{9}$ parameter combinations. We simulate each parameter instance for 20 s with time step $\Delta t=0.01$ (2000 steps). The control signal $u(t)=[\alpha(t),\delta(t)]^{T}$ is pre-defined and shared by all parameter instances; the corresponding state trajectories are generated by the RK4 integrator. For each trajectory, the first 1500 steps are used for training and the remaining 500 steps for testing (extrapolation).

FNODE learns the mapping from the full state and control inputs to the state time derivative, i.e., $(q,u,\Theta_{p})\mapsto\dot{q}$ as defined in Eqs. (38)–(40). All models are trained on the same trajectories and evaluated on the held-out extrapolation interval. The hyperparameters are summarized in Table 11.

Figure 25 compares the predicted trajectories of $q$ against the ground truth under identical control inputs. Tests 1-2 are drawn from the corners of the parameter grid, whereas tests 3-4 lie in the interior (Table 10) and therefore assess interpolation in parameter space. All methods track well within the training interval; discrepancies arise under extrapolation, where integration errors accumulate and parameter variations compound. FNODE remains stable across the grid and delivers comparatively accurate predictions over the test window, while the baselines exhibit larger drift and occasional divergence. Because the system’s high dimensionality makes the baseline datasets prohibitively large, we report only FNODE results here. Future work will investigate scalable baseline implementations for high-dimensional parameter spaces.

4 Discussion

Figure 26 summarizes the computational cost and predictive accuracy of the four models for the triple mass-spring-damper benchmark. Two observations emerge.

First, the training process of FNODE converges with lower computational cost than MBD-NODE. The computational cost is quantified by the training time of the model. As shown in Figure 26, the train loss of FNODE converges with lower training time compared to MBD-NODE. This gap arises because FNODE bypasses the adjoint sensitivity method and does not embed an ODE solver in the backpropagation loop. Compared to purely data-driven baselines such as FCNN, FNODE’s training times are of the same order of magnitude.

Second, FNODE achieves consistently lower prediction error than the black-box baselines and is competitive with, or superior to, MBD-NODE. In dissipative and chaotic systems (single/triple mass-spring-damper, double pendulum), FNODE outperforms all baselines by large margins, maintaining stable long-term predictions. For the slider-crank mechanism, FNODE achieves the best accuracy ($1.8\!\times\!10^{-3}$ MSE), significantly outperforming both MBD-NODE and the data-driven models. Only in the cart-pole case does MBD-NODE obtain a slightly lower error than FNODE. These results underscore the importance of accurate acceleration targets: FFT provides precision for smooth trajectories, while FD ensures robustness for non-periodic signals and trajectory endpoints.

FNODE provides a good efficiency vs. accuracy balance: it is orders of magnitude faster to train than MBD-NODE while delivering consistently higher fidelity than LSTM and FCNN, and near-parity or better accuracy than MBD-NODE across benchmarks. These results suggest that FNODE is computationally practical for higher-dimensional constrained multibody systems and can accommodate parameter variation, as illustrated by the added higher-DOF vehicle example and parameter-conditioned experiments.

In its current form, FNODE is best suited to systems with smooth or piecewise-smooth trajectories, where accelerations are well-defined and can be estimated reliably from sampled data (here via FFT/FD differentiation). Our results indicate that, despite the resulting non-smoothness, a Coulomb friction model remains amenable to FNODE, which can learn an accurate approximation of the associated acceleration field. In contrast, genuinely nonsmooth events, e.g., impacts (chattering) and complementarity-based stick-slip transitions, are deferred to future work, as discussed in the Conclusions.

5 Conclusions and Future Work

We introduced FNODE, a framework for data-driven modeling of multibody dynamics that learns accelerations directly rather than states. By reformulating the training objective and leveraging FFT- and finite-difference-based differentiation, FNODE eliminates the adjoint bottleneck of traditional Neural ODEs. Across benchmarks ranging from simple oscillators to chaotic and constrained mechanisms, FNODE trained one to two orders of magnitude faster than MBD-NODE while outperforming black-box baselines such as LSTM and FCNN. Key outcomes include:

• •

Good accuracy and generalization: FNODE achieved lower MSE than all baselines, performing well in both interpolation and extrapolation.

• •

Data efficiency: Accurate long-term predictions were obtained from short training trajectories, e.g., thousands of steps extrapolated in the slider-crank test.

• •

Multiscale robustness: FNODE captured dynamics across time scales, as shown in the triple-mass-spring-damper system.

• •

Reproducibility: Our open-source implementation and benchmark suite, which includes all NNs and dynamics models used in this study, are provided for reproducibility studies and further research.

FFT-based spectral differentiation yields high accuracy for smooth, uniformly sampled trajectories when periodicity (or approximate periodicity) is a reasonable assumption, since the derivative is obtained from a global Fourier representation. However, for finite, non-periodic trajectories the implicit periodic extension can introduce endpoint artifacts (spectral leakage/Gibbs-type effects), motivating the use of local finite-difference (FD) stencils near boundaries. The FD schemes employed in this work were of lower order (e.g., $\mathcal{O}(\Delta t^{2})$ in the interior and $\mathcal{O}(\Delta t)$ at endpoints for the stencils used here), but are robust for truncated data and do not rely on periodicity. Higher order FD, which have not been tested in this work, could be more accurate but may amplify noise. By employing the FFT and FD schemes in a complementary way, i.e., FFT where its assumptions are approximately satisfied and FD where boundary effects or non-periodicity dominate, we balance accuracy and robustness without relying on either method exclusively.

In the MSE comparison test in Section 3.1, we compared the MSE of Neural Network learned field mismatch and the error of numerical differentiation in the same setting. The result shows that the error of Neural Network learned field mismatch is dominated by the error of numerical differentiation. This suggests that with accurate acceleration calculated by numerical differentiation, FNODE is able to achieve good accuracy without numerical integrator in training process.

There are several limitations of the approach which we stopped short of addressing and thus are deferred to future work. First, FNODE relies on accurate acceleration estimates from trajectory data; both FFT and FD methods can suffer from noise and, in the case of FFT, Gibbs artifacts. Further work should be carried out to investigate the sensitivity of the approach to the choice of differentiation scheme. Second, integration error compounds over time, reducing long-term accuracy in chaotic systems such as the double pendulum. However, even when this is the case, the proposed approach still yields accurate predictions for short-term and intermediate-term horizons, a key requirement for controls applications.

Future work should investigate also noise-resilient methods for estimating accelerations from sparse or imperfect data. The goal is to reduce long-term error growth, for example by embedding physical constraints or conservation laws into training or integration, in line with the methodology embraced for the slider-crank example, which has kinematic constraints. An important direction for future work is a systematic numerical study of differentiation-induced errors, including evaluation of the stencil-disagreement indicator $\eta^{\prime}$ and the trajectory reconstruction residual $\eta_{z}$ across datasets and model variants. Finally, one rich direction of research is extending FNODE to problems with friction and contact, in the vein of the slider-crank example. Handling friction is expected to be relatively straightforward for smoothed, penalty-based formulations in which friction is represented through history-dependent ‘bristle’ dynamics, as in models such as LuGre and Dahl. The problem is more challenging for complementarity-based or Lagrange-multiplier formulations, where stick-slip transitions are enforced non-smoothly and can introduce abrupt regime changes. Beyond multibody dynamics, several avenues for future work include applying FNODE ideas to PDE-governed systems, integrating uncertainty quantification, and exploring closed-loop control applications.

The code and scripts used for this work are publicly available as open source for reproducibility studies and further research [44].

Acknowledgements

This work was carried out in part with support from National Science Foundation project CMMI2153855.

Declaration of Generative AI and AI-assisted technologies in the writing process

During the preparation of this work the authors used OpenAI ChatGPT in order to improve language and readability of the text. After using this service, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Algorithms for training the FNODE

Appendix B Double Pendulum Hamiltonian Formulation

For completeness, we outline the Hamiltonian derivation of the double pendulum dynamics used in Section 3. The Hamiltonian is obtained via a Legendre transform of the Lagrangian and is written as

[TABLE]

with kinetic and potential energy

[TABLE]

The canonical momenta are defined by $p_{\theta_{i}}=\tfrac{\partial T}{\partial\dot{\theta}_{i}}$. Using the expression for $T$ above yields

[TABLE]

Inverting these relations eliminates $\dot{\theta}_{1}$ and $\dot{\theta}_{2}$ in favor of $p_{\theta_{1}}$ and $p_{\theta_{2}}$, producing $H(\theta_{1},\theta_{2},p_{\theta_{1}},p_{\theta_{2}})$ and the explicit Hamiltonian equations below.

From Hamilton’s equations

[TABLE]

we obtain the explicit form of the governing equations:

[TABLE]

with auxiliary terms

[TABLE]

Appendix C A slider-crank mechanism, with rigid bodies

Based on the setup illustrated in Section 3.5, we formulate the equation of motion of slider-crank mechanism:

The mass matrix $M\in\mathbb{R}^{9\times 9}$ is:

[TABLE]

where:

[TABLE]

[TABLE]

[TABLE]

The states of the slider crank mechanism $(x_{1},y_{1},\theta_{1},x_{2},y_{2},\theta_{2},x_{3},y_{3},\theta_{3})$ follow the below constraints $\Phi:\mathbb{R}^{9}\rightarrow\mathbb{R}^{8}$ on the position:

[TABLE]

Then the constraints Jacobian assumes the following form:

[TABLE]

The external forces vector $F_{e}\in\mathbb{R}^{9}$ is formulated as:

[TABLE]

where:

[TABLE]

where $|\lambda_{7}|$ is the absolute value of the scalar Lagrange multiplier (reaction force) associated with the constraint $y_{3}=0$ in Eq. (53), and $\mu$ is the friction coefficient.

The rearranged constraint equations on the acceleration can be formulated from (53):

[TABLE]

Then we can get the $\gamma_{c}$ as follows:

[TABLE]

Appendix D Analytical Linearization for Cart-Pole

Based on the setup in Section 21, we formulate the cart-pole as a nonlinear continuous system:

[TABLE]

where $z=[\theta,x,\dot{\theta},\dot{x}]^{T}$ is the system state, $u$ is the control input, $\dot{z}=[\dot{\theta},\dot{x},\ddot{\theta},\ddot{x}]^{T}$ is the time derivative of $z$, and $f$ is the nonlinear function that defines the system dynamics.

To implement the MPC, we linearize the nonlinear system around the equilibrium point $(z_{0},u_{0})$ by first-order Taylor expansion:

[TABLE]

In equilibrium state, $z_{0}=0,u_{0}=0,f(z_{0},u_{0})=0$, and the deviations from equilibrium are $\Delta z=z-z_{0}=z$ and $\Delta u=u-u_{0}=u$. The linearized system in Eq. (60) can be expressed as:

[TABLE]

Then we define the linearized system matrices $A$ and $B$, where

[TABLE]

As it can be assumed that the pole remains near $\theta=0$, we can linearize the system around this equilibrium point to obtain a linearized model for MPC design.

[TABLE]

Then substitute Eq. (63) into the nonlinear system dynamics Eq. (35) to obtain the linearized model.

[TABLE]

By rearranging Eq. (64), the accelerations $\ddot{\theta}$ and $\ddot{x}$ can be expressed as:

[TABLE]

Thus, the linearized system matrices $A$ and $B$ can be derived as:

[TABLE]