API

System objects and types

The derivative of f w.r.t. x of degree deg

HarmonicBalance.DifferentialEquation Type

julia

mutable struct DifferentialEquation

Holds differential equation(s) of motion and a set of harmonics to expand each variable. This is the primary input for HarmonicBalance.jl. After inputting the equations, the harmonics ansatz needs to be specified using add_harmonic!.

Fields

equations::OrderedCollections.OrderedDict{Num, Equation}: Assigns to each variable an equation of motion.
harmonics::OrderedCollections.OrderedDict{Num, OrderedCollections.OrderedSet{Num}}: Assigns to each variable a set of harmonics.

Example

julia

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

# equivalent ways to enter the simple harmonic oscillator
julia> DifferentialEquation(d(x,t,2) + ω0^2 * x - F * cos(ω*t), x);
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.HarmonicVariable Type

julia

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.
type::String: Type of the variable (u or v for quadratures, a for a constant, Hopf for Hopf etc.)
ω::Num: The harmonic being described.
natural_variable::Num: The natural variable whose harmonic is being described.

source

HarmonicBalance.HarmonicEquation Type

julia

mutable struct HarmonicEquation

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

Fields

equations::Vector{Equation}: A set of equations governing the harmonics.
variables::Vector{HarmonicVariable}: A set of variables describing the harmonics.
parameters::Vector{Num}: The parameters of the equation set.
natural_equation::DifferentialEquation: The natural equation (before the harmonic ansatz was used).
jacobian::Matrix{Num}: The Jacobian of the natural equation.

source

HarmonicBalance.rearrange_standard Function

julia

rearrange_standard(
    eom::HarmonicEquation
) -> HarmonicEquation

Rearrange eom to the standard form, such that the derivatives of the variables are on one side.

source

HarmonicBalance.rearrange_standard! Function

julia

rearrange_standard!(eom::DifferentialEquation)
rearrange_standard!(eom::DifferentialEquation, degree)

Rearranges the differential equations in eom to standard form, where the highest derivative of each variable (specified by degree, default 2) appears isolated on the left-hand side. Modifies the equations in place.

source

julia

rearrange_standard!(
    eom::HarmonicEquation
) -> HarmonicEquation

Rearrange eom to the standard form, such that the derivatives of the variables are on one side.

source

HarmonicBalance.first_order_transform! Function

julia

first_order_transform!(diff_eom::DifferentialEquation, time)

Transforms a higher-order differential equation system into an equivalent first-order system by introducing additional variables. Modifies the system in place. The time parameter specifies the independent variable used for differentiation.

source

HarmonicBalance.is_rearranged_standard Function

julia

is_rearranged_standard(eom::DifferentialEquation) -> Any
is_rearranged_standard(
    eom::DifferentialEquation,
    degree
) -> Any

Checks if the differential equations in eom are arranged in standard form, where the highest derivative of each variable appears isolated on the left-hand side. The default degree is 2, corresponding to second-order differential equations.

source

HarmonicBalance.get_equations Function

julia

get_equations(eom::DifferentialEquation) -> Vector{Equation}

Return the equations of eom.

source

HarmonicBalance.get_harmonic_equations Function

julia

get_harmonic_equations(
    diff_eom::DifferentialEquation;
    fast_time,
    slow_time,
    degree,
    jacobian
) -> HarmonicEquation

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.

Example

julia

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

HarmonicBalance.KrylovBogoliubov.get_krylov_equations Function

julia

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.

Krylov-Bogoliubov averaging method can be applied up to order = 2.

Example

julia

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: ω, F, ω0

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

Harmonic equations:

