pyHamSys

pyHamSys (short for Python Hamiltonian Systems) is a Python library for scientific computing involving Hamiltonian systems. It provides tools to model, analyze, and simulate Hamiltonian systems. In particular, the library offers a streamlined and user-friendly environment for implementing and running symplectic-split integrators. These specialized numerical methods are designed to preserve the geometric structure of Hamiltonian systems, ensuring accurate and stable simulations of their dynamics over long time periods.

Installation

Installation within a Python virtual environment:

python3 -m pip install pyhamsys

For more information on creating a Python virtual environment, click here.

Symplectic Integrators

pyHamSys features a dedicated class, SymplecticIntegrator, which provides a comprehensive implementation of symplectic-split integrators. These integrators are designed specifically for the numerical integration of Hamiltonian systems, ensuring the preservation of the symplectic structure of phase space—a key property that underpins the stability and accuracy of long-term simulations of such systems. Symplectic-split integrators decompose the Hamiltonian into subcomponents that are individually solvable, allowing for efficient and accurate integration. This decomposition is particularly effective for complex or high-dimensional systems, as it minimizes numerical drift and conserves critical invariants like energy over extended time intervals. The SymplecticIntegrator class offers a variety of splitting methods, enabling users to select the most appropriate scheme for their specific Hamiltonian system and computational requirements. Each integrator is implemented to handle both autonomous and non-autonomous systems, supporting applications in classical mechanics, molecular dynamics, astrophysics, and quantum mechanics.

• Verlet (order 2, all purpose), also referred to as Strang or Störmer-Verlet splitting
• From Forest, Ruth, Physica D 43, 105 (1990):
  • FR (order 4, all purpose)
• From Yoshida, Phys. Lett. A 150, 262 (1990):
  • Yo#: # should be replaced by an even integer, e.g., Yo6 for 6th order symplectic integrator (all purpose)
  • Yos6: (order 6, all purpose) optimized symplectic integrator (solution A from Table 1)
• From McLachlan, SIAM J. Sci. Comp. 16, 151 (1995):
  • M2 (order 2, all purpose)
  • M4 (order 4, all purpose)
• From Omelyan, Mryglod, Folk, Comput. Phys. Commun. 146, 188 (2002):
  • EFRL (order 4) optimized for H = A + B
  • PEFRL and VEFRL (order 4) optimized for H = A(p) + B(q). For PEFRL, chi should be exp(h X_A)exp(h X_B). For VEFRL, chi should be exp(h X_B)exp(h X_A).
• From Blanes, Moan, J. Comput. Appl. Math. 142, 313 (2002):
  • BM4 (order 4, all purpose) refers to S₆
  • BM6 (order 6, all purpose) refers to S₁₀
  • RKN4b (order 4) refers to SRKN₆^b optimized for H = A(p) + B(q). Here chi should be exp(h X_B)exp(h X_A).
  • RKN6b (order 6) refers to SRKN₁₁^b optimized for H = A(p) + B(q). Here chi should be exp(h X_B)exp(h X_A).
  • RKN6a (order 6) refers to SRKN₁₄^a optimized for H = A(p) + B(q). Here chi should be exp(h X_A)exp(h X_B).
• From Blanes, Casas, Farrés, Laskar, Makazaga, Murua, Appl. Numer. Math. 68, 58 (2013):
  • ABA104 (order (10,4)) optimized for H = A + ε B. Here chi should be exp(h X_A)exp(h X_B).
  • ABA864 (order (8,6,4)) optimized for H = A + ε B. Here chi should be exp(h X_A)exp(h X_B).
  • ABA1064 (order (10,6,4)) optimized for H = A + ε B. Here chi should be exp(h X_A)exp(h X_B).

All purpose integrators are for any splitting of the Hamiltonian H=∑_k A_k in any order of the functions A_k. Otherwise, the order of the operators is specified for each integrator. These integrators are used in the functions solve_ivp_symp and solve_ivp_sympext by specifying the entry method (default is BM4).

---

HamSys class

