Extracting harmonic equations

Harmonic Balance method

Once a DifferentialEquation is defined and its harmonics specified, one can extract the harmonic equations using get_harmonic_equations, which itself is composed of the subroutines harmonic_ansatz, slow_flow, fourier_transform! and drop_powers.

The harmonic equations use an additional time variable specified as slow_time in get_harmonic_equations. This is essentially a label distinguishing the time dependence of the harmonic variables (expected to be slow) from that of the oscillating terms (expected to be fast). When the equations are Fourier-transformed to remove oscillating terms, slow_time is treated as a constant. Such an approach is exact when looking for steady states.

HarmonicBalance.get_harmonic_equations Function

get_harmonic_equations(
    diff_eom::DifferentialEquation;
    fast_time,
    slow_time,
    degree,
    explicit_jacobian
) -> HarmonicEquation{DifferentialEquation}

Apply the harmonic ansatz, followed by the slow-flow, Fourier transform and dropping higher-order derivatives to obtain a set of ODEs (the harmonic equations) governing the harmonics of diff_eom.

The harmonics evolve in slow_time, the oscillating terms themselves in fast_time. If no input is used, a variable T is defined for slow_time and fast_time is taken as the independent variable of diff_eom.

By default, all products of order > 1 of slow_time-derivatives are dropped, which means the equations are linear in the time-derivatives.

explicit_jacobian fills in the symbolic Jacobian of the rearranged system, see add_jacobian!. It defaults to false: stability analysis evaluates the Jacobian implicitly and never reads this field, while deriving it needs a symbolic mass-matrix inversion that dominates the cost of this function on larger systems. Pass explicit_jacobian=true to inspect the matrix, or call get_Jacobian on the result.

Example

julia> @variables t, x(t), ω0, ω, F;

# enter the simple harmonic oscillator
julia> diff_eom = DifferentialEquation( d(x,t,2) + ω0^2 * x ~ F *cos(ω*t), x);

# expand x in the harmonic ω
julia> add_harmonic!(diff_eom, x, ω);

# get equations for the harmonics evolving in the slow time T
julia> harmonic_eom = get_harmonic_equations(diff_eom)

A set of 2 harmonic equations
Variables: u1(T), v1(T)
Parameters: ω0, ω, F

Harmonic ansatz:
x(t) = u1*cos(ωt) + v1*sin(ωt)

Harmonic equations:

