Hamiltonian System API

Bases: object

Class for defining, integrating, and analyzing Hamiltonian systems.

This class represents Hamiltonian systems H(q, p), where q denotes the generalized coordinates and p the conjugate momenta. Two ways of specifying the system are supported:

Separable systems, H(q, p) = T(p) + V(q), specified via gradient functions for the kinetic and potential energies (grad_T, grad_V), with optional Hessians (hess_T, hess_V) for tangent-space computations. These are integrated with explicit symplectic methods (velocity Verlet, fourth-order Yoshida).
General (possibly non-separable) systems H(q, p), specified via the full equations of motion (eom) zdot = f(z), with z = (q, p) and the Hessian of H with respect to the combined state z (hess_H). These are integrated with the implicit midpoint method.

A system can be created either from one of the built-in models or from user-supplied functions of either kind above.

The class provides symplectic fixed-step integration routines together with tools for trajectory generation and nonlinear-dynamics analysis, including Poincaré sections, Lyapunov exponents, covariant Lyapunov vectors (CLVs), SALI, LDI, GALI, recurrence time entropy, and Hurst exponent estimation.

Parameters

modelstr or None, optional: Name of a built-in Hamiltonian model.
grad_Tcallable or None, optional: Gradient of the kinetic energy with respect to the momenta. Required (with grad_V) for separable systems.
grad_Vcallable or None, optional: Gradient of the potential energy with respect to the coordinates. Required (with grad_T) for separable systems.
hess_Tcallable or None, optional: Hessian of the kinetic energy with respect to the momenta.
hess_Vcallable or None, optional: Hessian of the potential energy with respect to the coordinates.
eomcallable or None, optional: Equations of motion of the system, with signature eom(q, p, parameters) -> (qdot, pdot), where qdot = dH/dp and pdot = -dH/dq. The return value must be a tuple in this exact order, (qdot, pdot), not (pdot, qdot). Required (with hess_H) for general, possibly non-separable systems integrated with the implicit midpoint method.
hess_Hcallable or None, optional: Full Hessian of the Hamiltonian H with respect to the combined state z = (q, p), with signature hess_H(q, p, parameters) returning an array of shape (2 * degrees_of_freedom, 2 * degrees_of_freedom). Required (with eom) for general, possibly non-separable systems integrated with the implicit midpoint method.
degrees_of_freedomint or None, optional: Number of degrees of freedom of the custom system.
parametersarray_like or None, optional: Parameter vector for the system.
number_of_parametersint or None, optional: Number of parameters expected by the custom system.

Notes

Custom systems must be specified either as a separable system via grad_T/grad_V (with optional hess_T/hess_V), or as a general system via eom/hess_H, together with degrees_of_freedom.
Lyapunov exponents, CLVs, SALI, LDI, and GALI require Hessian functions (hess_T/hess_V for separable systems).
Ensemble trajectory and reduced-map methods accept multiple initial conditions when supported by the corresponding wrapper.

See Also

ContinuousDynamicalSystem : Class for general continuous-time dynamical systems. DiscreteDynamicalSystem : Class for discrete-time maps.

classmethod available_models() → List[str][source]

List the available predefined Hamiltonian models.

Returns

list of str: Names of the supported models.

classmethod available_integrators() → List[str][source]

List the available predefined Hamiltonian models.

Returns

list of str: Names of the supported models.

property info: Dict[str, Any]

Information dictionary for the selected model.

Returns

dict: Dictionary containing metadata such as description, gradients, Hessians, degrees of freedom, and parameters.

Raises

ValueError: If no predefined model was used to initialize the system.

property integrator_info: Dict[str, Any]

Information dictionary for the current integrator.

Returns

dict: Dictionary containing the integrator description and associated step functions.

Set the symplectic integrator and integration time step.

Parameters

integratorstr: Name of the integrator. Available options are: - ‘svy4’: 4th-order Yoshida method - ‘vv2’: 2nd-order velocity-Verlet method
time_stepnumeric_t, optional: Integration time step. Must be a positive real number.
tolnumeric_t, optional: Newton convergence tolerance on the residual norm.
max_iterint, optional: Maximum Newton iterations per step.

Raises

TypeError: If integrator is not a string. If time_step, tol, or max_iter are not real numbers.
ValueError: If time_step, tol, or max_iter are not positive. If integrator is not implemented.