The HamSys class provides a robust framework for defining and integrating Hamiltonian systems. It allows users to specify the number of degrees of freedom, coordinate representations, and key attributes like the Hamiltonian and associated equations of motion.

Parameters

• ndof : integer or half-integer, optional
  The number of degrees of freedom in the Hamiltonian system. Half integers denote an explicit time dependence. Default is 1.
• complex : bool, optional
  If False, the dynamical variables (q, p) are real and canonically conjugate. If True, the dynamical variables are (ψ, ψ^*) where $\psi=(q + i p)/\sqrt{2}$. Default is False.

Parameters and Attributes

• y_dot : callable, optional
  A function of (t, y) which returns {y,H(t,y)} where y is the state vector and H is the Hamiltonian. In (real) canonical coordinates (used, e.g., in solve_ivp_sympext) where y = (q, p), this function returns (∂H/∂p, -∂H/∂q). In complex coordinate ψ, this function returns -i ∂H/∂ψ^*. For practical implementation, the state vector y should be represented as a one-dimensional array with a shape of (n,), where n denotes the total number of dynamical variables in the system. This ensures compatibility with numerical solvers and facilitates efficient computation of the system's evolution.
• k_dot : callable, optional
  A function of (t, y) which returns {k,H(t,y)} = -∂H/∂t where k is canonically conjugate to t and H is the Hamiltonian.
• hamiltonian : callable, optional
  A function of (t, y) which returns the Hamiltonian H(t,y) where y is the state vector.

Functions

• compute_vector_field : from a callable function (Hamiltonian in canonical coordinates) written with symbolic variables (SymPy), computes the vector fields, y_dot and k_dot.

  Determine Hamilton's equations of motion from a given scalar function –the Hamiltonian– H(q, p, t) where q and p are respectively positions and momenta. However, it is preferrable to code explicitly and optimize y_dot and k_dot.

  Parameters

  • hamiltonian : callable
    Function H(q, p, t) –the Hamiltonian expressed in symbolic variables–, expressed using SymPy functions.
  • output : bool, optional
    If True, displays the equations of motion. Default is False.

  The function compute_vector_field determines the HamSys function attributes y_dot and k_dot to be used in solve_ivp_sympext. The derivatives are computed symbolically using SymPy.

• compute_energy : callable
  A function of sol –a solution provided by solve_ivp_sympext– and a boolean maxerror. If maxerroris True the function compute_energy returns the maximum error in total energy; otherwise, it returns all the values of the total energy.

  Parameters

  • sol : OdeSolution
    Solution provided by solve_ivp_sympext.
  • maxerror : bool, optional
    Default is True.

---

solve_ivp_symp and solve_ivp_sympext

The functions solve_ivp_symp and solve_ivp_sympext solve an initial value problem for a Hamiltonian system using an element of the class SymplecticIntegrator, an explicit symplectic splitting scheme (see [1]). These functions numerically integrate a system of ordinary differential equations given an initial value:
  dy / dt = {y, H(t, y)}
  y(t₀) = y₀
Here t is a 1-D independent variable (time), y(t) is an N-D vector-valued function (state). A Hamiltonian H(t, y) and a Poisson bracket {. , .} determine the differential equations. The goal is to find y(t) approximately satisfying the differential equations, given an initial value y(t₀) = y₀.

The function solve_ivp_symp solves an initial value problem using an explicit symplectic integration. The Hamiltonian flow is defined by two functions chi and chi_star of (h, t, y) (see [2]). This function works for any set of coordinates, canonical or non-canonical, provided that the splitting H=∑_k A_k leads to facilitated expressions for the operators exp(h X_k) where X_k = {A_k , ·}.

The function solve_ivp_sympext solves an initial value problem using an explicit symplectic approximation obtained by an extension in phase space (see [3]). This symplectic approximation works for canonical Poisson brackets, and the state vector should be of the form y = (q, p).

Parameters:

• chi (for solve_ivp_symp) : callable
  Function of (h, t, y) returning exp(h X_n)...exp(h X₁) y at time t. If the selected integrator is not all purpose, refer to the list above for the specific ordering of the operators. The operator X_k is the Liouville operator associated with the function A_k, i.e., for Hamiltonian flows X_k = {A_k , ·} where {· , ·} is the Poisson bracket. chi must return an array of the same shape as y.
• chi_star (for solve_ivp_symp) : callable
  Function of (h, t, y) returning exp(h X₁)...exp(h X_n) y at time t. chi_star must return an array of the same shape as y.
• hs (for solve_ivp_sympext) : element of class HamSys
  The attribute y_dot of hs should be defined. If check_energy is True, the attribute hamiltonian and if the Hamiltonian system has an explicit time dependence (i.e., the parameter ndof of hs is a half-integer), the attribute k_dot of hs should be specified.
• t_span : 2-member sequence
  Interval of integration (t₀, t_f). The solver starts with t=t₀ and integrates until it reaches t=t_f. Both t₀ and t_f must be floats or values interpretable by the float conversion function.
• y0 : array_like
  Initial state. For solve_ivp_sympext, the vector y0 should be with shape (n,).
• step : float
  Step size.
• t_eval : array_like or None, optional
  Times at which to store the computed solution, must be sorted, and lie within t_span. If None (default), use points selected by the solver.
• method : string, optional
  Integration methods are listed on pyhamsys. Default is 'BM4'.
• omega (for solve_ivp_sympext) : float, optional
  Coupling parameter in the extended phase space (see [3]). Default is 10.
• command : void function of (t, y), optional
  Void function to be run at each step size (e.g., plotting an observable associated with the state vector y, modify global or mutable variables, or register specific events).
• check_energy (for solve_ivp_sympext) : bool, optional
  If True, the attribute hamiltonian of hs should be defined. Default is False.

Returns:

  Bunch object with the following fields defined:

• t : ndarray, shape (n_points,)
  Time points.
• y : ndarray, shape y0.shape + (n_points,)
  Values of the solution y at t.
• k (for solve_ivp_sympext) : ndarray, shape (n_points,)
  Values of k at t. Only for solve_ivp_sympext and if check_energy is True for a Hamiltonian system with an explicit time dependence (i.e., the parameter ndof of hs is half an integer).
• err (for solve_ivp_sympext) : float
  Error in the computation of the total energy. Only for solve_ivp_sympext and if check_energy is True.
• step : step size used in the computation.

Remarks:

• Use solve_ivp_symp if the Hamiltonian can be split and if each partial operator exp(h X_k) can be easily and explicitly expressed/computed. Otherwise use solve_ivp_sympext if your coordinates are canonical, i.e., in $(q,p)$ or $(\psi,\psi^*)$ variables.
• The step size is slightly readjusted so that the final time t_f corresponds to an integer number of step sizes. The step size used in the computation is recorded in the solution as sol.step.
• For integrating multiple trajectories at the same time, extend phase space and define a state vector y = (y₁, y₂,...y_N) where N is the number of trajectories. The Hamiltonian is given by $H(t,\mathbf{y})=\sum_{i=1}^N h(t, y_i)$.

References:

• [1] Hairer, Lubich, Wanner, 2003, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (Springer)
• [2] McLachlan, Tuning symplectic integrators is easy and worthwhile, Commun. Comput. Phys. 31, 987 (2022); arxiv:2104.10269
• [3] Tao, M., Explicit symplectic approximation of nonseparable Hamiltonians: Algorithm and long time performance, Phys. Rev. E 94, 043303 (2016)

Example

import numpy as xp
import sympy as sp
import matplotlib.pyplot as plt
from pyhamsys import HamSys, solve_ivp_sympext
hs = HamSys()
hamiltonian = lambda q, p, t: p**2 / 2 - sp.cos(q)
hs.compute_vector_field(hamiltonian, output=True)
sol = solve_ivp_sympext(hs, (0, 20), xp.asarray([3, 0]), step=1e-1, check_energy=True)
print(f"Error in energy : {sol.err}")
plt.plot(sol.y[0], sol.y[1])
plt.show()

For more examples, see the folder Examples

---

For more information: [email protected]