((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

HarmonicBalance.add_harmonic! Function

julia

add_harmonic!(diff_eom::DifferentialEquation, var::Num, ω)

Add the harmonic ω to the harmonic ansatz used to expand the variable var in diff_eom.

Example

define the simple harmonic oscillator and specify that x(t) oscillates with frequency ω

julia

julia> @variables t, x(t), y(t), ω0, ω, F, k;
julia> diff_eq = DifferentialEquation(d(x,t,2) + ω0^2 * x ~ F * cos(ω*t), x);
julia> add_harmonic!(diff_eq, x, ��) # expand x using ω

System of 1 differential equations
Variables:       x(t)
Harmonic ansatz: x(t) => ω;

(ω0^2)*x(t) + Differential(t)(Differential(t)(x(t))) ~ F*cos(t*ω)

source

HarmonicBalance.get_independent_variables Function

julia

get_independent_variables(
    diff_eom::DifferentialEquation
) -> Any

Return the independent dependent variables of diff_eom.

source

julia

get_independent_variables(
    eom::HarmonicEquation
) -> Vector{Num}

Return the independent variables (typically time) of eom.

source

Symbolics.get_variables Function

julia

get_variables(diff_eom::DifferentialEquation) -> Vector{Num}

Return the dependent variables of diff_eom.

source

julia

get_variables(eom::HarmonicEquation) -> Vector{Num}

Get the internal symbols of the independent variables of eom.

source

Solving and transforming solutions

HarmonicBalance.get_steady_states Function

julia

get_steady_states(problem::HarmonicEquation,
                    method::HarmonicBalanceMethod,
                    swept_parameters,
                    fixed_parameters;
                    show_progress=true,
                    sorting="nearest",
                    classify_default=true)

Solves problem with the method over the ranges specified by swept_parameters, keeping fixed_parameters constant. swept_parameters accepts pairs mapping symbolic variables to arrays or ranges. fixed_parameters accepts pairs mapping symbolic variables to numbers.

Keyword arguments

show_progress: Indicate whether a progress bar should be displayed.
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".
classify_default: If true, the solutions will be classified using the default classification method.

Example

solving a simple harmonic oscillator $m \ddot{x} + γ \dot{x} + ω_{0}^{2} x = F \cos (ω t)$ to obtain the response as a function of $ω$

julia

# having obtained a Problem object, let's find steady states
julia> range = (ω => range(0.8, 1.2, 100) ) # 100 parameter sets to solve
julia> fixed = ParameterList(m => 1, γ => 0.01, F => 0.5, ω_0 => 1)
julia> get_steady_states(problem, range, fixed)

A steady state result for 100 parameter points

    Solution branches:   1
       of which real:    1
       of which stable:  1

    Classes: stable, physical, Hopf, binary_labels

It is also possible to perform 2-dimensional sweeps.

julia

# The swept parameters take precedence over fixed -> use the same fixed
julia> range = (ω => range(0.8,1.2,100), F => range(0.1,1.0,10) )

# The swept parameters take precedence over fixed -> the F in fixed is now ignored
julia> get_steady_states(problem, range, fixed)

A steady state result for 1000 parameter points

    Solution branches:   1
       of which real:    1
       of which stable:  1

    Classes: stable, physical, Hopf, binary_labels

source

Methods

HarmonicBalance.WarmUp Type

julia

WarmUp

The Warm Up method prepares a warmup system with the Total Degree method using the parameter at index perturbed by perturbation_size. The warmup system is used to perform a homotopy using all other systems in the parameter sweep. It is very efficient for systems with minimal bifurcation in the parameter sweep. The Warm Up method should in theory guarantee to find all solutions, however, if the start_parameters is not proper (to close to the real line) it could miss some solutions.

SeeHomotopyContinuation.jl for more information.

Fields

warm_up_method::Union{Polyhedral{T}, TotalDegree{T}} where T: Method used for the warmup system.
start_parameters::Vector: Start parameters.
thread::Bool: Boolean indicating if threading is enabled.
tracker_options::HomotopyContinuation.TrackerOptions: Options for the tracker.
endgame_options::HomotopyContinuation.EndgameOptions: Options for the endgame.
compile::Union{Bool, Symbol}: Compilation options.
seed::UInt32: Seed for random number generation.

source

HarmonicBalance.TotalDegree Type

julia

TotalDegree

The Total Degree homotopy method performs a homotopy $H (x, t) = γ t G (x) + (1 - t) F (x)$ from the trivial polynomial system $F (x) = x ᵢ^{d ᵢ} + a ᵢ$ with the maximal degree $d ᵢ$ determined by the Bezout bound. The method guarantees to find all solutions, however, it comes with a high computational cost. See HomotopyContinuation.jl for more information.

Fields

gamma::Complex: Complex multiplying factor of the start system G(x) for the homotopy
thread::Bool: Boolean indicating if threading is enabled.
tracker_options::HomotopyContinuation.TrackerOptions: Options for the tracker.
endgame_options::HomotopyContinuation.EndgameOptions: Options for the endgame.
compile::Union{Bool, Symbol}: Compilation options.
seed::UInt32: Seed for random number generation.

source

HarmonicBalance.Polyhedral Type

julia

Polyhedral

The Polyhedral homotopy method constructs a homotopy based on the polyhedral structure of the polynomial system. It is more efficient than the Total Degree method for sparse systems, meaning most of the coefficients are zero. It can be especially useful if you don't need to find the zero solutions (only_non_zero = true), resulting in a speed up. See HomotopyContinuation.jl for more information.

Fields

only_non_zero::Bool: Boolean indicating if only non-zero solutions are considered.
thread::Bool: Boolean indicating if threading is enabled.
tracker_options::HomotopyContinuation.TrackerOptions: Options for the tracker.
endgame_options::HomotopyContinuation.EndgameOptions: Options for the endgame.
compile::Union{Bool, Symbol}: Compilation options.
seed::UInt32: Seed for random number generation.

source

Access solutions

HarmonicBalance.get_solutions Function

julia

get_solutions(
    res::Result, x::String;
    branches=1:branch_count(res), realify=false, class=["stable"], not_class=[]
    )
get_solutions(res::Result; branches=1:branch_count(res), class=["stable"], not_class=[])

Extract solution vectors from a Result object based on specified filtering criteria given by the class keywords. The first method allows extracting a specific solution component by name x. The second method returns complete solution vectors.

Keyword arguments

branches=1:branch_count(res): Range of branches to include in the output
realify=false: Whether to convert complex solutions to real form
class=["physical", "stable"]: Array of classification labels to include
not_class=[]: Array of classification labels to exclude

Returns

Filtered solution vectors matching the specified criteria

source

HarmonicBalance.attractors Function

julia

attractors(res::Result{D}; class="stable", not_class=[]) where D

Extract attractors from a Result object. Returns an array of dictionaries, where each dictionary maps branch identifier to the attractor. The attractors are filtered by their corresponding class.

Keyword arguments

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

Returns

Array{Dict,D}: Vector of dictionaries mapping branch indices to points satisfying the stability criteria at each parameter value

source

HarmonicBalance.phase_diagram Function

julia

phase_diagram(res::Result{D}; class="physical", not_class=[]) where {D}

Calculate the phase diagram from a Result object by summing over the number of states at each swept parameters.

Keyword arguments

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

Returns

Array{Int64,D}: Sum of states after applying the specified class masks

source

HarmonicBalance.get_single_solution Function

julia

get_single_solution(
    res::HarmonicBalance.Result{D, S, P, F} where F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}};
    branch,
    index
)

