FFT and harmonic balance

FFT and harmonic balance#

HBMEquation: Harmonic balance for nonautonomous ODEs#

class skhippr.cycles.hbm.HBMEquation(ode: AbstractODE, omega: float, fourier: Fourier, initial_guess: ndarray = None, period_k: float = 1, stability_method: AbstractStabilityHBM = None)#

HBMEquation is a subclass of AbstractEquation which encodes the harmonic balance equations to find a periodic solution of a non-autonomous AbstractODE with periodic excitation.

This class computes a nonlinear harmonic balance residual in the frequency domain using a Fourier object. The resulting nonlinear harmonic balance equations in the frequency domain are then solved using the methods of the parent NewtonProblem class.

As the HBMEquation is a subclass of AbstractCycleEquation, attribute calls and updates are delegated to the underlying AbstractODE whenever applicable.

Attributes:

  • ode: AbstractODE object representing the non-autonomous ODE of which a periodic solution is sought.

  • omega: Frequency of the system forcing.

  • Fourier object used to perform the FFT.

  • X: Vector of Fourier coefficients of the periodic solution.

  • period_k: The period time of the sought-after periodic solution is period_k times the forcing period given by omega.

  • stability_method: The stability method.

x_time() → ndarray#

Return the time series of the solution at the FFT sample points as a self.ode.n_dof x self.fourier.L_DFT array.

aft(X=None) → ndarray#

Compute the HBM residual and its derivatives using the AFT (Alternating Frequency/Time) method.

Parameters:

X (np.ndarray, optional) – State vector in the frequency domain. If None, uses self.X.

Returns:

R – Residual vector in the frequency domain.

Return type:

np.ndarray

residual_function()#

Returns self.aft(self.X) as the residual function of the HBM equations.

closed_form_derivative(variable)#

Return the closed-form derivative of the residual function w.r.t. the variable. by Fourier transforming the corresponding derivatives in the time domain,evaluated at the samples.

hill_matrix(real_formulation: bool = None) → ndarray#

Return the Hill matrix, which is the derivative of the HBM equations w.r.t. X.

Parameters:

real_formulation (bool, optional) – If True, returns the Hill matrix in real formulation, otherwise in complex formulation. If None, uses the value of self.fourier.real_formulation.

ode_coeffs() → ndarray#

Extract and return the Fourier coefficients of the Jacobian matrix df/dx(x(t)) from the Hill matrix.

Note

  • In the real formulation, coefficients are computed according to Bayer et al. 2024, Eq.s (82)-(83) (see references), with the first N cos / sin coefficients appearing in the first row block.

  • Otherwise, the most centered coefficients are taken from the central row block .

Returns:

A 3D array of shape (n_dof, n_dof, 2*N_HBM+1) containing the Fourier coefficients of df/dx(x(t))

Return type:

np.ndarray

References

Bayer et al. 2024, https://doi.org/10.1016/j.ijnonlinmec.2024.104894

ode_samples(fourier=None) → ndarray#

Generate samples of the Jacobian matrix J(t) in time from the Hill matrix.

Parameters:

fourier (Optional) – An object providing the matrix_inv_DFT method. If None, uses self.fourier.

Returns:

The df/dx(x(t)) at the FFT samples.

Return type:

np.ndarray

stability_criterion(eigenvalues) → bool#

The stability criterion for the HBM equations is that all Floquet multipliers must have a modulus less than 1. If the ode is autonomous, the Floquet multiplier corresponding to the phase freedom is removed from the stability criterion.

exponential_decay_parameters(threshold=5e-15)#

Estimate the exponential decay of the Fourier coefficients of the Jacobian matrix J(t).

This method computes parameters a and b such that the norm of the k-th Fourier coefficient of df/dx(x(t)) is bounded by max(threshold, a * exp(-b * |k|)).

Note

