diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 64119195e..b81793ded 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2023-12-27T16:59:18","documenter_version":"1.2.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.9.4","generation_timestamp":"2023-12-27T18:12:18","documenter_version":"1.2.1"}} \ No newline at end of file diff --git a/dev/background/harmonic_balance/index.html b/dev/background/harmonic_balance/index.html index dc77bbcd8..22d97de09 100644 --- a/dev/background/harmonic_balance/index.html +++ b/dev/background/harmonic_balance/index.html @@ -37,4 +37,4 @@ \\ \frac{dv_2}{dT} &= \frac{1}{6 \omega_d} \left[ \left(9\omega_d^2 - \omega_0^2\right) {u_2} - \frac{\alpha}{4} \left( u_1^3 + 3 u_2^3 + 6 u_1^2 u_2 - 3 v_1^2 u_1 + 3 v_2^2 u_2 + 6 v_1^2 u_2\right) \right] \:. \end{split} - \end{align}\]

In contrast to the single-frequency ansatz [Eqs. \eqref{eq:ansatz1}], we now have 4 equations of order 3, allowing up to $3^4=81$ solutions (the number of unique real ones is again generally far smaller). The larger number of solutions is explained by higher harmonics which cannot be captured perturbatively by the single-frequency ansatz. In particular, those where the $3 \omega_d$ component is significant. Such solutions appear, e.g., for $\omega_d \approx \omega_0 / 3$ where the generated $3 \omega_d$ harmonic is close to the natural resonant frequency. See the examples for numerical results.

+ \end{align}\]

In contrast to the single-frequency ansatz [Eqs. \eqref{eq:ansatz1}], we now have 4 equations of order 3, allowing up to $3^4=81$ solutions (the number of unique real ones is again generally far smaller). The larger number of solutions is explained by higher harmonics which cannot be captured perturbatively by the single-frequency ansatz. In particular, those where the $3 \omega_d$ component is significant. Such solutions appear, e.g., for $\omega_d \approx \omega_0 / 3$ where the generated $3 \omega_d$ harmonic is close to the natural resonant frequency. See the examples for numerical results.

diff --git a/dev/background/limit_cycles/index.html b/dev/background/limit_cycles/index.html index a9b10a710..d181d956c 100644 --- a/dev/background/limit_cycles/index.html +++ b/dev/background/limit_cycles/index.html @@ -1,2 +1,2 @@ -Limit cycles · HarmonicBalance.jl
GitHub
+Limit cycles · HarmonicBalance.jl
GitHub
diff --git a/dev/background/stability_response/index.html b/dev/background/stability_response/index.html index 2f2fc81d4..dc6f29bce 100644 --- a/dev/background/stability_response/index.html +++ b/dev/background/stability_response/index.html @@ -26,4 +26,4 @@ \end{multline}\]

Keeping in mind that $L(x)_{x_0, \gamma} = L(x + \Delta)_{x_0 + \Delta, \gamma}$ and the normalization $\delta \hat{u}_{i,j}^2 + \delta \hat{v}_{i,j}^2 = 1$, we can rewrite this as

\[\begin{equation} |\chi [\delta x _i](\tilde{\omega})|^2 = \sum_{j=1}^{M_i} \left( 1 + \alpha_{i,j} \right) L(\tilde{\omega})_{\omega_{i,j} - \text{Im}[\lambda], \text{Re}[\lambda]} + \left( 1 - \alpha_{i,j} \right) L(\tilde{\omega})_{\omega_{i,j} + \text{Im}[\lambda], \text{Re}[\lambda]} -\end{equation}\]

where

\[\alpha_{i,j} = 2\left( \text{Im}[\delta \hat{u}_{i,j}] \text{Re}[\delta \hat{v}_{i,j}] - \text{Re}[\delta \hat{u}_{i,j}] \text{Im}[\delta \hat{v}_{i,j}] \right)\]

The above solution applies to every eigenvalue $\lambda$ of the Jacobian. It is now clear that the linear response function $\chi [\delta x _i](\tilde{\omega})$ contains for each eigenvalue $\lambda_r$ and harmonic $\omega_{i,j}$ :

Sidenote: As $J$ a real matrix, there is an eigenvalue $\lambda_r^*$ for each $\lambda_r$. The maximum number of peaks in the linear response is thus equal to the dimensionality of $\mathbf{u}(T)$.

The linear response of the system in the state $\mathbf{u}_0$ is thus fully specified by the complex eigenvalues and eigenvectors of $J(\mathbf{u}_0)$. In HarmonicBalance.jl, the module LinearResponse creates a set of plottable Lorentzian objects to represent this.

Check out this example of the linear response module of HarmonicBalance.jl

+\end{equation}\]

where

\[\alpha_{i,j} = 2\left( \text{Im}[\delta \hat{u}_{i,j}] \text{Re}[\delta \hat{v}_{i,j}] - \text{Re}[\delta \hat{u}_{i,j}] \text{Im}[\delta \hat{v}_{i,j}] \right)\]

