Phase space integrity in neural network models of Hamiltonian dynamics: A Lagrangian descriptor approach
Authors: Abrari Noor Hasmi, Haralampos Hatzikirou, Hadi Susanto
Organizations: Department of Mathematics, Khalifa University of Science and Technology, PO Box 127788, Abu Dhabi, United Arab Emirates · Public Health and Epidemeology Department, Khalifa University of Science and Technology, PO Box 127788, Abu Dhabi, United Arab Emirates.
We propose Lagrangian Descriptors (LDs) as a diagnostic framework for evaluating neural network models of Hamiltonian systems beyond conventional trajectory-based metrics. Standard error measures quantify short-term predictive accuracy but provide little insight into global geometric structures such as orbits and separatrices. Existing evaluation tools in dissipative systems are inadequate for Hamiltonian dynamics due to fundamental differences in the systems. By constructing probability density functions weighted by LD values, we embed geometric information into a statistical framework suitable for information-theoretic comparison. We benchmark physically constrained architectures (SympNet, HénonNet, Generalized Hamiltonian Neural Networks) against data-driven Reservoir Computing across two canonical systems. For the Duffing oscillator, all models recover the homoclinic orbit geometry with modest data requirements, though their accuracy near critical structures varies. For the three-mode nonlinear Schrödinger equation, however, clear differences emerge: symplectic architectures preserve energy but distort phase-space topology, while Reservoir Computing, despite lacking explicit physical constraints, reproduces the homoclinic structure with high fidelity. These results demonstrate the value of LD-based diagnostics for assessing not only predictive performance but also the global dynamical integrity of learned Hamiltonian models.
Figures & tables
Figure 1: LDs and LD-weighted PDF for the Duffing oscillator. Top panels: 3D visualization of LD fields (left) and the corresponding PDF (right). Bottom panels: 2D projections of the same data. Singular features in the LD align with separatrices, while the weighted PDF highlights these regions as dynamically dominant. The LD functional in Eq. ( 2 ) is used with c=0.7 and an integration time of τ=4 . The PDF employs the weighting function g(x)=1/x .
NN type
Learned Hamiltonians
Layers
Neurons per Layer
Trainable parameters
#FLOPs (Inference)
SympNet
10
1
50
3000
6000
HénonNet
10
1
50
755
6020
GHNN
3
2
15
1710
6120
RC
–
–
400
800
6318
Table 1: NN configurations and FLOP counts for the Duffing equation.
Figure 2: Training (left) and validation (right) loss for the NN models in the Duffing case. SympNet (red dotted), HénonNet (green dash–dotted), GHNN (purple dash–dotted), and Reservoir Computing (RC, orange dashed). Note that RC training does not involve iterative optimization; the output weights are obtained by solving a regularized normal equation.
Figure 3: Forward and backward predictions of the Duffing equation starting from t0=0 . Top row: (q(0),p(0))=(−0.71,−0.51) (inside homoclinic orbit); Middle row: (q(0),p(0))=(−0.71,−0.61) (near homoclinic orbit); Bottom row: (q(0),p(0))=(−0.1,−0.71) (outside homoclinic orbit). From left to right, the columns show: the pointwise error in position q(t)−q^(t) , the pointwise error in momentum p(t)−p^(t) , the normalized prediction error (Eq. ( 26 )) as function of the integration time τ , and normalized LD (Eq. ( 27 )) as function of the integration time τ . Predictions are shown over the interval t∈[−10,10] . The reference solution (blue solid) is compared with SympNet (red dotted), HénonNet (green dash-dotted), GHNN (purple dash-dotted), and RC (orange dashed). All NNs were trained on the dataset corresponding to the 200-trajectory case.
Figure 4: Same as Fig. 3 , but with p(0)>0 . Top row: (q(0),p(0))=(−0.71,0.51) (inside homoclinic orbit); Middle row: (q(0),p(0))=(−0.71,0.61) (near homoclinic orbit); Bottom row: (q(0),p(0))=(−0.1,−0.71) (outside homoclinic orbit).
Figure 5: The prediction error e^ (left) and normalized prediction error e (right) of Duffing trajectories across initial conditions : (A-B) SympNet, (C-D) HénonNet, (E-F) GHNN, (G-H) RC. The white dots indicate the initial conditions of the trajectories shown in Figs. 3 and 4 .
Figure 6: Distribution of initial conditions q(0) in the training dataset for the 200-trajectory case. The red vertical lines indicate q(0)=±2 , corresponding to the intersections of the homoclinic orbit with the line p=0
Figure 7: LD-weighted PDFs for the Duffing equation. The left panels: (A) Numerical integration, (B) SympNet, (C) HénonNet, (D) GHNN, (E) RC. All NNs were trained on the dataset corresponding to the 200-trajectory case. High-density regions correspond to slow dynamics, revealing the fixed points near (q,p)=(±1,0) . The white lines indicate the phase-space cross sections shown in the right panels. Right panels: one-dimensional slices of the LD-weighted PDF across all methods at selected phase-space locations. From top to bottom: q=−1.0 , q=−0.7 , p=−0.4 , p=0.0 , and p=0.4 .
Figure 8: Comparison of homoclinic orbits predicted by NNs against the exact Duffing solution (Eq. ( 23 ), blue dashed). SympNet (red dotted), HénonNet (green dash-dotted), GHNN (purple dash-dotted), and RC (brown dashed). LDs reveal deviations in the homoclinic structure: HénonNet’s predicted orbit diverges significantly from the reference solution, while the other NNs align more closely. All models were trained on the dataset corresponding to the 200-trajectory case.
NN models
Number of training trajectories N
10
20
50
100
200
300
400
500
SympNet
2.17e-01
2.12e-01
2.19e-01
2.13e-01
1.32e-03
3.24e-04
1.01e-03
6.87e-04
HénonNet
1.76e-01
1.80e-01
1.84e-02
1.07e-02
8.03e-03
4.64e-03
4.16e-04
8.04e-03
GHNN
2.17e-01
7.11e-03
6.30e-04
1.02e-03
2.75e-03
2.71e-04
1.74e-04
1.28e-04
RC
1.07e-01
6.16e-01
3.59e-01
1.02e-01
6.20e-05
1.10e-04
8.03e-03
8.40e-04
Table 2: KL divergence between model-predicted and reference LD-weighted PDFs for the Duffing equation, evaluated across training dataset sizes N={10,20,50,100,200,300,400,500} . Lower values indicate better reconstruction of the global phase-space structure.
Figure 9: Diagnostic visualization of model failures with insufficient training data. Left column (A–C): normalized prediction error distributions across phase space. Right column (D–F): corresponding LD-weighted PDFs. Top row (A,D): SympNet (100 trajectories) misses the figure-eight topology, learning instead a single center structure at (0,0) . Middle row (B,E): HénonNet (50 trajectories) produces a distorted homoclinic orbit with asymmetric lobes and incorrect curvature. Bottom row (C,F): RC (100 trajectories) captures the global homoclinic boundary but fails to reproduce the two center fixed points inside the orbit. LD analysis reveals specific structural errors that remain hidden in conventional error metrics.
NN type
Learned Hamiltonians
Layers
Neurons per layer
Trainable parameters
FLOPs at inference
SympNet
10
1
50
4000
10000
HénonNet
10
1
50
1010
10040
GHNN
5
2
15
10500
10800
RC
–
–
400
2000
9996
Table 3: Hyperparameters of the NN architectures for the three-mode NLS system.
Figure 10: Training (left) and validation (right) loss for the NN models in the three-mode NLS case. The legend is similar with Fig. 2 .
Figure 11: Forward and backward forecasts of the three-mode NLS system, initialized near the homoclinic orbit at u0=(η(0),ϕ(0))=(0.4,0.57) . Top rows: mode-0 dynamics ( q0 , left; p0 , right) and mode-1 dynamics ( q1 , left; p1 , right). Bottom row: pointwise prediction error (left) and normalized LD/ τ (right) versus integration time τ . Model comparisons follow the color/marker scheme of Fig. 3 . All NNs were trained on the 500-trajectory dataset.
Figure 12: Same as Fig. 11 , but for initial conditions (η(0),ϕ(0))=(0.4,0.4) (inside the homoclinic orbit).
Figure 13: Same as Fig. 11 , but for initial conditions (η(0),ϕ(0))=(0.4,1.57) (outside the homoclinic orbit).
Figure 14: Same as Fig. 11 , but for initial conditions (η(0),ϕ(0))=(0.99,1.57) (near the phase-space boundary).
Figure 15: Spatial distribution of the normalized prediction error in the (η,ϕ) plane for the three-mode NLS system, as defined in Eq. ( 40 ). Errors are evaluated at τ=10 and compared with the reference solution: (A) SympNet, (B) HénonNet, (C) GHNN, (D) RC.
Figure 16: LD-weighted PDFs for the three-mode NLS system in polar coordinates (η,ϕ) . Panel A: reference solution from numerical integration. Panels B–E: predictions from SympNet, HénonNet, GHNN, and RC, respectively. RC most accurately reproduces the reference structure, resolving the homoclinic orbit and fixed points. SympNet, GHNN, and HénonNet display distorted density distributions, with displaced or missing fixed points inside the homoclinic orbit.
Figure 17: Homoclinic orbit predictions in polar coordinates (η,ϕ) for NN models compared with the exact orbit (Eq. ( 39 ), blue dashed). Color scheme as in Fig. 8 . RC reproduces the orbit with high fidelity, while the symplectic models preserve the overall figure-eight structure but underestimate its diameter and curvature. All models trained on 500 trajectories.
Training distribution
NN model
Number of training trajectories N
125
250
500
1000
Uniform (inside + outside)
SympNet
2.14e-03
2.81e-03
2.36e-03
3.39e-03
HénonNet
3.59e-03
3.60e-03
3.73e-03
1.18e-02
GHNN
3.90e-03
1.87e-03
3.47e-03
1.86e-03
RC
1.72e-04
8.11e-05
1.04e-05
3.60e-04
Inside only
SympNet
3.35e-03
5.07e-03
4.45e-03
4.76e-03
Table 4: KL divergence between model-predicted and reference LD-weighted PDFs for the three-mode NLS, evaluated across training dataset sizes N={125,250,500,1000} . Lower values indicate better reconstruction of the global phase-space structure. Bold entries mark the best-performing model for each condition. RC achieves superior performance (by 1–2 orders of magnitude) with uniform and outside-only training, while all models perform similarly when trained exclusively inside the homoclinic orbit.
Figure 18: Normalized prediction error in the (η,ϕ) phase plane for GHNN (top row) and RC (bottom row) trained on different data distributions: inside-homoclinic only (left), uniform (middle), and outside-homoclinic only (right). Training with biased data yields complementary error patterns: inside-only (outside-only) training produces large errors in outside (inside) regions. Even with uniform data, all models show increased errors along the homoclinic boundary, though RC maintains substantially lower errors overall.
Figure 19: LD-weighted PDFs for cases correspond to Fig. 18 . Training distribution dramatically affects orbit reconstruction: inside-only training produces contracted orbits, while outside-only training yields blurred orbits near the saddle point at ϕ=0 . RC better preserves global topology across all cases compared to the symplectic architectures.
Figure 20: Panels show cross-sections of computed LD values (left column) and corresponding LD-weighted PDFs ρ at p=0.4 (right column) for varying integration time τ (top row), LD exponent c (middle row), and weighting function g (bottom row). The baseline parameters are τ=4 , c=0.7 , and g(x)=1/x (shown as red curves). Gray shaded regions indicate areas inside the homoclinic orbit. Different parameter choices systematically affect the relative importance assigned to various phase-space regions.
Figure 21: Similar to Fig. 20 , but for p=0 .
Section
Parameter
SympNet
HénonNet
GHNN
RC
τ
2.0
1.28e-03 (2)
9.11e-03 (4)
2.95e-03 (3)
3.09e-04 (1)
4.0*
1.32e-03 (2)
8.03e-03 (4)
2.75e-03 (3)
6.20e-05 (1)
6.0
1.74e-03 (2)
8.50e-03 (4)
3.24e-03 (3)
3.22e-05 (1)
10.0
1.30e-03 (2)
6.90e-03 (4)
2.55e-03 (3)
2.51e-05 (1)
Δt
0.1*
1.32e-03 (2)
8.03e-03 (4)
2.75e-03 (3)
6.20e-05 (1)
0.2
1.31e-03 (2)
8.02e-03 (4)
2.74e-03 (3)
6.25e-05 (1)
Table 5: Sensitivity analysis of the LD-weighted PDF framework for the Duffing oscillator system using 200 training trajectories. The table summarizes KL divergence values between model-predicted and reference LD-weighted PDFs under variations in integration time τ , time discretization Δt , exponent c in the LD definition, choice of weighting function g . Rankings (in parentheses) indicate relative model performance, with 1 being the best. RC consistently outperforms symplectic architectures across all parameter settings, demonstrating the robustness of our comparative framework.
Appendix figures & tables1 asset
Supplementary material from the paper’s appendix.
Appendix
Figure 22: LD for the Duffing equation correspond to Fig. 7 . The left panels: (A) Numerical integration, (B) SympNet, (C) HénonNet, (D) GHNN, (E) RC. All NNs were trained on the dataset corresponding to the 200-trajectory case. High-density regions correspond to slow dynamics, revealing the fixed points near (q,p)=(±1,0) . ight panels: one-dimensional slices of the LD-weighted PDF across all methods at selected phase-space locations. From top to bottom: q=−1.0 , q=−0.7 , p=−0.4 , p=0.0 , and p=0.4 .
We investigate the applicability of Lagrangian and Hamiltonian Neural Network models to a dissipative system that has explicit time dependence in its Lagrangian, Hamiltonian, and total energy. To do so we consider these neural network models for simulated systems of a harmonic one-dimensional, one-component oscillator with damping, as well as without damping for comparison. We find that both the Lagrangian and Hamiltonian approaches are able to predict the empirical physical behavior of the damped oscillator systems and to effectively ``learn'' to varying degrees the underlying Lagrangians and Hamiltonians, as has previously been shown to be the case with undamped oscillator systems. These investigations elucidate important properties of Lagrangian and Hamiltonian mechanics, including properties that are not manifest when considering systems without explicit time dependence.
We develop a geometric framework for learning deterministic and stochastic forced Hamiltonian systems with neural networks. Motivated by the Lagrange-d'Alembert principle and the theory of variational integrators, we introduce the notion of a Lagrange-d'Alembert map and establish a Cr convergence theorem for first-order one-step methods. Building on these results, we propose Generalized Forced Hamiltonian Neural Networks (GFHNNs), a class of structure-preserving neural networks obtained by concatenating Lagrange-d'Alembert-Euler maps, and prove a universal approximation theorem for this architecture. We further extend the framework to parameter-dependent systems, leading to Parametric Generalized Forced Hamiltonian Neural Networks (PGFHNNs). By interpreting the multiple Stratonovich integrals appearing in the Stratonovich-Taylor expansion as parameters, the same framework can be applied to stochastic forced Hamiltonian systems whenever information about the underlying Wiener process is available. Our numerical experiments demonstrate that the proposed geometric architectures provide significantly improved long-time stability and accuracy compared to non-geometric residual neural networks, while requiring substantially less training data to achieve a comparable level of performance.
Benedikt Brantner, Tomasz Tyranowski
Max-Planck-Institut für Plasmaphysik Boltzmannstraße 2, 85748 Garching, Germany · Technische Universität München, Zentrum Mathematik Boltzmannstraße 3, 85748 Garching, Germany · University of Twente, Department of Applied Mathematics PO Box 217, 7500AE Enschede, The Netherlands
Hamiltonian Neural Networks (HNNs) parameterize conservative dynamics through a learned scalar Hamiltonian, providing an architectural prior that is absent from generic vector-field neural networks. We evaluate this prior under a controlled protocol in which an HNN and a parameter-matched feedforward baseline are trained on the same RK4-generated trajectories, use the same central-difference derivative targets and optimization settings, and are integrated at inference with the same RK4 scheme. Results are reported over five independent training seeds. On the nonlinear pendulum, the HNN reduces mean energy drift by 42-fold and mean trajectory MSE by 15.8-fold at T = 100, approximately 16 pendulum periods. Its energy drift also remains bounded and exhibits substantially lower seed-to-seed variability than the standard-network baseline. An energy-stratified analysis shows that the difference becomes more pronounced as trajectories explore more nonlinear regions of phase space. As an additional diagnostic, we examine an explicit Störmer--Verlet-style rollout of the learned HNN. Because the learned Hamiltonian is not constrained to the separable form H(q,p) = T(p) + V(q), the standard symplecticity guarantee of velocity Verlet does not directly apply. We further apply the same matched-integrator protocol to the three-dimensional Kepler two-body problem. The HNN again exhibits lower trajectory, energy, and angular-momentum drift than the parameter-matched baseline. These experiments provide a controlled study of how Hamiltonian parameterization affects long-horizon prediction and physical consistency across two conservative dynamical systems.