Examples

>>> from pynamicalsys import HamiltonianSystem
>>> HamiltonianSystem.available_integrators()
['svy4', 'vv2', 'imp']
>>> ds = HamiltonianSystem(model="henon heiles")
>>> ds.integrator("svy4", time_step=0.001)
>>> ds.integrator("vv2", time_step=0.001)
>>> ds.integrator("imp", time_step=0.01, tol=1e-14, max_iter=100)

set_parameters(parameters: ndarray[tuple[int, ...], dtype[float64]] | Sequence[float] | float) → None[source]

Set the parameter vector of the dynamical system.

This method validates and stores the model parameters. The input can be a scalar, a sequence of floats, or a NumPy array. It is internally converted into a float64 NumPy array of the appropriate size.

Parameters

parametersfloat or sequence of float or ndarray of shape (P,): The parameter set to be used by the system.

Returns

None

get_parameters() → ndarray[tuple[int, ...], dtype[float64]] | None[source]

Return the current parameter vector of the dynamical system.

Returns

ndarray of float64, shape (P,): The parameter vector currently stored in the system.

Advance the Hamiltonian system by one integration step.

Parameters

qnumeric_like_t: Generalized coordinates. Must define either a 1D array of shape (dof,) or a 2D array of shape (num_ic, dof).
pnumeric_like_t: Generalized momenta. Must have the same shape as q.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.

Returns

tuple[NDArray[np.float64], NDArray[np.float64]]

Updated coordinates and momenta after one integration step.

If q and p are 1D, returns (q_new, p_new) with shape (dof,).
If q and p are 2D, returns (q_new, p_new) with shape (num_ic, dof).

Raises

ValueError: If q and p do not have the same shape. If the number of parameters does not match the expected number.
TypeError: If q or p cannot be interpreted as valid initial conditions.

Generate a trajectory for the Hamiltonian system.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define either a 1D array of shape (dof,) or a 2D array of shape (num_ic, dof).
pnumeric_like_t: Initial generalized momenta. Must have the same shape as q.
total_timenumeric_t: Total integration time.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.

Returns

NDArray[np.float64]

Trajectory data.

If q and p are 1D, returns an array of shape (num_steps + 1, 2 * dof + 1), where the first column is time, the next dof columns are the coordinates, and the final dof columns are the momenta.
If q and p are 2D, returns an array of shape (num_ic, num_steps + 1, 2 * dof + 1).

Raises

ValueError: If q and p do not have the same shape. If total_time is not positive. If the number of parameters does not match the expected number.
TypeError: If total_time is not a real number. If q or p cannot be interpreted as valid initial conditions.

Compute a Poincaré section of the Hamiltonian trajectory.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define either a 1D array of shape (dof,) or a 2D array of shape (num_ic, dof).
pnumeric_like_t: Initial generalized momenta. Must have the same shape as q.
num_intersectionsint_t: Number of section crossings to record.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
section_indexint_t, optional: Index of the coordinate used to define the section.
section_valuenumeric_t, optional: Value of the selected coordinate at which the section is taken.
crossingint_t, optional: Crossing rule: - -1 for downward crossings - 0 for all crossings - 1 for upward crossings
periodic_section_coordinatebool, optional: If True, treats q[section_index] as a periodic coordinate on S¹ and performs crossing detection using modulo arithmetic. If False, uses standard Euclidean crossing detection.
periodnumeric_t, optional: Period of the angular coordinate when periodic_section_coordinate=True. Typically 2π for action-angle systems.
max_workersint_t, optional: The maximum number of processes that can be used to execute the given calls. If -1 or not given then as many worker processes will be created as the machine has processors. Only used when passing an ensemble of initial conditions.

Returns

NDArray[np.float64]

Poincaré-section points.

If q and p are 1D, returns an array of shape (num_intersections, 2 * dof + 1), where the first column is the crossing time, the next dof columns are the coordinates, and the final dof columns are the momenta.
If q and p are 2D, returns an array of shape (num_ic, num_intersections, 2 * dof + 1).

Raises

ValueError: If q and p do not have the same shape. If num_intersections is negative. If section_index is outside [0, dof). If crossing is not one of -1, 0, or 1. If the number of parameters does not match the expected number.
TypeError: If num_intersections or section_index is not an integer. If section_value is not a valid real number. If periodic_section_coordinate is not a boolean.

