pyHamSys

pyHamSys (short for Python Hamiltonian Systems) is an open-source 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. For a summary with the main steps, 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.

Pre-defined integrators are:

• 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.
• btype : str, optional Information on the Poisson bracket used in the equations of motion. For btype='pq', a canonical Poisson bracket in (p,q) is used, i.e., the dynamical variables (q, p) are real and canonically conjugate. If btype='psi', the dynamical variables are (ψ, ψ^*) where $\psi=(q + i p)/\sqrt{2}$. Default is 'pq'. For other btypes, the function coupling should be specified for the element of the class HamSys.

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 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.
• chi and chi_star : callable, optional
  Functions of (h, t, y) which returns, respectively, exp(h X_N)...exp(h X₁)y and exp(h X₁)...exp(h X_N)y at time t for its use in symplectic-split integration (without phase space extension). See [2] for more details.
• hamiltonian : callable, optional
  A function of (t, y) which returns the Hamiltonian H(t,y) where y is the state vector.
• coupling (for exended phases space with the restraint as in [3]) : callable, optional
  A function of (h, y, ω) which advances y from time t to t+h for the coupling Hamiltonian $\omega (y - \bar{y})^2/2$. This function is already computed for the types btype='pq' and 'psi'. For any other type, it should be provided.

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. For systems with n degrees of freedom, the Hamiltonian should be of the form $$H(q_1,\ldots,q_n,p_1,\ldots,p_n,t)$$ (see Example 1).

  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.

• integrate : callable
  Integrate the Hamiltonian system using either a pre-defined symplectic solver (see above for a complete list) or a standard IVP solver ('RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA'). Supports optional symplectic extension and energy conservation checks.

  • z0 (array_like)
    Initial condition(s) of the system.
  • t_eval (array_like)
    Times at which the solution is evaluated and stored.
  • params (Parameters)
    Parameters for the integration of the Hamiltonian system. Step must be provided.
  • 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).

  Pa