UEPI: Universal  Energy-Behavior-Preserving Integrators for Energy Conservative/Dissipative  Differential Equations
Authors: Elena Celledoni
Brynjulf Owren
Chong Shen
Baige Xu
Takaharu Yaguchi
Venue: NeurIPS 2025 poster ()
Keywords: energy-preserving integrator
energy-dissipating integrator
Hamiltonian system
numerical integrator
discrete gradient method
Primary area: machine_learning_for_sciences
TL;DR: We propose a method for learning energy-preserving or -dissipating numerical schemes. Our approach admits the universal approximation property.

Abstract: Physical phenomena in the real world are often described by energy-based modeling theories, such as Hamiltonian mechanics or the Landau theory. It is known that physical phenomena based on these theories have an energy conservation law or a dissipation law. Therefore, in the simulations of such physical phenomena, numerical methods that preserve the energy-conservation or dissipation laws are desirable. However, because various energy-behavior-preserving numerical methods have been proposed, it is difficult to discover the best one. In this study, we propose a method for learning highly accurate energy-behavior-preserving integrators from data. Numerical results show that our approach certainly learns energy-behavior-preserving numerical methods that are more accurate than existing numerical methods for various differential equations, including chaotic Hamiltonian systems, dissipative systems, and a nonlinear partial differential equation. We also provide universal approximation theorems for the proposed approach.

UEPI: Universal Energy-Bchavior-Preervin Integrators for Energy Conservative/Disipative Differential Equations
Elena Celledoni
Department of Mathematical Sciences
Norwegian University of Science and Technology
Trondheim, Norway
elena. wwwweeeeeeeeeeeeee www.eeniiiinnu.uo..mmmmm no
Brynjulf Owren
Chong Shen
Department of Mathematical Sciences
Graduate School Science
Norwegian University of Science and Technology
Kobe University
Trondheim, Norway brynjulf.co wwwenentnu..m
Kobe, Japan wttp:o2660stu..obeeeeoo ac. jp
Baige Xu Graduate School of Science
Takaharu Yaguchi
Graduate School of Science
Kobe University Kobe, Japan
Kobe University RIKEN AIP Kobe, Japan
saigexu@stu.kobe-teo ac.
yaguchi@pear) kobe-u. u. ac.
Abstract
Physical phenomena in the real world are often described by energy-based mod- eling theories, such as Hamiltonian mechanics or the Landau theory. It is known that physical phenomena based on these theories have an energy conservation law or a dissipation law. Therefore, in the simulations of such physical phenom- ena, numerical methods that preserve the are desirable. However, because various mnergybeehavior methods have been proposed, it is difficult to discover the best one. In this study, we propose method for learning highly accurate energy beergy-behavior.. htape-bbbavvoooreeevinnee integrators from data. Numerical results show that our approach certainly learns numerical methods that more accurate than existing numerical methods for various differential equations, including chaotic Hamilto- nian systems, dissipative systems, and a nonlinear partial differential equation. We also provide universal approximation theorems for the proposed approach.
or dissipation laws www--bbaho-rreeeevveeeee preserving numerical
mergg-bechaioor-peseeee
Introduction
Differential equations serve as fundamental tools modeling the dynamics various physical systems. However, most real-world differential equations cannot be solved analytically, necessitating the use numerical methods to obtain approximate solutions. critical challenge this context is ensuring both accuracy and stability over long simulation times, especially for systems governed by complex or nonlinear dynamics.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Table 1:
Comparison of the proposed method, discrete gradients, SympNets Bi]. Because the
proposed method is a numerical integrator, the approximation accuracy is always guaranteed to certain extent, regardless of the results of training.
naive Runge-Kutta
discrete gradient
SympNets
Proposed
numerical integrators
energy conservation energy dissipation learnable from data
In recent years, advances in machine learning have inspired a new class of data-driven numerical methods 112 21, where neural networks used optimize the coefficients or design high-order schemes for classical integrators such as Runge-Kutte methods. These methods offer remarkable adaptability and high short-term accuracy when tailored to specific problems. However, significant limitation is their lack of preservation of key invariants, such as energy conservation or symplecticity, which are essential for long-term simulations. In contrast, geometric numerical integrators [24 are specifically designed to preserve such intrinsic structures, thereby ensuring better performance over long time intervals. Among them, energy prereeving methods such as the discrete gradient methods 47 have gained attention for their ability to conserve the system's 's energy evolution precisely. To the best of our knowledge, existing literature has not yet explored the potential neural networks in optimizing energy- preegy-preserving geometric integrators.
In this study, we propose method for learning numerical integrators that preserve the energy behaviors of target differential equations. More specifically, the focus of this paper is on dynamical systems characterized by their energy-conserving dissipating properties. We present general form of the discrete gradient method and propose method for learning highly accurate mergy-behavion preserving integrators from data by employing neural networks. In addition, we show universal approximation theorems for the proposed approach.
The main contributions of this paper are as follows.
We characterize and then parameterize the energy-behavior-pessseee numerical integrators.
We propose neural network models with universal approximation properties for the energy- behavior- parevivo-presrvvnggg. numerical integrators. Numerical experiments have confirmed that the method can be applied wide range of phenomena, including chaotic conservative systems, dissipative systems, and nonlinear partial differential equations.
Numerical experiments also show that the proposed method, in fact, finds highly accurate computational methods, which are difficult to discover theoretically. In particular, it is numerically shown that numerical schemes derived by the proposed method are highly accurate and stable even for chaotic systems.
2
Related Work
Neural networks for improving numerical methods represent a hybrid paradigm. One line of re- search focuses on optimizing classical numerical methods through data-driven parameterization rather than approximating dynamics directly. Anastassi [I developed artificial neural networks(ANN) to generate the coefficients of two-stage Runge- Kutta methods specifically optimized for solving the two-body problem. Dehghanpour et al. 112 utilized ANNs to automatically compute the optimal coefficients of third-order Runge-Kutta methods. Similarly, Guo [21] explored machine learning methods based on the Runge-Kutta-in Kutta-igrtttt architecture, capable of automatically generating high-order integrators for specific families of ODEs. Universal Numerical Integrators (UNIs) 60, 61 proposed framework, which combines neural networks as universal approximators with numerical integrators. second distinct line of work investigates learning discretization schemes from data, Bar-Sinai et al. (2) introduced a data-driven discretization approach that employs neural networks to optimize the computation of spatial derivatives. Maddu et al. 44 developed STENCIL-NET, which can adaptively learn local discretization operators based on solution data, addressing long-term prediction challenges for unknown nonlinear PDEs. Ranade et (55) introduced DiscretizationNet. generative baN-based approach to model discretization schemes. Our approach follows this hybrid

---

paradigm but shifts the focus from general integrator optimization to preserving the energy behavior of dynamical systems. We achieve this by combining neural networks with the structure-prest prreerire-reservving discrete gradient method, rather than modifying traditional integration integrators.
Geometric numerical integration focuses on designing numerical methods that preserve certain qualitative characteristics of dynamical systems. Unlike traditional methods such as the explicit Euler method or the classical Runge-Kutta method, which may lead long-term errors, geometric integrators maintain properties that make them particularly suitable for long-time simulations of physical systems.
Geometric numerical integrators are generally divided into symplectic integrators and energy- preserving integrators, as Zhong [67] and Chartier et al. [6] have shown that it is impossible for numerical integrator to preserve both properties. Symplectic integrators (116) (66) are integrators methods aim to conserve the total energy of the system. Among energy-pree peegyy-preerrving methods, the discrete gradient method 37 47l) is particularly successful and has also been extended to partial differential equations (15 46) ) Eidnes [15] develop order theory for discrete gradient methods. Projection
that preserve the symplectic structure of Hamiltonian systems. Emergy-preserviie
2 can also be employed to preserve energy by projecting the numerical solution back onto the manifold that satisfies a certain energy. Norton et al. (50) show that linear projection methods are a subset of discrete gradient methods. Miyatake and Butcher [48] provided a condition to be mneryy-preserye numerical integrators, see
methods(
for continuous stage Runge- Kutta methods. For other geometric 57)) for details. In our study, we place particular emphasis on discrete gradient methods, as they inherently preserve the energy behaviors, making them highly relevant for simulations involving Hamiltonian dynamics and dissipative systems.
Neural network-Ban Based Physical Modelings have been widely applied to the study of differential equations by combining data-driven learning with physical principles. Chen al. 7 introduced Neural ODEs, which model uattnuuussttime dynamics via parameterized vector fields. This idea has been extended in various directions 3 3 63]. Another influential framework is that of Physics- Informed Neural Networks (PINNs) [54] which incorporates differential equation residuals into the loss function to enforce physical consistency during training. Numerous extensions have been developed to enhance their accuracy, robustness, efficiency, and scalability [30
32
49
54
[58).
Operator learning 38 42 aims to learn mappings between infinitediimensional function spaces. Dynamic Mode Decomposition (DMD) and the Koopman operator theory have emerged as powerful data-driven techniques for approximating nonlinear dynamical systems with linear representations neural networks have been proposed to better reflect the physical properties of dynamical systems. Hamiltonian Neural Networks (HNNs) (20 model the Hamiltonian to enforce energy conservation. Extensions of HNNs
34 43 (59 63). Beyond generic data-driven models, grructure-preserving
45 have incorporated flexibility and accuracy improvements. Lagrangian Neural Networks (LNNs) [9] take similar approach via Lagrangian mechanics. Many studies also focus the symplecticity 88 31551 62). For muergy-dissipativ systems, some studics have introduced neural networks grounded in thermodynamic principles such as the GENERIC framework to ensure stability, interpretability, and physical consistency in modeling irreversible dynamics 26 27 66). Despite their effectiveness in capturing qualitative physical properties, these methods are not numerical integrators. They focus on learning dynamics or flow maps, and their approximation quality can be difficult to control or quantify in classical numerical analysis sense. This limits their applicability in contexts where rigorous accuracy guarantees or long-term stability are required.
11
25
In summary, our method bridges geometric numerical integration and neural networks by learning integrators. This offers improved long-term stability and accuracy in simulating wide range of systems.
3 Target Energy-Basse -Baee Dynamical Systems
We focus differential equations written in the form
= L=S(u)VH()), u u R4.
(I)
Here H E C² R) and S E C(Rd, Rdxd). Differential equations in this form describe vari- ous physical systems modeled by ordinary differential equations (ODEs) and discretized partial differential equations (PDEs) [ 52). Two classes of problems are usually considered:
3

---

Conservative case: S(u)² T = =-((u) (skew-symin memmetyymmeric.)
Dissipative case: S(u) - O (negative semi-definite.)
Note that S O means that the matrix S negative semi-definite, i.e., u Su S 0 for any vector u. Also, if matrix S is skew-symmetric, then u Su for any vector u.
The systems described by the above equations have the following properties.
Theorem 3.1. The system has the energy dissipation law if S & O and the energy conservation law S is skew-symmetric.
See Appendix A for the proof. Such differential equations include various equations that appear in physics. For example, if S is skew semmeeyiittiic and invertible and also satisfies certain condition called the Jacobi identity, the equation describes Hamiltonian systems. As the Hamilton equation is generalization of Newton's 's equation of motion, mechanical systems that be described by classical mechanics, such as the motion of a pendulum. As an example outside of physics, the Lotka-Volterr predator-prey model is also Hamiltonian system. On the other hand, this equation describes dissipative systems when S negative-definite A typical example of dissipative systems the equation of motion for mechanical systems with friction. Thermodynamic systems can often be expressed as dissipative systems as well. Also, this equation has a geometric background; see Appendix for details.
The system I also describes semi-discretized Hamiltonian or dissipative PDEs in space. For example, the PDEs of the form
du a 8H
Öt Ör 8u
are an example of Hamiltonian PDEs, where Huu, u,) is the energy density and 8 6H/óu denotes the variational derivative of H1 which is defined by otf 8u == 0H a aH This equation admits
Du. 8a du
the energy conservation law: di J Hdr under, for example, the periodic boundary condition. Similarly, the PDEs
du 82 8H
at 2x 8u
admits the energy dissipation property: 4 ʃ Hdx S under the periodic boundary condition. See Appendix (fr details on the semi-discretizaae
Hamiltonian PDEs include the Maxwell equation and the shallow water equations such as the KdV equation, the advection equation, and the Burgers equation. Dissipative PDEs express physical systems derived from the Landau free-energy minimization including the Cahn Cahn-Hilliard equation and the phase-fiet field models for phase transitions and pattern formulations. Other target equations include the equations with complex state variables, such as the Schrödinger equation and the Ginzburg- Landau equation. See, e.g., 17 for details.
4 Methods
4.1 Discrete Gradient Methods as Snergy-Behavior-Prreeriigg Integrators
A popular way of discretizing TT is by a class of numerical schemes called discrete gradient methods defined as
( (n)
(n)()), ( (m+
(H() (()+))
h
where h is the step size. In this discrete counterpart to ① the skew-symmetrix matrix S(u) has been replaced by the matrix S(u, v) depending
arguments, but still required to be skew-symmetric. VH RRd x Rdd + Rd is the discrete grudient defined as follows.
Definition 4.1. A discrete gradient for a smooth function H : .R + R is a continuous mapping H Rd Rd + Rd that satisfies the following properties:
H(v) Hu)) = Hu, v), u) for all u, E Rdd,
(5)
H u) (u) for all e d
(6)

---