Compute Lyapunov exponents for a Hamiltonian system.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
total_timenumeric_t: Total integration time.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
num_exponentsint_t | None, optional: Number of Lyapunov exponents to compute. If None, the full spectrum of size 2 * degrees_of_freedom is computed.
return_historybool, optional: If True, return the time evolution of the exponents.
seedint_t, optional: Random seed used to initialize the deviation vectors.
log_basenumeric_t, optional: Base of the logarithm used in the exponent calculation.
qr_intervalint_t, optional: Number of integration steps between successive QR reorthonormalizations.
methodstr, optional: QR decomposition method: - “QR”: internal reduced modified Gram-Schmidt QR - “QR_HH”: numpy.linalg.qr based on Householder reflections

Returns

np.float64 | NDArray[np.float64]

If return_history=True, returns a 2D array whose first column contains the sampled times and whose remaining columns contain the Lyapunov exponents at those times.
If return_history=False and num_exponents == 1, returns the largest Lyapunov exponent as a scalar.
If return_history=False and num_exponents > 1, returns a 1D array containing the final Lyapunov exponents.

Raises

ValueError: If Hessian functions are missing. If q and p do not have the same shape. If total_time is not positive. If num_exponents is not in [1, 2 * degrees_of_freedom]. If log_base == 1. If method is not “QR” or “QR_HH”.
TypeError: If method is not a string. If seed, qr_interval, or num_exponents are not integers. If total_time or log_base is not a valid real number.

Compute covariant Lyapunov vectors (CLVs) for a Hamiltonian system.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
total_timenumeric_t: Total integration time over which CLVs are computed.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
num_clvsint_t | None, optional: Number of CLVs to compute. If None, computes the full set of size 2 * degrees_of_freedom.
warmup_timenumeric_t, optional: Forward warmup time used to drive the orthonormal tangent basis toward the backward Lyapunov vectors.
tail_timenumeric_t, optional: Additional forward integration time after the storage window, used to initialize the backward recursion.
qr_time_stepnumeric_t | None, optional: Time interval between successive QR factorizations. If None, defaults to the integrator time step.
seedint_t, optional: Random seed used in the tangent-basis initialization and backward recursion.
poincare_sectionbool, optional: If True, return the sampled trajectory restricted to a Poincaré section.
section_indexint_t | None, optional: Index of the coordinate defining the Poincaré section.
section_valuenumeric_t | None, optional: Value of the section coordinate.
crossingint_t | None, optional: Crossing rule: - -1 for downward crossings - 0 for all crossings - 1 for upward crossings
periodic_section_coordinatebool, optional: If True, treats q[:, section_index] as a periodic coordinate on S¹ with the given period, accumulating unbounded across samples (never re-wrapped). Crossing detection shifts the wrapped offset using delta arithmetic, mirroring generate_poincare_section. If False, uses standard Euclidean crossing detection.
periodnumeric_t, optional: Period of the angular coordinate when periodic_section_coordinate=True. Typically 2π.
methodstr, optional: QR decomposition method: - “QR”: internal reduced modified Gram-Schmidt QR - “QR_HH”: numpy.linalg.qr based on Householder reflections

Returns

tuple[NDArray[np.float64], NDArray[np.float64]]

clvs: array of shape (T, 2 * dof, num_clvs)
traj: sampled trajectory array with columns [t, q…, p…]

Raises

ValueError: If q and p do not have the same shape. If total_time is not positive. If warmup_time or tail_time is negative. If num_clvs is not in [1, 2 * degrees_of_freedom]. If qr_time_step is smaller than the integration time step. If Poincaré-section parameters are inconsistent. If method is not “QR” or “QR_HH”.
TypeError: If poincare_section is not boolean. If method is not a string. If time-like arguments are not valid real numbers.

