Linear response

This module currently has two goals. One is calculating the first-order Jacobian, used to obtain stability and approximate (but inexpensive) the linear response of steady states. The other is calculating the full response matrix as a function of frequency; this is more accurate but more expensive.

The methodology used is explained in Jan Kosata phd thesis.

Stability

The Jacobian is used to evaluate stability of the solutions. It can be shown explicitly,

HarmonicBalance.get_Jacobian Function

get_Jacobian(eom)

Obtain the symbolic Jacobian matrix of eom. This is the linearised left-hand side of F(u) = du/dT.

source

Obtain a Jacobian from a DifferentialEquation by first converting it into a HarmonicEquation.

source

Get the Jacobian of a set of equations eqs with respect to the variables vars.

source

Linear response

The response to white noise can be shown with plot_linear_response. Depending on the order argument, different methods are used.

First order

The simplest way to extract the linear response of a steady state is to evaluate the Jacobian of the harmonic equations. Each of its eigenvalues $\lambda$ describes a Lorentzian peak in the response; $\text{Re}[\lambda ]$ gives its center and $\text{Im}[\lambda ]$ its width. Transforming the harmonic variables into the non-rotating frame (that is, inverting the harmonic ansatz) then gives the response as it would be observed in an experiment.

The advantage of this method is that for a given parameter set, only one matrix diagonalization is needed to fully describe the response spectrum. However, the method is inaccurate for response frequencies far from the frequencies used in the harmonic ansatz (it relies on the response oscillating slowly in the rotating frame).

Behind the scenes, the spectra are stored using the dedicated structs Lorentzian and JacobianSpectrum.

HarmonicBalance.LinearResponse.get_jacobian_response Function

get_jacobian_response(
    res::HarmonicBalance.Result{D, S, P, F} where F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}},
    nat_var::Num,
    Ω_range,
    branch::Int64;
    show_progress
) -> Matrix

Calculate the Jacobian response spectrum for a given system. Computes the magnitude of the Jacobian response for stable solutions across specified frequency ranges.

Arguments

• res::Result: Result object containing the system's solutions

• nat_var::Num: Natural variable to evaluate in the response

• Ω_range: Range of frequencies to evaluate

• branch::Int or followed_branches::Vector{Int}: Branch number(s) to analyze

• show_progress=true: Whether to show a progress bar

• force=false: Force recalculation of spectrum even if already exists

Returns

• Array{P,2}: Complex response matrix where rows correspond to frequencies and columns to solutions

source

HarmonicBalance.LinearResponse.JacobianSpectrum Type

mutable struct JacobianSpectrum{T<:Real}

Holds a set of Lorentzian objects belonging to a variable.

Fields

• peaks::Array{HarmonicBalance.LinearResponse.Lorentzian{T}, 1} where T<:Real

Constructor

JacobianSpectrum(res::Result; index::Int, branch::Int)

source

HarmonicBalance.LinearResponse.Lorentzian Type

struct Lorentzian{T<:Real}

Holds the three parameters of a Lorentzian peak, defined as A / sqrt((ω-ω0)² + Γ²).

Fields

• ω0::Real

• Γ::Real

• A::Real

source

Higher orders

Setting order > 1 increases the accuracy of the response spectra. However, unlike for the Jacobian, here we must perform a matrix inversion for each response frequency.

HarmonicBalance.LinearResponse.get_linear_response Function

get_linear_response(
    res::HarmonicBalance.Result{D, S, P, F} where F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}},
    nat_var::Num,
    Ω_range,
    branch::Int64;
    order,
    show_progress
)

Calculate the linear response of the system for a given branch. Evaluates the linear response by solving the linear response ODE for each stable solution and input frequency in the given range.

Arguments

• res: Result object containing the system's solutions

• nat_var::Num: Natural variable to evaluate in the response

• Ω_range: Range of frequencies to evaluate

• branch::Int: Branch number to analyze

• order: Order of the response to calculate

• show_progress=true: Whether to show a progress bar

Returns

• Array{P,2}: Response matrix where rows correspond to frequencies and columns to stable solutions

source

HarmonicBalance.LinearResponse.ResponseMatrix Type

struct ResponseMatrix

Holds the compiled response matrix of a system.

Fields

• matrix::Matrix{Function}: The response matrix (compiled).

• symbols::Vector{Num}: Any symbolic variables in matrix to be substituted at evaluation.

• variables::Vector{HarmonicVariable}: The frequencies of the harmonic variables underlying matrix. These are needed to transform the harmonic variables to the non-rotating frame.

source

HarmonicBalance.LinearResponse.get_response Function

get_response(
    rmat::HarmonicBalance.LinearResponse.ResponseMatrix,
    s::OrderedCollections.OrderedDict,
    Ω
) -> Any

For rmat and a solution dictionary s, calculate the total response to a perturbative force at frequency Ω.

source

HarmonicBalance.LinearResponse.get_response_matrix Function

get_response_matrix(
    diff_eq::DifferentialEquation,
    freq::Num;
    order
) -> Matrix

Obtain the symbolic linear response matrix of a diff_eq corresponding to a perturbation frequency freq. This routine cannot accept a HarmonicEquation since there, some time-derivatives are already dropped. order denotes the highest differential order to be considered.

source

Rotating frame

HarmonicBalance.LinearResponse.eigenvalues Function

eigenvalues(res::Result, branch; class=["physical"])

Calculate the eigenvalues of the Jacobian matrix of the harmonic equations of a branch for a one dimensional sweep in the Result struct.