Return an ordered dictionary specifying all variables and parameters of the solution in result on branch at the position index.

source

HarmonicBalance.transform_solutions Function

julia

transform_solutions(
    res::HarmonicBalance.Result{D, S, ParType, F} where {ParType<:Number, F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}}},
    func;
    branches,
    realify
) -> Vector

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

Classify

HarmonicBalance.classify_solutions! Function

julia

classify_solutions!(
    res::HarmonicBalance.Result,
    func::Union{Function, String},
    name::String;
    physical
) -> Any

Creates a solution class in res using the function func (parsed into Symbolics.jl input). The new class is labeled with name and stored under res.classes[name]. By default, only physical (real) solutions are classified, and false is returned for the rest. To also classify complex solutions, set physical=false.

Example

julia

# solve a previously-defined problem
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

HarmonicBalance.get_class Function

julia

get_class(
    res::HarmonicBalance.Result,
    branch::Int64,
    class::String
) -> Any

Returns an array of booleans classifying branch in the solutions in res according to class.

source

julia

get_class(
    soln::HarmonicBalance.Result,
    class::String
) -> Vector

Returns an array of booleans classifying each branch in the solutions in res according to class.

source

Plotting

RecipesBase.plot Function

julia

plot(
    res::HarmonicBalance.Result{D, S, P, F} where F<:FunctionWrappers.FunctionWrapper{Array{S, 2}, Tuple{Array{S, 1}}},
    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().

Limit cycles

HarmonicBalance.LimitCycles.get_limit_cycles Function

julia

get_limit_cycles(
    eom::HarmonicEquation, method::HarmonicBalanceMethod, swept, fixed, ω_lc; kwargs...)

Variant of get_steady_states for a limit cycle problem characterised by a Hopf frequency (usually called ω_lc)

Solutions with ω_lc = 0 are labelled unphysical since this contradicts the assumption of distinct harmonic variables corresponding to distinct harmonics.

source

HarmonicBalance.LimitCycles.get_cycle_variables Function

julia

get_cycle_variables(
    eom::HarmonicEquation,
    ω_lc::Num
) -> Vector{HarmonicVariable}

Return the harmonic variables which participate in the limit cycle labelled by ω_lc.

source

HarmonicBalance.LimitCycles.add_pairs! Function

julia

add_pairs!(eom::DifferentialEquation; ω_lc::Num, n=1)

Add a limit cycle harmonic ω_lc to the system Equivalent to adding n pairs of harmonics ω +- ω_lc for each existing ω.

source

Linear Response

HarmonicBalance.LinearResponse.eigenvalues Method

julia

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 Method

julia

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_jacobian_response Method

julia

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.get_linear_response Method

julia

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.get_rotframe_jacobian_response Method

julia

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

HarmonicBalance.get_Jacobian Function

julia

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

Extensions

OrdinaryDiffEq

HarmonicBalance.AdiabaticSweep 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{Num, Function}: Maps each swept parameter to a function.

Examples

julia

# create a sweep of parameter a from 0 to 1 over time 0 -> 100
julia> @variables a,b;
julia> sweep = AdiabaticSweep(a => [0., 1.], (0, 100));
julia> sweep[a](50)
0.5
julia> sweep[a](200)
1.0

# do the same, varying two parameters simultaneously
julia> sweep = AdiabaticSweep([a => [0.,1.], b => [0., 1.]], (0,100))

Successive sweeps can be combined,

julia

sweep1 = AdiabaticSweep(ω => [0.95, 1.0], (0, 2e4))
sweep2 = AdiabaticSweep(λ => [0.05, 0.01], (2e4, 4e4))
sweep = sweep1 + sweep2

multiple parameters can be swept simultaneously,

julia

sweep = AdiabaticSweep([ω => [0.95;1.0], λ => [5e-2;1e-2]], (0, 2e4))

and custom sweep functions may be used.

julia

ωfunc(t) = cos(t)
sweep = AdiabaticSweep(ω => ωfunc)

source

HarmonicBalance.follow_branch Function

julia

follow_branch(
    starting_branch::Int64,
    res::HarmonicBalance.Result;
    y,
    sweep,
    tf,
    ϵ
) -> Tuple{Any, Any}

Return the indexes and values following stable branches along a 1D sweep. When a no stable solutions are found (e.g. in a bifurcation), the next stable solution is calculated by time evolving the previous solution (quench).

Keyword arguments

y: Dependent variable expression (parsed into Symbolics.jl) to evaluate the followed solution branches on .
sweep: Direction for the sweeping of solutions. A right (left) sweep proceeds from the first (last) solution, ordered as the sweeping parameter.
tf: time to reach steady
ϵ: small random perturbation applied to quenched solution, in a bifurcation in order to favour convergence in cases where multiple solutions are identically accessible (e.g. symmetry breaking into two equal amplitude states)

source

HarmonicBalance.plot_1D_solutions_branch Function

julia

plot_1D_solutions_branch(
    starting_branch::Int64,
    res::HarmonicBalance.Result;
    x,
    y,
    sweep,
    tf,
    ϵ,
    class,
    not_class,
    kwargs...
)

Plot a bifurcation diagram from a continuation sweep starting from starting_branch using the Result struct res. Time integration is used to determine what follow up branch in the continuation is.

Keyword arguments

x::String: Expression for the x-axis variable
y::String: Expression for the y-axis variable
sweep::String="right": Direction to follow the branch ("right" or "left")
tf::Real=10000: Final time for time integration
ϵ::Real=1e-4: Tolerance for branch following
kwargs...: Additional plotting arguments passed to Plots.jl
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

Returns

A Plots.jl plot object containing the bifurcation diagram with the followed branch

Description

This function creates a bifurcation diagram using follow_branch. The followed branch is plotted as a dashed gray line.

source

HarmonicBalance.follow_branch Method

julia

follow_branch(
    starting_branch::Int64,
    res::HarmonicBalance.Result;
    y,
    sweep,
    tf,
    ϵ
) -> Tuple{Any, Any}

Keyword arguments

y: Dependent variable expression (parsed into Symbolics.jl) to evaluate the followed solution branches on .
sweep: Direction for the sweeping of solutions. A right (left) sweep proceeds from the first (last) solution, ordered as the sweeping parameter.
tf: time to reach steady
ϵ: small random perturbation applied to quenched solution, in a bifurcation in order to favour convergence in cases where multiple solutions are identically accessible (e.g. symmetry breaking into two equal amplitude states)

source

SteadyStateSweep

HarmonicBalance.steady_state_sweep Function

julia

steady_state_sweep(prob::SteadyStateProblem, alg::DynamicSS; varied::Pair, kwargs...)

Sweeps through a range of parameter values using a dynamic steady state solver DynamicSS of the SteadyStateDiffEq.jl package. Given a steady state problem and a parameter to vary, computes the steady state solution for each value in the sweep range. The solutions are returned as a vector where each element corresponds to the steady state found at that parameter value.

source

julia

steady_state_sweep(prob_np::NonlinearProblem, prob_ss::SteadyStateProblem,
                  alg_np, alg_ss::DynamicSS; varied::Pair, kwargs...)

Performs a parameter sweep by combining nonlinear root alg_np and steady state solvers alg_ss. For each parameter value, it first attempts a direct nonlinear root solver and checks its stability. If the solution is unstable or not found, it switches to a dynamic steady state solver. This hybrid approach is much faster then only using a steady state solver.

source

ModelingToolkit

SciMLBase.ODEProblem Type

julia

ODEProblem(
    eom::HarmonicEquation,
    fixed_parameters;
    sweep,
    u0,
    timespan,
    perturb_initial,
    kwargs...
)

Creates an ODEProblem object used by OrdinaryDiffEqTsit5.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 u0 is specified, it is used as an initial condition; otherwise the values from fixed_parameters are used.

source

julia

ODEProblem(
    eom::Union{DifferentialEquation, HarmonicEquation},
    u0,
    tspan::Tuple,
    p::AbstractDict;
    in_place,
    kwargs...
) -> Any

Creates and ModelingToolkit.ODEProblem from a DifferentialEquation or HarmonicEquation.

Example

julia

using ModelingToolkit, StaticArrays

@variables α ω ω0 F γ t x(t)
diff_eq = DifferentialEquation(
    d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x
)
add_harmonic!(diff_eq, x, ω) #
harmonic_eq = get_harmonic_equations(diff_eq)

# in place (most performant for large systems)
ODEProblem(harmonic_eq, [1.0, 0.0], (0, 100), param)

# out of place (most performant for small systems with StaticArrays)
ODEProblem(
    harmonic_eq, [1.0, 0.0], (0, 100), param;
    in_place=false, u0_constructor=x -> SVector(x...)
)

source

ModelingToolkit.ODESystem Type

julia

ODESystem(eom::HarmonicEquation) -> Any

Creates and ModelingToolkit.ODESystem from a HarmonicEquation.

Example

julia

using ModelingToolkit

@variables α ω ω0 F γ t x(t)
diff_eq = DifferentialEquation(
    d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x
)
add_harmonic!(diff_eq, x, ω) #
harmonic_eq = get_harmonic_equations(diff_eq)

sys = ODESystem(harmonic_eq)
param = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01, ω => 1.1)
ODEProblem(sys, [1.0, 0.0], (0, 100), param)

