Stability methods

The stability methods are implemented as subclasses of abstract AbstractStabilityMethod. Each such class defines a method which returns stability-defining eigenvalues from the derivative information of a given AbstractEquation, and a method to assert stability based on these eigenvalues.

Abstract parent class

class skhippr.stability.AbstractStabilityMethod.AbstractStabilityMethod(label: str, tol: float)

Abstract base class for stability analysis methods. This class defines the commmon interface for stability methods in SKHiPPR.

Attributes:

labelstr

A descriptive label for the stability method.

tolfloat

Tolerance value used in stability determination.

abstractmethod determine_eigenvalues(equation: AbstractEquation) → ndarray

Determine the eigenvalues that govern the stability from a given NewtonProblem.

Parameters:

equation (AbstractEquation) – The problem instance containing the (converged) system for which eigenvalues are to be computed.

Returns:

A 1-D array of eigenvalues associated with the stability of the system.

Return type:

np.ndarray

Notes

Each concrete subclass of this abstract base class can have a different interpretation of the eigenvalues, e.g. Floquet multipliers.

abstractmethod determine_stability(eigenvalues: ndarray) → bool

Determine the stability based on a set of eigenvalues as returned by self.determine_eigenvalues().

Parameters:

eigenvalues (np.ndarray) – Array of eigenvalues.

Returns:

True if the system is stable, False otherwise.

Return type:

bool

Notes

The specific stability criterion depends on the implementation. For instance, for equilibria, stability is indicated if all eigenvalues have negative real parts; for periodic solutions of nonautonomous systems, all eigenvalues (Floquet multipliers) should have magnitudes less than one.

Stability of equilibria

class skhippr.stability.AbstractStabilityMethod.StabilityEquilibrium(n_dof: int, tol: float = 0)

Class for assessing the stability of an equilibrium of a an AbstractODE.

Parameters:

• n_dof (int) – Number of degrees of freedom of the problem.

• tol (float, optional) – Tolerance for stability assessment. Default is 0.

determine_eigenvalues(ode) → ndarray

Stability is governed immediately by the eigenvalues of the Jacobian matrix, which are returned by this function.

determine_stability(eigenvalues) → bool

returns True if all eigenvalues have negative real part (smaller than self.tol).

Abstract parent class for Hill matrix methods

class skhippr.stability.AbstractStabilityHBM.AbstractStabilityHBM(label: str, fourier: Fourier, tol: float, autonomous: bool = False)

Abstract base class for stability methods for periodic solutions, applicable to HBM problems.

In addition to the methods inherited from AbstractStabilityMethod, it requires a method to compute the fundamental solution matrix of the periodic solution. With this, the Floquet multipliers (eigenvalues), which govern stability based on their magnitude, can be determined.

Attributes:

fourierFourier

Fourier object representing the FFT configuration (in particular: N_HBM and real/complex formulation) of the considered problem class.

tolfloat

Tolerance for stability checks.

autonomousbool, optional

Indicates if the system is autonomous and has one Floquet multiplier at 1 due to the freedom of phase. Default: False.

abstractmethod fundamental_matrix(t_over_period: float, hbm: HBMEquation) → ndarray

Compute the fundamental matrix at a given normalized time for a specified periodic solution.

Parameters:

• t_over_period (float) – Normalized time over the period (typically between 0 and 1).

• hbm (HBMEquation) – The (solved) HBMEquation of whose solution the fundamental matrix is sought.

Returns:

The computed fundamental matrix as a NumPy array.

Return type:

np.ndarray

determine_eigenvalues(hbm: HBMEquation) → ndarray

Determine the Floquet multipliers of the periodic solution (eigenvalues of the monodromy matrix) for the periodic solution encoded in the given hbmProblem.

determine_stability(eigenvalues) → bool

Determines the stability of a periodic solution based on the Floquet multipliers (eigenvalues).

The periodic solution is asserted stable if all Floquet multipliers lie within the unit circle.

