Symplectic Integrators

The symplectic module provides high-order explicit symplectic integrators for polynomial Hamiltonian systems.

Main Classes

class hiten.algorithms.integrators.symplectic._ExtendedSymplectic[source]

Bases: _Integrator

High-order explicit Tao symplectic integrator for polynomial Hamiltonian systems.

Parameters:

  • order (int, optional) – Even order of the underlying scheme (>= 2). Default is 6.

  • c_omega_heuristic (float, optional) – Scaling coefficient used in the heuristic _get_tao_omega() frequency. Default is 20.0.

  • **options – Additional keyword options stored verbatim in options.

name

Human-readable identifier, e.g. "Symplectic6".

Type:

str

order

Same as the order constructor argument.

Type:

int

c_omega_heuristic

Same as the constructor argument.

Type:

float

Examples

>>> from hiten.algorithms.integrators.symplectic import _ExtendedSymplectic
>>> integrator = _ExtendedSymplectic(order=8, c_omega_heuristic=25.0)
>>> sol = integrator.integrate(hamiltonian_system, y0, t_vals)

Notes

The target system must expose three public attributes:

  • jac_H - Jacobian of the Hamiltonian given as a nested list of

    polynomial coefficient arrays compatible with _polynomial_evaluate().

  • clmo_H - co-efficient layout mapping objects for the same

    polynomials.

  • n_dof - number of degrees of freedom (must equal 3 for this

    implementation).

property order: int

Order of accuracy of the symplectic method.

validate_system(system)[source]

Validate that the system is compatible with symplectic integration.

Parameters:

system (_HamiltonianSystemProtocol) – The system to validate

Raises:

ValueError – If the system doesn’t provide the required Hamiltonian structure

Return type:

None

integrate(system, y0, t_vals, *, event_fn=None, event_cfg=None, **kwargs)[source]

Integrate the Hamiltonian system using symplectic method with event detection.

Parameters:

  • system (_HamiltonianSystemProtocol) – The Hamiltonian system to integrate (must provide polynomial structure)

  • y0 (numpy.ndarray) – Initial state vector [Q, P], shape (2*n_dof,)

  • t_vals (numpy.ndarray) – Array of time points at which to evaluate the solution

  • event_fn (Callable[[float, numpy.ndarray], float], optional) – Scalar event function evaluated as g(t, y). A zero crossing may terminate integration or mark an event.

  • event_cfg (_EventConfig | None) – Configuration controlling event directionality, terminal behavior, and tolerances.

  • **kwargs – Additional integration options (currently unused)

Returns:

Integration results containing times and states

Return type:

_Solution

Notes

This method delegates to the existing _integrate_symplectic function, which expects the system to provide: - jac_H: Jacobian polynomial coefficients - clmo_H: Coefficient layout mapping objects - n_dof: Number of degrees of freedom

The integration preserves the symplectic structure exactly, making it ideal for long-term integration of Hamiltonian systems where energy conservation is critical.

When event detection is enabled, the main trajectory remains symplectic. Event refinement uses Hermite interpolation which is accurate but not symplectic within the final step.

Factory Classes

class hiten.algorithms.integrators.symplectic.ExtendedSymplectic[source]

Bases: object

Implement a factory for extended symplectic integrators.

This factory class provides a convenient way to create symplectic integrators of different orders without directly instantiating the underlying implementation classes.

Parameters:

  • order (int, default 6) – Order of the symplectic method. Must be 2, 4, 6, or 8.

  • **opts – Additional options passed to the integrator constructor.

Returns:

A symplectic integrator instance of the specified order.

Return type:

_ExtendedSymplectic

Raises:

ValueError – If the specified order is not supported.

Examples

>>> integrator = ExtendedSymplectic(order=6)
>>> solution = integrator.integrate(hamiltonian_system, y0, t_vals)
static __new__(cls, order=6, **opts)[source]

Create a symplectic integrator of the specified order.

Parameters:

  • order (int, default 6) – Order of the symplectic method. Must be 2, 4, 6, or 8.

  • **opts – Additional options passed to the integrator constructor.

Returns:

A symplectic integrator instance of the specified order.

Return type:

_ExtendedSymplectic

Raises:

ValueError – If the specified order is not supported.

Functions

hiten.algorithms.integrators.symplectic._get_tao_omega()[source]

Calculate the frequency parameter for the symplectic integrator.

Parameters:

  • delta (float) – Time step size

  • order (int) – Order of the symplectic integrator

  • c (float, optional) – Scaling parameter, default is 10.0

Returns:

Frequency parameter tau*omega for the symplectic scheme

Return type:

float

Notes

The calculated value scales with (c*delta)^(-order) to ensure numerical stability for larger time steps.

hiten.algorithms.integrators.symplectic._construct_6d_eval_point()[source]

Construct a 6D evaluation point from N-DOF position and momentum vectors. Assumes N_SYMPLECTIC_DOF is 3 for this specific 6D polynomial evaluation context.