(ω0^2)*u1(T) + (2//1)*ω*Differential(T)(v1(T)) - (ω^2)*u1(T) ~ F

(ω0^2)*v1(T) - (ω^2)*v1(T) - (2//1)*ω*Differential(T)(u1(T)) ~ 0

source

HarmonicVariable and HarmonicEquation types

The equations governing the harmonics are stored using the two following structs. When going from the original to the harmonic equations, the harmonic ansatz $x_i(t)=\sum _{j=1}^M{u}_{i,j}(T)\cos ({\omega }_{i,j}t)+v_{i,j}(T)\sin ({\omega }_{i,j}t)$ is used. Internally, each pair $(u_{i,j},v_{i,j})$ is stored as a HarmonicVariable. This includes the identification of ${\omega }_{i,j}$ and $x_i(t)$, which is needed to later reconstruct $x_i(t)$.

QuestBase.HarmonicVariable Type

mutable struct HarmonicVariable

Holds a variable stored under symbol describing the harmonic ω of natural_variable.

Fields

• symbol::Num: Symbol of the variable in the HarmonicBalance namespace.

• name::String: Human-readable labels of the variable, used for plotting.

• kind::QuestBase.HarmonicKind: Role of the variable in the ansatz.

• ω::Num: The harmonic being described.

• natural_variable::Num: The natural variable whose harmonic is being described.

source

When the full set of equations of motion is expanded using the harmonic ansatz, the result is stored as a HarmonicEquation. For an initial equation of motion consisting of $M$ variables, each expanded in $N$ harmonics, the resulting HarmonicEquation holds $2NM$ equations of $2NM$ variables. Each symbol not corresponding to a variable is identified as a parameter.

A HarmonicEquation can be either parsed into a steady-state HarmonicSteadyState.HomotopyContinuationProblem or solved using a dynamical ODE solver.

QuestBase.HarmonicEquation Type

mutable struct HarmonicEquation{T}

Holds a set of algebraic equations governing the harmonics of a system of equations of motion. The system it was derived from is retained in source_equations and is retrievable with source; its type is the type parameter T and is retrievable with source_type.

Typically T is a DifferentialEquation, but a HarmonicEquation can also be derived from a QuantumCumulants.MeanfieldEquations.

Fields

• equations::Vector{Equation}: A set of equations governing the harmonics.

• variables::Vector{QuestBase.HarmonicVariable}: A set of variables describing the harmonics.

• parameters::Vector{Num}: The parameters of the equation set.

• jacobian::Matrix{Num}: The Jacobian of the harmonic equations.

• source_equations::Any: The system of equations of motion the harmonic equations were derived from.

source

Krylov-Bogoliubov Averaging Method

The Krylov-Bogoliubov averaging method is an alternative high-frequency expansion technique used to analyze dynamical systems. Unlike the Harmonic Balance method, which is detailed in the background section, the Krylov-Bogoliubov method excels in computing higher orders in $1/\omega$, enabling the capture of faster dynamics within a system.

HarmonicBalance.get_krylov_equations Function

get_krylov_equations(
    diff_eom::DifferentialEquation;
    order,
    fast_time,
    slow_time
)

Apply the Krylov-Bogoliubov averaging method to a specific order to obtain a set of ODEs (the slow-flow equations) governing the harmonics of diff_eom.

The van der Pol transformation turns diff_eom into du/dT = f(u, t), with f periodic in the fast time t. A near-identity change of variables then removes the time dependence order by order, leaving the autonomous slow flow Σₙ gₙ / n!. Each order solves a homological equation ∂ₜWₙ = gₙ - Kₙ for the generator Wₙ of the Lie transform, where Kₙ collects the lower orders through Deprit's recursion (Hori-Deprit averaging): gₙ is the average of Kₙ and Wₙ the zero-mean antiderivative of its oscillating part. Both are read off the Fourier series of Kₙ in t.

The first two orders coincide with the classical Krylov-Bogoliubov expansion. Higher orders may differ from other averaging schemes by a near-identity change of the slow variables, under which steady states and their stability are unchanged.

The harmonics evolve in slow_time, the oscillating terms themselves in fast_time. If no input is used, a variable T is defined for slow_time and fast_time is taken as the independent variable of diff_eom.

Example

julia> @variables t, x(t), ω0, ω, F;

# enter the simple harmonic oscillator
julia> diff_eom = DifferentialEquation( d(x,t,2) + ω0^2 * x ~ F *cos(ω*t), x);

# expand x in the harmonic ω
julia> add_harmonic!(diff_eom, x, ω);

# get equations for the harmonics evolving in the slow time T to first order
julia> harmonic_eom = get_krylov_equations(diff_eom, order = 1)

A set of 2 harmonic equations
Variables: u1(T), v1(T)
Parameters: ω, ω0, F

Harmonic ansatz:
x(t) = u1(T)*cos(ωt) + v1(T)*sin(ωt)

Harmonic equations:

(-(1//2)*v1(T)*(ω^2) + (1//2)*v1(T)*(ω0^2)) / ω ~ Differential(T, 1)(u1(T))

((1//2)*F + (1//2)*u1(T)*(ω^2) - (1//2)*u1(T)*(ω0^2)) / ω ~ Differential(T, 1)(v1(T))

source

Purpose and Advantages

The primary advantage of the Krylov-Bogoliubov method lies in its ability to delve deeper into high-frequency components, allowing a more comprehensive understanding of fast dynamical behaviors. By leveraging this technique, one can obtain higher-order approximations that shed light on intricate system dynamics.

Any order is supported. The averaging is carried out as a Lie transform (Hori-Deprit averaging), which reduces to the classical Krylov-Bogoliubov expansion at orders 1 and 2. To first order, the slow flow coincides with the harmonic equations of get_harmonic_equations for the same ansatz $x=u\cos (\omega t)+v\sin (\omega t)$.

However, it's essential to note a limitation: this method cannot handle multiple harmonics within a single variable, unlike some other high-frequency expansion methods.

For further information and a detailed understanding of this method, refer to Krylov-Bogoliubov averaging method on Wikipedia.