The above solution applies to every eigenvalue $\lambda$ of the Jacobian. It is now clear that the linear response function $\chi [\delta x _i](\tilde{\omega})$ contains for each eigenvalue $\lambda_r$ and harmonic $\omega_{i,j}$ :

Sidenote: As $J$ a real matrix, there is an eigenvalue $\lambda_r^*$ for each $\lambda_r$. The maximum number of peaks in the linear response is thus equal to the dimensionality of $\mathbf{u}(T)$.

The linear response of the system in the state $\mathbf{u}_0$ is thus fully specified by the complex eigenvalues and eigenvectors of $J(\mathbf{u}_0)$. In HarmonicBalance.jl, the module LinearResponse creates a set of plottable Lorentzian objects to represent this.

Check out this example of the linear response module of HarmonicBalance.jl

diff --git a/dev/examples/limit_cycles/index.html b/dev/examples/limit_cycles/index.html index b885158b7..bcbca9cb0 100644 --- a/dev/examples/limit_cycles/index.html +++ b/dev/examples/limit_cycles/index.html @@ -8,4 +8,4 @@ Solution branches: 100 of which real: 4 - of which stable: 4

The results show a fourfold degeneracy of solutions. The automatically created solution class unique_cycle filters the degeneracy out,

plot(result, ω_lc)

fig1

plot(result, ω_lc, class="unique_cycle")

fig2

Driven system - coupled Duffings

Under construction, see Chapter 6.2.2 of Jan's thesis

+ of which stable: 4

The results show a fourfold degeneracy of solutions. The automatically created solution class unique_cycle filters the degeneracy out,

plot(result, ω_lc)

fig1

plot(result, ω_lc, class="unique_cycle")

fig2

Driven system - coupled Duffings

Under construction, see Chapter 6.2.2 of Jan's thesis

diff --git a/dev/examples/linear_response/index.html b/dev/examples/linear_response/index.html index 2a69342af..b7002c7d2 100644 --- a/dev/examples/linear_response/index.html +++ b/dev/examples/linear_response/index.html @@ -29,4 +29,4 @@ plot(result, "sqrt(u1^2 + v1^2)", xscale=:log) plot_linear_response(result, x, branch=1, - Ω_range=LinRange(0.9,1.1,300), order=1, logscale=true, xscale=:log)

We see that for low $F$, quasi-linear behaviour with a single Lorentzian response occurs, while for larger $F$, two peaks form in the noise response. The two peaks are strongly unequal in magnitude, which is an example of internal squeezing.

+ Ω_range=LinRange(0.9,1.1,300), order=1, logscale=true, xscale=:log)

We see that for low $F$, quasi-linear behaviour with a single Lorentzian response occurs, while for larger $F$, two peaks form in the noise response. The two peaks are strongly unequal in magnitude, which is an example of internal squeezing.

diff --git a/dev/examples/parametron/index.html b/dev/examples/parametron/index.html index 4e9c7ebbb..929d3be63 100644 --- a/dev/examples/parametron/index.html +++ b/dev/examples/parametron/index.html @@ -11,4 +11,4 @@ plot!(result, "sqrt(u1^2 + v1^2)", not_class="large")

This overrides the default plotting style (stable=full, unstable=dashed) to give

Alternatively, we may visualise all underlying solutions, including complex ones,

plot(result, "sqrt(u1^2 + v1^2)", class="all")

⠀

2D parameters

The parametrically driven oscillator boasts a stability diagram called "Arnold's tongues" delineating zones where the oscillator is stable from those where it is exponentially unstable (if the nonlinearity was absence). We can retrieve this diagram by calculating the steady states as a function of external detuning $\delta=\omega_L-\omega_0$ and the parametric drive strength $\lambda$.

To perform a 2D sweep over driving frequency $\omega$ and parametric drive strength $\lambda$, we keep fixed from before but include 2 variables in varied

varied = (ω => range(0.8, 1.2, 50), λ => range(0.001, 0.6, 50))
 result_2D = get_steady_states(harmonic_eq, varied, fixed);

Now, we count the number of solutions for each point and represent the corresponding phase diagram in parameter space. This is done using plot_phase_diagram. Only counting stable solutions,

plot_phase_diagram(result_2D, class="stable")

⠀

In addition to phase diagrams, we can plot functions of the solution. The syntax is identical to 1D plotting. Let us overlay 2 branches into a single plot,

# overlay branches with different colors
 plot(result_2D, "sqrt(u1^2 + v1^2)", branch=1, class="stable", camera=(60,-40))
-plot!(result_2D, "sqrt(u1^2 + v1^2)", branch=2, class="stable", color=:red)

⠀

Note that solutions are ordered in parameter space according to their closest neighbors. Plots can again be limited to a given class (e.g stable solutions only) through the keyword argument class.

+plot!(result_2D, "sqrt(u1^2 + v1^2)", branch=2, class="stable", color=:red)

⠀

Note that solutions are ordered in parameter space according to their closest neighbors. Plots can again be limited to a given class (e.g stable solutions only) through the keyword argument class.