Compute CLV-based angle diagnostics for a Hamiltonian system.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
total_timenumeric_t: Total integration time used for the angle diagnostics.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
subspacestuple[tuple[tuple[int, …], tuple[int, …]], …] | None, optional: Requested subspace-angle diagnostics. Each entry is a pair (A, B), where A and B are tuples of CLV indices defining two subspaces.
pairstuple[tuple[int, int], …] | None, optional: Requested pairwise CLV angles. Each entry (i, j) returns the angle between CLV i and CLV j.
warmup_timenumeric_t, optional: Forward QR warmup duration passed to the CLV computation.
tail_timenumeric_t, optional: Backward-recursion convergence duration passed to the CLV computation.
qr_time_stepnumeric_t | None, optional: Time spacing between successive QR factorizations. If None, defaults to the integrator time step.
seedint_t, optional: Seed used in the CLV computation.
poincare_sectionbool, optional: If True, return the trajectory restricted to a Poincaré section.
section_indexint_t, optional: Index of the coordinate defining the section.
section_valuenumeric_t, optional: Value of the section coordinate.
crossingint_t, optional: Crossing rule: - -1 for downward crossings - 0 for all crossings - 1 for upward crossings
periodic_section_coordinatebool, optional: If True, treats q[:, section_index] as a periodic coordinate on S¹ with the given period, accumulating unbounded across samples (never re-wrapped). Crossing detection shifts the wrapped offset using delta arithmetic, mirroring generate_poincare_section. If False, uses standard Euclidean crossing detection.
periodnumeric_t, optional: Period of the angular coordinate when periodic_section_coordinate=True. Typically 2π.
methodstr, optional: QR decomposition method: - “QR”: internal reduced modified Gram-Schmidt QR - “QR_HH”: numpy.linalg.qr based on Householder reflections

Returns

tuple[NDArray[np.float64], NDArray[np.float64]]

angles: array of shape (T, M) containing the requested angle time series
traj: sampled trajectory array with columns [t, q…, p…]

Raises

ValueError: If Hessian functions are missing. If q and p do not have the same shape. If total_time is not positive. If warmup_time or tail_time is negative. If qr_time_step is smaller than the integration time step. If both subspaces and pairs are missing or empty. If any subspace or pair specification is invalid. If Poincaré-section parameters are inconsistent. If method is not “QR” or “QR_HH”.
TypeError: If poincare_section is not boolean. If method is not a string.

Compute the Smaller Alignment Index (SALI).

SALI distinguishes between chaotic and regular motion by evolving two deviation vectors and monitoring their alignment over time. In chaotic motion, SALI tends exponentially to zero, while in regular motion it remains bounded away from zero.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
total_timenumeric_t: Total integration time.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
return_historybool, optional: If True, return the time evolution of SALI.
seedint_t, optional: Random seed used to initialize the deviation vectors.
thresholdnumeric_t, optional: Early stopping threshold. Integration stops when SALI <= threshold.

Returns

NDArray[np.float64]

If return_history=True, returns an array of shape (N, 2) with columns [time, SALI].
If return_history=False, returns an array of shape (2,) containing the final [time, SALI].

Raises

ValueError: If Hessian functions are missing. If q and p do not have the same shape. If total_time is not positive. If threshold is negative. If the number of parameters does not match the expected number.
TypeError: If return_history is not boolean. If seed is not an integer. If total_time or threshold is not a valid real number.

Compute the Linear Dependence Index (LDI).

LDI measures the linear dependence among k deviation vectors evolved along a trajectory. It is computed from the product of singular values of the deviation matrix. In chaotic motion, LDI tends rapidly to zero, while in regular motion it remains bounded away from zero.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
total_timenumeric_t: Total integration time.
kint_t: Number of deviation vectors to evolve.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
return_historybool, optional: If True, return the time evolution of LDI.
seedint_t, optional: Random seed used to initialize the deviation vectors.
thresholdnumeric_t, optional: Early stopping threshold. Integration stops when LDI <= threshold.

Returns

NDArray[np.float64]

If return_history=True, returns an array of shape (N, 2) with columns [time, LDI].
If return_history=False, returns an array of shape (2,) containing the final [time, LDI].

Raises

ValueError: If Hessian functions are missing. If q and p do not have the same shape. If total_time is not positive. If k is not in [2, 2 * degrees_of_freedom]. If threshold is negative. If the number of parameters does not match the expected number.
TypeError: If return_history is not boolean. If k or seed is not an integer. If total_time or threshold is not a valid real number.

Compute the Generalized Alignment Index (GALI).