Parameters:

  • Q_current_ndof (numpy.ndarray) – Position vector (dimension N_SYMPLECTIC_DOF, e.g., [q1, q2, q3])

  • P_current_ndof (numpy.ndarray) – Momentum vector (dimension N_SYMPLECTIC_DOF, e.g., [p1, p2, p3])

Returns:

6D evaluation point for polynomial evaluation, ordered [q1,q2,q3,p1,p2,p3]

Return type:

numpy.ndarray

Notes

This function maps N-DOF coordinates to a 6D vector suitable for the polynomial evaluation, which expects variables in a specific order.

hiten.algorithms.integrators.symplectic._eval_dH_dQ()[source]

Evaluate derivatives of Hamiltonian with respect to generalized position variables.

Parameters:

  • Q_eval_ndof (numpy.ndarray) – Position vector ([q1,q2,q3]) at which to evaluate derivatives

  • P_eval_ndof (numpy.ndarray) – Momentum vector ([p1,p2,p3]) at which to evaluate derivatives

  • jac_H (List[List[numpy.ndarray]]) – Jacobian of Hamiltonian as list of polynomial coefficients for 6 variables

  • clmo_H (List[numpy.ndarray]) – List of coefficient layout mapping objects for the polynomials

Returns:

Vector of partial derivatives dH/dQ (e.g., [dH/dq1, dH/dq2, dH/dq3])

Return type:

numpy.ndarray

hiten.algorithms.integrators.symplectic._eval_dH_dP()[source]

Evaluate derivatives of Hamiltonian with respect to generalized momentum variables.

Parameters:

  • Q_eval_ndof (numpy.ndarray) – Position vector ([q1,q2,q3]) at which to evaluate derivatives

  • P_eval_ndof (numpy.ndarray) – Momentum vector ([p1,p2,p3]) at which to evaluate derivatives

  • jac_H (List[List[numpy.ndarray]]) – Jacobian of Hamiltonian as list of polynomial coefficients for 6 variables

  • clmo_H (List[numpy.ndarray]) – List of coefficient layout mapping objects for the polynomials

Returns:

Vector of partial derivatives dH/dP (e.g., [dH/dp1, dH/dp2, dH/dp3])

Return type:

numpy.ndarray

hiten.algorithms.integrators.symplectic._phi_H_a_update_poly()[source]

Apply the first Hamiltonian splitting operator (phi_a) in the symplectic scheme.

Parameters:

  • q_ext (numpy.ndarray) – Extended state vector [Q, P, X, Y] to be updated in-place

  • delta (float) – Time step size

  • jac_H (List[List[numpy.ndarray]]) – Jacobian of Hamiltonian as list of polynomial coefficients

  • clmo_H (List[numpy.ndarray]) – List of coefficient layout mapping objects for the polynomials

Notes

Implements the symplectic update step: - P <- P - delta * dH/dQ(Q,Y) - X <- X + delta * dH/dP(Q,Y)

This modifies q_ext in-place through views/slices. Q, P, X, Y are now N_SYMPLECTIC_DOF dimensional.

hiten.algorithms.integrators.symplectic._phi_H_b_update_poly()[source]

Apply the second Hamiltonian splitting operator (phi_b) in the symplectic scheme.

Parameters:

  • q_ext (numpy.ndarray) – Extended state vector [Q, P, X, Y] to be updated in-place

  • delta (float) – Time step size

  • jac_H (List[List[numpy.ndarray]]) – Jacobian of Hamiltonian as list of polynomial coefficients

  • clmo_H (List[numpy.ndarray]) – List of coefficient layout mapping objects for the polynomials

Notes

Implements the symplectic update step: - Q <- Q + delta * dH/dP(X,P) - Y <- Y - delta * dH/dQ(X,P)

This modifies q_ext in-place through views/slices. Q, P, X, Y are now N_SYMPLECTIC_DOF dimensional.

hiten.algorithms.integrators.symplectic._phi_omega_H_c_update_poly()[source]

Apply the rotation operator (phi_c) in the symplectic scheme.

Parameters:

  • q_ext (numpy.ndarray) – Extended state vector [Q, P, X, Y] to be updated in-place

  • delta (float) – Time step size

  • omega (float) – Frequency parameter for the rotation

Notes

Implements a rotation in the extended phase space with mixing of coordinates. The transformation is implemented using trigonometric functions and temporary variables to ensure numerical stability.

This step is crucial for high-order symplectic integration methods with the extended phase-space technique. Q, P, X, Y are now N_SYMPLECTIC_DOF dimensional.

hiten.algorithms.integrators.symplectic._recursive_update_poly()[source]

Apply recursive symplectic update of specified order.