The first condition corresponds to the discrete chain-rule dH (Óu; u) (H( ou for the Fréchet derivative (::u) of u. where du an infinitesimal change of u. The second condition ensures that a discrete gradient V H is certainly typically the standard Hermitian inner product for ODEs and the discrete I,a inner product for discretized PDEs. For details on discrete gradient methods, see, e.g. 15].
approximation of the gradient VH. The inner product
The discrete gradient schemes preserve the energy behavior of the target differential equations in the following sense.
Theorem 4.2. The discrete gradient scheme admits the discrete energy conservation law
H( H(l!)))
i is skew-symmetric and the discrete energy dissipation law
H( S (u!
if & O. In particular, if the system is dissipative, the amount of energy dissipation an approxima- tion of that of the original differential equation.
Proof. The proof is exactly the same as that for Theorem 3.1 See Appendix A for details.
There several popular choices for discrete gradients. The average vector field (AVF) discrete gradient [53] is defined as
(H(u, VH((1 - (c)
The Gonzales discrete gradient
19
is defined by
u +
H(v) H(u) (VH (x
u)
(H( = H
+
(v u)
2
{u u. u))
Another example is the Itoh-Abe discrete gradient 29
H Uj Meguon H tta
21 H ug.....
Pg
ItohhAb H(,
H H
H( 1884
4.2 Universal Ehergg-Behavvir-PPreeeving Integrators
The above three discrete gradients are just examples from the family of maps H satisfying (5) and (6) above. this paper, we fully characterize the set of all possible discrete gradients as follows.
Let V H be particular choice of discrete gradient, henceforth denoted the reference gradient. Then, from it follows that any other discrete gradient HH, (u, v) satisfies
H(, v) Ho, v), v- u) =
Thus, the set of all discrete gradients be described, and hence. parameterized for learning as follows.
Theorem 4.3 (Characterization of discrete gradients). For fixed function H, the of all discrete gradient operators V Rd Ed + IRd has the following parametrization
VH(,,)) VoHl (u, v) + A(u, v)(v -u),
(8)
where A(u, w) is continuous and skew symmetric for all E TRd, and where Vo v) is fixed reference gradient satisfying and
Proof. For the proof, see Appendix A
In fact, we can prove a stronger statement: all energy-conserving numerical methods for the Hamilton equations can be written in this form.
5

---

Theorem 4,4. For a fixed function H, any energy-preserving integrators for the Hamilton equation
du dt = (v)
can be written as a discrete gradient scheme with a discrete gradient of the following ( v) = (u, v) + A(u, v)(u - -u),
where A(u, v) is skew symmetric for all u, v E I& and VoH v) is a fixed reference gradient satisfying (5) and (6)
Proof. See Appendix A for the proof.
Note that we cannot obtain similar theorem for dissipative systems in straightforward way because, in the proof of the above theorem, we use the energy-conservation law of the conservative systems and also the non-degeneracy of S.
Proposed Method Although the discrete gradient method preserves the energy behaviors, it still introduces numerical errors that accumulate and grow over time. In this study, we propose method for learning highly accu- rate mnergy-behavior-preseevv can be summarized as Figure [ This paper aims to op- timize the discrete gradient based on the above charac- terization of discrete gradients. with the aim of not only preserving the energy but also yielding numerical solutions closer to the true solution.
usn ((++)))
integrators. The method
(H (
MNN
The proposed method is formulated follows:
V xuH xNH(( 21 HH(u!)) aoa- MMx (u(an))) (000)
NX H (u(m), u(x+))) = (u(( w(a+1))
+ MNN (n
- (n))), (9)
MNNH(ul"), wentiy
where MN)(u(n), ())) skew-symmetric matrix learned by neural network that takes (n) and (n+1)) as the input. The universal approximation properties the direct consequences of the above theorems.
Figure 1: Outline of proposed method
Theorem 4.5 (Universal approximation theorem for mergy-preserree integrators for Hamiltonian systems). Suppose that the phase space is compact and MxN modeled by using a multilayer perceptron (MLP) with sufficiently smooth activation function so that the MLPs admit the uni- versal approximation property for the continuous skew symmetric matrices M (u, v) ss. Then the proposed method admits the universal approximation property for energy-preserving integrators for Hamiltonian systems.
Theorem 4.6 (Universal approximation theorem for energy-dissiptee dissiptive discrete discrete-gradient integrators for dissipative systems). Suppose that the phase space is compact and MsN is modeled by using MLP with a sufficiently smooth activation function so that the MLPs admit the universal approximation property for the continuous skew symmetric matrices Mu, v) 's. Then the proposed method admits the universal appreximation property for energy-dissipee dissyddiss discrete discrete-grudien integrators for dissipative systems.
Given a time series data (n)'s at time ((m)ss 's the proposed model is trained to minimize the squared emror between the left- and right-hand sides of (9
(n+1) (n)
minimize Σ
SV NYNH (w un
(10)
h(m)
where (m)) ((++1) ((m).
If the matrix MnN satisfies MNN(U, v) MNx u), then the proposed integrator ( at least has the second-order accuracy. Otherwise, the proposed integrator g is the first-order method (24]

---

Experiment
We used standard multilayer perceptrons with the hyperbolic tangent (tanh) as the activation function to learn the matrix MNN. The hyperparameters for training are provided in Appendix D.I. The optimizer is set Adam B3 As for the loss function, we chose MSELoss (Mean Squared Error Loss). Final loss values are provided in Appendix D.2] All experiments were performed on an Intel 15-13500H CPU, and computations were done in double precision.
In our experiments, we employ the AVF and Itoh-Abe methods as reference discrete gradient. As demonstrated in the experiments provided in the Appendix [D.2] for the choice of the reference discrete gradient, we prefer to use the AVF method for its accuracy when it is computationally feasible. Otherwise, the Itoh-Abe method offers a simpler approach, making it more suitable for general equations. Moreover, to improve the performance of the Itoh-Abe method, we employ its time-symmetric version as the reference method all subsequent experiments follows:
V VitohAbesyas (u, v) (VitohAbe (u, v) + Hto u)) /2.
Although the proposed method can theoretically guarantee at most second-order accuracy, we conduct comparisons with reference numerical methods using half step size and double steps, as well with fourth-order AVF and Itoh-Abe methods constructed by the composition method 24
The implementation code is available in supplementary material on OpenReview:
https://
5.1
AVF Diserete-Graddent-BBas Method for Energy-Conservatiy Nonlinear System
Hénon-Heiles System stellar motion restricted to plane. shows chaotic motion and has a separable and polynomial Hamiltonian given by:
The Hénon-Heiles model (Figure 2 describes the
1 (0f p3) u0 +
Hp. g) =
Figure 2: Hénon- Heiles
Since this is a polynomial Hamiltonian equation, the AVF method can be easily applied. Due to the chaotic characteristics of the system, even small errors can evolve into large deviations over long iterations. Although the AVF method conserves energy, can still lead to significant errors for this system.
We use two sets of data: the flow data and the irregular data for the Hénon-Heiles equation. The training flow data {(u(x-1), 1), (())
) n 1, ..., NN consists of N = 100 pairs of points obtained from the trajectory calculated using a high-order integrator sampled with a time step h = 0.3, starting from
(0) = (01") (0) (0) (0) ) (0.3, 0.3, -0.2, -0.2). After training, we use the trained integrator to compute the flow starting u for 300 steps.
02 P1 PP2
Ground truit
Propored metted
Stepe
Kal onders
So moners
GINNN WTT
L23
055
055
-25
6.5
-55
06
05
00
-02
45
-55
60
45
80
Figure 3: Results for the Hénon-Heiles eq. (predicted trajectory and global error.) The error plot uses the same color as the trajectory plots.
We also consider the computational time because the proposed numerical scheme is computationally more expensive than the AVF method due to the use of the neural network. When using a 10-layer, 100-width network for the Hénon-Heiles experiment, the computation time for the trajectory of 300 points falls between the AVF method and the AVF (half step size), which uses half the time step size and doubles the number of steps. For example, The computational times were: Proposed: 5.522 [sec], AVF: 3.608 [sec], AVF (half step size):: 6.616 [sec]. In fact, the proposed method outperforms the fourth-order method, which offer higher accuracy while require more computation time than AVF(half step size). We also compared the trajectory computed using 10-layer, 100-width G-sympnets trained for 100, 000 iterations with learning rate 0.001 and the same training dataset by modifying the officially released code [31]. The predicted trajectory and the corresponding global
7

---

error ([lu(*) - 222 comparison are illustrated in Figure from which we can deduce that the proposed method yields trajectory closer to the ground truth, and for details see Appendix [D4
We fürther conducted experiments using irregular data. The training dataset {((u((-1), (((u((t1),,ii (,()),,,))) u( i = - 1,..., /}
consists of / = 100 one-step input-output pairs, where the elements of uli are randomly sampled from the domain [-1, 11, and the corresponding step sizes h( are randomly chosen
from [0, 0.5]. Each (9) represents the one-step solution computed from (i-1)) using step size h()),
After training. we simulate the flow starting from the initial state (0)) (9) 90 Pi P2 ( = (0.3, 0.3, -02 -022) using fixed step size h 0.3 for 100 time steps. The numerical results indicate that the proposed method yields significantly improved accuracy over baseline methods. For detailed comparisons, please refer to Appendix D.4
(0) (0)
We also conducted experiments on the pendulum system to validate the generality of our proposed method. Detailed experimental settings and results are provided in Appendix D.3
5.2 Itoh-Abe Seecett--rradent
Method for Seergy-Conservative Nonlinear System
Two-body Problem In the two-body problem, the motion of two bodies that attract each other considered. We assume that one of the bodies is fixed at the center of the coordinate system and the position of the other one is given as q = (q), (2) (Figure The Hamiltonian is given as follows:
91
(q1,, 42 PI, P1, = 2 + P
Figure 4: Two- body
V V + 0
The training flow data {(u(n ())))
1, N consists of N 30 pairs with fixed time (0)
step h = 0.1, starting from (0) = (4) trained integrator to compute the flow starting at u (N) for 300 steps.
10) = ( (- 0.8, 0,,,0,0 After training, we use the
92
Ground buch
Propoued method
koh-Abe vrlution
Boh-Abe half stepsize
Booh-Abe oecent
Giball emorr
2.5
0.0
0.
9
2.0
-15
0.5
133
42
Figure 5: Results of the Two-body (predicted trajectory and global error.) The error plot uses the same color the trajectory plots.
Figure compares the trajectories predicted by the proposed method with Itoh-Abe methods. While the Abe methods yield solutions that generally follow the elliptical shape of the ground-truth, they exhibit noticeable deviations from the true trajectory. For details, see Appendix D.4
Double Pendulum
The double pendulum is known to exhibit complex behaviors with strong sensitivity to initial conditions, which makes the computation very challenging.
This system has a non non-polynomial Hamiltonian as follows:
((g Q2 P1. P2) pi 2pp 2p1p2 cos(1) q2)
2 cos Q1 COS q2-
2(1 + sin (q1 -q2))
The training flow data [n] (m)) n = 1,... N} consists of = 100 pairs with fixed time step h 0.3, starting from u (0)) = (g 92 PI P2z
Figure 6: 2-
(0) (0) (0) = (g Mr. 0)..
pendulum
After training, use the trained integrator to compute the flow starting at u (N) for 300 steps.
Although the chaotic nature of the double pendulum makes learning difficult, the proposed method provides a solution that is close to the true trajectory (Figure 7). For more details, see Appendix D.4
5.3
AVF Siscrete-Gradient-Baeeee Method for Dissipative Nonlinear System
The discrete gradient methods also be applied to dissipative systems
d dt
0 = 0
1
P
-1
-8
VH
8

---

0.75
SAUcion
0.75
Rohwibe helf stepsize
6.75
Iroh-Abe andeed
eno
0.52
03.3
0.23
0.25
0:23
0.00
8.00
-0.2
-0.25
-0.75
-1.25
2:50
-0.0
-0.35
-0.75
-0.
Figure 7: Results of the double pendulum (predicted trajectory and global error.) The error plot uses the same color as the trajectory plots.
to construct dirrypiiisppatigg integrators, where denotes the dissipative coefficient. Therefore, our proposed approach can be employed to optimize them.
Dissipative Pendulum The Hamiltonian for the dissipative pendulum is the same that for the pendulum (g. p) = 222² - cos(g))
Global errr
Energy
3.5
Gouund
Ahh cttos artert
Ax>
-325
The dissipation coefficient was set to 8 = 0.01. We train the network with the flow data ((u(R-1), u!") I n 1, ..., N) consists of N = 100 pairs with fixed time step
330
0.00
10
22 20
25 36
30
15
20
25
30
0.3,
Figure 8: Results of the dissipative pendulum
starting from (0)) (q(0) (0)) = (2, 0). After training, the trained integrator was used to compute the flow starting at u (N) for 100 steps. As seen in Figure & the accuracy of the integrator was greatly
improved again by the proposed method while preserving the dissipation property.
Duffing Oscillator The Duffing equation is a non-lineat differential equation that models certain damped oscillators with the Hamiltonian:
1 1 aaa H p) 2 + ==a²²= 09
where we set the linear stiffness coefficient to a 1, the non-linearity coefficient to B -5 and the damping coefficient to 8 = 0.05. We train the flow data {(u(n-1), ) 1, 1,,, N} consists of
N = 200 pairs with fixed time step h 0.1, starting from (g) (0)) (4, 0). After training, we computed the flow starting at (N) for 300 steps.
Figure D shows that the proposed method again improves the accuracy; however, compared to energy-conssee seeggyconserrving systems, the numerical error is slightly larger. This is because the discrete gradient method only guarantees the energy dis- sipation property but does not ensure the accu- racy of energy for energy-dissee energy-dissipative systems.
Global enoor
Encrgy
Amppoee methoe
bal Hespis Aander
20
25
20
15
Figure 9: Results of the Duffing oscillator
5.4 AVF Discrete-Gradient-Be
Method for Energy-Conservatte Nonlinear PDEs
As the discrete gradient method is applicable to semi-discretized PDEs, our approach is also applicable. We applied the proposed method to a discretized Hamiltonian PDE, the Korteweg-de Vries (KdV) equation:
Ou
6u Ou 8au
Öt
=
Ox
3z
The semi-discrete Hamiltonian is defined as H(u)
Global emrrr
( TA (us+1 n2)² uz) where the index denotes spatial dis- cretization points at fixed time level with the spatial step size Ar.
6.000
0.001
Propoced mittool Aue
0.002
This equation a partial differential equation that describes shallow water waves. We consider this equation on the interval [0, 20) under the periodic boundary condition. We used the training data with the time length of = 30 and the time step size of At 0.01. The spatial step size was Az 0.5. For testing,
Figure Results of KdV
10:

---

the computation was carried out from the time length T = 30, with the starting being at the final state of the training dataset.
The numerical errors shown in Figure 10 confirm that the proposed method yields more accurate solutions compared to the AVF method. For more details, Appendix D.4
5.5 Supplementary Experiments
Impact of Training Data We also conducted comparative experiments to investigate how varying the size of irregularly training data influences the performance of the learned integrator. The experimental results demonstrate that increasing the amount of training data beneficial for improving the training performance of the learned integrator. See Appendix D.2 for details
Comparison with Other Methods For completeness of exposition, we compare our approach with other methods, Our method focuses on learning structure-preserving numerical integrators, whereas many existing approaches (e.g.. PINNs) aim directly approximate continuous solutions. Because these objectives differ, such comparisons are inherently difficult and should be interpreted with care. We therefore report supplementary results with PINNs and TRS-ODEN on standard dynamical systems. The protocol and complete results are given in Appendix D.5 Overall, our method shows smaller errors and more stable long-horizon trajectories.
Computational Time One limitation of the method the increased computational time associated with neural networks during training and inference. In particular, the proposed method requires time for data preparation and network training. To see the actual increase in total computation time, we conducted additional numerical experiments. Specifically, we discuss the computation time in detail, especially for the case of the KdV equation, which is the largest system in the experiments. For other experiments, Appendix D.6 Also, GPUs can be used for training to accelerate calculations. In the experiment, we used an Intel Xeon 6900P CPU and an AMD Instinct MI300A GPU for computation.
We first set the number of nodes for spatial discretization of the KdV equation to 100. The computation time for generating the training data was 28. 709 sec using a CPU and the scipy odeint function. The training time was 6.407 sec using MI300A GPU. The integration time during testing was 52.738 (proposed), 22.785 (2nd order AVF) and 127.280 (4th order AVF). Thus the total computational time of the proposed method was 28. 709 + 6.407 + 528 87.85 87.854, which is less than the 127.28 sec required by the 4th order AVF method. Similarly, we set the number the computational time for each method was 50.957 (data generation) 8.311 (training) 76.460 (integration) = 135.728 (proposed), 33.034 (2nd order AVF), 200.003 (4th order AVF). The total computational time of the proposed method was again shorter than that of the 4th-order AVF method. Because our proposed approach typically outperforms the 4th-order AVF method, the proposed approach can compute more accurate numerical solutions in shorter amount of time, even when the training time and the time for generating data are included.
the nodes to 120,
A closer look at the execution time of the above experiments shows that the computation of the neural network itself is actually not so large. Rather, the matrix-vecte vector multiplication of the skew-symmetric matrix affects the computational complexity. However, the effect of this matrix multiplication does not seem to be so signíficant, since the multiplication with the sem-ssymmttiic or begative-definie matrix placed in front of the gradient of the energy takes about the same amount of computation. It is also possible to reduce computational cost by low-rank approximation of the learned matrices. In addition, our method does not necessarily need to be used at all time steps. By using it only when higher computational accuracy is desired, computational time can be further reduced.
Conclusion
In this work, we proposed method for learning numerical integrators that preserve the energy conservation or dissipation laws with universal approximation properties. We validated the proposed approach across a diverse set of systems, including chaotic nonlinear Hamiltonian systems, dissi- pative systems, and nonlinear partial differential equation. In all scenarios, the method exhibited advantages, particularly in complex or chaotic regimes. With further optimization, we believe this method will play important role in application areas in the future.
10

---

Acknowledgments and Disclosure of Funding
Funding in direct support of this work: JST CREST Grant Number JPMICR1914 and JPMICR24Q5. ISPS KAKENHI Grant Number 25K15148, JST ASPIRE JPMIAP2329 and NIFS Collaborative Research NIFS25KISC015. This work is also supported by the Horizon Europe, MSCA-SE project 101131557 (REMODEL).
References
[1] Angelos Anastassi. Constructing Runge- Kutta methods with the use of artificial neural networks. Neural Computing and Applications, 25:229-25 236, 2014.
[2] Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P Brenner. Learning data-driven discretizations for partial differential equations. Proc. Natl. Acad. Sei. U. S. A., 116(31): 15344-15344 July 2019.
[3] Marin Biloš, Johanna Sommer, Syama Sundar Rangapuram, Tim Januschowski, and Stephan Günnemann. Neural flows: Efficient alternative to neural ODEs. Advances in neural information processing systems, 34:21325-21333 21337, 2021.
[4] C. J. Budd and M. D. Piggott. Geometric integration and its applications. 2003. ISSN 15708659,
[5] Elena Celledoni, Volker Grimm, Robert McLachlan, David McLaren, Dion O' Neale, Brynjulf Owren, and GRW Quispel. Preserving energy resp. dissipation in numerical PDEs using the "Average Vector Field" method. Journal of Computational Physics, 2377776777, 2012.
[6] Philippe Chartier, Erwan Faou, and Ander Murua. An algebraic approach invariant preserving integators: the case of quadratic and Hamiltonian invariants. Numerische Mathematik, 103: 575-590, 2006.
[7] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
[8] Zhengdao Chen, Jianyu Zhang, Martin Arjovsky, and Léon Bottou. Symplectic recurrent neural networks, 2020. URL httppp:://aa....o///aaaaaa//////
13334
[9] Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho. Lagrangian neural networks. arXiv preprint arXiv:2003. 04630, 2020.
[10] Morten Dahlby, Brynjulf Owren, and Takaharu Yaguchi. Preserving multiple first integrals by discrete gradients. Journal of Physics A: Mathematical and Theoretical, 44(30)330505,, 2011.
[11] Marco David and Florian Méhats. Symplectic learning for Hamiltonian neural networks. Journal of Computational Physics, 494:1 12495, 2023.
[12] M Dehghanpour, A Rahati, and E Dehghanian. ANN-based modeling of third order Runge Kutta method. J. Adv. Comput. Sei. Technol., 4(1): April 2015.
[13] Emilien Dupont, Amaad Doucet, and Yee whye Teh. Augmented neural ODEs. Advances in neural information processing systems, 32, 2019.
[14] Edda Eich. Convergence results for coordinate projection method applied to mechanical systems with algebraic constraints. SIAM Journal on Numerical Analysis, 3065):1477-1482 1993.
[15] Sølve Eidnes. Order theory for discrete gradient methods. BIT, 6204):1207-11555 December 2022.
[16] Kang Feng and Mengzhao Qin. Symplectic Geometric Algorithms for Hamiltonian Systems. Springer- Verlag, Berlin, 2010. doi: (0010077/978-33-42--17777
Epringer-Verlag,
[17] Daisuke Furihata. Finite Difference Schemes for 0u/al Conservation or Dissipation Property. Journal of Computational Physics, 156(1):181-20) noy 1999. ISSN 19219991. doi: 10. 00.1066/jch..1999666
(0/0x)) =0 8G/ōu That Inherit Energy

---

[18] Daisuke Furihata and Takayasu Matsuo. Discrete Variational Derivative Method: A Structure- Preserving Numerical Method for Partial Differential Equations. Chapman and Hal//CRC, dec 2010. ISBN 08804291484661
doi: 10. 20.1200//10387
[19] O. Gonzalez. Time integration and discrete Hamiltonian systems. Journal of Nonlinear Science,
6, 1996. ISSN 09388974. doi: 10. 10.1007/BFF0244112
[20] Samuel Greydanus, Misko Dzamba, and Jason Yosinski. Hamiltonian neural networks, Advances neurul information processing systems, 32, 2019.
[21] Yue Guo, Felix Dietrich, Tom Bertalan, Danimir T Doncevic, Manuel Dahmen, Ioannis G Kevrekidis, and Qianxiao Li. Personalized algorithm generation: case study in learning ODE integrators. arXiv [math.NA], May 2021.
[22] Ernst Hairer. Symmetric projection methods for differential equations manifolds. BIT Numerical Mathematics, 40:726-734, 2000.
[23] Ernst Hairer. Geometric integration of ordinary differential equations on manifolds. BIT Numerical Mathematics, 11:996-1007, 2001.
[24] Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric numerical integration, vol- ume 31 of Springer Series Computational Mathematics. Springer- Verlag, Berlin, second edition, 2006. ISBN 3-540-30663-3; 978-3-540-30063-4, Structure-preserving algorithms for ordinary differential equations.
[25] Chen-Di Han, Bryan Glaz, Mulugeta Haile, and Ying-Chene Lai. Adaptable Hamiltonian neural networks. Physical Review Research, 3(2):22355,, 2021.
[26] Quereus Hernández, Alberto Badías, David González, Francisco Chinesta, and Elías Cueto. Structure-preservvigg neural networks. Journal of Computational Physics, 126:1099900 2021.
[27] Quercus Hernández, Alberto Badías, Francisco Chinesta, and Elías Cueto. Thermodynamics- informed graph neural networks. IEEE Transactions on Artificial Intelligence, 5(3):967-976, 2022.
[28] In Huh, Eunho Yang, Sung Ju Hwang, and Jinwoo Shin. T ODE network. Advances in Neural Information Processing Systems, 33:19016-19027, 2020.
[29] Toshiaki Itoh and Kanji Abe. Hamittonian-conserving discrete canonical equations based variational difference quotients. Journal of Computational Physics, 76, 1988. ISSN 10902716.
doi:
[30] Ameya D Jagtap and George Em Karniadakis. Extended physics-informet informeed neural networks (XPINNs): A generalized space-time domain decomposition based deep learning framework for nonlinear partial differential equations. Communications in Computational Physics, 28(5),
2020.
[31] Pengzhan Jin, Zhen Zhang, Aiqing Zhu, Yifa Tang, and George Em Karniadakis. SympNets: Intrinsic structure-preset symplectic networks for identifying Hamiltonian systems. Neural Networks, 132, 2020. ISSN 18792782. doi: n0.10066j.ncumet22000 neumet.2020.08..
[32] Ehsan Kharazmi, Zhongqiang Zhang, and George Em Karniadakis. Variational physics-informed neural networks for solving partial differential equations. arXiv preprint 192.00391, 2019.
[33] Diederik P. Kingma and Jimmy Lei Ba. Adam: A method for stochastic optimization. In 3rd International Conference on Learning Representations, ICLR 2015 Conference Track Proceedings, 2015.
[34] Nathan Kutz, Steven L. Brunton, Bingni W. Brunton, and Joshua Proctor. Dynamic Mode Decompoxition. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2016. doi: 10. 1 137/1 97816119745000 URL. https: ://pubs.//pbp/////// siam. www..oogg//o/oo/// aas///b//////// 1137/1
9781611974508
12

---

[35] Robert A LaBudde and Donald Greenspan. Discrete mechanics a general treatment. Journal of Computational Physics, 15(2):134-167, 1974.
[36] Robert LaBudde and Donald Greenspan. Energy and momentum conserving methods of arbitrary order for the numerical integration of equations of motion: I. motion of a single particle. Numerische Mathematik, 25:323-346, 1975.
[37] Robert LaBudde and Donald Greenspan. Energy and momentum conserving methods of arbitrary order for the numerical integration of equations of motion: II. motion of system of particles. Numerische Mathematik, 26:1-16, 1976.
[38] Haolin Li, Yuyang Miao, Zahra Sharif Khodaei, and MH Aliabadi. An architectural analysis of DeepOnet and a general extension of the physics-informed DeepOnet model on solving nonlinear parametric partial differential equations. Neurocomputiigg 61 1:1675,,,,, 2025.
[39] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differ- ential equations. In ICLR 2021 9th International Conference on Learning Representations, 2021.
[40] Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for PDEs on general geometries. Journal of Machine Learning Research, 24(388):1 2023.
[41] Zongyi Li, Hongkai Zheng, Nikola Kovachki, David Jin, Haoxuan Chen, Burigede Liu, Kamyar Azizzadenesheli, and Anima Anandkumar. Physics-informed neural operator for learning partial differential equations. ACM/JMS Journal of Data Science, 1(3):1-27, 2024.
[42] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 3, 2021. ISSN 25225839. doi: 01036//422256002-0000023
[43] Bethany Lusch, J Nathan Kutz, and Steven Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature communications, 9(1):4950, 2018
[44] Suryanarayana Maddu, Dominik Sturm, Bevan L Cheeseman, Christian Müller, and Ivo Sbalzarini. STENCIL-NET: Data-driven solution-adaptive discretization of partial differential equations. arXiv [math.NA], January 2021.
[45] Marios Mattheakis, David Sondak, Akshunna S Dogra, and Pavlos Protopapas. Hamiltonian neural networks for solving equations of motion. Physical Review E, 105(6))065305, 2022.
[46] Robert I. McLachlan and G. R.W. Quispel. Discrete gradient methods have an energy conserva- tion law. Discrete and Continuous Dynamical Systems- Series A, 34, 2014. ISSN 10780947. doi: 10.3934/dcds.20144 22200.444.1999
[47] Robert I McLachlan, G R W Quispel, and Nicolas Robidoux. Geometric integration using discrete gradients. Philos. Trans. A Math. Phys. Eng. Sci., 357(1754):1011-1045 April 1999.
[48] Yuto Miyatake and John C Butcher. characterization of energy-preserving methods and the construction of parallel integrators for Hamiltonian systems. SIAM Journal Numerical Analysis, 1993-20133, 2016.
[49] Ben Moseley, Andrew Markham, and Tarje Nissen-Meyer. Finite basis physics-informed neural networks (FBPINNs): scalable domain decomposition approach for solving differential equations. Advances in Computational Mathematics, 49(4):62, 2023.
[50] Richard A. Norton, David I. McLaren, G. R. Quispel, Ari Stern, and Antonella Zanna. Projection methods and discrete gradient methods for preserving first integrals of ODEs. Discrete and Continuous Dynamical Systems- Series A, 35, 2015. ISSN 15535231. doi: 10.3934/dcds 2015.35.2079.
[51] Christian Offen and Sina Ober-Blöbaum. Symplectic integration of learned Hamiltonian systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(1), 2022.
13

---

[52] G. R.W. Quispel and H. W. Capel. Solving ODEs numerically while preserving a first integral. Physics Letters, Section A: General, Atomic and Solid State Physics, 218(3-6):223-22)) 1996. ISSN 03759601. doi: 10. (016600775-660196))))))
[53] G. R.W. Quispel and D. I. McLaren. A new class of energy- preygevne numerical integration methods. Journal of Physics A: Mathematical and Theoretical, 41, 2008. ISSN 17518113. doi: 10.1088/1751-88 13/41/4/045206.
[54] M. Raissi, P. Perdikaris, and G. E. Kamiaadkii. Physics-informen neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378, 2019. ISSN 10902716. doi: 10.1016/j..p.20018 10.00.4
[55] Rishikesh Ranade, Chris Hill, and Jay Pathak. DiscretizationNet: A machine-learning based solver for Savier-Stoked equations using finite volume discretization. Comput. Methods Appl. Mech. Eng., (78(113722):
May 2021.
Jesus M Sanz-Serna. Symplectic integrators for Hamiltonian problems: overview. Acta [56]
numerica, 1 1:243-286, 1992.
[57] Jesus M Sanz-Serna and Mari-Paz Calvo. Numerical Hamiltonian problems, volume 7. Courier Dover Publications, 2018.
[58] Khemraj Shukla, Ameya D Jagtap, and George Em Karniadakis. Parallel physics-informed neural networks via domain decomposition. Journal of Computational Physics, 447:110683, 2021.
[59] Naoya Takeishi, Yoshinobu Kawahara, and Takehisa Yairi. Learning Koopman invariant subspaces for dynamic mode decomposition. Advances in neural information processing systems, 30, 2017.
[60] Paulo M Tasinaffo, Luiz A Dias, and Adilson M da Cunha. A survey about universal numerical integrators (UNIs): Part II or quantitative approach. Hum-Cent Intell Syst, November 2024.
[61] Paulo M Tasinaffo, Luiz A V Dias, and Adilson M da Cunha. A qualitative approach to universal numerical integrators (UNIs) with computational application. Hum--ent Intell Syst, (4):571-598, November 2024.
[62] Yunjin Tong, Shiying Xiong, Xingzhe He, Guanghan Pan, and Bo Zhu Symplectic neural networks Taylor series form for Hamiltonian systems. Journal of Computational Physies, 437:110325, 2021.
[63] Enoch Yeung, Soumya Kundu, and Nathan Hodas. Learning deep neural network representations for Koopman operators of nonlinear dynamical systems. In 2019 American Control Conference (ACC), pages 4832-4839. IEEE, 2019.
[64] Huaiqian You, Quinn Zhang, Colton J Ross, Chung-Hao Lee, and Yue Yu. Learning deep implicit Fourier neural operators (IFNOs) with applications to heterogeneous material modeling. Computer Methods in Applied Mechanics and Engineering. 398: 115296, 2022.
[65] Chengxi Zang and Fei Wang. Neural dynamics on complex networks. In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pages 892-902, 2020.
[66] Zhen Zhang, Yeonjong Shin, and George Em Karniadakis. GFINNs: GENERIC formalism informed neural networks for deterministic and stochastic dynamical systems. Philosophical Transactions of the Royal Society A, 38002229):22210007,, 2022.
[67] G Zhong. Lie-Poisson integrators and Lie-Poisson Hamilton-Jacol theory. Phys. Lett. A, 133: 134-139, 1988.
14

---

A Proofs
Proof of Theorem 3.7 By using the chain rule, holds that
dh
du
VH
dt
dc
Substitution of the equation yields
dH
=
(if S is skew symmetric)
VH
SVH
dt
< 0
(if O.)
Proof of Theorem 4.2 The proof is discrete analogue of the proof of Theorem 3.1 discrete chain rule, holds that
By using the
H H(u(+1)) - H((ln)) AH
(((+1) (())
h
Substitution of the discrete gradient scheme yields
H H(ln+1)) Hul()) = =VH SVH
(if 3 is skew symmetric)
=
h
I 0
(if S < O.)
Proof of Theorem 4.3 First, note that a vector w is perpendicular to another vector v - u, then there exists skew symmetric matrix A such that = A(v u). Hence, the proof is completed if we show that all discrete gradient operators have the following representation:
VH(,, u) Hu, u) + w(u, v), where w(u, v)
u for all u, Rdd
Let V VH(uu v) be of the form
(8) for some w(u, v). and we check that it satisfies (5) and (6
(VH(u, v), v-u) = (Voll(u, v), v- + v), v- = (V, v), = H(v)-H(u))
By definition, AH (u, u) VoH (u. u) + w(u, u) VH(w) (u). Conversely, if HHHHH, v) is a discrete gradient, then define w(u. v) = VH(u.v)
(u, w) which is continuous. By 7 w(u, 1 u. and w(u, u) (H u) V.H (u, u) by since both of the discrete gradients are consistent.
Proof of Theurem
Suppose that numerical scheme
(n+1) (n)
h f(u(a+1),
(11)
preserves the energy
H(u()(1) H(ul))).
Let (n+1) and (m)) be solutions to (II). is sufficient to show that f can be written as fu, S (V, w(u, v)), where w(u,v) 1. with reference discrete gradient VoH (u, v) of H.
From the property of the discrete gradient and the energy conse-conseervation property, we get
H (u(a+1)) - (u())) Voll(u((t)) uni T ((++1) (1)
Voll(u (n) T f(u(nt2), (n))..
h
Hence / must be perpendicular to the discrete gradient V.H. Because we assume that the target system is a Hamiltonian equation, the matrix S is invertible. Therefore, we have
LSVol (m (n)))) T S-1 f(u(x+1) ()) Voll(u(nt1) ( . ST (()((t))) (x Vo Hu (x+11) ((n))) ((u(1)1), vn(t)) 0
(nl)
15

---

Thus, f is in the orthogonal complement M of SV Because S is snee-syymmetiic VoH is perpendicular SV which means that VoH in the orthogonal complement M. Therefore. we can obtain orthogonal hasis (e1, ,,.-.., 6N) of M, where VoH and N = dim M. By using this basis, f can be written as
S S = = VoH +e
where 6j Ik and span(ee, ... EN.. We want to show that is perpendicular to wln+1) u(n) hf(laat1) (. (x))) The inner product of and f is - -u(1) -u( () = h(f.e) h(f, - o1 = h(f,5-1s) - cahls, VoH) =0
because S-1 is skew-symmetric, ( = Hence (S. O. Thus e is perpen-
dicular to u(++1) u(m), so there exists skew-symmetric matrix A(u(at1), ((()) such that A(u(n+1), (u())))( (
Finally, we show that Ci = 1. To this end, we take the limit h - in
s = c VoH +e V.H A(u
(")).
Since the scheme is consistent, f / as h - 0. On the other hand, the right-hand side tends èi H because of the consistency condition ( of discrete gradients. Hence we get V 41 H and hence = 1.
B Geometric Description of Target Differential Equations
The target equations in this study are characterized as gradient flows on symplectic Riemannian manifolds. Let (M,w) be a symplectic manifold and (M,g) be a Riemannian manifold. For symplectic manifolds, w is a symplectic 2-form. In Riemannian manifolds, denotes the inner product. These define a bilinear function Wa TuM x T.M - IR or 9u T.M x T.M - R for each M. It is assumed that and g are nondegenerate in the sense that the matrix M, representing these bilinear functions is nondegenerate for any E M.
For given energy function H M IR, we can define differential equations
ύ= X, -
on symplectic manifolds and
== x, ((X,)) -dH(()
on Riemannian manifolds, where dH TM R is the differential of H. For symplectic manifolds, this equation is the Hamilton equation, and for Riemannian manifolds, it is the gradient flow.
These equations have energy conservation and the energy dissipation laws, respectively. In fact, in Hamiltonian systems, it holds that
dH dt dH(X) = w(X
because wi is skew-symmetric due to the property of symplectic forms. Also, for Riemannian manifolds,
dH dX( -xxx 50
di
By using the relation between the gradient and the derivative
dX(X) = ( XV
these equations can be rewritten as vector-matrix representations:
du
M, H
di
for symplectic manifolds and
du
- M, VH
di
for Riemannian manifolds, respectively. These equations the target equations of this study.
16

---

C Target Partial Differential Equations and Its Semi-Discretization
The proposed method in this paper also be applied to semi-discretized partial differential equations in the spatial direction. Typical examples of the target partial differential equations are equations of the following form
Ou Ōt
8H
(12)
dx
5u
where is the energy density, which is function of u, Ur. Maz, and so on. uh denotes Ou and Uaz dx
denotes 8Pu Da S. f/u is called the variational derivative of H. For example, if I/ is function of and Uz then this is defined by
811 du
01
0H
8u
dx ou,
Equations of the form 12) include many partial differential equations, such as the advection equation, the KdV equation, the Burgers equation, the Allen-Cahn equation, and the Cahn-HHiliard equation. Equations that are slightly different but can be expressed in a similar form include the Maxwell equation, the wave equation, the Klein-Gordon equation, and the nonlinear Schrodinger equation. Although the above equations are defined in one-dimensional space, they can be extended to multiple dimensions in a straightforward way. See, e.g.. Furihata and Matsuo [8 for details.
is known that 112) has the global energy conservation law or the energy dissipation or increasing law under certain boundary conditions, such the periodic boundary condition. Suppose that the equation is defined on the interval [0, L). For simplicity, we assume the periodic boundary condition. In this case, the time differentiation of the total energy yields
J di J
d L Hdr = S0
L
OH Ou
J/"
OO Ou
δ 0H
OH Uu
L
Ou
)
Unt
dx
Ox Ou.
u.dæ
Ou.
Because we assume the periodic boundary condition, the boundary term must vanish. Therefore, we have
a % Hdz -
di [Hude
( O OH
u,da.
&u
Oe Ou.
Substituting the definition of the variational derivatives and also the equation into the above equality, we obtain
d
L' Hdr (-
(OO du
[H
θ OH
6"
8H du
8H
u.dr
dr.
dt
dr du.
dr
8u
If is odd and expressed we get
2n with an integer n, then by repeating the integration by parts,
1 5H
à 2n
SH
dr = (-1)" 1"
8H
5H
da
du
dx
du
Ex
8u
dx
du
= (-11ntt f"
n+1
8H
0
8H
dx.
Ox
8u
Ör
du
Therefore, we have
(-1)a /
8H
a+1
8H
#+1
8H
0
8H
da
- (-1)
dx
0x
5u
Ör
ou
Ür
8u
Fx
8u
and hence the energy conservation law follows:
5H
n+1
8H
2(-1)"
u =
dr
du
Ox
du
Similarly, when a = 2n, we have
/
1
6H
2n+1
8H
// 1
[
"8
urd = (-1)a
udr.
8u
dx
8u
dz
8u
17

---

Thus, the energy decreasing or increasing law follows. Thus. it can be seen that these equations have the same energy behavior the target ordinary differential equations of this study.
Next, we semi-discretize these partial differential equations only in the spatial direction using the finite difference method while preserving the energy behaviors. We devide the interval [0, Z.] into M equal parts with the spatial step size a. The approximation of utt, j(rr is denoted by uy(t), and we donote the vector (uo(t), ..., (t)) as Let D1. D4, and D be difference matrices representing the central, forward, and backward differences, respectively:
1
1 Di 2Ar
1
D+
D
2Ar
2Ar
9
First, the energy function can be approximated by approximating us in the energy function by using the finite difference
J
Hu,,uy)r u2)dr Ha D D4 D =: H.)
We also discretize the differential operator
(13)
Ox
in the equation. For = 2n, we can use n copies of D and D_ to form a matrix
D D D_ D_
to approximate (3).. This matrix positive or negative semi-definite, since D = -D_. Thus, the semi-discretized equation
du dt = D......DD.. ... ...D....D.VAA. (H)
is an ordinary differential equation with the energy dissipating or increasing property. Similarly, if a 2n + 1, then the equation is semi-discretized into an ordinary differential equation with an energy conservation law:
dū dt
D4 ... D
Note that D4... D,D ... D_ is skew-symmetric because D2 is skew-symmetric. The proposed method yields numerical integrators that preserve the energy behavior for each of these semi- discretized equations.
D
Supplementary experiments
D.1 hyperparameters
Table 2 shows the hyperparameters for training.
D.2 Experiment loss
Loss under Different Reference Methods with Varying Network Depths and Widths the Hénon-Heiles model:
We use
1 H(p. - (p) p2)
(of +4 o2
-
The training flow data (1((1-1), (())) n 1,..., N) consists of N 100 pairs of points obtained from the trajectory calculated using a high-order integrator sampled with time step h = 0.3, starting from (0) (41) 02 "P,", .p,) (0.3, 0.3, -0.2, -0.2). After training, we use the flow starting
(0)
at (N) for 300 steps test.
18

---

Table 2: hyperparameters
PROBLEM
LR
EPOCHS
LAYER
WIDTH
PENDULUM
-FLOW
0.0001
10000
50
-IRREGULAR
0.0001
10000
50
-DISSPATIVE
0.0001
10000
50
HÉNON-HEILES
-FLOW
0.00001
30000
10
100
-IRREGULAR
0.00001
30000
10
100
2-BODY
0.0001
10000
50
2- 2-PENDULUM
0.0001
10000
50
DUFFING
0.0001
10000
50
KDV
0.0001
10000
200
width
0.000001
50
100
AVF2layer
AVF2-test
AVFSlayer
AVFS-test
0,0000001
AVF1Clayer
AVF10-test
ITHH2layer
ITOH2-lest
loss
1E-08
ITHHSlayer
-.. ITOHS-test
IOH10llseer
ITOH10-test
SYMITOH2layer
SYMITOH2-test
1E-09
SYMITOHS-test
SYMITOH10lajer SYMITOHID-test
1E-10
Figure 11: Loss for Different Reference Methods with Varying Network Depths and Widths
We conducted experiments using both AVF, Itoh-Abe method and time-symmetric Itoh-Abe method as baseline schemes. The learning rate was set to 0.00001, and the number of epochs was 30,000, The network architecture was varied with 5, and 10 layers, and width of 10, 50, and 100 per layer. For each configuration, we performed five independent training runs and reported the final loss as results, as shown in the Figure 11.
As expected, the AVF-based method, being second-order scheme, generally outperformed the Itoh-Abe method.
Based on the above experimental results and observations, we conclude that in the initial selection of the baseline method, the AVF scheme should be prioritized whenever it is computationally feasible, due to its higher accuracy and time symmetry. In cases where the original Hamiltonian system is too complex for the AVF discrete gradient to be computed, the Itoh-Abe method can serve as practical alternative. Furthermore, it is worth exploring the possibility of constructing a time-symmetric version of the Itoh-Abe method enhance its accuracy, potentially achieving second-order scheme.
For this reason, we adopt time-symmetric construction Itoh-Abe method to achieve second-order. And the results in the Figure 11 indicate that the time-symmetric Itoh-Abe method outperforms the standard one.
V HtchAbeeyyHHlu,, (VtohhbbeHl VithAbe u)
19

---

Experiment loss
The training loss was computed by the mean-squared-erroo between Mirain Shaann
Prain
and SV NNNN H. (rraii " J, and the test loss computed by the mean-squared-errro between and NY (u flest" () Table 3 shows the losses in the form of mean + standard deviation measured from 5 independent experiments. Note that all computations were performed in double precision. It can be seen that the models were trained very accurately.
(n)
(16)
Rigus
Table 3: The values of the train and test loss.
PROBLEM
TRAIN LOSS
TEST LOSS
PENDULUM
-FLOW
7.75 x 10 10 + 7.46 x 10 10-10
1.95 x 10 + 2.21 X 10
-IRREGULAR
1.33 10- 4 1.07 10-?
6.08 x 10 + 3.95 10-7
-DISSPATIVE
4.69 x 1 2.06 x 10
1.70 x 10 9 x 10
HÉNON-HEILES
-FLOW
1.19 10-10 10 11
1.72 10
10-10
IRREGULAR
1.40 x 10 x 10-9
2.14 x 10-6 x 10-9
2-BODY
1.14 10-8 10
2.92 10-
10
2- 2-PENDULUM
6.34 x 10 1..06 10
8.69 x 100 x 10
DUFFING
9.24 100=0100--110 1.95
10
7.73 x 1 5.91 10
KDV
1.09 x 1 6 = 1.25
x 10 10 10
1.69 x 10 = 1.61 x 10 10
Experiment loss with different data and initial conditions We conducted additional experiments on both the Ménon-Heiles, and the double pendulum. The aim was to evaluate the impact of training data distribution, data size, and initial conditions on the performance of the learned integrator.
Hénon-Heiles system: We used irregular training data (((m((11), (,()),, (0) | .,,,, I}, where the elements of (ff-1) are randomly sampled from the domain [-1, 1] and each step size h6 is drawn independently from [0, 0.5). Each ((()) represents the one-step solution computed from (i- using step size (+) with high-order solver.
We performed experiments using 100, 1000, 5000 training samples, respectively. The neural network consists of 10 layers, each with 100 width. We set the learning rate 0.00001 and trained the model for 30000 iterations. For testing, we selected random initial conditions from [-0.4, 0.4]4, For each test trajectory, a fixed step size h was randomly chosen from [0, 0.5) and applied over 100 time steps.
Double pendulum: We generated training data as follows. First, we randomly sampled initial angles (gl. q2) from the interval [-0.8, 0.8] and set the initial momenta (p1, p2) to zero. For each initial condition, we randomly chose fixed step size from the interval [0, 0.5]. and integrated the system forward for 100 time steps using high-order solver. This process was repeated to generate 500 trajectories, each of 100 steps. From the resulting collection of one-step pairs, we randomly shuffled and selected 100, 1000, 5000 samples to construct training datasets. The neural network consists of 10 layers, each with 100 width. We set the learning rate to 0.0001 and trained the model for 30000 iterations.
During testing, we followed the same protocol as in training: each test trajectory starts from randomly chosen initial condition, with angles in [-0.8, 0.8] and momenta set to zero. A fixed step size h E [0, 0.5) was randomly selected per test and kept constant over 100 steps.
Table shows the average training and test loss with standard deviations for different training dataset sizes. Fach result is computed based 5 independent training runs with different data. For both systems, the training loss consistently converges to a similar small value and the test loss improves as the number of training samples increases.
D.3 Pendulum Experiment
We conduct experiments on the pendulum (Figure 12.. The mathematical pendulum (mass m string length 1 = 1, gravity = 1 is system with the Hamiltonian H(g. p) = 1/22 -cos(q).
20

---

Table 4: The values of the train and test loss with different data and initial conditions.
PROBLEM
TRAIN LOSS
TEST LOSS
HENON-HEILES
-100 TRAINING DATA
1.64 x 10 6.09 x 10-3
3.43 x 10 1.98 x 10-6
-1000 TRAINING DATA
1.83 10 2.66 10
5,29 x 10 4 x 10
-5000 TRAINING DATA
1.90 x 10 x 10-9
1.58 x 10-9 +1.6 x 10-9
2-PENDULUM
-100 TRAINING DATA
2.09 10 + 10-7
1.22
10
1000 TRAINING DATA
3.44 x 10-1 x 10-7
6,43 x 10 1 10
5000 TRAINING DATA
3.09 10 + 10
4.24 10 3.52 x 10
The position coordinate q denotes the angle from the vertical suspension point and p = denotes the momentum. Although the trigonometric function exists in the Hamiltonian, the integral in the AVF method is analytically computable, and hence the AVF method is available.
Figure 12: The Pen- dulum
In this experiment, we used two sets of data: flow dataset and irregular dataset like the example shown in the Figure [3
The training flow data consists of N = 15 points sampled from trajectory for each fixed time step h = 0.5, calculated using a high-order integrator starting from (( p())) = 0). The training flow data are pairs of numerical solutions at certain time step along with the one at the next time step: {n((==), (n))) = N.. After training, we use the trained gradient to compute the flow starting at (NN for 200 steps.
The training irregular data ((u1(111) M.), ( 11=1.... 1) consists of T 100 pairs of points. (f-1) are randomly sampled from [-3,3] x [-3, 3] and the time steps (() are randomly chosen from [O, 0.8). x(9) are numerical solutions at the next step of (i-1) with the step size (0) After training, we use the trained gradient compute the flow starting (0) (2, 0) with 0.5 for 200 steps.
Figure 14, we show the additional for the pendulum experiment. The plot in the top-left corner shows the irregular training data. The top-right plot displays the test trajectories of the trained network and numerical integrators during the test. The two middle plots show the values of and over time. The bottom-left plot shows the error evolutions over time, while the bottom-right plot illustrates the energy conservation.
Figure 13: Flow and irregular data with 15 pairs of points
In Figure 15, we show the additional results for the pendulum experiment with the irregular data. The plot in the top-left corner shows the irregular training data. The top-right plot displays the test trajectories the trained network and numerical integrators during the test. The two middle plots show the values of y and p over time. The bottom-left plot shows the error evolutions over time, while the bottom-right plot illustrates the energy conservation.
D.4 Supplementary figures of the experimental results
Results for the Hénon-Heiles system In Figure 16, we show the additional results for the Hénon- Heiles experiment. The plot in the top-left corner shows the positions of the training flow data, while the top-middle displays the values of test 91. The plot in the top-right shows the change in error over time, and the center-right shows energy conservation. The remaining plots display the test trajectories of the ground truth, the solution by the AVF method, that by the AVF method with a half step size, order, and the solution by the proposed method.
In Figure 17, we show the additional results for the Hénon- Heills experiment with the irregular data. The plot in the tef-left corner shows the positions of the randomly selected irregular training data, while the top-middle displays the values of test qı- The plot in the top-right shows the change in error over time, and the center-right shows energy conservation. The remaining plots display the test trajectories of the ground truth, the solution by the AVF method, that by the AVF method with a half step size, 4 order, and the solution by the proposed method.
21

---

Traindata
Pendlum
Graundrruhh
Prosesed methos
Avl
1-
Ground bruth
Grond South
Troposed meshod
Proposed methon Ail
0-
Al
AUT chopa
Art hal soppstes ordure
Energy
stopsize
0.4164
0.4162
0.2
0.415
0.0
40
200
40
200
Figure 14: Results for the pendulum experiment
Traindata
Pendium
Propoeed
Propered
Aue
Art hall starize Avt
Avd hut saassize ArF ordere
200
Energy
Proposed merred
0.41.6
Gound tuth
Proposed methos
stepstice
0.4
0.4162
0.41.0
0.415
Figure 15: Results for the pendulum experiment with the irregular data
Results for 2body In Figure 18, in addition to the numerical errors and the evolutions of the energy mentioned the main text, the top-left corner shows the training flow data, while the top-middle display the values of test Q1- The plot in the center-right shows energy conservation.
Results for Double pendulum In Figure 19, the plot in the top-left corner shows the positions of training flow data, while the two plots at the center-right show the evolutions of the emoor the energy over time. The bottom plots display the trajectories of the ground truth, the solution by the AVF method, the solution by the AVF method with a half step size, and the solution by the proposed method.
22

---

Results for the dissipative pendulum In Figure 20, in addition to the numerical errors and the evolutions of the energy mentioned in the main text, the top-left corner shows the training flow data, the top-right corner displays the test trajectory after training, the middle display the values of test q. and the two plots the bottom correspond the values the coordinates over time.
Results for the Duffing oscillator In Figure 21, we show the results for the Duffing oscillator. In addition to the numerical errors and the evolutions of the energy mentioned in the main text, the top-left corner shows the training flow data, the top-right corner displays the test trajectory after training, the middle display the values of test q; p. and the two plots the bottom correspond the values of the coordinates over time.
consider this equation on the interval [0, Lj, = 20 under the periodic boundary condition. We used the training data with the time length of T = 30 and the time step size of Át = 0.01. The spatial step size was Ax = 0.5. The initial condition data is constructed by superimposing two solitary wave solutions of the KdV equation: uz(2) Va 21, c1) Ueol (z 22, c2). where (I, c) f sech ( vō x)) with 0.75, 21 0.33 C2 = 0.4, 21 = 0.65L. The time evolution is computed by solving the KdV equation using the solver with high accuracy, where spatial derivatives were approximated using finite difference methods. For testing, the computation was carried out from the time length T 30, with the start point being at the final state of the training dataset.
Results for KDV
the training
In Figure 22, we show the results for the Kdv. The plot in the top-left corner shows the positions of the training flow data. The plot in the bottom-left shows the change in error over time, and the pottom-middie shows energy conservation. The remaining plots display the test trajectories of the ground truth, the solution by the AVF method, and the solution by the proposed method.
Traindata
Q1
Global error
0.2
0.75
0.0
2.50
0.75
-0.4
9.5
0.0
0.5
-10
50
75
5.00
50
75
AVI half stepsioe
AVT ordend
Energy
Ground buth
9.59
0.50
Propesed mothod
AVT
0.25
0.34860
0.60
0.00
-3.5
0.0
75
22
Ground truth
Proposed method
AVF solution
0.53
0.53
3.50
0.50
0.03
5.00
-0.59
0.50
-0.50
61
Figure 16: Results for the Hénon-Heiles experiment
23

---

Trannaaa
Global error
0.3
MW
-0.5
-1.0
-1.
AVF half stepsice
orderd
Energy
0.59
0.59
0.14803
Proposed
AVT
0.14802
6.00
0.03
-0.50
-0.53
0.14793
Ground
Predicted
0.75
0.75
0.50
0.52
0.52
2.2
0.25
0.25
0.00
0.02
0.02
-0.5
-0.50
-0.53
Figure 17: Results for the experiment with the irregular data
Taindata
Global error
0.50
05-
0.25
0.00
0.0
3.50
-0.5
0.0
0.5
1.0
10
20
10
ItekwAbe half stepsize
Roh-Abe orderd
Energy
9.5
05
-0, 74095
method
(buh Abc
0
-0.75003
-2.5
-05
Ground truth
Proposed method
Itch-Abe solution
0.5
05.-
0.5
-0.5
-05
-0.5
-3.5
0.5
Figure 18: Results for 2body
24

---

Traindata
Global error
0.5
0.3
0.9
-10
75
Energy
Itot-Abe half stepsize
Roh-Abe ordend
RohiAze
3.
-0.25
-0.00
0.5
-0.35
-0.75
Ground
Proposed method
solution
0.75
0.25
0.75
-0.50
-055
-0.00
-5.75
-0.75
-6.75
Figure 19: Results for the double pendulum
Traindata
Dissloative Rendum
Eropored
W
Ground buth
Cround buth tropesed metsed Aye AvF nalf steas AVF onde
Prosesed mnchol
AlF
zoF Fal stepsiee
AVT oroent
-1-
Energy
3.25
0.10
3.00
-3.05
0.00
15 20
Figure 20: Results for the dissipative pendulum
25

---

Traindara
Test
mmmmio Prososed methad JVF
Prpposed methhe
HAT ha
AVT har
HNT Sorder
AVT dcrter
Energy
Propesed
Ground
1.5
0.0
Figure 21: Results for the Duffing oscillator
Traindata
533
5.
10
20
Ground truth
Predicted salution
AVT salution
23
23
02
0.0
Energy
+0.3200-
mothod
0.035
-0.3105
Proposed mottled
-0.3700
-0.3355
0.000
Figure 22: Results for the KdV
26

---

D.5 Comparison with other methods
We implemented the PINNs as [54] with energy conservation/diisipation regularization terms, using the DeepXDE library, and applied it to three benchmark systems: Hénon Hénon-Heiles, double pendulum, and Duffing oscillator. Both our method and the PINNs use the same network architecture (10 layers, width 100), optimizer (Adam, learning rate 0.0001), and 30000 training iterations. We sampled 2000 random collocation points over the spatiotemporal domain.
We also implemented the Time-Reversal Symmetric ODE Network (TRS-ODEN) 28 For the time-reverse symmetry loss term, we conducted experiments with three different weight values:) 10, 0.5, which are employed in the reference paper. The network architecture (10 layers, width 100), optimizer (Adam, learning rate 0.0001), and 30000 training iterations are kept the same.
To evaluate the performance, we generated ground-trute reference trajectory for each system using high-order solver with 100 time steps, starting from the same initial conditions. We then simulated trajectories using both the trained PINNs, TRS-ODEN and our proposed method under the same initial condition. We report the maximum error over the time domain (0, T), defined as the maximum deviation from the ground-truth trajectory. Furthermore, since PINNs tends to degrade over long time horizons, we computed the maximum error accumulated up to each checkpoint (T' = 5, 10, 50, 100), i.e. the largest deviation from the reference trajectory from step 0 up to the given time step T'
Hénon-Heiles System We set the initial condition as (0) (0.3,0.. 0.3, -.2, -0.2), and performed simulations over 100 steps with a fixed time step size = 0.3.
(
(0)
(0)
42
P1
P2
Double Pendulum We set the initial condition as u = (4i
.p2 = (0/6, 0/0.0)
92
and performed simulations over 100 steps with fixed time step size h 0.3.
Duffing Oscillator We set the initial condition as u (0) = (g p0)) (2.95, -3.08) and per-
formed simulations over 100 steps with fixed time step size h = 0.1.
The maximum errors of both models, measured against the reference solution, are summarized in table
These results in all three systems suggest that our structure-preservingg method offers improved accu- racy and stability over the other two methods, especially in long-term simulations. This improvement attributed its formulation as a numerical integrator designed to respect the system's 's physical structure.
Table 5: Comparison of Maximum Errors
TRS ODEN
TIME STEP
PROPOSED METHOD
PINNs
À
==1
0.5
HÉNON-HEIL
5
1.43 x 10
2.86 X 10
9.11 x 10
3.72 x 10
2.37 x 10
10
1.54 x 10
9.77 x 10
2.16 x 10
1.24 x 10
4.35 x 10
50
1.10 x 10
1.48
7.40 x 10
3.16 x 10
5.21 x 10
100
9.53 x 10-3
1.50
1.32
3.26 x 10
1.43 x 10
PENDUULUM
5
1.72 x 10
1.92 x 10
1.32 x 10
1.08 x 10
9.58 x 10
10
8.28 x 10
5.39 x 10
1.85 x 10
1.18 x 10
9.58 x 10
50
3.50 x 1
1.98 x 10
L.1 X 10
3.85 X 10
4.37 x 10
100
6.18 x 10
2.38 x 10
2.64 X 10
7.59 x 10
1.24 X 10
DUFFING
5
8.40 x 10
4.20 x 1
1.4 X 10
3.90 x 10
2.96 X 10
10
9.38 x 10
6.99 x 10
2.06 x 10
3.90 X 10
3.55 x 10
50
5.18 x 10
4.50
2.86
5.22
3.33
100
1.12 X 10
4.50
6.40
5.38
5.12
27

---

D.6 Computational Time
We performed additional experiments re-evaluate the computation time across different systems. In the experiments, we used an Intel Xeon 6900P CPU.
The results, the size of the neural networks used and the number of time steps in the numerical calculations are shown in 6 The computation time is generally similar to that of the Hénon- Heiles problem reported in the main body. In particular, the proposed method achieves better accuracy with less computational time than the fourth-order methods. For the KdV equation experiment, the network was made as small as possible without sacrificing accuracy. The network used is quite small, yet the test loss was 1.706 x 10 6, which is not much different from the test loss of 1.695 x 10-6 when larger network with 5 layers with width of 200 is used. For other experiments, the computation time can be further reduced by making the network smaller.
Table 6: Comparison of Computational Time (sec)
Problem
Layer
Width
Steps
Proposed
2nd Order Methods
4th Order Methods
Pendulum
200
0.579
0.503
1.137
Hénon-Heiles
10
100
300
2.090
1.437
3,533
2-Body
5
50
300
2.862
2.443
7.693
2-Pendulum
5
50
300
5.768
5.197
18.043
Duffing
5
50
300
1.020
0.685
1.084
Kdv
2
30
3000
6.566
1.634
8.781
28

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claims in the abstract and introduction clearly reflect the key contributions and findings of the paper, and are well supported by both the theoretical analysis and experimental results.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: One limitation of the method lies in the increased computational cost due to the use of neural networks.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification. asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
29

---

Answer: [Yes]
Justification: The assumptions required for our theoretical results are described in the main paper, and the proofs are provided in the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification:
We provide detailed experimental settings and the code.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
30

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
We have provided a ZIP file containing the experimental code.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips //ips.cc//pill wttpcccpppppppppeeeeeeeeeeeeeeee httpa///giiieee/eeeeebiiisssssss /Cc////oddiiisssseeeeeee ssinsoliccyyyyyc.ccc for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification:
We provide detailed experimental settings and the code.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We provide detailed experimental losses in the form of mean + standard deviation measured from independent experiments the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for
example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
31

---

The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (c.g. negative emror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We include description the hardware.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't 't make it into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer:
[Yes]
Justification:
We adhere to the NeurlPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: There is no socictal impact of the work performed.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
32

---

The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, scraped datasets)?
Question:
Answer: [NA]
Justification: The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We used the official code of SympNet. Although the information about the license is not available online, properly acknowledged the code.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the Datasets.papersstt paperswi thcode. cww/ datasets package should be provided. For popular datasets,
has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
33

---

If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not introduce any new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question:
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and expect authors to adhere to the NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
34

---

Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that if the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology. scientific rigorousness, or originality the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does not involve LLMs any important, original, or non-standard components.
Guidelines:
The answer NA
that the method development in this research does involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy for what should or should be described.
35


Reviews:
[{"review_id": "xcziMwHX8h", "reviewer": "Reviewer_gBGZ", "summary": "This paper introduces a method for learning numerical integrators from data for deterministic energy-conserving or dissipative systems. The key idea is to use the neural network to parametrize the discrete gradients in the classical implicit discrete-gradient methods. The learned discrete gradient is adaptive to the time-stepping size of the training data; thus, it can achieve large time-stepping while keeping the accuracy. Numerical results show that the proposed method can learn accurate large time-stepping integrators for various differential equations, including chaotic Hamiltonian systems (the Henon-Heiles Model and the two-body problem, the double pendulum), dissipative pendulum, and the KdV equation.", "questions": "- What is the matrix $\\bar{S}$ in Eq.(5)? Is it the same as the matrix $S$ in Figure 1? How to derive or approximate $S$ or $\\bar{S}$ in applications? \n \n- Line 263: is $u^{(i)}$ a typo for $u^{(i-1)}$? Is there any reason for selecting such random training data?\n\n- Line 306+: the natation x is a typo of q. What if the system is stiff?\n\n- KdV example: What is the training data, and how is it generated? How does the training data affect the performance of the proposed model?", "limitations": "yes", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths: \n- The paper proposes a robust method for learning a data-driven numerical integrator that incorporates energy conservation laws into the design of neural network approximations.   \n- The paper validates the proposed approach in a diverse set of systems, including nonlinear Hamiltonian systems and dissipative systems. \n- The paper is well-organized, and the main idea is clearly presented. Additionally, the paper provides a detailed description of the experiments and the code, ensuring the results are reproducible. \n \nWeakness: \nThe paper is solid overall. The following are minor weakness points: \n- Lack of statistical robustness tests. Since the proposed integrator is learned from data, it is important to discuss its dependence on training data, including data size and the methods used to generate training data (single or multiple trajectories, which can be interpreted as the distribution of the training data). Additionally, it is crucial to test the learned integrator with different initial conditions, not just at the end of the training trajectory. Then, report the mean and standard deviation of the errors (in time, as shown in the last plot of Figure 3) for each example.\n- Missing computational time for the examples, except the first example. Computation efficiency is a major motivation of the proposed method over classical methods (which use a smaller time step size). However, the complex NN approximation will take more time to compute the implicit integrator. It is unclear if the proposed method is more efficient. Therefore, it is important to report computation time for each example.  \n- Theorem 3.1 and Theorem 4.2 are more like definitions or simple facts. They are not sufficient to be called a theorem. In particular, the last sentence of Theorem 4.2 is not clear enough to be a statement in a theorem. \n- The KdV example lacks the details, and the results are rough.", "was_revised": true, "final_justification": "This is a technically solid paper, presenting a robust method for learning a data-driven numerical integrator that incorporates energy conservation laws into the design of neural network approximations. The additional tests and comparison with the baselines further improve the paper. Thus, I keep my rating.", "rebuttal": "We sincerely thank the reviewer for the valuable comments and suggestions. We have carefully addressed the concerns raised regarding statistical robustness, computational efficiency, the clarity and rigor of the theorems, and the completeness of the KdV example. We believe the revisions and clarifications provided below improve the overall quality and rigor of the paper. We respond to each point in detail as follows.\n\n- **Lack of statistical robustness tests. Additionally, it is crucial to test the learned integrator with different initial conditions, not just at the end of the training trajectory. Then, report the mean and standard deviation of the errors for each example.**\n\nTo address the concern, we conducted additional experiments on both the Hénon-Heiles and the double pendulum. The aim was to evaluate the impact of training data distribution, data size, and initial conditions on the performance of the learned integrator.\n\n**Hénon-Heiles system**: We used irregular training data {$((u^{(i-1)},h^{(i)}),u^{(i)}) \\mid i=1,...,I$ }, where the elements of $u^{(i-1)}$ are randomly sampled from the domain $[-1, 1]$ and each step size $h^{(i)}$ is drawn independently from $[0, 0.5]$. Each $u^{(i)}$ represents the one-step solution computed from $u^{(i-1)}$ using step size $h^{(i)}$ with a high-order solver. \n\nWe performed experiments using $100$, $1000$, and $5000$ training samples, respectively. The neural network consists of $10$ layers, each with $100$ width. We set the learning rate to $0.00001$ and trained the model for $30000$ iterations.\nFor testing, we selected random initial conditions from $[-0.4, 0.4]^4$, including high-energy cases with total energy exceeding $1/6$. For each test trajectory, a fixed step size $h$ was randomly chosen from $[0, 0.5]$ and applied over 100 time steps.\n\n**double pendulum**: We generated training data as follows. First, we randomly sampled initial angles $(q_1, q_2)$ from the interval $[-0.8, 0.8]$, and set the initial momenta $(p_1, p_2)$ to zero. For each initial condition, we randomly chose a fixed step size $h$ from the interval $[0, 0.5]$, and integrated the system forward for 100 time steps using a high-order solver. This process was repeated to generate $500$ trajectories, each of $100$ steps. From the resulting collection of one-step pairs, we randomly shuffled and selected $100$, $1000$, and $5000$ samples to construct training datasets. The neural network consists of $10$ layers, each with $100$ width. We set the learning rate to $0.0001$ and trained the model for $30000$ iterations.\n\nDuring testing, we followed the same protocol as in training: each test trajectory starts from a randomly chosen initial condition, with angles in $[-0.8, 0.8]$ and momenta set to zero. A fixed step size $h \\in [0, 0.5]$ was randomly selected per test and kept constant over 100 steps.\n\n| Problem       | 100 train data                              | 1000 train data                             | 5000 train data                             |\n|---------------|---------------------------------------------|---------------------------------------------|---------------------------------------------|\n| **Hénon-Heiles** |                                             |                                             |                                             |\n| - training loss   | $1.64 \\times 10^{-8} \\pm 6.09 \\times 10^{-9}$ | $1.83 \\times 10^{-8} \\pm 2.66 \\times 10^{-9}$ | $1.90 \\times 10^{-8} \\pm 1.30 \\times 10^{-9}$ |\n| - test loss    | $3.43 \\times 10^{-8} \\pm 1.98 \\times 10^{-8}$ | $5.29 \\times 10^{-9} \\pm 4.45 \\times 10^{-9}$ | $1.58 \\times 10^{-9} \\pm 1.60 \\times 10^{-9}$ |\n| **2-Pendulum** |                                             |                                             |                                             |\n| - training loss   | $2.09 \\times 10^{-7} \\pm 2.37 \\times 10^{-7}$ | $3.44 \\times 10^{-7} \\pm 1.11 \\times 10^{-7}$ | $3.09 \\times 10^{-7} \\pm 3.05 \\times 10^{-8}$ |\n| - test loss    | $1.22 \\times 10^{-5} \\pm 1.64 \\times 10^{-5}$ | $6.43 \\times 10^{-8} \\pm 1.78 \\times 10^{-8}$ | $4.24 \\times 10^{-8} \\pm 3.52 \\times 10^{-8}$ |\n\nThe above table shows the average training and test losses with standard deviations for different training dataset sizes. Each result is computed based on 5 independent training runs with different random data. For both systems, the training loss consistently converges to a similar small value. However, the test loss improves as the number of training samples increases. \n\n- **Missing computational time for the examples, except the first example**. \n\nWe conducted a re-evaluation of the computation time across different systems. The computation time is generally similar to that of the Hénon-Heiles problem reported in the submitted paper. In particular, the proposed method achieves better accuracy with less computational time than the fourth-order methods. For the KdV equation experiment, the network was made as small as possible without sacrificing accuracy. The network used is quite small, yet the test loss was $1.706 \\times 10^{-6}$, which is not much different from the test loss of $1.695 \\times 10^{-6}$ when a larger network with 5 layers with width of 200 is used. For other experiments, the computation time can be further reduced by making the network smaller. We will include these results in the final version.\n\n| Problem        | Layer | Width | Steps | Proposed (sec) | 2 Order Methods (sec) | 4 Order Methods (sec) |\n|----------------|:-----:|:-----:|:-----:|----------------:|-------------------------:|-------------------------:|\n| Pendulum       |   5   |  50   |  200  |       0.579          |           0.503               |           1.137               |\n| Hénon-Heiles   |  10   | 100   |  300  |          2.090  |                   1.437  |                   3.533  |\n| 2-Body         |   5   |  50   |  300  |          2.862  |                   2.443 |                   7.693  |\n| 2-Pendulum     |   5   |  50   |  300  |          5.768  |                   5.197  |                  18.043  |\n| Duffing        |   5   |  50   |  300  |          1.020  |                    0.685  |                   1.084  |\n| KdV            |   2  | 30   | 3000  |        6.566  |                   1.634  |                  8.781  |\n\n- **Theorem 3.1 and Theorem 4.2 are more like definitions or simple facts. They are not sufficient to be called a theorem. In particular, the last sentence of Theorem 4.2 is not clear enough to be a statement in a theorem**.\n\nThank you very much. We will correct the statements of these theorems appropriately. \n\n- **The KdV example lacks the details, and the results are rough**. \n\nWe apologize for the lack of details. We consider this equation on the interval $[0,L], L=20$ under the periodic boundary condition. We used the training data with the time length of $T=30$ and the time step size of $\\Delta t=0.01$. The spatial step size was $\\Delta x=0.5$. The initial condition of the training data is constructed by superimposing two solitary wave solutions of the KdV equation: $u_0(x) = u_{sol}(x - x_1, c_1) + u_{sol}(x - x_2, c_2)$, where $u_{sol}(x, c) = \\frac{1}{2}c \\cdot \\mathrm{sech}^2\\left(\\frac{1}{2}\\sqrt{c} \\cdot x \\right)$ with $c_1=0.75, x_1=0.33L,c_2=0.4, x_2=0.65L$. The time evolution is computed by solving the KdV equation using the odeint solver with high accuracy, where spatial derivatives were approximated using finite difference methods. For testing, the computation was carried out from the time length $T=30$, with the start point being at the final state of the training dataset. We will include these details in the final version.\n\n- **What is the matrix $\\bar{S}$ in Eq.(5)? Is it the same as the matrix $S$ in Figure 1? How to derive or approximate $S$ or $\\bar{S}$ in applications?**\n\n$\\overline{S}$ can be $S$ itself or an approximation of $S$ as far as $\\overline{S}$ is skew-symmetric or negative-definite. A typical example where $S$ must be approximated is the KdV equation. In this equation, $S$ is the differential operator $\\partial/\\partial x$. Therefore, it is necessary to approximate this operator by, e.g., the finite difference method. Since the forward or backward difference does not result in a skew-symmetric matrix, the central difference method is usually used.\n\n- **Line 263: is $u^{(i)}$ a typo for $u^{(i-1)}$? Is there any reason for selecting such random training data?**\n\nThank you very much. $u^{(i-1)}$ is indeed a typo. Yes, there is a reason behind selecting the training data randomly within the range [−1,1]. Our primary goal is to capture the characteristic chaotic behavior of the Hénon–Heiles system within the bounded region. This choice allows the training set to cover the region where trajectories remain bounded and exhibit complex, chaotic motion. Some of the data can have higher energy, providing information about the system's behavior near the boundary between bounded and unbounded motion.\n\n- **Line 306+: the notation x is a typo of q. What if the system is stiff?**\n\nThank you very much. We will correct the typo. If the system is stiff, we should use the discrete gradient method for stiff problems as the base method. For example, the exponential AVF method (e.g., Wu et al., In Recent Developments in Structure-Preserving Algorithms for Oscillatory Differential Equations, Springer, 2018. Gu et al., arXiv:2110.04092, 2021)  is combined with the proposed approach."}, {"review_id": "7EOZrK1IVd", "reviewer": "Reviewer_U6fP", "summary": "The paper addresses the challenge of accurately simulating physical systems while preserving their intrinsic energy laws. Traditional numerical methods often fail to maintain these invariants, leading to instability. UEPI proposes a novel hybrid approach, integrating neural networks with discrete gradient methods to learn highly accurate, structure-preserving integrators from data. Physical phenomena are often described by energy-based theories, requiring numerical methods that preserve energy conservation or dissipation for accuracy and long-term stability. The paper demonstrates superior accuracy over  diverse systems like chaotic Hamiltonian, dissipative, and nonlinear PDEs.", "questions": "- What is the precision at which the experiments are run? At low-precision (FP32) common in machine learning, it is almost impossible for any solver (learned or otherwise) to capture chaotic behavior [1]\n- Since the paper proposes better numerical solvers (and not a better discrete gradient method), additional benchmarks using DAEs, symplectic integrators, etc. (at sufficient precision) should be demonstrated.\n\n[1] https://www.stochasticlifestyle.com/how-chaotic-is-chaos-how-some-ai-for-science-sciml-papers-are-overstating-accuracy-claims/", "limitations": "No, see the questions section.", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "### Strengths\n\n- The UEPI methodology is built upon discrete gradient methods, numerical schemes for differential equations with energy conservation or dissipation laws. The formulation ensures that energy preservation / dissipation is built-in into the formulation.\n\n### Weaknesses\n\n- Currently looking at the results, it seems like the proposed method is a better discrete gradient method. However, it is not clear where it stands compared to other numerical integrators. A simple example would be for the chaotic Henon-Heiles model [1] satisfies invariants even more consistently than the proposed method.\n- Henon-Heiles is very moderately chaotic in the tested regime (E < 1 / 6) [2]. Claiming that this method is general enough and works for chaotic systems is not validated with proper experiments.\n- Several examples (double pendulum / two body problem) can be solved with a DAE solver or with a ODE Mass Matrix solver. Additionally, comparison with Sympnets has been provided, but other data-driven methods that can conserve invariants (often outperforming Sympnets) like COMET [3], Stabilized Neural ODEs [4], etc. are omitted from discussion.\n\n\n[1] https://docs.sciml.ai/DiffEqDocs/stable/examples/classical_physics/#H%C3%A9non-Heiles-System\n\n[2] https://adsabs.harvard.edu/full/1964AJ.....69...73H\n\n[3] https://arxiv.org/abs/2208.10387\n\n[4] https://arxiv.org/abs/2306.09739", "was_revised": true, "final_justification": "I am updating my score since the core issue stemmed from a mistake in calculation on my end with no fault of the authors. I have updated my score (and reduced confidence for this regard) accordingly.", "rebuttal": "We sincerely thank the reviewer for the valuable comments and suggestions. We have carefully addressed the concerns raised by the reviewer. We believe the revisions and clarifications provided below improve the overall quality of the paper. We respond to each point in detail as follows.\n\n- **Currently looking at the results, it seems like the proposed method is a better discrete gradient method. However, it is not clear where it stands compared to other numerical integrators. A simple example would be for the chaotic Henon-Heiles model [1] satisfies invariants even more consistently than the proposed method.**\n- **Henon-Heiles is very moderately chaotic in the tested regime (E < 1 / 6) [2]. Claiming that this method is general enough and works for chaotic systems is not validated with proper experiments.**\n\nTo address the reviewer's concern, we conducted additional experiments on the Hénon-Heiles. As is well known (e.g., Costin et al. arXiv:2411.16071, 2024), when the energy exceeds 1/6, the state often diverges to infinity. For this reason, we performed experiments with energies marginally above 1/6. \n\nIn addition, as suggested by another reviewer, to address the concern regarding the statistical robustness and the dependence of the proposed method on training data, we evaluate the impact of training data distribution, data size, and initial conditions on the performance of the learned integrator.  **In particular, the initial conditions for testing include the cases with energies above 1/6**.  \n\nWe used irregular training data {$ ((u^{(i-1)},h^{(i)}),u^{(i)}) \\mid i=1,...,I $}, where the elements of $u^{(i-1)}$ are randomly sampled from the domain $[-1, 1]$ and each step size $h^{(i)}$ is drawn independently from $[0, 0.5]$. Each $u^{(i)}$ represents the one-step solution computed from $u^{(i-1)}$ using step size $h^{(i)}$ with a high-order solver. Each sample corresponds to a one-step transition computed from $u^{(i-1)}$ using a high-order solver. We performed experiments using $100$, $1000$, and $5000$ training samples, respectively. The neural network consists of $10$ layers, each with $100$ width. We set the learning rate to $0.00001$ and trained the model for $30000$ iterations.\n\nFor testing, we selected random initial conditions from $[-0.4, 0.4]^4$, including both general cases and high-energy cases with total energy exceeding $1/6$. For each test trajectory, a fixed step size $h$ was randomly chosen from $[0, 0.5]$ and applied over 100 time steps.\n\n| Problem       | 100 train data                              | 1000 train data                             | 5000 train data                             |\n|---------------|---------------------------------------------|---------------------------------------------|---------------------------------------------|\n| **Hénon-Heiles** |                                             |                                             |                                             |\n| - train loss   | $1.64 \\times 10^{-8} \\pm 6.09 \\times 10^{-9}$ | $1.83 \\times 10^{-8} \\pm 2.66 \\times 10^{-9}$ | $1.90 \\times 10^{-8} \\pm 1.30 \\times 10^{-9}$ |\n| - test loss    | $3.43 \\times 10^{-8} \\pm 1.98 \\times 10^{-8}$ | $5.29 \\times 10^{-9} \\pm 4.45 \\times 10^{-9}$ | $1.58 \\times 10^{-9} \\pm 1.60 \\times 10^{-9}$ |\n\nThe above table shows the average training and test losses (with standard deviations) for different training dataset sizes (100, 1000, and 5000) on the Hénon-Heiles. Each result is computed based on 5 independent training runs with different random data and test trajectories. The training loss consistently converges to a similar small value; however, the test loss improves as the number of training samples increases. This result shows that the proposed method also works well when including the cases with energies greater than 1/6.\n\n- **Comparison with Sympnets has been provided, but other data-driven methods that can conserve invariants (often outperforming Sympnets) like COMET [3], Stabilized Neural ODEs [4], etc. are omitted from discussion**.\n\nTo enrich our comparison with data-driven baselines that can conserve invariants, we additionally evaluated with a representative Physics-Informed Neural Network (PINN) baseline under aligned training settings. We chose PINNs for the additional comparison because this is suggested by another reviewer, and as far as we understood, COMET and Stabilized Neural ODEs are methods for finding differential equations. Additionally, as for COMET, a method combined with the discrete gradient method has been proposed (Matsubara et al., ICLR2023). Therefore, our method can be combined with COMET in the same way as Matsubara et al. We will include a discussion of the relationship between our method and these approaches in the final version.\n\nSpecifically, we implemented **the PINN baseline** as in [Raissi et al., 2019] **with energy conservation/dissipation regularization terms**, using the DeepXDE library, and applied it to three benchmark systems: Hénon–Heiles, double pendulum, and Duffing oscillator. Both our method and the PINN baseline used the same network architecture (10 layers, width 100), optimizer (Adam, learning rate 1e-4), and number of training iterations (30,000). We sampled 2,000 random collocation points over the spatiotemporal domain.\n\nTo evaluate the performance, we generated a ground-truth reference trajectory for each system using `scipy.integrate.odeint` with 100 time steps, starting from the same initial conditions. We then simulated trajectories using both the trained PINN and our trained model under the same initial condition. We report the maximum error over the time domain $[0, T]$, defined as the maximum deviation from the ground-truth trajectory. Furthermore, since PINN tends to degrade over long time horizons, we computed the maximum error accumulated up to each checkpoint (time steps $T=5, 10, 50$, and $100$), i.e., the largest deviation from the reference trajectory from step 0 up to the given time step $T$.\n\n**Experiment 1: Hénon--Heiles system** We set the initial condition as $(q_1, q_2, p_1, p_2) = (0.3, 0.3, -0.2, -0.2)$, and performed simulations over 100 steps with a fixed time step size $h = 0.3$. The maximum errors of both models, measured against the reference solution, are summarized in the table below.\n\n| Time step | Max error (ours) | Max error (PINN) |\n| --- | --- | --- |\n| 5 | 1.4300e-04 | 2.8557e-02 |\n| 10 | 1.5437e-04 | 9.7690e-02 |\n| 50 | 1.0973e-03 | 1.4770e+00 |\n| 100 | 9.5319e-03 | 1.5013e+00 |\n\n**Experiment 2: double pendulum system** We set the initial condition as $(q_1, q_2, p_1, p_2) = (\\pi/6, \\pi/8, 0, 0)$, and performed simulations over 100 steps with a fixed time step size $h = 0.3$.\n\n| Time step | Max error (ours) | Max error (PINN) |\n| --- | --- | --- |\n| 5 | 1.7222e-04 | 1.9158e-02 |\n| 10 | 8.2759e-04 | 5.3900e-02 |\n| 50 | 3.4987e-03 | 1.9816e-01 |\n| 100 | 6.1786e-03 | 2.3785e-01 |\n\n**Experiment 3: Duffing oscillator system** We set the initial condition as $(q, p) = (2.95, -3.08)$, and performed simulations over 100 steps with a fixed time step size $h = 0.1$.\n\n| Time step | Max error (ours) | Max error (PINN) |\n| --- | --- | --- |\n| 5 | 8.4021e-04 | 4.1957e-01 |\n| 10 | 9.3785e-04 | 6.9862e-01 |\n| 50 | 5.1848e-02 | 4.4981e+00 |\n| 100 | 1.1232e-01 | 4.4981e+00 |\n\nThe results in all three systems suggest that our method offers improved accuracy and stability over the tested PINN baseline, especially in long-term simulations.\n\n- **What is the precision at which the experiments are run? At low-precision (FP32) common in machine learning, it is almost impossible for any solver (learned or otherwise) to capture chaotic behavior**.\n\nAll experiments are performed with double precision (FP64). In fact, the neural networks used to learn the skew-symmetric matrices can be small and hence can be trained even on a CPU. \n\n- **Since the paper proposes better numerical solvers (and not a better discrete gradient method), additional benchmarks using DAEs, symplectic integrators, etc. (at sufficient precision) should be demonstrated**.\n\nWe indeed propose a better discrete gradient method. In the research field of structure-preserving numerical integrators, it is known as the Ge-Marsden theorem (Ge and Marsden, Phys. Lett. A, 133, 134, 1988) that numerical integrators that are both symplectic and energy-preserving cannot be designed. Whether symplectic or energy-preserving integrators are better depends on the problem. Therefore, the study of energy-preserving integrators, particularly discrete gradient methods, is important. \n\nA straightforward method that preserves conservation laws is generally the projection method, not DAE solvers. The projection method is a method of projecting numerical solutions computed by a numerical integrator onto the manifold that is defined by the target conservation law. This method is used to design energy-preserving integrators, but it is known that numerical solutions by the projection method are typically poorer than those by the discrete gradient method.\n\nComparisons between these methods and the discrete gradient method have long been made in the field of geometric numerical integration methods. We will include the typical comparisons between these methods in the final version of the paper as an appendix."}, {"review_id": "Cg4PkrE0AQ", "reviewer": "Reviewer_woku", "summary": "This paper proposes to learn an energy-preserving numerical integrator for energy-based dynamical systems (e.g., Hamiltonian mechanics or dissipative systems), building upon the the discrete gradient method. The key idea is to augment the discrete gradient with a learnable skew-symmetric term, modeled by a neural network. The authors claim that the proposed method generalizes the traditional discrete gradient method and is theoretically capable of universally approximating energy-preserving integrators. The proposed method is evaluated on various dynamical systems and outperforms conventional methods such as the AVF and Ito-Abe integrators.", "questions": "- Could the authors discuss the actual computational cost of the proposed method and compare it with that of the high-order method used to generate the training trajectories?\n\n- Could the authors provide a more detailed explanation of how the proposed method handles dissipative systems? For example, how is Theorem 4.6 derived? Additionally, why can the MLP-based **skew-symmetric** matrices approximate general **dissipative** integrators?\n\n- Is the learned neural skew-symmetric matrix generalizable across different integrator orders? For instance, if it is trained using a 4th-order discrete gradient scheme, can it still yield good results when tested with a 1st-order integrator, or vice versa?\n\n- There are some typos, e.g., in line 118, \"...if a matrix $S$ is skew-symmetric, then $u^T \\mathbf{G} u = 0$…\", ...", "limitations": "Yes, the authors acknowledge the computational cost as a limitation of the proposed method. However, I believe this issue is more critical than currently addressed. As noted in the Weaknesses section, a more thorough discussion is needed to support their efficiency claims.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\n- The proposed method appears to work well empirically, as the authors validate their approach across a range of physical systems, including the Henon–Heiles ODE, the two-body problem, the double pendulum, the dissipative pendulum, the Duffing oscillator, and the KdV PDEs.\n\n**Weaknesses**\n\n- The proposed approach requires learning from trajectory data to approximate a proper skew-symmetric matrix using a neural network. In other words, it necessitates obtaining a sufficient amount of high-quality trajectory data generated by a high-order integrator. Moreover, the performance of the learned numerical integrator heavily depends on the quality of the provided data. Therefore, the actual overall computational cost of the proposed method should include the time needed to generate the training data, the training time of the neural network, and the integration time during testing. However, current discussions of computational cost appear to focus solely on post-training integration time, overlooking these significant additional costs.\n\n- While the paper is generally well written, it is not entirely clear how the proposed method can universally approximate dissipative discrete gradient integrators. In line 204, the authors claim that the universal approximation properties (Theorem 4.5 for energy-preserving systems and Theorem 4.6 for dissipative systems) follow directly from Theorem 4.4. However, since Theorem 4.4 is valid only for conservative systems, and the authors state in lines 188–190 that it is impossible to derive an analogous result for dissipative systems, it remains unclear how Theorem 4.6 is justified.", "was_revised": true, "final_justification": "The authors provided a strong rebuttal, including an evaluation of computational complexities and benchmarking against other competitors. Based on this, I have increased my score to 4.", "rebuttal": "We are very grateful for the insightful comments provided by the reviewer. To address the reviewer's concerns, we respond to each point in detail as follows.\n\n- **The proposed approach requires learning from trajectory data to approximate a proper skew-symmetric matrix using a neural network. In other words, it necessitates obtaining a sufficient amount of high-quality trajectory data generated by a high-order integrator. Therefore, the actual overall computational cost of the proposed method should include the time needed to generate the training data, the training time of the neural network, and the integration time during testing**. \n\nAs many reviewers are concerned about computation time, we discuss the computation time in detail, especially for the case of the KdV equation, which is the largest system in the experiments. Firstly, GPUs can be used for training to accelerate calculations. In the following experiment, we used an Intel Xeon 6900P CPU and an AMD Instinct MI300A GPU for computation. \n\nWe first set the number of nodes for spatial discretization of the KdV equation to 100. The computation time for generating the training data was 28.709 sec using a CPU and the scipy odeint function. The training time was 6.407 sec using a GPU. The integration time during testing was 52.738 (proposed), 22.785 (2nd order AVF) and 127.280 (4th order AVF). **The total computational time of the proposed method was 28.709 + 6.407 + 52.738 = 87.854,  which is less than the 127.280 sec required by the 4th order AVF method.** \n\nSimilarly, if we set the number of the nodes to 120, the computational time for each method was 50.957 (data generation) + 8.311 (training) + 76.460 (integration) =135.728 (proposed), 33.034 (2nd order AVF), 200.003 (4th order AVF).  **The total computational time of the proposed method was again shorter than that of the 4th-order AVF method.** \n\nPlease note that, as shown in the submitted paper, our proposed approach typically outperforms the 4th-order AVF method.  Thus, the proposed approach can compute more accurate numerical solutions in a shorter amount of time, even when the training time and the time for generating data are included.\n\nA closer look at the execution time of the above experiments shows that the computation of the neural network itself is actually not so large. Rather, the matrix-vector multiplication of the skew-symmetric matrix takes $O(N^2)$ for the system with $N$ variables, which affects the computational complexity because the computation of a gradient typically requires $O(N)$ computations. However, the effect of this matrix multiplication does not seem to be so significant, since the multiplication with the skew-symmetric or negative-definite matrix placed in front of the gradient of the energy takes about the same amount of computation. It is also possible to reduce computational cost by low-rank approximation of the learned matrices.\n\nAnother possibility to reduce the computational time is to use the recently proposed **Kolmogorov-Arnold networks** (KANs, Liu et al., arXiv:2404.19756, 2024). This method is expected to approximate functions with fewer parameters than MLPs, thereby reducing computational time. Importantly, it is also known to be able to apply **symbolic regression to learned functions**. If a numerical method learned using the proposed method with KANs can be regressed symbolically, **the learned scheme becomes a standard numerical integration method expressed in mathematical form** rather than a black box method by machine learning. Also, the resultant mathematical expression may give clues to manually derive high-precision numerical integrators by hand.\n\n- **Could the authors discuss the actual computational cost of the proposed method and compare it with that of the high-order method used to generate the training trajectories?**\n\nIn the above experiment, the time periods computed to create the training data and the time periods for the test data are identical. Therefore, we can compare the computation times of the high-order method for creating the data and the proposed method for the integration for testing. It can be seen that the proposed method takes the same amount of time or less than twice as long as the high-order method; however, this time is implementation-dependent, and it is possible to reduce the computation time by improving the implementation. At least, adaptive time stepping must be introduced. Also, Reviewer tjF9 suggested using the neural network not every time stepping, but only when needed. This is expected to increase accuracy without a significant increase in computation time. Additionally, the high-order method does not satisfy the energy conservation/dissipation laws, and hence it is not suitable for long-term prediction. Therefore, it should be noted that it is not necessarily appropriate to compare them only in terms of computation speed.\n\n- **While the paper is generally well written, it is not entirely clear how the proposed method can universally approximate dissipative discrete gradient integrators**. \n- **Could the authors provide a more detailed explanation of how the proposed method handles dissipative systems? For example, how is Theorem 4.6 derived? Additionally, why can the MLP-based skew-symmetric matrices approximate general dissipative integrators?**\n\nWe apologize for the confusion. Please note that there is a small but important difference between Theorem 4.5 and Theorem 4.6; Theorem 4.5 is the universal approximation property for all energy-preserving numerical schemes and Theorem 4.6 is for all energy-dissipative discrete gradient schemes. \n\nFirstly, the proposed method can universally approximate any discrete gradient schemes, that is, the numerical scheme in the following form:\n\n$$\n\\frac{u^{(n+1)}-u^{(n)}}{\\Delta t } =S \\overline{\\nabla} H(u^{(n+1)}, u^{(n)})\n$$\n\nwhere $\\overline{\\nabla} H(u^{(n+1)}, u^{(n)})$ is a discrete gradient. These numerical schemes are dissipative if the matrix $S$ is negative-semidefinite. The MLP-based skew-symmetric matrices are used to universally approximate the discrete gradient $\\overline{\\nabla} H(u^{(n+1)}, u^{(n)})$. \n\nSecondly, any energy-conserving numerical scheme is shown to be a discrete gradient scheme; any energy-conservative numerical scheme must be of the above form. Therefore, our method has the universal approximation property for energy-preserving numerical schemes. \n\nHowever, we have not yet shown that all energy-dissipative numerical methods can be written in the above form. Therefore, for dissipative numerical methods, our method only has the universal approximation property for dissipative numerical schemes that can be expressed as a discrete gradient scheme, which is the statement of Theorem 4.6. We will clarify it in the final version of the manuscript by adding a remark.\n\n- **Is the learned neural skew-symmetric matrix generalizable across different integrator orders? For instance, if it is trained using a 4th-order discrete gradient scheme, can it still yield good results when tested with a 1st-order integrator, or vice versa?**\n\nThank you very much indeed for your insightful comments. In general, the learned skew-symmetric matrix cannot be used across different integrators. However, a matrix trained using a lower-order integrator can be used to derive a higher-order scheme by the composition method (Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A., Vol. 150, 262-268, 1990). The composition method designs a higher-order numerical integrator by applying a low-order method multiple (typically three) times with different step sizes. Thus, a low-order scheme improved by the proposed method can be further enhanced in this way. \n\n- **There are some typos, e.g., in line 118, \"...if a matrix $S$ is skew-symmetric, then $u^T \\mathbf{G} u = 0$…\"**, ...\n\nMany thanks. We will correct the typos."}, {"review_id": "N1mdbEvhA2", "reviewer": "Reviewer_tjF9", "summary": "This paper proposes a novel method that uses a neural network to learn highly accurate numerical solvers for physical systems. It builds upon discrete gradient methods, which inherently guarantee the conservation or dissipation of energy. The key innovation is using a neural network to optimize a corrective term within this physically-consistent framework. This approach ensures the learned solver respects fundamental physical laws by design, while training on data makes it significantly more accurate. Experiments show it excels in long-term simulations of complex and chaotic systems, outperforming standard geometric integrators.", "questions": "See above. I am happy to raise my score if my questions are well-addressed in the weakness section.", "limitations": "See above", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Pros:\n\n1. The paper is well-written in general.\n\n2. The proposed UEPI is technically sound, with detailed derivation towards its guaranteed physical properties.\n\n3. The experiments covers a diverse set of physical systems, showing the expressive power of the proposed approach.\n\n\nCons:\n\n1. There are some PINN work that also studies energy-conservative/disppastive systems, such as [1,2]. The authors only compared with  AVF, Itoh-Abe method and time-symmetric Itoh-Abe method. A lot of PINN methods should also be discussed and considered. \n\n2. Using a neural network at each step of the integration process is more computationally expensive than traditional methods. While the paper argues this cost is justified by the significant accuracy gains, it remains a notable overhead, especially for high-dimensional systems. Is there any comprehensive time comparison analysis regarding the computational overhead, in terms of prediction length, complexity of physical systems, etc? Have the authors study scheduling policy that only applys the neural network at certain integration steps, similar to [3]?\n\n[1] Physics-Informed Regularization for Domain-Agnostic Dynamical System Modeling.\n\n[2] Time-Reversal Symmetric ODE Network.\n\n[3] PhysDiff: Physics-Guided Human Motion Diffusion Model", "was_revised": false, "final_justification": "", "rebuttal": "We are really grateful for the insightful comments provided by the reviewer. To address the reviewer's concerns, we respond to each point in detail as follows.\n\n- **There are some PINN work that also studies energy-conservative/disppastive systems, such as [1,2]. The authors only compared with AVF, Itoh-Abe method and time-symmetric Itoh-Abe method. A lot of PINN methods should also be discussed and considered**.\n\nThank you for your thoughtful comment. PINN is a class of methods that solve differential equations by minimizing physics-informed loss functions. Our method is fundamentally different in that it constructs a numerical integrator rather than learning the solution directly; however, we understand the importance of empirical comparison and thus designed experiments as best as possible.\n\nSpecifically, we implemented **the PINN baseline** as in [Raissi et al., 2019] **with energy conservation/dissipation regularization terms**, using the DeepXDE library, and applied it to three benchmark systems: Hénon–Heiles, double pendulum, and Duffing oscillator. Both our method and the PINN baseline used the same network architecture (10 layers, width 100), optimizer (Adam, learning rate 1e-4), and number of training iterations (30,000). We sampled 2,000 random collocation points over the spatiotemporal domain.\n\nTo evaluate the performance, we generated a ground-truth reference trajectory for each system using `scipy.integrate.odeint` with 100 time steps, starting from the same initial conditions. We then simulated trajectories using both the trained PINN and our trained model under the same initial condition. We report the maximum error over the time domain $[0,T]$, defined as the maximum deviation from the ground-truth trajectory. Furthermore, since PINN tends to degrade over long time horizons, we computed the maximum error accumulated up to each checkpoint ($T = 5$, $10$, $50$, and $100$), i.e., the largest deviation from the reference trajectory from step 0 up to the given time step $T$.\n\n**Experiment 1: Hénon--Heiles system** We set the initial condition as $(q_1, q_2, p_1, p_2) = (0.3, 0.3, -0.2, -0.2)$, and performed simulations over 100 steps with a fixed time step size $h = 0.3$. The maximum errors of both models, measured against the reference solution, are summarized in the table below.\n\n| Time step | Max error (ours) | Max error (PINN) |\n| --- | --- | --- |\n| 5 | 1.4300e-04 | 2.8557e-02 |\n| 10 | 1.5437e-04 | 9.7690e-02 |\n| 50 | 1.0973e-03 | 1.4770e+00 |\n| 100 | 9.5319e-03 | 1.5013e+00 |\n\n**Experiment 2: double pendulum system** We set the initial condition as $(q_1, q_2, p_1, p_2) = (\\pi/6, \\pi/8, 0, 0)$, and performed simulations over 100 steps with a fixed time step size $h = 0.3$.\n\n| Time step | Max error (ours) | Max error (PINN) |\n| --- | --- | --- |\n| 5 | 1.7222e-04 | 1.9158e-02 |\n| 10 | 8.2759e-04 | 5.3900e-02 |\n| 50 | 3.4987e-03 | 1.9816e-01 |\n| 100 | 6.1786e-03 | 2.3785e-01 |\n\n**Experiment 3: Duffing oscillator system** We set the initial condition as $(q, p) = (2.95, -3.08)$, and performed simulations over 100 steps with a fixed time step size $h = 0.1$.\n\n| Time step | Max error (ours) | Max error (PINN) |\n| --- | --- | --- |\n| 5 | 8.4021e-04 | 4.1957e-01 |\n| 10 | 9.3785e-04 | 6.9862e-01 |\n| 50 | 5.1848e-02 | 4.4981e+00 |\n| 100 | 1.1232e-01 | 4.4981e+00 |\n\nThese results in all three systems suggest that our structure-preserving method offers improved accuracy and stability over the tested PINN baseline, especially in long-term simulations. This improvement is attributed to its formulation as a numerical integrator designed to respect the system's physical structure.\n\n- **Using a neural network at each step of the integration process is more computationally expensive than traditional methods. While the paper argues this cost is justified by the significant accuracy gains, it remains a notable overhead, especially for high-dimensional systems. Is there any comprehensive time comparison analysis regarding the computational overhead, in terms of prediction length, complexity of physical systems, etc? Have the authors study scheduling policy that only applys the neural network at certain integration steps, similar to [3]?**\n\nThank you very much indeed for the very insightful comments. We have not considered the scheduling policy similar to [3], and your suggested approach is in fact helpful for designing efficient numerical integrators. Indeed, by combining with local error estimation techniques, we can apply the neural network only when the estimated errors are very large. This would greatly reduce the computational complexity while improving accuracy, particularly for long-term predictions. We will discuss this approach in the final version of the paper.\n\nTo address the concerns about computation time, we performed additional experiments to check the computation time in detail, especially for the case of the KdV equation, which is the largest system in the experiments. Firstly, GPUs can be used to accelerate calculations for training. In the following experiment, we used an Intel Xeon 6900P and an AMD Instinct MI300A GPU for computation. \n\nWe first set the number of nodes for spatial discretization of the KdV equation to 100. The computation time for generating the training data was 28.709 sec using a CPU and the scipy odeint function. The training time was 6.407 sec using a GPU. The integration time during testing was 52.738 (proposed), 22.785 (2nd order AVF) and 127.280 (4th order AVF). The total computational time of the proposed method was 28.709 + 6.407 + 52.738 = 87.854,  which is less than the 127.280 sec required by the 4th order AVF method.\n\nSimilarly, if we set the number of the nodes to 120, the computational time for each method was 50.957 (data generation) + 8.311 (training) + 76.460 (integration) =135.728 (proposed), 33.034 (2nd order AVF), 200.003 (4th order AVF).  The total computational time of the proposed method was again shorter than that of the 4th-order AVF method.\n\nAs seen above, as the number of nodes increases and the complexity of the equation rises, all numerical integrators become time-consuming. In particular, computation time with the 4th-order AVF method increases significantly. Since the computation time of the proposed method is not so large compared to the 4th-order AVF method, the proposed method may be more effective when the complexity increases and the prediction period is long, particularly with the help of your suggested approach.\n\nA closer look at the execution time of the above experiments shows that the computation of the neural network itself is actually not so large. Rather, the matrix-vector multiplication of the skew-symmetric matrix takes $O(N^2)$ for the system with $N$ variables, which affects the computational complexity because the computation of a gradient typically requires $O(N)$ computations. However, the effect of this matrix multiplication does not seem to be so significant, since the multiplication with the skew-symmetric or negative-definite matrix placed in front of the gradient of the energy takes about the same amount of computation. It is also possible to reduce computational cost by low-rank approximation of the learned matrices.\n\nAnother possibility to reduce the computational time is to use the recently proposed **Kolmogorov-Arnold networks** (KANs, Liu et al., arXiv:2404.19756, 2024). This method is expected to approximate functions with fewer parameters than MLPs, thereby reducing computational time. Importantly, it is also known to be able to apply **symbolic regression to learned functions**. If a numerical method learned using the proposed method with KANs can be regressed symbolically, **the learned scheme becomes a standard numerical integration method expressed in mathematical form** rather than a black box method by machine learning. Also, the resultant mathematical expression may give clues to manually derive high-precision numerical integrators by hand."}]

Author rebuttal:

Meta-review: This paper introduces a method for learning numerical integrators that preserve or dissipate energy, extending discrete gradient methods with neural parametrizations. The approach is theoretically grounded, with universal approximation results, and is validated on a diverse set of Hamiltonian, dissipative, and PDE systems. Reviewers noted that the paper is well written, reproducible, and technically sound, with potential impact in machine learning for the sciences.

The main concerns centered on computational cost and the scope of baseline comparisons. While some reviewers initially questioned efficiency claims and the relevance of comparisons to methods such as PINNs, COMET, and Stabilized Neural ODEs, the authors provided detailed timing analyses and additional experiments with PINN and time-reversal ODE baselines, demonstrating favorable trade-offs and strong performance. Clarifications were also given for the theoretical results, dissipative cases, and robustness to training data.

After rebuttal, two reviewers increased their scores, one maintained a positive evaluation, and one remained borderline reject, largely due to preferences for broader baseline coverage. Overall, I find that the strengths (including sound methodology, thorough experiments, and substantial author clarifications) outweigh the remaining concerns. Therefore, I recommend acceptance.
Decision: Accept (poster)