Linear response (WIP)

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.LinearResponse.get_Jacobian — Function.

get_Jacobian(eom)

Obtain the symbolic Jacobian matrix of eom (either a HarmonicEquation or a DifferentialEquation). 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.

# HarmonicBalance.LinearResponse.plot_linear_response — Function.

plot_linear_response(res::Result, nat_var::Num; Ω_range, branch::Int, order=1, logscale=false, show_progress=true, kwargs...)

Plot the linear response to white noise of the variable nat_var for Result res on branch for input frequencies Ω_range. Slow-time derivatives up to order are kept in the process.

Any kwargs are fed to Plots' gr().

Solutions not belonging to the physical class are ignored.

source

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.JacobianSpectrum — Type.

mutable struct JacobianSpectrum

Holds a set of Lorentzian objects belonging to a variable.

Fields

• peaks::Vector{HarmonicBalance.LinearResponse.Lorentzian}

Constructor

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

source

# HarmonicBalance.LinearResponse.Lorentzian — Type.

struct Lorentzian

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

Fields

• ω0::Float64

• Γ::Float64

• A::Float64

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.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{Num, ComplexF64},
    Ω
) -> 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=2)

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