In the autonomous case, one Floquet multiplier at (numerically) 1 is expected, irrespective of the stability properties (due to the freedom of phase). Therefore, if self.autonomous is True, the method identifies and excludes the Floquet multiplier closest to 1. If this multiplier deviates from 1 by more than self.tol, a warning is issued.

Warns:

UserWarning – In the autonomous case, if the Floquet multiplier associated with the freedom of phase is not sufficiently close to 1. This indicates that not enough harmonics were used to accurately cover the nonlinearities of the system.

error_bound(t, a, b)

Compute an error bound for the fundamental solution matrix based on the exponential decay of Fourier coefficient matrices as returned by exponential_decay_parameters.

a and b bound the norm of the Fourier coefficients by:

|| J_k || <= a*np.exp(-b*np.abs(k)).

Parameters:

• t (float) – The time at which the error bound is evaluated.

• a (float) – Factor in front of the exponential decay

• b (float) – Exponential decay

Returns:

The computed error bound for the fundamental solution matrix.

Return type:

float

Raises:

AttributeError – If the stability method does not provide an error bound.

Koopman-Hill projection methods

class skhippr.stability.KoopmanHillProjection.KoopmanHillProjection(fourier: Fourier, tol: float = 0, autonomous=False)

Direct Koopman Hill projection method for stability analysis of periodic solutions.

This subclass of AbstractStabilityHBM implements the abstract method fundamental_matrix() of its parent class using the direct Koopman-Hill projection formula:

fundamental_matrix = C @ D_time(t) @ np.expm(hill_matrix*t) @ W

Upon initialization, the projection matrices self.C and self.W are constructed once-and-for-all depending on the formulation (real or complex), and stored as attributes.

Notes

It is possible to overwrite the direct projection matrix self.C manually with a nontrivial choice (which must respect the normalization constraint, see Bayer & Leine 2023) after instantiation of the object. However, all such nontrivial choices show reduced convergence (see Bayer & Leine 2025).

Parameters:

• fourier (Fourier) – The Fourier object containing harmonic balance settings and transformation matrices.

• tol (float, optional) – Tolerance for stability computations (default is 0).

• autonomous (bool, optional) – Whether the system is autonomous (default is False).

Attributes:

Cnp.ndarray

Direct projection matrix for the Koopman Hill method. Is zero almost everywhere. In the real formulation, an identity matrix is in the first block. In the complex formulation, an identity matrix is in the central block.

Wnp.ndarray

Initial condition matrix for the Koopman Hill method. In the complex formulation, it is a stack of identity matrices. In the real formulation, it consists of an identity matrix, then a stack of identity matrices multiplied by 2, then a stack of zeros.

References

• Complex formulation: Bayer and Leine (2023): Sorting-free Hill-based stability analysis of periodic solutions through Koopman analysis. Nonlinear Dyn 111, 8439–8466, https://doi.org/10.1007/s11071-023-08247-7.

• Real formulation: Bayer et al. (2024): Koopman-Hill Stability Computation of Periodic Orbits in Polynomial Dynamical Systems Using a Real-Valued Quadratic Harmonic Balance Formulation. International Journal of Non-Linear Mechanics, 167, 104894, https://doi.org/10.1016/j.ijnonlinmec.2024.104894.

• Convergence guarantee: Bayer and Leine (2025, preprint): Explicit error bounds and guaranteed convergence of the Koopman-Hill projection stability method for linear time-periodic dynamics, https://arxiv.org/abs/2503.21318

fundamental_matrix(t_over_period: float, hbm: HBMEquation) → ndarray

Compute the fundamental solution matrix for the given periodic solution and normalized time using direct Koopman-Hill projection.

The fundamental matrix is calculated as:

self.C @ D_time(t) @ np.expm(hill_matrix * t) @ self.W

where:

• C and W are projection matrices that were computed once-and-for-all at initialization.

