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

1. Overview

Data-driven modeling of constrained multibody dynamics faces two persistent obstacles: (i) the training cost of Neural ODEs, which typically backpropagate through an ODE solver, and (ii) error accumulation over long rollout predictions. We introduce FNODE, a Flow-Matching Neural ODE framework that learns the acceleration mapping directly from trajectory data. Because supervision targets accelerations rather than integrated states, training reduces to plain supervised regression — the ODE-adjoint/solver backpropagation bottleneck disappears. Acceleration targets are obtained efficiently by 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, and 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, comparing against MBD-NODE, LSTM, and fully connected baselines. On the single mass–spring–damper we also decompose the rollout error into the contributions of numerical differentiation and of the learned vector field: the network’s approximation error dominates, confirming that numerical differentiation of the training targets is robust. Across all benchmarks, FNODE improves both prediction accuracy and training and runtime efficiency, while the partitioning procedure maintains constraint satisfaction. Our code and scripts are released as open source to support reproducibility and follow-on research.

Keywords Multibody dynamics · Neural ODE · Constrained dynamics · Friction dynamics · Scientific machine learning · Flow-Matching

Contributions

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.
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.
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.
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. Method

2.1 Supervise accelerations, not states

A conventional Neural ODE obtains future states by solving an initial value problem, so every training step must differentiate through — or adjoint-solve alongside — the integrator. FNODE instead models the acceleration directly on the augmented state Z = (z^⊤, ż^⊤)^⊤: z̈(t, μ) = f_Θ(Z(t), μ), with μ an optional parameter input (e.g., a friction coefficient). Training is plain regression on the vector field,

ℒ(Θ) = (1/N) ∑_i ‖ f(z_i, ż_i; Θ) − z̈_i ‖₂²,

with no ODE solve anywhere in the loop. The solver enters only at inference — and can be swapped freely: adaptive steppers for stiff phases, or symplectic leapfrog when energy behavior matters.

2.2 Hybrid FFT–FD acceleration targets

For smooth, near-periodic signals, spectral differentiation is highly accurate. A finite sampling window, however, breaks periodicity, and the resulting Gibbs error scales as ∼|Δ_boundary|/(2Δt) — independent of sample count, and growing as the step shrinks. The FFT pipeline therefore detrends, mirror-reflects with a cosine taper, applies a Tukey window (α = 0.2), and low-pass filters in the frequency domain (σ ≈ N/20).

Non-periodic data use finite differences: three-point central stencils in the interior (O(Δt²)) and one-sided stencils at endpoints (O(Δt)) — deliberately conservative, since wider stencils amplify noise. A Huber loss can absorb occasional label spikes without changing the architecture.

2.3 Constraints by coordinate partitioning

Constrained mechanisms evolve on the manifold of their kinematic constraint equations; learning accelerations for all redundant coordinates would let rollouts drift off it. Following the generalized coordinate partitioning of Wehage and Haug, FNODE learns the acceleration field only for the independent generalized coordinates; the dependent coordinates are recovered at every step by solving the position-level constraint equations — constraint satisfaction by construction.

Concretely, the slider-crank’s 18-dimensional redundant state collapses to the two learned quantities (θ₁, θ̇₁); the kinematics return everything else.

2.4 An error budget that behaves

Rollout error separates into integrator truncation, learned-field mismatch, and differentiation-induced label error. An empirical decomposition on the single mass–spring–damper (15 step sizes, Δt ∈ [5×10⁻⁴, 10⁻²]) shows that the learned-field mismatch dominates — numerical differentiation of the training labels is not the bottleneck, whether the model is trained on analytical or finite-difference targets.

Every stage converges at its theoretical order: measured log–log slopes ≈ 3 for boundary FD, ≈ 4 for interior FD, and ≈ 8 for the RK4 trajectory MSE (O(Δt^2p), p = 4).

3. Benchmarks and results

FNODE is evaluated on six mechanical systems of increasing difficulty against MBD-NODE, LSTM, and FCNN baselines — and, on the conservative oscillator, against Hamiltonian and Lagrangian Neural Networks as well. Figure and table numbers follow the published article; click any figure to enlarge it.

Rendered schematic of three masses connected in series to a wall by springs and dampers. — Fig. 7. Triple mass-spring-damper system. The setup is similar to the single mass-spring-damper system, except for the addition of two more masses, springs, and dampers. The system admits an analytical solution that was used for accuracy studies.