source

julia

ODESystem(diff_eq::DifferentialEquation) -> Any

Creates and ModelingToolkit.ODESystem from a DifferentialEquation.

Example

julia

using ModelingToolkit

@variables α ω ω0 F γ t x(t)
diff_eq = DifferentialEquation(
    d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x
)
sys = ODESystem(diff_eq)

param = (α => 1.0, ω0 => 1.1, F => 0.01, γ => 0.01, ω => 1.1)

ODEProblem(sys, [1.0, 0.0], (0, 100), param)

source

SciMLBase.SteadyStateProblem Type

julia

SteadyStateProblem(
    eom::HarmonicEquation,
    u0,
    p::AbstractDict;
    in_place,
    kwargs...
) -> Any

Creates and ModelingToolkit.SteadyStateProblem from a HarmonicEquation.

Example

julia

using ModelingToolkit, StaticArrays

@variables α ω ω0 F γ t x(t)
diff_eq = DifferentialEquation(
    d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x
)
add_harmonic!(diff_eq, x, ω) #
harmonic_eq = get_harmonic_equations(diff_eq)


SteadyStateProblem(harmonic_eq, [1.0, 0.0], param)

source

SciMLBase.NonlinearProblem Type

julia

NonlinearProblem(
    eom::HarmonicEquation,
    u0,
    p::AbstractDict;
    in_place,
    kwargs...
) -> Any

Creates and ModelingToolkit.NonlinearProblem from a HarmonicEquation.

Example

julia

using ModelingToolkit, StaticArrays

@variables α ω ω0 F γ t x(t)
diff_eq = DifferentialEquation(
    d(x, t, 2) + ω0^2 * x + α * x^3 + γ * d(x, t) ~ F * cos(ω * t), x
)
add_harmonic!(diff_eq, x, ω) #
harmonic_eq = get_harmonic_equations(diff_eq)


NonlinearProblem(harmonic_eq, [1.0, 0.0], param)

source

API ​

Table of contents ​

System objects and types ​

Solving and transforming solutions ​

Methods ​

Access solutions ​

Classify ​

Plotting ​

Limit cycles ​

Linear Response ​

Extensions ​

OrdinaryDiffEq ​

SteadyStateSweep ​

ModelingToolkit ​

API

Table of contents

System objects and types

Solving and transforming solutions

Methods

Access solutions

Classify

Plotting

Limit cycles

Linear Response

Extensions

OrdinaryDiffEq

SteadyStateSweep

ModelingToolkit