Parameters:

  • q_ext (numpy.ndarray) – Extended state vector [Q, P, X, Y] to be updated in-place

  • timestep (float) – Time step size

  • order (int) – Order of the symplectic integrator (must be even and >= 2)

  • omega (float) – Frequency parameter for the rotation

  • jac_H (List[List[numpy.ndarray]]) – Jacobian of Hamiltonian as list of polynomial coefficients

  • clmo_H (List[numpy.ndarray]) – List of coefficient layout mapping objects for the polynomials

Notes

For order=2, applies the basic second-order symplectic scheme:

phi_a(delta/2) o phi_b(delta/2) o phi_c(delta) o phi_b(delta/2) o phi_a(delta/2)

For higher orders, applies a recursive composition method with carefully chosen substeps to achieve the desired order of accuracy.

hiten.algorithms.integrators.symplectic._integrate_symplectic()[source]

Integrate Hamilton’s equations using a high-order symplectic integrator for a system with N_SYMPLECTIC_DOF degrees of freedom (e.g., 3 DOF for a 6D phase space).

Parameters:

  • initial_state_6d (numpy.ndarray) – Initial state vector [Q, P] (e.g., [q1,q2,q3,p1,p2,p3]) (shape: 2*N_SYMPLECTIC_DOF)

  • t_values (numpy.ndarray) – Array of time points at which to compute the solution

  • jac_H (List[List[np.ndarray]]) – Jacobian of Hamiltonian as a list of polynomial coefficients (for 6 variables)

  • clmo_H (List[numpy.ndarray]) – List of coefficient layout mapping objects for the polynomials

  • order (int) – Order of the symplectic integrator (must be even and >= 2)

  • c_omega_heuristic (float, optional) – Scaling parameter for the frequency calculation, default is 20.0

Returns:

Trajectory array of shape (len(t_values), 2*N_SYMPLECTIC_DOF)

Return type:

numpy.ndarray

Notes

Uses an extended phase-space technique to implement a high-order symplectic integration method for the polynomial Hamiltonian.

The method is particularly suitable for center manifold dynamics where energy conservation over long time integration is crucial.

The algorithm: 1. Creates an extended phase space with auxiliary variables [Q, P, X, Y] 2. Recursively applies composition of basic symplectic steps 3. Returns trajectory only for the physical variables [Q, P]

For optimal energy conservation, higher c_omega_heuristic values may be used at the cost of potentially smaller effective timesteps.

Event Functions

hiten.algorithms.integrators.symplectic._hermite_refine_event_symplectic()[source]

Refine event location using bisection on Hermite interpolant.

Performs bisection search on the interval [0, 1] to locate the precise time when the event function crosses zero.

Parameters:

  • event_fn (callable) – Compiled event function g(t, y) -> float

  • t0 (float) – Time at start of step

  • y0 (numpy.ndarray) – State at start of step

  • f0 (numpy.ndarray) – Derivative at start of step

  • t1 (float) – Time at end of step

  • y1 (numpy.ndarray) – State at end of step

  • f1 (numpy.ndarray) – Derivative at end of step

  • h (float) – Step size (t1 - t0)

  • direction (int) – Event direction (-1, 0, +1)

  • xtol (float) – Tolerance on time location

  • gtol (float) – Tolerance on event function value

Returns:

  • t_hit (float) – Time of event

  • y_hit (numpy.ndarray) – State at event

Notes

The bisection is performed on the normalized interval [0, 1], then mapped back to [t0, t1]. The Hermite interpolation provides smooth, differentiable approximation but is NOT symplectic.

hiten.algorithms.integrators.symplectic._integrate_symplectic_until_event()[source]

Integrate with symplectic method until an event is detected.

This function performs symplectic integration and monitors an event function for zero crossings. When a crossing is detected, it uses Hermite interpolation to refine the event location within the step.

Parameters:

  • initial_state_6d (numpy.ndarray) – Initial state vector [Q, P] (shape: 2*N_SYMPLECTIC_DOF)

  • t_values (numpy.ndarray) – Array of time points at which to check for events

  • jac_H (List[List[numpy.ndarray]]) – Jacobian of Hamiltonian as list of polynomial coefficients

  • clmo_H (List[numpy.ndarray]) – List of coefficient layout mapping objects

  • order (int) – Order of the symplectic integrator

  • event_fn (callable) – Compiled event function g(t, y) -> float

  • direction (int) – Event crossing direction: -1 (negative), 0 (any), +1 (positive)

  • xtol (float) – Tolerance on event time location

  • gtol (float) – Tolerance on event function value

  • c_omega_heuristic (float, optional) – Scaling parameter for frequency calculation

Returns:

  • hit (bool) – True if event was detected, False otherwise

  • t_hit (float) – Time of event (or final time if no event)

  • y_hit (numpy.ndarray) – State at event (or final state if no event)

  • trajectory (numpy.ndarray) – Trajectory up to (but not including) the event

Notes

The main integration remains fully symplectic. Event refinement uses Hermite interpolation which is NOT symplectic, but provides accurate event location within the last step.

The symplectic structure is preserved throughout the trajectory except for the final interpolated point when an event is detected.