Schematic of a double pendulum with link lengths l1, l2 and angles θ1, θ2 measured from the vertical. — Fig. 10. Schematic of the double pendulum.

Rendered schematic of a slider-crank mechanism: crank, connecting rod, and slider with spring, damper, and friction force. — Fig. 13. Slider-crank mechanism with motor torque `τ`₁, rotational dampers at all joints, and slider friction. The friction force is shown in red, see Appendix C, Eq. (C.5c), for its expression. It acts between body 3 and ground. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Line-drawing schematic of a cart-pole: cart of mass M under force F carrying a pole of length l at angle θ. — Fig. 21. Cart-pole system.

Four panels comparing ground truth, HNN, LNN, FNODE with RK4, FNODE with leapfrog, and RK4: phase portrait, total energy versus time, and position and velocity time series over 30 seconds. — Fig. 2. Phase-space trajectories and energy evolution for the single mass–spring system. The left panel shows the phase portrait (`v` vs. `x`), while the right panel shows the total energy. The energy drift ratios are: FNODE_RK4 (15.70%), FNODE_LF (0.46%), RK4 (13.39%), HNN (0.98%), and LNN (−0.24%). The mse for each model is provided: FNODE_RK4 (`ε` = 7.50e-2), FNODE_LF (`ε` = 9.66e-5), RK4 (`ε` = 2.02e-3), HNN (`ε` = 5.51e-2), and LNN (`ε` = 2.27e-4). The ground truth is provided by the analytical solution of the system.

Energy behavior over long horizons

The model is trained on 3 s of the undamped oscillator and rolled out for 30 s. Because the learned field is independent of the solver, it can be paired with any integrator after training: under a symplectic leapfrog scheme (FNODE_LF) it holds the energy drift to 0.46% and reaches the lowest MSE of all compared methods. It matches the structure-preserving HNN and LNN — in 450 training epochs where HNN needs 30,000. The same field under RK4 drifts like the classical RK4 baseline: the drift is an integrator property, not a modeling failure.

Model	MSE	Energy drift
FNODE_LF	9.66×10⁻⁵	0.46%
LNN	2.27×10⁻⁴	−0.24%
RK4 (exact EOM)	2.02×10⁻³	13.39%
HNN	5.51×10⁻²	0.98%
FNODE_RK4	7.50×10⁻²	15.70%

Eight panels of displacement and velocity versus time for the single mass-spring-damper under FNODE, MBD-NODE, FCNN, and LSTM: only FNODE stays on the ground truth beyond the training window. — Fig. 3. Temporal evolution of state variables for the single-mass-spring-damper system. ID stands for “in distribution” and OOD stands for “out of distribution”. The left and right columns depict the displacement (`x`) and velocity (`ẋ`), respectively. Dashed lines indicate model predictions over the training interval (`t` ∈ [0, 3]), while dotted lines represent extrapolated predictions on the test interval. The mean squared error (MSE) for each model is provided: FNODE (`ε` = 7.6e-5), MBD-NODE (`ε` = 1.4e-3), LSTM (`ε` = 1.5e-2), and FCNN (`ε` = 1.6e-1). The ground truth is provided by the analytical solution of the system.

Dissipative dynamics

Every model fits the 3-second training window of the damped oscillator; only FNODE stays on the analytical solution beyond it. LSTM fails to capture the dissipation, and FCNN diverges entirely. FNODE also trains in roughly half MBD-NODE’s wall-clock time (112 s vs. 202 s).

Model	Rollout MSE
FNODE	7.6×10⁻⁵
MBD-NODE	1.4×10⁻³
LSTM	1.5×10⁻²
FCNN	1.6×10⁻¹

Displacement and velocity of all three masses in the triple mass-spring-damper for four models: FNODE stays on ground truth in the test interval while the baselines fan away. — Fig. 8. Predicted temporal evolution for the triple-mass-spring-damper benchmark. The left column presents the displacement (`x`) and the right column presents the velocity (`ẋ`) for each of the three masses. Dashed and dotted lines denote predictions on the training (`t` ∈ [0, 3]) and test sets, respectively. The associated MSEs are: FNODE (`ε` = 6.1e-3), MBD-NODE (`ε` = 5.4e-2), LSTM (`ε` = 2.5e-1), and FCNN (`ε` = 8.9e-2).

Multiscale dynamics

Masses of 100, 10, and 1 kg span two orders of magnitude, so the light mass oscillates rapidly while the heavy mass responds slowly. FNODE preserves accuracy across the separated time scales. MBD-NODE accumulates solver error and struggles on the slow, low-variance heavy mass; LSTM replays its training patterns; FCNN diverges.