diff --git a/dev/examples/simple_Duffing/index.html b/dev/examples/simple_Duffing/index.html index 49a1d23c2..1be9b71e5 100644 --- a/dev/examples/simple_Duffing/index.html +++ b/dev/examples/simple_Duffing/index.html @@ -52,4 +52,4 @@ Classes: stable, physical, Hopf, binary_labels

Although 9 branches were found in total, only 3 remain physical (real-valued). Let us visualise the amplitudes corresponding to the two harmonics, $\sqrt{U_1^2 + V_1^2}$ and $\sqrt{U_2^2 + V_2^2}$ :

p1 = plot(result, "sqrt(u1^2 + v1^2)", legend=false)
 p2 = plot(result, "sqrt(u2^2 + v2^2)")
-plot(p1, p2)

fig3

The contributions of $\omega$ and $3\omega$ are now comparable and the system shows some fairly complex behaviour! This demonstrates how an exact solution within an extended Fourier subspace [Eq. (3)] goes beyond a perturbative treatment.

+plot(p1, p2)

fig3

The contributions of $\omega$ and $3\omega$ are now comparable and the system shows some fairly complex behaviour! This demonstrates how an exact solution within an extended Fourier subspace [Eq. (3)] goes beyond a perturbative treatment.

diff --git a/dev/examples/time_dependent/index.html b/dev/examples/time_dependent/index.html index 305221020..8b6175ae4 100644 --- a/dev/examples/time_dependent/index.html +++ b/dev/examples/time_dependent/index.html @@ -67,4 +67,4 @@ time_problem = ODEProblem(harmonic_eq, initial_state, sweep=sweep, timespan=(0,2*T)) time_evo = solve(time_problem, saveat=100);

Inspecting the amplitude as a function of time,

plot(time_evo, "sqrt(u1^2 + v1^2)", harmonic_eq)

we see that initially the sweep is adiabatic as it proceeds along the steady-state branch 1. At around $T = 2E6$, an instability occurs and $u_1(T)$ starts to rapidly oscillate. At that point, the sweep is stopped. Under free time evolution, the system then settles into a limit-cycle solution where the coordinates move along closed trajectories.

By plotting the $u$ and $v$ variables against each other, we observe the limit cycle shapes in phase space,

p1 = plot(time_evo, ["u1", "v1"], harmonic_eq)
 p2 = plot(time_evo, ["u2", "v2"], harmonic_eq)
-plot(p1, p2)
+plot(p1, p2) diff --git a/dev/index.html b/dev/index.html index e5b3aa1c2..8511100f2 100644 --- a/dev/index.html +++ b/dev/index.html @@ -26,4 +26,4 @@ of which real: 3 of which stable: 2 -Classes: stable, physical, Hopf, binary_labels
plot(result, "sqrt(u1^2 + v1^2)")

⠀

Documentation

This documentation was built using Documenter.jl.

Documentation built 2023-12-27T16:59:09.733 with Julia 1.9.4
+Classes: stable, physical, Hopf, binary_labels
plot(result, "sqrt(u1^2 + v1^2)")

⠀

Documentation

This documentation was built using Documenter.jl.

