On the algebraic properties of difference approximations of Hamiltonian systems

Full Text

1. Introduction In the realm of computational mathematics, the precise simulation of complex dynamical systems is paramount for enhancing scientists’ understanding of both natural and engineered phenomena. Hamiltonian described systems have many applications across diverse fields, including astrophysics, quantum mechanics, and both mechanical and electrical engineering. One of the most fundamental models in celestial mechanics and theoretical physics is Hamiltonian equations of motion system. These systems delineate the evolution of mechanical and physical phenomena over time and are essential for comprehending the behavior of a wide range of physical systems: from simple pendulums to intricate planetary orbits. In numerous instances, the Hamiltonian of these systems is represented as a polynomial or algebraic function of canonical variables ( 1,..., ) and momenta ( 1,..., ). The way these equations are structured helps to keep important quantities like energy and momentum constant, it ensures that the obtained solutions are consistent with physical principles and empirical data [1]. © 2025 Lapshenkova, L. O., Malykh, M. D., Matyukhina, E. N. This work is licensed under a Creative Commons “Attribution-NonCommercial 4.0 International” license. One of the main numerical study problem lies in the development of difference approximations that accurately capture the qualitative behavior of the system. Traditional methods, such as the widely utilized Runge-Kutta schemes, provide high accuracy but may fail to conserve certain geometric properties intrinsic to Hamiltonian systems, including phase-space volume and symplectic structure[2]. Over extended timescales, these numerical methods can lead to significant deviations in the system’s trajectory from its true path, particularly in chaotic or highly sensitive regimes [3]. The concept of “ordered chaos” finds profound resonance in both natural systems and technological applications, often analyzed through the lens of Hamiltonian mechanics. In nature, Hamiltonian systems describe the evolution of dynamic systems like planetary motion or fluid dynamics, where apparent chaos emerges from deterministic laws. Similarly, in technology, Hamiltonian principles underlie optimization in quantum computing and energy conservation in algorithms, demonstrating structured complexity. Whether modeling weather patterns or designing artificial intelligence systems, ordered chaos reveals how intricate behaviors can arise from foundational rules, bridging the gap between unpredictability and order in both realms. Consequently, alternative approaches that preserve the qualitative integrity of Hamiltonian dynamics are worth developing [4, 5]. 2. Numerical approaches In recent years, researchers have turned their attention to symplectic integration methods, specifically designed to preserve the symplectic structure of Hamiltonian flows [6]. These methods aim to conserve quantities such as energy and phase-space volume, providing a more faithful representation of the true motion of the system [7]. However, symplectic methods are often insufficient to fully capture the complexity of polynomial Hamiltonian systems, which may exhibit rich and diverse behaviors, including chaotic trajectories and resonance phenomena [8, 9]. Chaotic dynamics are utilized in optimization algorithms, particularly in exploring large, complex solution spaces. The unpredictable yet deterministic nature of chaotic trajectories helps avoid becoming stuck in local minima, improving the efficiency of finding global solutions in systems like neural networks, genetic algorithms, and quantum computation. The study presented here focuses on difference approximations for Hamiltonian systems with polynomial Hamiltonians, a subset that includes many physically significant cases. We explore the idea of birational correspondences, wherein each step of the approximation represents a rational transformation between the initial and final states of the system [10]. This approach allows for exact preservation of certain structural properties and often leads to more accurate long-term simulations than traditional methods. The finite difference method proposes replacing the system of differential equations = ( 1,…, ), = 1,…, , or, for short, = ( ), (1) with a system of algebraic equations ( , , ) = 0,̂ = 1,…, , (2) relating the value of the solution at some moment in time with the value ̂ of the solution at the moment in time + . The system of the algebraic equation (2) itself will be called a difference scheme for a system of the differential equation (1). Discrete models are essential in connecting analytical solutions with numerical simulations in dynamical systems, particularly for Hamiltonian systems with polynomial Hamiltonians. Kahan’s method is notable for preserving energy integrals and approximating elliptic oscillators, while symplectic schemes, like the midpoint method, ensure energy conservation in quadratic Hamiltonians, facilitating closed trajectories. Appelroth’s quadratization technique simplifies complex systems by transforming polynomial right-hand sides into quadratic forms [11, 12]. Geometric integrators are designed to maintain the analytical properties of the original systems, ensuring qualitative behaviors are preserved. Additionally, combining Appelroth’s quadratization with Kahan’s method enhances the analytic properties of difference approximations for higher-genus curves. Together, these discrete models provide effective tools for simulating dynamical systems and deepen the understanding of the interplay between algebraic structures and physical phenomena. In mechanics the quantity has often been treated as a finite increment, and it was implied that Newton’s equations were actually difference equations [13]. For example, the explicit Euler scheme ⃗ - ⃗ =̂ ( )⃗ ⃗ for linear oscillator preserves the energy = 2 + 2 only at → 0. The problem is that classical difference schemes (explicit Runge-Kutta schemes) are not rich in algebraic properties. We do describe properties of discrete models by looking back at continuous models. Kahan’s method is intrinsically linked to the development of discrete models for dynamical systems, particularly in the context of Hamiltonian mechanics. Discrete models serve as numerical approximations of continuous systems, allowing for the simulation of complex dynamics over time. Kahan’s discretization specifically addresses the challenges associated with preserving the geometric and analytical properties of Hamiltonian systems when transitioning from continuous to discrete representations [14, 15]. The method is designed to preserve the energy integral, which is crucial for maintaining the qualitative behavior of the system over time. In the context of a Hamiltonian system, the energy integral can be expressed as: ( , ) = , where is a constant representing the total energy of the system. Kahan’s method allows for the discretization of the Hamiltonian equations, leading to a difference scheme that can be written in the form of a quadrature: +1 = ∫ , where represents an elliptic integral of the first kind. This formulation ensures that the points of the approximate solution lie on an elliptic curve, thus inheriting the geometric properties of the original Hamiltonian system. Moreover, Kahan’s discretization approximates the solution by defining a birational transformation on the integral curve ( , ) = , which does not extend to a birational transformation of the entire phase space . This method allows for the correction of the integral curve while preserving its genus, enabling the simulation of an elliptic oscillator. Although the symplectic structure is not preserved exactly in case of the methods or conditions of the system fail to respect the intrinsic geometric properties of Hamiltonian dynamic, it is inherited to a significant extent, making Kahan’s method a powerful tool for imitating the dynamics of systems governed by cubic Hamiltonians [16]. In 1990s, the concept of geometric integrators was introduced, marking a significant advancement in the numerical analysis of differential equations. These integrators are designed to construct numerical schemes that preserve specific algebraic and geometric properties intrinsic to the original continuous model. It leads to maintaining the qualitative features of the dynamical system, such as symplectic structure, conservation laws, and invariants, which are often lost in traditional numerical methods. This approach enhances the fidelity of simulations and ensures that the long-term behavior of the system is accurately represented, thereby providing a robust framework for the analysis of complex dynamical phenomena. Due to the following expression + ̂ - = ̂ ( ) (3) 2 we can outline that it is -symmetric and symplectic. Also, (3) preserves quadratic integrals (Cooper’s theorem). The midpoint scheme is a numerical method used for solving ordinary differential equations (ODEs), particularly in the context of Hamiltonian systems. It is a type of symplectic integrator, which means it is designed to preserve the geometric properties of the original continuous system, such as energy conservation and symplectic structure. 3. Quadratic Hamiltonian A quadratic Hamiltonian refers to a Hamiltonian function that is a quadratic polynomial in the phase space variables, typically the coordinates and momenta . In mathematical terms, a quadratic Hamiltonian can be expressed in the form: 1 2 + 1 2, ( , ) = 2 2 where is the mass of the particle, is a constant related to the spring constant in the case of harmonic oscillators, is the momentum, and is the position. Quadratic Hamiltonians are significant in classical mechanics because they describe simple harmonic motion, such as that of a mass on a spring or a pendulum for small angles. The dynamics of systems with quadratic Hamiltonians can be analyzed using symplectic integrators, which preserve the structure of Hamiltonian mechanics, ensuring that energy is conserved over time. This makes them particularly useful in numerical simulations, as they provide accurate representations of the system’s behavior [17]. The midpoint scheme perfectly imitates a Hamiltonian system = , = - (4) with a quadratic Hamiltonian , for example, a harmonic oscillator with Hamiltonian = 2 + 2. - According to Cooper’s theorem, the energy integral is preserved exactly on the scheme, and the approximate solution itself is a sequence of points = ( , ) of the circle 2 + 2 = . - Each step of the approximate solution is a rotation by an angle +1 = ∫ , √ - 2 which does not depend on . 4. Cubic Hamiltonian Cubic Hamiltonians extend Hamiltonian mechanics beyond quadratic forms, expressed as: ( , ) = 3 + 2 + , with constants , , and . The cubic term introduces nonlinear dynamics, potentially leading to chaotic behavior. It provides insight into complex transitions, enhances the realism of models, and enables the study of phenomena that cannot be explained by linear systems. Traditional symplectic integrators may struggle with these complexities, but methods like Kahan’s discretization effectively approximate these systems, preserving integral properties and simulating relationships between initial and final states. This enhances numerical accuracy and provides greater insight into cubic Hamiltonian systems in both theoretical and applied contexts. In the 1990s, it was anticipated that preserving a symplectic structure would enable the imitation of a continuous model in the nonlinear case. The midpoint scheme is a symplectic Runge-Kutta scheme, i.e. ̂ ∧ = ∧ .̂ However, now the energy integral is not conserved, but is inherited in a very tricky formulation. Theorem 1 (J. M. Sanz-Serna and M.P. Calvo, 1994). For any ∈ ℕ, there exists a polynomial ( , ) such that goes to at to 0, ( ,̂ , ) = ̂ ( , , ) + ( ). Thus, in computer experiments, it seems that approximate solution lies on closed curve ( , ) = at sufficiently large . Let’s get back to Kahan’s method. Kahan’s method for Hamiltonian systems: 1. Geometric properties of Kahan’s method restricted to quadratic vector fields. 2. For the systems with cubic hamiltonian, Kahan’s method conserved the modified Hamiltonian -1 + ( - ⃗) .⃗ ∇ 3 2 ⃗ 3. For the systems with cubic Hamiltonian, Kahan’s method preserves the measure 1 ∧ 2 ⋯ ∧ . ⃗ det ( - ) 2 ⃗ Kahan’s scheme perfectly imitates a Hamiltonian system with a cubic Hamiltonian , for example, an elliptic ℘-oscillator. 1. According to 1st Celledoni’s theorem, the symplectic structure is inherited, i.e. ̂ ∧ = (1 + ( )) ∧ .̂ 2. According to 2nd Celledoni’s theorem, the energy integral is inherited, thus the approximate solution itself is a sequence of points = ( , ) of an elliptic curve ( , , ) = . Consider more closely the narrowing of Cremona map to the invariant curve ( , , ) = . Using constructions from Picard’s theorem, it follows that the difference scheme can be again represented using quadrature ⃗ ̂ ∫ ( , , ) = ( ), ⃗ where 1 is an elliptic integral of the 1st kind on invariant curve and, of course, → (at → 0). 5. Polinomial Hamiltonian If the Hamiltonian is a polynomial of degree > 3, then the exact solution to the continuous model lies on an algebraic curve ( , ) = , whose genus is greater than 1. Thus the quadrature = + on the curve is Abelian integral of the 1st kind. The integral cannot be inverted, and the functions ( ), ( ) (Jacobi problem). Formally, our method is suitable only for dynamical systems with quadratic right-hand sides, as the nonlinearity introduced by higher-order terms complicates the inversion process and the analytic behavior of the solutions. For more general systems, alternative methods, such as perturbation techniques or numerical approaches, are required to obtain meaningful results. Theorem 2 (Appelroth, 1902). Any dynamical system with polynomial right-side can be rewritten as dynamical systems with quadratic right-side in new variables. In the 21st century, the process of reducing a dynamical system characterized by a greater number of new variables than initial variables to one with a quadratic right-hand side is termed quadratization. The method developed by Appelroth in the early 20th century facilitates this transformation specifically for systems with polynomial right-hand sides, thereby enabling the application of symplectic Runge-Kutta schemes, which preserve integrals of the system, though not all properties of the original formulation may be retained. Consequently, any polynomial Hamiltonian system can be effectively integrated using a reversible scheme that consists of two key steps: first, the quadratization as proposed by Appelroth, followed by the discretization method established by Kahan [18]. Although Newton’s equations must define a one-to-one correspondence between the initial and final positions of a dynamical system, these equations do not actually define a one-to-one correspondence between the initial and final states of the system [10]. We can impose violently this property on the difference scheme. Currently, there are several implementations of quadratization algorythm [19]: Qbee by A. Bychkov and BioCham by M. Hemery et al. The example was produced via Sage using Qbee library made by Bychkov [20]. For more information, visit the https://github.com/AndreyBychkov/ QBee repository on GitHub. Example 1. Let’s figure out how the following system can be solved: ⎧ ̇1 = 2, ⎨⎩ ̇2 = 4. Consider the system with Hamiltonian of this system is 13 25 = - . 3 5 The solution is described by the quadrature 2 ′. ∫= + + 3 The particular solution of the system (4) with the initial conditions 1 = 2 = 1 at = 0 has a branching point ≈ 0.52. In the process of solving the given system of equations, we utilized Sage along with the Qbee library made by Bychkov. We started by defining the system, where the first equation is ̇1 = 2 and the second equation is ̇2 = 4. x1, x2 = functions("x1, x2") system = [ (x1, eval("x2**4")), (x2, eval("x1**2")) ] res = quadratize(system) Using the capabilities of the Qbee library, we proceeded to perform the necessary calculations. The main problem we faced with was the Qbee output difference in comparison with Cremona in Sage which, in turn, is used to find the approximate solution of the system. After quadratization implemented by Bychkov and restructured by us we’ve got the following system: ⎧ ̇1 = 0 ⋅ 2, ⎪ ⎪⎪ 2, ̇2 = 1 ⎪ ⎨ ̇0 = 2 ⋅ 1 ⋅ 1, ⎪ ⎪⎪ ̇1 = 02 + 2 ⋅ 2 ⋅ 2, ⎪ ⎩ ̇2 = 3 ⋅ 12, which can be rewritten into the following form to be solved using Cremona in Sage: ['x1', 'x2', 'w0', 'w1', 'w2'] ['w0*x2', 'x1**2', '3*w1*x1', 'w0**2 + 2*w2*x2', '3*w1**2'] The following problem is solved in the field of real numbers: a1,a2,a3,a4,a5=1 pr2 = Initial_problem([eval(dn) for dn in derivative_names], [eval(ed) for ed in equations_derivatives], [a1,a2,a3,a4,a5], 0.6) sol2=cremona_scheme(pr2, N=100, field=RR) The approximate solution passes through the branch point as a pole. Before the branch point the exact and approximate solutions coincide. After this point the exact solution in imaginary, but the approximate is real. The integral that connect old and additional variables is not preserved by reversible schemes. Before the branch point the expression 0 = 23 is equal to 0 on exact solution and is small on approximate solution, but its value at the branch point is very large (about 1029). Figure 1. Approximate solution - 23 + 0 on dependency The appropriateness of applying the combined Appelroth-Kahan approach to natural phenomena depends on which properties of these phenomena are important for research and which can be sacrificed [21]. This approach is good if the one-to-one correspondence between the initial and final positions of the system is most important [22]. 6. Discussion The design of difference schemes that accurately replicate systems with polynomial Hamiltonians involves a complex interplay between algebraic properties and physical principles. A critical question arises regarding whether the correspondence between initial and final positions should be one-toone, particularly in Hamiltonian dynamics, where preserving certain properties is vital for model integrity [23]. Historically, the quadrature (∫ = ) has been fundamental in understanding Hamiltonian dynamics. However, it is known that this representation fails to express ( ) as a single-valued analytic function of ( ) when the genus of the curve defined by ( ( , ) = ) exceeds 1, underscoring the complexities of nonlinear systems and the challenges in simulating their behavior. Our exploration reveals that while traditional symplectic methods, such as the midpoint scheme, effectively preserve energy in quadratic Hamiltonians, they may not extend this preservation to systems with cubic or higher-order Hamiltonians. This raises questions about whether preserving symplecticity guarantees that numerical solutions accurately reflect the system’s dynamics [24]. Introducing additional variables, as seen in the many-body problem, presents a potential solution. Reformulating the system to include distances and reciprocal distances allows for a new set of differential equations that, despite losing Hamiltonian structure, preserve classical integrals of motion [25]. This approach highlights the importance of integrals in governing dynamical behavior and suggests that focusing on these may yield more accurate numerical approximations. Furthermore, combining methods from Appelroth and Kahan offers a promising avenue for approximating solutions to polynomial Hamiltonian systems. Utilizing Cremona transformations can simplify the analytic properties of difference approximations compared to the original Hamiltonian model, raising the question of whether these simplified models are superior to their continuous counterparts [26]. Ultimately, the design of difference schemes transcends technical challenges, raising important questions about the nature of dynamical systems. The delicate balance between preserving essential properties and achieving computational feasibility significantly impacts the fidelity of numerical solutions. 7. Conclusion In conclusion, investigating difference approximations for Hamiltonian systems necessitates further exploration of the governing principles. As we refine our understanding, we must consider both algebraic and physical implications of our numerical methods. The pursuit of accurate simulations of dynamical systems continues, with insights from this research contributing to advancements in theoretical and applied mathematics.

About the authors

Lyubov O. Lapshenkova

RUDN University

Email: lapshenkova_lo@pfur.ru
ORCID iD: 0000-0002-1053-4925

PhD student of the chair of Mathematical Modeling and Artificial Intelligence

6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation

Mikhail D. Malykh

RUDN University; Joint Institute for Nuclear Research

Author for correspondence.
Email: malykh_md@pfur.ru
ORCID iD: 0000-0001-6541-6603
Scopus Author ID: 6602318510
ResearcherId: P-8123-2016

Doctor of Physical and Mathematical Sciences, Head of the department of Mathematical Modeling and Artificial Intelligence of RUDN University and research fellow of LIT JINR

6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation; 6 Joliot-Curie St, Dubna, 141980, Russian Federation

Elena N. Matyukhina

RUDN University

Email: matykhina_en@pfur.ru
Senior lecturer of the chair of Mathematical Modeling and Artificial Intelligence 6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation

References

Supplementary files