Arguments

• res::Result: Result object containing solutions and jacobian information

• branch: Index of the solution branch to analyze

• class=["physical"]: Filter for solution classes to include, defaults to physical solutions

Returns

• Vector of filtered eigenvalues along the solution branch

Notes

• Currently only supports 1-dimensional parameter sweeps (D=1)

• Will throw an error if branch contains NaN values

• Eigenvalues are filtered based on the specified solution classes

source

HarmonicBalance.LinearResponse.eigenvectors Function

eigenvectors(res::Result, branch; class=["physical"])

Calculate the eigenvectors of the Jacobian matrix of the harmonic equations of a branch for a one dimensional sweep in the Result struct.

Arguments

• res::Result: Result object containing solutions and jacobian information

• branch: Index of the solution branch to analyze

• class=["physical"]: Filter for solution classes to include, defaults to physical solutions

Returns

• Vector of filtered eigenvectors along the solution branch

Notes

• Currently only supports 1-dimensional parameter sweeps (D=1)

• Will throw an error if branch contains NaN values

• Eigenvectors are filtered based on the specified solution classes

source

HarmonicBalance.LinearResponse.get_rotframe_jacobian_response Function

get_rotframe_jacobian_response(
    res::HarmonicBalance.Result{D, S, P, F} where F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}},
    Ω_range,
    branch::Int64;
    show_progress,
    damping_mod
)

Calculate the rotating frame Jacobian response for a given branch. Computes the rotating frame Jacobian response by evaluating eigenvalues of the numerical Jacobian and calculating the response magnitude for each frequency in the range.

Arguments

• res::Result: Result object containing the system's solutions

• Ω_range: Range of frequencies to evaluate

• branch::Int: Branch number to analyze

• show_progress=true: Whether to show a progress bar

• damping_mod: Damping modification parameter

Returns

• Array{P,2}: Response matrix in the rotating frame

source

Plotting

HarmonicBalance.plot_linear_response Function

plot_linear_response(
    res::HarmonicBalance.Result{D, SolType} where SolType<:Number,
    nat_var::Num,
    branch::Int64;
    Ω_range,
    order,
    logscale,
    show_progress,
    kwargs...
)

Plot the linear response to white noise of the variable nat_var for Result res on branch identifier.

Keyword arguments

• Ω_range: Range of frequency of the noise probe

• order: Order of slow-time derivatives to keep (default: 1)

• logscale: Whether to plot response in log scale (default: false)

• show_progress: Show progress bar during computation (default: true)

• kwargs...: Additional arguments passed to Plots.heatmap

Returns

A Plots.jl heatmap showing the linear response magnitude across parameter and frequency space.

source

plot_linear_response(
    res::HarmonicBalance.Result,
    nat_var::Num,
    followed_branches::Vector{Int64};
    Ω_range,
    logscale,
    show_progress,
    switch_axis,
    force,
    kwargs...
)

Plot the linear response to white noise of the variable nat_var for Result res on the followed_branches identifiers with the size of Ω_range.

Keyword arguments

• Ω_range: Range of frequency of the noise probe

• order: Order of slow-time derivatives to keep (default: 1)

• logscale: Whether to plot response in log scale (default: false)

• show_progress: Show progress bar during computation (default: true)

• kwargs...: Additional arguments passed to Plots.heatmap

Returns

A Plots.jl heatmap showing the linear response magnitude across parameter and frequency space.

source

HarmonicBalance.plot_rotframe_jacobian_response Function

plot_rotframe_jacobian_response(
    res::HarmonicBalance.Result{D, S, P, F} where F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}},
    branch::Int64;
    Ω_range,
    logscale,
    damping_mod,
    show_progress,
    kwargs...
)

Plot the linear response to white noise in the rotating frame defined the harmonic ansatz for Result res on branch identifier.

Keyword arguments

• Ω_range: Range of frequencies to analyze

• logscale: Whether to plot response in log scale (default: true)

• damping_mod: Multiplier for the real part of Jacobian eigenvalues (default: 1.0)

• show_progress: Show progress bar during computation (default: true)

• kwargs...: Additional arguments passed to Plots.heatmap

Returns

A Plots.jl heatmap showing the response magnitude in the rotating frame.

Notes

• Setting damping_mod < 1 can help distinguish between peaks with similar frequencies

• Solutions not belonging to the physical class are ignored

source

HarmonicBalance.plot_eigenvalues Function

plot_eigenvalues(
    res::HarmonicBalance.Result{D, S, P, F} where F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}},
    branch::Int64;
    class,
    type,
    projection,
    cscheme,
    kwargs...
) -> Any

Visualize the eigenvalues of the Jacobian in the rotating frame for branch identifier in the Result res.

Keyword arguments

• class: Array of solution classes to include (default: ["physical"])

• type: Which part of eigenvalues to plot (:real or :imag, default: :imag)

• projection: Function mapping eigenvectors to colors (default: v->1)

• cscheme: Color scheme for plotting (:default or custom scheme)

• kwargs...: Additional arguments passed to Plots.scatter

Returns

A scatter plot of eigenvalues colored by the projection of their eigenvectors.

Example

# Plot imaginary parts of eigenvalues
plot_eigenvalues(result, branch=1)

# Plot real parts with custom coloring based on the norm of eigenvectors of the first harmonic
plot_eigenvalues(result, branch=1, type=:real, projection=v->sqrt(v[1]^2+v[2]^2))

source