GALI extends SALI by considering the evolution of k deviation vectors. It measures the volume spanned by the normalized deviation vectors. In chaotic motion, GALI tends rapidly to zero, while in regular motion it typically follows a slower decay or remains bounded.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
total_timenumeric_t: Total integration time.
kint_t: Number of deviation vectors to evolve.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
return_historybool, optional: If True, return the time evolution of GALI.
seedint_t, optional: Random seed used to initialize the deviation vectors.
thresholdnumeric_t, optional: Early stopping threshold. Integration stops when GALI <= threshold.
methodstr, optional: Method used to compute GALI: - “DET”: determinant of the Gram matrix - “QR”: product of diagonal entries from the internal QR routine - “QR_HH”: product of diagonal entries from numpy.linalg.qr

Returns

NDArray[np.float64]

If return_history=True, returns an array of shape (N, 2) with columns [time, GALI].
If return_history=False, returns an array of shape (2,) containing the final [time, GALI].

Raises

ValueError: If Hessian functions are missing. If q and p do not have the same shape. If total_time is not positive. If k is not in [2, 2 * degrees_of_freedom]. If threshold is negative. If method is not “DET”, “QR”, or “QR_HH”. If the number of parameters does not match the expected number.
TypeError: If return_history is not boolean. If k or seed is not an integer. If method is not a string. If total_time or threshold is not a valid real number.

Compute the recurrence time entropy (RTE) for a Hamiltonian system.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
num_intersectionsint_t: Number of Poincaré-section crossings used in the recurrence analysis.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
section_indexint_t, optional: Index of the coordinate used to define the Poincaré section.
section_valuenumeric_t, optional: Value of the section coordinate.
crossingint_t, optional: Crossing rule: - -1 for downward crossings - 0 for all crossings - 1 for upward crossings
periodic_section_coordinatebool, optional: If True, treats q[section_index] as a periodic coordinate on S¹ and performs crossing detection using modulo arithmetic. If False, uses standard Euclidean crossing detection.
periodnumeric_t, optional: Period of the angular coordinate when periodic_section_coordinate=True. Typically 2π for action-angle systems.
**kwargsAny: Additional keyword arguments passed to RTEConfig, including: - metric - std_metric - threshold - threshold_mode - threshold_std - lmin - return_final_state - return_recmat - return_p

Returns

float or tuple: The RTE value, optionally followed by: - the final Poincaré-section point without time - the recurrence matrix - the white-vertical-line distribution

Raises

ValueError: If q and p do not have the same shape. If num_intersections is negative. If section_index is outside [0, degrees_of_freedom). If crossing is not one of -1, 0, or 1. If the number of parameters does not match the expected number.
TypeError: If section_value is not a valid real number. If section_index, crossing, or num_intersections is not an integer.

Estimate the Hurst exponent from a Hamiltonian Poincaré section.

Parameters

qnumeric_like_t: Initial generalized coordinates. Must define a 1D array of shape (dof,).
pnumeric_like_t: Initial generalized momenta. Must define a 1D array of shape (dof,) and have the same shape as q.
num_intersectionsint_t: Number of Poincaré-section crossings used in the analysis.
parametersnumeric_like_t | None, optional: System parameters. If None, the parameters stored in the instance are used.
wminint_t, optional: Minimum window size used in the rescaled-range calculation.
section_indexint_t, optional: Index of the coordinate used to define the Poincaré section.
section_valuenumeric_t, optional: Value of the section coordinate.
crossingint_t, optional: Crossing rule: - -1 for downward crossings - 0 for all crossings - 1 for upward crossings
periodic_section_coordinatebool, optional: If True, treats q[section_index] as a periodic coordinate on S¹ and performs crossing detection using modulo arithmetic. If False, uses standard Euclidean crossing detection.
periodnumeric_t, optional: Period of the angular coordinate when periodic_section_coordinate=True. Typically 2π for action-angle systems.

Returns

NDArray[np.float64]: Estimated Hurst exponent values for the reduced Poincaré-section coordinates.

Raises

ValueError: If q and p do not have the same shape. If num_intersections is negative. If section_index is outside [0, degrees_of_freedom). If crossing is not one of -1, 0, or 1. If wmin < 2 or wmin >= num_intersections // 2. If the number of parameters does not match the expected number.
TypeError: If section_value is not a valid real number. If wmin, section_index, crossing, or num_intersections is not an integer.