These decay parameters play a central role in the error bound of Bayer and Leine, 2025 (https://arxiv.org/pdf/2503.21318).

Parameters:

threshold (float, optional) – The minimum threshold for considering the norm of the coefficients. Defaults to 5e-15.

Returns:

  • np.ndarray – an array of shape (m, 2) where each row contains a pair [a, b] for an exponential envelope which is exact at two Fourier coefficients and upperbounds the others.

  • Warns –

    UserWarning: If the norm of the coefficients does not decay below the threshold within the

    available harmonic range (N_HBM), a warning is issued suggesting to increase N_HBM.

error_bound_fundamental_matrix(t: float | ndarray = None, _as=None, bs=None)#

Attempts to compute an error bound for the fundamental solution matrix at specified time(s) t. For each time instant, the lowest bound produced by a combination of (a, b) is returned.

Notes

  • If both _as and bs are provided, exponential decay parameters are not computed from the Hill matrix and the given parameters are used directly.

  • If only one of them is provided, all applicable exponential

decay parameter combinations are computed and the ones closest to the provided values are used. * If neither is provided, all possible combinations of decay parameters are considered. At each time instant, the combination of (a, b) that produces the lowest bound is used.

Parameters:

  • t (float or np.ndarray, optional) – Time or array of times at which to evaluate the error bound. If None, uses the solution’s period.

  • _as (array_like, optional) – (Target) parameters for the exponential decay.

  • bs (array_like, optional) – (Target) parameters for the exponential decay.

Returns:

E_bound – Array of error bounds evaluated at each time in t.

Return type:

np.ndarray

HBMSystem: Treating autonomous and nonautonomous ODEs#

class skhippr.cycles.hbm.HBMSystem(ode, omega, fourier, initial_guess: ndarray = None, period_k: float = 1, stability_method: AbstractStabilityHBM = None, harmo_anchor: int = 1, dof_anchor: int = 0)#

This subclass of EquationSystem instantiates a HBMEquation and considers it as the first equation. The Fourier coefficient vector X is the first unknown.

If the underlying ODE is autonomous, the frequency omega of the periodic solution is not known in advance and is appended to the unknowns. Correspondingly, a HBMPhaseAnchor equation is appended to the equations.

class skhippr.cycles.hbm.HBMPhaseAnchor(fourier, X, harmo, dof)#

This class implements an anchor equation for the harmonic balance method (HBM) in autonomous systems to ensure that the phase of a specified degree of freedom and harmonic of the periodic solution does not change during the HBM solution procedure.

  • Complex formulation:

    exp(i*phi) = X+/X- = const

  • Real formulation:

    -tan(phi) = c_k/s_k = const

Hereby, harmo and dof specify the harmonic and degree of freedom for which the phase is anchored.

residual_function()#

Always returns zero.

closed_form_derivative(variable)#

Return the anchor as derivative w.r.t X and zero otherwise.

Fourier: Fast Fourier transform methods#

class skhippr.Fourier.Fourier(N_HBM: int, L_DFT: int, n_dof: int, real_formulation: bool = True)#

Methods for (Fast) Fourier Transform (FFT) and related operations in real or complex formulation.

Attributes:

  • n_dof: Number of states (degrees of freedom) of the considered signals.

  • N_HBM: Largest harmonic to be considered in the Fourier series.

  • L_DFT: Number of samples for the fast Fourier transform (FFT). Must be >= 2 * (N_HBM + 1) to avoid aliasing.

  • real_formulation: Whether Fourier coefficients are returned in real or complex formulation.

Useful methods:

Notes

This class supports both real-valued and complex-valued formulations.

  • For the complex-valued formulation, Fourier coefficients are represented in increasing frequency order (from -N_HBM to N_HBM).

  • For the real-valued formulation, the 0-th Fourier coefficient is stored first, then all cosine coefficients in increasing frequency order, then all sine coefficients in increasing frequency order.

time_samples(omega: float, periods: float = 1) → ndarray#

Generate a uniformly spaced vector of time samples starting at 0 with L_DFT samples per period.

  • If periods is integer, the last sample is T*periods-dt.

  • If periods is non-integer, it is rounded down to obtain an integer number of samples that are equally spaced at the DFT sampling frequency. The last sample lies in [T*periods-dt, T*periods).

Parameters:

  • omega (float) – The angular frequency (in radians per unit time).

  • periods (float, optional) – Number of periods to sample. Defaults to 1. Can be non-integer.

Returns:

1-D array of time samples.

Return type:

np.ndarray

DFT(x_samp: ndarray) → ndarray#

Compute the Fourier coefficient vector of the input signal using scipy.fft.fft.

Parameters:

x_samp (np.ndarray) – Input signal array of shape (n_dof, …, L_DFT)

Returns:

Array containing the stacked Fourier coefficients with shape (n_dof * (2* N_HBM + 1 , …).

Return type:

np.ndarray

Notes

If real_formulation is True, the real-valued Fourier coefficients are returned; otherwise, the complex-valued Fourier coefficients.

inv_DFT(X: ndarray) → ndarray#

Compute the inverse Discrete Fourier Transform (DFT) of the given Fourier coefficient vector.

This method reconstructs the time-domain samples from the provided Fourier coefficients.

Parameters:

X (np.ndarray) – Fourier coefficient vector (as generated by DFT() ). First dimension must have n_dof * (2 * N_HBM + 1) entries. May have more dimensions.

Returns:

Time-domain samples reconstructed from the Fourier coefficients, returned as an array of shape [n_dof, …, L_DFT].

Return type:

np.ndarray

matrix_DFT(A: ndarray) → ndarray#

Compute the matrix discrete Fourier transform of samples of a square matrix, as required for the HBM Jacobian, using efficient scaling and FFT.

The result of this function is equivalent to DFT_matrix @ A_block @ iDFT_matrix, where A_block is a block-diagonal matrix formed from the samples of A, but the computation is carried out efficiently using scipy.fft.fft.

Parameters:

A (np.ndarray) – Input array of shape (n_dof, n_dof, L_DFT), collecting samples of the time-periodic matrix.

Returns:

The result of the matrix DFT operation, a square array of length n_dof * (2 * N_HBM + 1).

Return type:

np.ndarray

Notes

  • The k-th column of the result is computed as the FFT of the k-th column of A_block @ iDFT_matrix.

  • Due to the block-diagonal structure, this is achieved by scaling A[:, :, l] by iDFT_small [k, l] for each l.

  • The resulting matrix has block-Toeplitz structure.

matrix_inv_DFT(A: ndarray) → ndarray#

Pseudo-inverse operation of matrix_DFT() (up to aliasing).

Parameters:

A (np.ndarray) – The result of a matrix DFT operation, a square array of length n_dof * (2 * N_HBM + 1) with block-Toeplitz structure.

Returns:

array of shape (n_dof, n_dof, L_DFT), collecting samples of the corresponding time-periodic matrix.

Return type:

np.ndarray

derivative_coeffs(X: ndarray, omega: float = 1) → ndarray#

Compute the Fourier coefficients of the derivative of a signal which is represented by Fourier coefficients.

Parameters:

  • X (np.ndarray) – A 1-D array of shape n_dof * (2 * N_HBM + 1), containing the Fourier coefficients of the original signal.

  • omega (float, optional) – The angular frequency to scale the derivative coefficients. Default is 1.

Returns:

A 1-D array of the same shape as X, containing the Fourier coefficients of the derivative of the signal.

Return type:

np.ndarray

differentiate(x_samp: ndarray, omega: float) → ndarray#

Compute the the derivative of a periodic signal in the frequency domain.

Parameters:

  • x_samp (np.ndarray) – A 2-D array of shape (n_dof, L_DFT), containing the samples of the original signal.

  • omega (float, optional) – The angular frequency of the signal. Default is 1.

Returns:

A 1-D array of the same shape as x_samp, containing the dfferentiated signal.

Return type:

np.ndarray

__replace__(**changes)#

Return a new instance of the Fourier class with updated attributes.

Parameters:

**changes (dict) – Keyword arguments (arguments of Fourier) specifying attribute values to replace in the new instance. If an attribute is not provided in changes, its value from the current instance is used.

Returns:

A new instance of the Fourier class with the specified changes applied.

Return type:

Fourier

Contents