Documentation built 2023-12-27T18:12:09.870 with Julia 1.9.4
diff --git a/dev/manual/Krylov-Bogoliubov_method/index.html b/dev/manual/Krylov-Bogoliubov_method/index.html index 809d5696a..03754fcc7 100644 --- a/dev/manual/Krylov-Bogoliubov_method/index.html +++ b/dev/manual/Krylov-Bogoliubov_method/index.html @@ -28,4 +28,4 @@ ((1//2)*(ω^2)*v1(T) - (1//2)*(ω0^2)*v1(T)) / ω ~ Differential(T)(u1(T)) -((1//2)*(ω0^2)*u1(T) - (1//2)*F - (1//2)*(ω^2)*u1(T)) / ω ~ Differential(T)(v1(T))source

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

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

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

diff --git a/dev/manual/entering_eom/index.html b/dev/manual/entering_eom/index.html index c12371ccc..a24b8a244 100644 --- a/dev/manual/entering_eom/index.html +++ b/dev/manual/entering_eom/index.html @@ -6,7 +6,7 @@ julia> DifferentialEquation(d(x,t,2) + ω0^2 * x ~ F * cos(ω*t), x); # two coupled oscillators, one of them driven -julia> DifferentialEquation([d(x,t,2) + ω0^2 * x - k*y, d(y,t,2) + ω0^2 * y - k*x] .~ [F * cos(ω*t), 0], [x,y]);source
HarmonicBalance.add_harmonic! — Function
add_harmonic!(
+julia> DifferentialEquation([d(x,t,2) + ω0^2 * x - k*y, d(y,t,2) + ω0^2 * y - k*x] .~ [F * cos(ω*t), 0], [x,y]);
source
HarmonicBalance.add_harmonic! — Function
add_harmonic!(
     diff_eom::DifferentialEquation,
     var::Symbolics.Num,
     ω
@@ -19,10 +19,10 @@
 Variables:       x(t)
 Harmonic ansatz: x(t) => ω;
 
-(ω0^2)*x(t) + Differential(t)(Differential(t)(x(t))) ~ F*cos(t*ω)
source
Symbolics.get_variables — Method
get_variables(
+(ω0^2)*x(t) + Differential(t)(Differential(t)(x(t))) ~ F*cos(t*ω)
source
Symbolics.get_variables — Method
get_variables(
     diff_eom::DifferentialEquation
 ) -> Vector{Symbolics.Num}
-

Return the dependent variables of diff_eom.

source
HarmonicBalance.get_independent_variables — Method
get_independent_variables(
+

Return the dependent variables of diff_eom.

source
HarmonicBalance.get_independent_variables — Method
get_independent_variables(
     diff_eom::DifferentialEquation
 ) -> Any
-

Return the independent dependent variables of diff_eom.

source
+

Return the independent dependent variables of diff_eom.

source diff --git a/dev/manual/extracting_harmonics/index.html b/dev/manual/extracting_harmonics/index.html index e44be2013..6336b1e82 100644 --- a/dev/manual/extracting_harmonics/index.html +++ b/dev/manual/extracting_harmonics/index.html @@ -21,15 +21,15 @@ (ω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)) ~ 0source
HarmonicBalance.harmonic_ansatz — Function
harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian")

Expand each variable of diff_eom using the harmonics assigned to it with time as the time variable. For each harmonic of each variable, instance(s) of HarmonicVariable are automatically created and named.

source
HarmonicBalance.slow_flow — Function
slow_flow(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)

Removes all derivatives w.r.t fast_time (and their products) in eom of power degree. In the remaining derivatives, fast_time is replaced by slow_time.

source
HarmonicBalance.fourier_transform — Function
fourier_transform(
+(ω0^2)*v1(T) - (ω^2)*v1(T) - (2//1)*ω*Differential(T)(u1(T)) ~ 0
source
HarmonicBalance.harmonic_ansatz — Function
harmonic_ansatz(eom::DifferentialEquation, time::Num; coordinates="Cartesian")

Expand each variable of diff_eom using the harmonics assigned to it with time as the time variable. For each harmonic of each variable, instance(s) of HarmonicVariable are automatically created and named.

source
HarmonicBalance.slow_flow — Function
slow_flow(eom::HarmonicEquation; fast_time::Num, slow_time::Num, degree=2)

Removes all derivatives w.r.t fast_time (and their products) in eom of power degree. In the remaining derivatives, fast_time is replaced by slow_time.

source
HarmonicBalance.fourier_transform — Function
fourier_transform(
     eom::HarmonicEquation,
     time::Symbolics.Num
 ) -> HarmonicEquation
-

Extract the Fourier components of eom corresponding to the harmonics specified in eom.variables. For each non-zero harmonic of each variable, 2 equations are generated (cos and sin Fourier coefficients). For each zero (constant) harmonic, 1 equation is generated time does not appear in the resulting equations anymore.

Underlying assumption: all time-dependences are harmonic.

source
HarmonicBalance.drop_powers — Function
drop_powers(expr, vars, deg)
+

Extract the Fourier components of eom corresponding to the harmonics specified in eom.variables. For each non-zero harmonic of each variable, 2 equations are generated (cos and sin Fourier coefficients). For each zero (constant) harmonic, 1 equation is generated time does not appear in the resulting equations anymore.

Underlying assumption: all time-dependences are harmonic.

source
HarmonicBalance.drop_powers — Function
drop_powers(expr, vars, deg)

Remove parts of expr where the combined power of vars is => deg.

Example

julia> @variables x,y;
 julia>drop_powers((x+y)^2, x, 2)
 y^2 + 2*x*y
 julia>drop_powers((x+y)^2, [x,y], 2)
 0
 julia>drop_powers((x+y)^2 + (x+y)^3, [x,y], 3)
-x^2 + y^2 + 2*x*y
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)$.

HarmonicBalance.HarmonicVariable — Type
mutable struct HarmonicVariable

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

Fields

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

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

  • type::String: Type of the variable (u or v for quadratures, a for a constant, Hopf for Hopf etc.)

  • ω::Symbolics.Num: The harmonic being described.

  • natural_variable::Symbolics.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 Problem or solved using a dynamical ODE solver.

HarmonicBalance.HarmonicEquation — Type
mutable struct HarmonicEquation

Holds a set of algebraic equations governing the harmonics of a DifferentialEquation.

Fields

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

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

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

  • natural_equation::DifferentialEquation: The natural equation (before the harmonic ansatz was used).

source
+x^2 + y^2 + 2*x*ysource

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)$.

HarmonicBalance.HarmonicVariable — Type
mutable struct HarmonicVariable

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

Fields

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

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

  • type::String: Type of the variable (u or v for quadratures, a for a constant, Hopf for Hopf etc.)

  • ω::Symbolics.Num: The harmonic being described.

  • natural_variable::Symbolics.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 Problem or solved using a dynamical ODE solver.

HarmonicBalance.HarmonicEquation — Type
mutable struct HarmonicEquation

Holds a set of algebraic equations governing the harmonics of a DifferentialEquation.

Fields

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

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

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

  • natural_equation::DifferentialEquation: The natural equation (before the harmonic ansatz was used).

source
diff --git a/dev/manual/linear_response/index.html b/dev/manual/linear_response/index.html index 765c6d352..f240959bb 100644 --- a/dev/manual/linear_response/index.html +++ b/dev/manual/linear_response/index.html @@ -1,8 +1,8 @@ Linear response (WIP) · HarmonicBalance.jl
GitHub

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

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{Symbolics.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(
+

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

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{Symbolics.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{Symbolics.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
+

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
diff --git a/dev/manual/plotting/index.html b/dev/manual/plotting/index.html index 213fb2ff2..9882408e1 100644 --- a/dev/manual/plotting/index.html +++ b/dev/manual/plotting/index.html @@ -6,10 +6,10 @@ rules, target_type ) -> Any -

Takes a Result object and a string f representing a Symbolics.jl expression. Returns an array with the values of f evaluated for the respective solutions. Additional substitution rules can be specified in rules in the format ("a" => val) or (a => val)

source

Plotting solutions

The function plot is multiple-dispatched to plot 1D and 2D datasets. In 1D, the solutions are colour-coded according to the branches obtained by sort_solutions.

RecipesBase.plot — Method
plot(res::Result, varargs...; cut, kwargs...) -> Plots.Plot
+

Takes a Result object and a string f representing a Symbolics.jl expression. Returns an array with the values of f evaluated for the respective solutions. Additional substitution rules can be specified in rules in the format ("a" => val) or (a => val)

source

Plotting solutions

The function plot is multiple-dispatched to plot 1D and 2D datasets. In 1D, the solutions are colour-coded according to the branches obtained by sort_solutions.

RecipesBase.plot — Method
plot(res::Result, varargs...; cut, kwargs...) -> Plots.Plot

Plot a Result object.

Class selection done by passing String or Vector{String} as kwarg:

class       :   only plot solutions in this class(es) ("all" --> plot everything)
 not_class   :   do not plot solutions in this class(es)

Other kwargs are passed onto Plots.gr().

See also plot!

The x,y,z arguments are Strings compatible with Symbolics.jl, e.g., y=2*sqrt(u1^2+v1^2) plots the amplitude of the first quadratures multiplied by 2.

1D plots

plot(res::Result; x::String, y::String, class="default", not_class=[], kwargs...)
-plot(res::Result, y::String; kwargs...) # take x automatically from Result

Default behaviour is to plot stable solutions as full lines, unstable as dashed.

If a sweep in two parameters were done, i.e., dim(res)==2, a one dimensional cut can be plotted by using the keyword cut were it takes a Pair{Num, Float64} type entry. For example, plot(res, y="sqrt(u1^2+v1^2), cut=(λ => 0.2)) plots a cut at λ = 0.2.

2D plots

plot(res::Result; z::String, branch::Int64, class="physical", not_class=[], kwargs...)

To make the 2d plot less chaotic it is required to specify the specific branch to plot, labeled by a Int64.

The x and y axes are taken automatically from res

source

Plotting phase diagrams

In many problems, rather than in any property of the solutions themselves, we are interested in the phase diagrams, encoding the number of (stable) solutions in different regions of the parameter space. plot_phase_diagram handles this for 1D and 2D datasets.

HarmonicBalance.plot_phase_diagram — Function
plot_phase_diagram(res::Result; kwargs...) -> Plots.Plot
+plot(res::Result, y::String; kwargs...) # take x automatically from Result

Default behaviour is to plot stable solutions as full lines, unstable as dashed.

If a sweep in two parameters were done, i.e., dim(res)==2, a one dimensional cut can be plotted by using the keyword cut were it takes a Pair{Num, Float64} type entry. For example, plot(res, y="sqrt(u1^2+v1^2), cut=(λ => 0.2)) plots a cut at λ = 0.2.

2D plots

plot(res::Result; z::String, branch::Int64, class="physical", not_class=[], kwargs...)

To make the 2d plot less chaotic it is required to specify the specific branch to plot, labeled by a Int64.

The x and y axes are taken automatically from res

source

Plotting phase diagrams

In many problems, rather than in any property of the solutions themselves, we are interested in the phase diagrams, encoding the number of (stable) solutions in different regions of the parameter space. plot_phase_diagram handles this for 1D and 2D datasets.

HarmonicBalance.plot_phase_diagram — Function
plot_phase_diagram(res::Result; kwargs...) -> Plots.Plot

Plot the number of solutions in a Result object as a function of the parameters. Works with 1D and 2D datasets.

Class selection done by passing String or Vector{String} as kwarg:

class::String       :   only count solutions in this class ("all" --> plot everything)
-not_class::String   :   do not count solutions in this class

Other kwargs are passed onto Plots.gr()

source

Plot spaghetti plot

Sometimes, it is useful to plot the quadratures of the steady states (u, v) in function of a swept parameter. This is done with plot_spaghetti.

HarmonicBalance.plot_spaghetti — Function
plot_spaghetti(res::Result; x, y, z, kwargs...)

Plot a three dimension line plot of a Result object as a function of the parameters. Works with 1D and 2D datasets.

Class selection done by passing String or Vector{String} as kwarg:

class::String       :   only count solutions in this class ("all" --> plot everything)
-not_class::String   :   do not count solutions in this class

Other kwargs are passed onto Plots.gr()

source
+not_class::String : do not count solutions in this class

Other kwargs are passed onto Plots.gr()

source

Plot spaghetti plot

Sometimes, it is useful to plot the quadratures of the steady states (u, v) in function of a swept parameter. This is done with plot_spaghetti.

HarmonicBalance.plot_spaghetti — Function
plot_spaghetti(res::Result; x, y, z, kwargs...)

Plot a three dimension line plot of a Result object as a function of the parameters. Works with 1D and 2D datasets.

Class selection done by passing String or Vector{String} as kwarg:

class::String       :   only count solutions in this class ("all" --> plot everything)
+not_class::String   :   do not count solutions in this class

Other kwargs are passed onto Plots.gr()

source
diff --git a/dev/manual/saving/index.html b/dev/manual/saving/index.html index 45e8a0bd9..9668ad140 100644 --- a/dev/manual/saving/index.html +++ b/dev/manual/saving/index.html @@ -1,5 +1,5 @@ Saving and loading · HarmonicBalance.jl
GitHub

Saving and loading

All of the types native to HarmonicBalance.jl can be saved into a .jld2 file using save and loaded using load. Most of the saving/loading is performed using the package JLD2.jl, with the addition of reinstating the symbolic variables in the HarmonicBalance namespace (needed to parse expressions used in the plotting functions) and recompiling stored functions (needed to evaluate Jacobians). As a consequence, composite objects such as Result can be saved and loaded with no loss of information.

The function export_csv saves a .csv file which can be plot elsewhere.

HarmonicBalance.save — Function
save(filename, object)
-

Saves object into .jld2 file filename (the suffix is added automatically if not entered). The resulting file contains a dictionary with a single entry.

source
HarmonicBalance.load — Function
load(filename)
-

Loads an object from filename. For objects containing symbolic expressions such as HarmonicEquation, the symbolic variables are reinstated in the HarmonicBalance namespace.

source
+

Saves object into .jld2 file filename (the suffix is added automatically if not entered). The resulting file contains a dictionary with a single entry.

source
HarmonicBalance.load — Function
load(filename)
+

Loads an object from filename. For objects containing symbolic expressions such as HarmonicEquation, the symbolic variables are reinstated in the HarmonicBalance namespace.

source
HarmonicBalance.export_csv — Function
export_csv(filename, res, branch)
+

Saves into filename a specified solution branch of the Result res.

source
diff --git a/dev/manual/solving_harmonics/index.html b/dev/manual/solving_harmonics/index.html index 142d9fae9..8d965772e 100644 --- a/dev/manual/solving_harmonics/index.html +++ b/dev/manual/solving_harmonics/index.html @@ -2,7 +2,7 @@ Solving harmonic equations · HarmonicBalance.jl
GitHub

Solving harmonic equations

Once a differential equation of motion has been defined in DifferentialEquation and converted to a HarmonicEquation, we may use the homotopy continuation method (as implemented in HomotopyContinuation.jl) to find steady states. This means that, having called get_harmonic_equations, we need to set all time-derivatives to zero and parse the resulting algebraic equations into a Problem.

Problem holds the steady-state equations, and (optionally) the symbolic Jacobian which is needed for stability / linear response calculations.

Once defined, a Problem can be solved for a set of input parameters using get_steady_states to obtain Result.

HarmonicBalance.Problem — Type
mutable struct Problem

Holds a set of algebraic equations describing the steady state of a system.

Fields

  • variables::Vector{Symbolics.Num}: The harmonic variables to be solved for.

  • parameters::Vector{Symbolics.Num}: All symbols which are not the harmonic variables.

  • system::HomotopyContinuation.ModelKit.System: The input object for HomotopyContinuation.jl solver methods.

  • jacobian::Any: The Jacobian matrix (possibly symbolic). If false, the Jacobian is ignored (may be calculated implicitly after solving).

  • eom::HarmonicEquation: The HarmonicEquation object used to generate this Problem.

Constructors

Problem(eom::HarmonicEquation; Jacobian=true) # find and store the symbolic Jacobian
 Problem(eom::HarmonicEquation; Jacobian="implicit") # ignore the Jacobian for now, compute implicitly later
 Problem(eom::HarmonicEquation; Jacobian=J) # use J as the Jacobian (a function that takes a Dict)
-Problem(eom::HarmonicEquation; Jacobian=false) # ignore the Jacobian
source
HarmonicBalance.get_steady_states — Function
get_steady_states(prob::Problem,
                     swept_parameters::ParameterRange,
                     fixed_parameters::ParameterList;
                     random_warmup=true,
@@ -32,7 +32,7 @@
        of which real:    1
        of which stable:  1
 
-    Classes: stable, physical, Hopf, binary_labels
source
HarmonicBalance.Result — Type
mutable struct Result

Stores the steady states of a HarmonicEquation.

Fields

  • solutions::Array{Vector{Vector{ComplexF64}}}: The variable values of steady-state solutions.

  • swept_parameters::OrderedCollections.OrderedDict{Symbolics.Num, Vector{Union{Float64, ComplexF64}}}: Values of all parameters for all solutions.

  • fixed_parameters::OrderedCollections.OrderedDict{Symbolics.Num, Float64}: The parameters fixed throughout the solutions.

  • problem::Problem: The Problem used to generate this.

  • classes::Dict{String, Array}: Maps strings such as "stable", "physical" etc to arrays of values, classifying the solutions (see method classify_solutions!).

  • jacobian::Union{Int64, Function}: The Jacobian with fixed_parameters already substituted. Accepts a dictionary specifying the solution. If problem.jacobian is a symbolic matrix, this holds a compiled function. If problem.jacobian was false, this holds a function that rearranges the equations to find J only after numerical values are inserted (preferable in cases where the symbolic J would be very large).

source

Classifying solutions

The solutions in Result are accompanied by similarly-sized boolean arrays stored in the dictionary Result.classes. The classes can be used by the plotting functions to show/hide/label certain solutions.

By default, classes "physical", "stable" and "binary_labels" are created. User-defined classification is possible with classify_solutions!.

HarmonicBalance.Result — Type
mutable struct Result

Stores the steady states of a HarmonicEquation.

Fields

  • solutions::Array{Vector{Vector{ComplexF64}}}: The variable values of steady-state solutions.

  • swept_parameters::OrderedCollections.OrderedDict{Symbolics.Num, Vector{Union{Float64, ComplexF64}}}: Values of all parameters for all solutions.

  • fixed_parameters::OrderedCollections.OrderedDict{Symbolics.Num, Float64}: The parameters fixed throughout the solutions.

  • problem::Problem: The Problem used to generate this.

  • classes::Dict{String, Array}: Maps strings such as "stable", "physical" etc to arrays of values, classifying the solutions (see method classify_solutions!).

  • jacobian::Union{Int64, Function}: The Jacobian with fixed_parameters already substituted. Accepts a dictionary specifying the solution. If problem.jacobian is a symbolic matrix, this holds a compiled function. If problem.jacobian was false, this holds a function that rearranges the equations to find J only after numerical values are inserted (preferable in cases where the symbolic J would be very large).

source

Classifying solutions

The solutions in Result are accompanied by similarly-sized boolean arrays stored in the dictionary Result.classes. The classes can be used by the plotting functions to show/hide/label certain solutions.

By default, classes "physical", "stable" and "binary_labels" are created. User-defined classification is possible with classify_solutions!.

HarmonicBalance.classify_solutions! — Function
classify_solutions!(
     res::Result,
     condition::String,
     name::String;
@@ -42,9 +42,9 @@
 res = get_steady_states(problem, swept_parameters, fixed_parameters)
 
 # classify, store in result.classes["large_amplitude"]
-classify_solutions!(res, "sqrt(u1^2 + v1^2) > 1.0" , "large_amplitude")
source

Sorting solutions

Solving a steady-state problem over a range of parameters returns a solution set for each parameter. For a continuous change of parameters, each solution in a set usually also changes continuously; it is said to form a ''solution branch''. For an example, see the three colour-coded branches for the Duffing oscillator in Example 1.

For stable states, the branches describe a system's behaviour under adiabatic parameter changes.

Therefore, after solving for a parameter range, we want to order each solution set such that the solutions' order reflects the branches.

The function sort_solutions goes over the the raw output of get_steady_states and sorts each entry such that neighboring solution sets minimize Euclidean distance.

Currently, sort_solutions is compatible with 1D and 2D arrays of solution sets.

Sorting solutions

Solving a steady-state problem over a range of parameters returns a solution set for each parameter. For a continuous change of parameters, each solution in a set usually also changes continuously; it is said to form a ''solution branch''. For an example, see the three colour-coded branches for the Duffing oscillator in Example 1.

For stable states, the branches describe a system's behaviour under adiabatic parameter changes.

Therefore, after solving for a parameter range, we want to order each solution set such that the solutions' order reflects the branches.

The function sort_solutions goes over the the raw output of get_steady_states and sorts each entry such that neighboring solution sets minimize Euclidean distance.

Currently, sort_solutions is compatible with 1D and 2D arrays of solution sets.

HarmonicBalance.sort_solutions — Function
sort_solutions(
     solutions::Array;
     sorting,
     show_progress
 ) -> Any
-

Sorts solutions into branches according to the method sorting.

solutions is an n-dimensional array of Vector{Vector}. Each element describes a set of solutions for a given parameter set. The output is a similar array, with each solution set rearranged such that neighboring solution sets have the smallest Euclidean distance.

Keyword arguments

  • sorting: the method used by sort_solutions to get continuous solutions branches. The current options are "hilbert" (1D sorting along a Hilbert curve), "nearest" (nearest-neighbor sorting) and "none".
  • show_progress: Indicate whether a progress bar should be displayed.
source
+

Sorts solutions into branches according to the method sorting.

solutions is an n-dimensional array of Vector{Vector}. Each element describes a set of solutions for a given parameter set. The output is a similar array, with each solution set rearranged such that neighboring solution sets have the smallest Euclidean distance.

Keyword arguments

source diff --git a/dev/manual/time_dependent/index.html b/dev/manual/time_dependent/index.html index dd6f7e228..58797aed2 100644 --- a/dev/manual/time_dependent/index.html +++ b/dev/manual/time_dependent/index.html @@ -5,7 +5,7 @@ x0::Vector, sweep::ParameterSweep, timespan::Tuple - )

Creates an ODEProblem object used by OrdinaryDiffEq.jl from the equations in eom to simulate time-evolution within timespan. fixed_parameters must be a dictionary mapping parameters+variables to numbers (possible to use a solution index, e.g. solutions[x][y] for branch y of solution x). If x0 is specified, it is used as an initial condition; otherwise the values from fixed_parameters are used.

source
HarmonicBalance.TimeEvolution.ParameterSweep — Type

Represents a sweep of one or more parameters of a HarmonicEquation. During a sweep, the selected parameters vary linearly over some timespan and are constant elsewhere.

Sweeps of different variables can be combined using +.

Fields

  • functions::Dict{Symbolics.Num, Function}: Maps each swept parameter to a function.

Examples

# create a sweep of parameter a from 0 to 1 over time 0 -> 100
+        )

Creates an ODEProblem object used by OrdinaryDiffEq.jl from the equations in eom to simulate time-evolution within timespan. fixed_parameters must be a dictionary mapping parameters+variables to numbers (possible to use a solution index, e.g. solutions[x][y] for branch y of solution x). If x0 is specified, it is used as an initial condition; otherwise the values from fixed_parameters are used.

source
HarmonicBalance.TimeEvolution.ParameterSweep — Type

Represents a sweep of one or more parameters of a HarmonicEquation. During a sweep, the selected parameters vary linearly over some timespan and are constant elsewhere.

Sweeps of different variables can be combined using +.

Fields

  • functions::Dict{Symbolics.Num, Function}: Maps each swept parameter to a function.

Examples

# create a sweep of parameter a from 0 to 1 over time 0 -> 100
 julia> @variables a,b;
 julia> sweep = ParameterSweep(a => [0., 1.], (0, 100));
 julia> sweep[a](50)
@@ -14,11 +14,11 @@
 1.0
 
 # do the same, varying two parameters simultaneously
-julia> sweep = ParameterSweep([a => [0.,1.], b => [0., 1.]], (0,100))
source

Plotting

RecipesBase.plot — Method
plot(soln::ODESolution, f::String, harm_eq::HarmonicEquation; kwargs...)

Plot a function f of a time-dependent solution soln of harm_eq.

As a function of time

plot(soln::ODESolution, f::String, harm_eq::HarmonicEquation; kwargs...)

f is parsed by Symbolics.jl

parametric plots

plot(soln::ODESolution, f::Vector{String}, harm_eq::HarmonicEquation; kwargs...)

Parametric plot of f[1] against f[2]

Also callable as plot!

source

Miscellaneous

Using a time-dependent simulation can verify solution stability in cases where the Jacobian is too expensive to compute.

HarmonicBalance.TimeEvolution.is_stable — Function
is_stable(
+julia> sweep = ParameterSweep([a => [0.,1.], b => [0., 1.]], (0,100))
source

Plotting

RecipesBase.plot — Method
plot(soln::ODESolution, f::String, harm_eq::HarmonicEquation; kwargs...)

Plot a function f of a time-dependent solution soln of harm_eq.

As a function of time

plot(soln::ODESolution, f::String, harm_eq::HarmonicEquation; kwargs...)

f is parsed by Symbolics.jl

parametric plots

plot(soln::ODESolution, f::Vector{String}, harm_eq::HarmonicEquation; kwargs...)

Parametric plot of f[1] against f[2]

Also callable as plot!

source

Miscellaneous

Using a time-dependent simulation can verify solution stability in cases where the Jacobian is too expensive to compute.

HarmonicBalance.TimeEvolution.is_stable — Function
is_stable(
     soln::OrderedCollections.OrderedDict{Symbolics.Num, ComplexF64},
     eom::HarmonicEquation;
     timespan,
     tol,
     perturb_initial
 )
-

Numerically investigate the stability of a solution soln of eom within timespan. The initial condition is displaced by perturb_initial.

Return true the solution evolves within tol of the initial value (interpreted as stable).

source
+

Numerically investigate the stability of a solution soln of eom within timespan. The initial condition is displaced by perturb_initial.

Return true the solution evolves within tol of the initial value (interpreted as stable).

source