Model	MSE (Table 13)
FNODE	6.1×10⁻³
MBD-NODE	6.8×10⁻³
FCNN	9.0×10⁻²
LSTM	2.5×10⁻¹

Angular position and velocity of the double pendulum for four models: FNODE follows the ground truth longest in extrapolation, LSTM drifts out of phase, FCNN flatlines. — Fig. 11. Dynamic response predictions for the double pendulum system. The left and right columns show the temporal evolution of angular position (`θ`) and angular velocity (`θ̇`), respectively. Dashed lines represent the fit to the training data (`t` ∈ [0, 3]), while dotted lines show the extrapolation on test data. The reported MSE values are: FNODE (`ε` = 1.4e-1), MBD-NODE (`ε` = 2.3e-1), LSTM (`ε` = 7.6e-1), and FCNN (`ε` = 1.9e0).

A chaotic system

On the double pendulum, every model fits the training window — and every model eventually diverges, as any pointwise prediction must once errors amplify at the rate set by the largest Lyapunov exponent. FNODE tracks the true trajectory longest; FCNN’s rollout collapses to a constant. The paper frames these results explicitly as finite-horizon predictions relative to the Lyapunov time — the regime that matters for control applications.

Model	Rollout MSE
FNODE	1.4×10⁻¹
MBD-NODE	2.3×10⁻¹
LSTM	7.6×10⁻¹
FCNN	1.9

Crank angle and angular velocity over a 45-second rollout for four models: FNODE stays locked to the damped oscillation while baselines drift or collapse. — Fig. 14. Long-term prediction results for the frictionless slider-crank’s independent coordinate (`θ`₁). The left column presents angular position versus time; the right column presents angular velocity versus time. Model performance on the training interval (`t` ∈ [0, 15]) is shown with dashed lines, while extrapolated performance on the test interval is shown with dotted lines. The corresponding MSEs are: FNODE (`ε` = 1.7e-1), MBD-NODE (`ε` = 8.7e-1), LSTM (`ε` = 1.4e0), and FCNN (`ε` = 3.1e0).

Constrained dynamics via partitioning

The central constrained benchmark: nine generalized coordinates, eight position constraints, one degree of freedom. FNODE learns only the crank acceleration; trained on 15 s and rolled out to 45 s, it holds amplitude and phase throughout while the baselines drift or flatline. Solving the position-level constraint equations then reconstructs every dependent position and velocity, so the constraints are satisfied by construction.

Model	Rollout MSE
FNODE	1.7×10⁻¹
MBD-NODE	8.7×10⁻¹
LSTM	1.4
FCNN	3.1

Reconstructed dependent states (Fig. 15)

Twelve panels of the slider-crank's dependent coordinates and velocities reconstructed from FNODE's predicted minimal coordinates, tracking ground truth for 45 seconds. — Fig. 15. Temporal evolution of the dependent state variables for the slider-crank mechanism, as reconstructed from the FNODE model’s prediction of the minimal coordinates.

Slider-crank trajectories at four friction coefficients: FNODE's dashed predictions overlay the ground truth for angle, angular velocity, phase portraits, and accelerations. — Fig. 19. Trajectories of the slider-crank mechanism under varying friction coefficients, as predicted by the FNODE model. The model was trained on friction coefficients sampled uniformly from [0.0, 0.6] and tested on interpolated values. In this figure, `c` = 0.0 and `c` = 0.6 are from the training set, while `c` = 0.15, 0.35 are interpolated test cases.

Generalizing across a non-smooth, parameterized field

With Coulomb friction at the slider, the friction coefficient enters FNODE as an extra input. A single model trained on trajectories with μ sampled uniformly from [0, 0.6] accurately predicts angle, angular velocity, phase portraits, and the discontinuous acceleration jumps at frictional transitions — at coefficients never seen in training. Capturing this piecewise-smooth field benefits from network depth (Table 7):

Hidden layers	Acceleration MSE (`μ` = 0.6)
2	3.70×10⁻²
3	1.72×10⁻²
4	1.47×10⁻²
5	1.44×10⁻²
6	1.41×10⁻²

4. Accuracy and efficiency

Across the five shared benchmarks, FNODE posts the lowest rollout MSE in the paper’s summary table — on most systems by an order of magnitude or more against the discrete-time baselines. Its training cost, meanwhile, stays on par with a plain FCNN, because no ODE solver sits in the backpropagation loop.