• hill_matrix() is obtained from the Jacobian matrix of the hbmProblem.

• D_time() scales the projection matrix – Only relevant if t_over_period is non-integer and self.C has been manually modified to be nontrivial.

Parameters:

• t_over_period (float) – The time normalized over the period (i.e., t/T, where T is the period).

• problem (HBMProblem) – A solved hbmProblem, encoding the periodic solution.

Returns:

The computed fundamental matrix as a NumPy array.

Return type:

np.ndarray

C_time(t_over_period: float) → ndarray

Compute the time-dependent, scaled projection matrix C at a given, nontrivial normalized time.

• If t_over_period is integer, this function always returns self.C.

• If self.C has not been modified manually after initialization, this function has no effect and returns self.C.

Parameters:

t_over_period (float) – The normalized time over the period (e.g., t / T), where T is the period.

Returns:

The scaled projection matrix C at the specified normalized time.

Return type:

np.ndarray

D_time(t_over_period: float) → ndarray

Constructs the time-dependent transformation matrix D(t) for the projection matrix at arbitrary (non-integer) times t.

For the complex-valued formulation, the k-th block column of self.C is scaled by np.exp(k * omega * t) and the resulting matrix is diagonal (see Bayer&Leine, 2023). The real-valued formulation follows by multiplication with the transformation matrix HBMProblem.T_to_cplx_from_real.

Parameters:

t_over_period (float) – The normalized time, expressed as a fraction of the period (t/T).

Returns:

The time-dependent transformation matrix D(t), with shape (n_dof * n_coeff, n_dof * n_coeff), where n_coeff depends on the Fourier formulation.

Return type:

np.ndarray

References

Bayer and Leine (2023): Sorting-free Hill-based stability analysis of periodic solutions through Koopman analysis. Nonlinear Dyn 111, 8439–8466, https://doi.org/10.1007/s11071-023-08247-7.

error_bound(t, a, b)

Compute the theoretical error bound for the fundamental solution matrix based on the exponential decay of Fourier coefficient matrices as returned by exponential_decay_parameters.

According to (Bayer & Leine, 2025), the error bound is given by:

|| E(t) || < (2*exp(-b))**N_HBM * exp(4*a*t)

Parameters:

• t (float) – The time at which the error bound is evaluated.

• a (float) – Factor in front of the exponential decay

• b (float) – Exponential decay

Returns:

The computed error bound for the fundamental solution matrix.

Return type:

float

class skhippr.stability.KoopmanHillProjection.KoopmanHillSubharmonic(fourier: Fourier, tol: float = 0, autonomous=False)

Subharmonic Koopman Hill projection method for stability analysis of periodic solutions.

This class modifies its parent class KoopmanHillProjection to implement the subharmonic Koopman-Hill projection formula:

fundamental_matrix = C @ D_time(t) @ np.expm(hill_matrix*t) @ W + C_subh @ D_sub(t) @ np.expm(hill_sub*t) @ W_subh

The subharmonic formulation is more accurate (error bound decays twice as fast) at the cost of approximately twice the computation time of the direct method.

Upon initialization, the projection matrices self.C and self.W as well as self.C_subh and self.W_subh are constructed once-and-for-all depending on the formulation (real or complex), and stored as attributes.

Parameters:

• fourier (Fourier) – The Fourier object containing harmonic balance settings and transformation matrices.

• tol (float, optional) – Tolerance for stability computations (default is 0).

• autonomous (bool, optional) – Whether the system is autonomous (default is False).

Attributes:

Cnp.ndarray

Projection matrix for the non-subharmonic components of the subharmonic Koopman Hill method.

C_subhndarray

Projection matrix for the subharmonic components of the subharmonic Koopman Hill method. Is the negative of self.C, with the block corresponding to the 0-th frequency removed.

Wnp.ndarray

Initial condition matrix for the non-subharmonic components of the subharmonic Koopman Hill method. Is the same as in KoopmanHillProjection.

