Skip to main content

Some tools for Hamiltonian systems

Project description

pyHamSys

pyHamSys is a Python package for scientific computing involving Hamiltonian systems

PyPI License

Installation:

pip install pyhamsys

Symplectic Integrators

pyHamSys includes a class SymplecticIntegrator containing the following symplectic splitting integrators:

All purpose integrators are for any splitting of the Hamiltonian H=∑k Ak in any order of the functions Ak. 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

Parameters

  • ndof : number of degrees of freedom of the Hamiltonian system 'ndof' should be an integer or half an integer. Half integers denote an explicit time dependence.

Attributes

  • hamiltonian : callable
    A function of (t, y) which returns the Hamiltonian H(t,y) where y is the state vector.
  • y_dot : callable
    A function of (t, y) which returns {y,H(t,y)} where y is the state vector and H is the Hamiltonian. In canonical coordinates (used, e.g., in solve_ivp_sympext) where y = (q, p), this function returns (∂H/∂p, -∂H/∂q).
  • k_dot : callable
    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.

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.

    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 maxerror, a boolean indicating whether the maximum error in total energy is given (if True) or all the values of the total energy (if False).

    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(t0) = y0
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(t0) = y0.

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 Ak leads to facilitated expressions for the operators exp(h Xk) where Xk = {Ak , ·}.

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 Xn)...exp(h X1) 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 Xk is the Liouville operator associated with the function Ak, i.e., for Hamiltonian flows Xk = {Ak , ·} 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 X1)...exp(h Xn) 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 attributes y_dot of hs should be defined. If check_energy is True. It the Hamiltonian system has an explicit time dependence (i.e., the parameter ndof of hs is half an integer), the attribute k_dot of hs should be specified.
  • t_span : 2-member sequence
    Interval of integration (t0, tf). The solver starts with t=t0 and integrates until it reaches t=tf. Both t0 and tf must be floats or values interpretable by the float conversion function.
  • y0 : array_like, shape (n,)
    Initial state.
  • step : float
    Step size.
  • t_eval : array_like or None, optional
    Times at which to store the computed solution, must be sorted and equally spaced, and lie within t_span. If None (default), use points selected by the solver.
  • method : string, optional
    Integration methods are listed on pyhamsys.
    'BM4' is the default.
  • omega (for solve_ivp_sympext) : float, optional
    Coupling parameter in the extended phase space (see [3]). Default = 10.
  • command : function of (t, y)
    Function to be run at each step size (e.g., plotting an observable associated with the state vector y, 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 (n, n_points)
    Values of the solution y at t.
  • k (for solve_ivp_sympext) : ndarray, shape (n//2, 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 is the Hamiltonian can be split and if each partial operator exp(h Xk) can be easily expressed/computed. Otherwise use solve_ivp_sympext if your coordinates are canonical.
  • If t_eval is a linearly spaced list or array, or if t_eval is None (default), the step size is slightly readjusted so that the output times contain the values in t_eval, or the final time tf corresponds to an integer number of step sizes. The step size used in the computation is recorded in the solution as sol.step.

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 information: cristel.chandre@cnrs.fr

Project details


Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.

Source Distribution

pyhamsys-0.31.tar.gz (13.3 kB view hashes)

Uploaded Source

Built Distribution

pyhamsys-0.31-py3-none-any.whl (13.8 kB view hashes)

Uploaded Python 3

Supported by

AWS AWS Cloud computing and Security Sponsor Datadog Datadog Monitoring Fastly Fastly CDN Google Google Download Analytics Microsoft Microsoft PSF Sponsor Pingdom Pingdom Monitoring Sentry Sentry Error logging StatusPage StatusPage Status page