Table 13. Summary of numerical MSE (lower is better).
Test case	FNODE	MBD-NODE	LSTM	FCNN
Single mass–spring–damper	7.6×10⁻⁵	1.4×10⁻³	1.5×10⁻²	1.6×10⁻¹
Triple mass–spring–damper	6.1×10⁻³	6.8×10⁻³	2.5×10⁻¹	9.0×10⁻²
Double pendulum	1.4×10⁻¹	2.3×10⁻¹	7.6×10⁻¹	1.9
Slider-crank	1.7×10⁻¹	8.7×10⁻¹	1.4	3.1
Cart-pole	6.7×10⁻⁴	1.1×10⁻³	1.6×10⁻¹	3.5×10⁻¹

Table 12 (condensed). Wall-clock training time in seconds, Nvidia GeForce 5070 Ti (16 GB). MBD-NODE embeds an RK4 solver in its training loop; the vehicle baselines were not trained (dataset size prohibitive).
Test case	FNODE	MBD-NODE	LSTM	FCNN
Single mass–spring–damper	112.46	201.90	171.99	103.28
Triple mass–spring–damper	108.23	244.80	172.05	102.78
Double pendulum	89.33	212.09	148.41	88.04
Slider-crank	565.01	1185.64	813.11	453.50
Cart-pole	59.46	144.51	96.08	59.36
4-DOF vehicle	3304.32	–	–	–

The accuracy–compute trade-off is summarized on the triple mass–spring–damper. In the test-loss-versus-training-time protocol of Fig. 26, FNODE’s error drops below 10⁻³ within roughly 50 s of cumulative training and settles at 7×10⁻⁴. MBD-NODE trains for over 220 s to reach its final 6.8×10⁻³, and the discrete-time baselines plateau roughly two orders of magnitude above FNODE’s settled error. (The rollout MSEs of Table 13 are computed under a separate evaluation protocol.)

Training setup common to the benchmarks

All models use 256 nodes per hidden layer, MSE loss, the Adam optimizer, and exponential learning-rate decay (0.98; 0.99 for FNODE on the slider crank); FNODE, MBD-NODE, and FCNN are fully connected with tanh activations, while the LSTM uses its sigmoid/tanh gating. Depth is tuned per system: two hidden layers for the mass–spring(–damper), double-pendulum, and cart-pole benchmarks, three for the slider crank, six for the 4-DOF vehicle (trained with AdamW at an initial learning rate of 10⁻⁴). FNODE and FCNN train for up to 450 epochs, MBD-NODE for 300–400, LSTM for 300; on the conservative single mass–spring benchmark, HNN requires 30,000 epochs. Complete per-system hyperparameter tables are given in the paper (Tables 2–6, 8, and 11).

Log-scale plot of test MSE versus cumulative training time for four models: FNODE converges fastest and to the lowest error; MBD-NODE, FCNN, and LSTM plateau at successively higher errors. — Fig. 26. The test MSE vs training time comparison of the all models on the triple mass-spring-damper benchmark.

5. Code and reproducibility

The complete implementation is released as open source in the Simulation-Based Engineering Lab reproducibility repository. It includes the FNODE training and evaluation pipelines; the MBD-NODE, FCNN, and LSTM baselines and the HNN/LNN references; ground-truth data generators and integrators for all six benchmark systems; and the friction-parameterized slider-crank, MPC cart-pole, and parameterized 4-DOF vehicle studies:

github.com/uwsbel/sbel-reproducibility / 2025/FNODE

6. Citation

Hongyu Wang, Jingquan Wang, Dan Negrut, “FNODE: Flow-matching for data-driven simulation of constrained multibody systems,” Computer Methods in Applied Mechanics and Engineering, Volume 455, 2026, 118912, ISSN 0045-7825. doi:10.1016/j.cma.2026.118912

A preprint is available as arXiv:2509.00183. When using the code, please also cite the repository: H. Wang, J. Wang, D. Negrut, Models, scripts, and meta-data: FNODE, github.com/uwsbel/sbel-reproducibility, 2025.

@article{Wang2026FNODE,
  title   = {{FNODE}: Flow-matching for data-driven simulation
             of constrained multibody systems},
  author  = {Hongyu Wang and Jingquan Wang and Dan Negrut},
  journal = {Computer Methods in Applied Mechanics and Engineering},
  volume  = {455},
  pages   = {118912},
  year    = {2026},
  issn    = {0045-7825},
  doi     = {10.1016/j.cma.2026.118912}
}