W_subh: np.ndarray

Initial condition matrix for the subharmonic components of the subharmonic Koopman Hill method. Is equal to self.W, with the block corresponding to the 0-th frequency removed.

References

• Bayer and Leine (2025, preprint): Explicit error bounds and guaranteed convergence of the Koopman-Hill projection stability method for linear time-periodic dynamics, https://arxiv.org/abs/2503.21318

fundamental_matrix(t_over_period: float, hbm: HBMEquation) → ndarray

Compute the fundamental solution matrix for the given periodic solution and normalized time using subharmonic Koopman-Hill projection.

The fundamental matrix is calculated as:

self.C @ D_time(t) @ np.expm(hill_matrix * t) @ self.W + self.C_subh_time(t) @ np.expm(hill_subh * t) @ self.W_subh

where:

• C, W, W_subh are projection matrices that were computed once-and-for-all at initialization.

• C_subh_time() is computed by scaling C_subh for non-integer times.

• hill_matrix() is obtained from the Jacobian matrix of the hbmProblem.

• hill_subh() is obtained from the Hill matrix by eliminating the constant row and column and shifting the frequencies by 0.5.

Parameters:

• t_over_period (float) – The time normalized over the period (i.e., t/T, where T is the period).

• problem (HBMProblem) – A solved hbmProblem, encoding the periodic solution.

Returns:

The computed fundamental matrix as a NumPy array.

Return type:

np.ndarray

C_subh_time(t_over_period: float) → ndarray

Computes the scaled subharmonic projection matrix at a given normalized time.

Caution

In contrast to C_time(), the scaling must be considered even for integer t_over_period.

Parameters:

t_over_period (float) – The normalized time over the period.

Returns:

The scaled subharmonic projection matrix evaluated at the specified normalized time.

Return type:

np.ndarray

hill_subh(equ: HBMEquation) → ndarray

Constructs the subharmonic Hill matrix for the given HBM problem.

As described in (Bayer and Leine, 2025), the subharmonic Hill matrix is given by the even row and column blocks of the Hill matrix evaluated with the halved frequency.

However, practically, the subharmonic Hill matrix is constructed immediately from the hill_matrix() by removing the 0-frequency row and column and shifting all frequency terms by omega/2.

Parameters:

problem (HBMProblem) – The harmonic balance method (HBM) problem instance containing the periodic solution and the Hill matrix.

Returns:

The subharmonic Hill matrix as a NumPy array. It has n_dof fewer rows / columns than the Hill matrix itself.

Return type:

np.ndarray

References

• Bayer and Leine, 2025: Subharmonic formulation

• Bayer et al., 2024, Appendix: Details on the block structure real-valued formulation.

error_bound(t, a, b)

Compute the theoretical error bound for the fundamental solution matrix based on the exponential decay of Fourier coefficient matrices as returned by exponential_decay_parameters.

According to (Bayer & Leine, 2025), the error bound in the subharmonic formulation is given by:

|| E(t) || < (2*exp(-b))**(2*N_HBM) * exp(4*a*t)

Parameters:

• t (float) – The time at which the error bound is evaluated.

• a (float) – Factor in front of the exponential decay

• b (float) – Exponential decay

Returns:

The computed error bound for the fundamental solution matrix.

Return type:

float

Classical sorting-based Hill methods

class skhippr.stability.ClassicalHill.ClassicalHill(fourier: Fourier, sorting_method: str, tol: float = 0, autonomous: bool = False)

Stability analysis for periodic solutions by solving the Hill eigenvalue problem for the Floquet exponents with subsequent sorting.

Caution

In accordance with the other methods, the Floquet multipliers (not the Floquet exponents) are returned by determine_eigenvalues().

fundamental_matrix(t_over_period: float, hbm: HBMEquation)

Compute the fundamental matrix at a given normalized time for a specified periodic solution.

Parameters:

• t_over_period (float) – Normalized time over the period (typically between 0 and 1).

• hbm (HBMEquation) – The (solved) HBMEquation of whose solution the fundamental matrix is sought.

Returns:

The computed fundamental matrix as a NumPy array.

Return type:

np.ndarray

hill_EVP(hbm: HBMEquation, visualize: bool = False) → tuple[ndarray, ndarray]

Solves the eigenvalue problem for the Hill matrix and performs sorting to identify the Floquet exponents.

Computes the eigenvalues of hbm.hill_matrix and then chooses those that minimize self.sorting_criterion.

The Floquet exponents can optionally be visualized in the complex plane.

Parameters:

• problem (HBMProblem) – encodes the solution whose stability is sought.

• visualize (bool, optional) – If True, plots the real and imaginary parts of all eigenvalues and the selected Floquet exponents (default is False).

Returns:

• floquet_exponents (np.ndarray) – Array of computed Floquet exponents (eigenvalues).

• eigenvectors (np.ndarray) – Array of eigenvectors corresponding to the selected Floquet exponents.

Caution

This function returns the Floquet exponents, not the Floquet multipliers. The conversion to Floquet multipliers happens in determine_eigenvalues().

Notes

All sorting methods can exhibit numerical issues for negative real Floquet multipliers, which require extra care. The function currently does not explicitly handle the case for negative real Floquet multipliers specially, being error-prone in this case.

determine_eigenvalues(hbm: HBMEquation) → ndarray

Determine the eigenvalues (Floquet multipliers) for the given periodic solution.

This method computes the Floquet exponents using the Hill eigenvalue problem and converts them to Floquet multipliers.

Parameters:

problem (HBMProblem) – The periodic solution to be investigated

Returns:

Array of Floquet multipliers corresponding to the computed Floquet multipliers, converted from Floquet exponents.

Return type:

np.ndarray

Explicit Runge Kutta single-pass time integration method

class skhippr.stability.SinglePass.SinglePassRK(fourier: Fourier, A: ndarray, b: ndarray, c: ndarray, label: str, stepsize: float = None, tol: float = 0)

SinglePassRK implements single-pass fixed-step explicit Runge-Kutta methods for stability analysis in the context of Harmonic Balance Methods (HBM).

Parameters:

fourierFourier

Fourier object containing discretization and transformation information.

Anp.ndarray

Butcher tableau A (stage weights). Must be strictly lower triangular, i.e., explicit Runge Kutta method.

bnp.ndarray

Butcher tableau b (quadrature weights).

cnp.ndarray

Butcher tableau c (normalized evaluation points), must be an array of Fractions with a low common denominator.

stepsizefloat, optional

Desired step size for the integration. If not provided, it is chosen such that no additional FFTs are required.

tolfloat, optional

Tolerance for the stability computation (default is 0).

fundamental_matrix(t_over_period: float, hbm: HBMEquation) → ndarray

Computes the fundamental matrix for a given normalized time.

This method fixed-step-integrates the variational equation up to t.

Parameters:

• t_over_period (float) – The time, expressed as a multiple of the period over which to integrate.

• problem (HBMProblem) – The problem instance containing system parameters, Fourier information, and Jacobian samples.

Returns:

The fundamental matrix (state transition matrix) after integrating over the specified time.

Return type:

np.ndarray

class skhippr.stability.SinglePass.SinglePassRK4(fourier: Fourier, stepsize: float = None, tol: float = 0)

Explicit Single-Pass method using the RK4 tableau.

This method only needs one intermediate sample between each integration step, allowing for a relatively low L_DFT at small time step.

References

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Examples (accessed on 06/17/2025)

class skhippr.stability.SinglePass.SinglePassRK38(fourier: Fourier, stepsize: float = None, tol: float = 0)

Explicit Single-Pass method using the 3/8-rule.

This method needs two intermediate samples per time step.

References

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Examples (accessed on 06/17/2025)