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{Num}: The harmonic variables to be solved for.

• parameters::Vector{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;
                    method=:warmup,
                    threading = Threads.nthreads() > 1,
                    show_progress=true,
                    sorting="nearest")

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

Keyword arguments

• method: If :warmup (default), a problem similar to prob but with random complex parameters is first solved to find all non-singular paths. The subsequent tracking to find results for all swept_parameters is then much faster than the initial solving. If method=:total_degree, each parameter point is solved separately by tracking the maximum number of paths (employs a total degree homotopy).

This takes far longer but can be more reliable.

• threading: If true, multithreaded support is activated. The number of available threads is set by the environment variable JULIA_NUM_THREADS.

• 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.

Example: solving a simple harmonic oscillator $m\ddot{x}+\gamma \dot{x}+{\omega }_0^{2}x=F\cos (\omega t)$ to obtain the response as a function of $\omega$

# having obtained a Problem object, let's find steady states
julia> range = ParameterRange(ω => LinRange(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 create multi-dimensional solutions plots.

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

# 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

# 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{Num, Vector{Union{Float64, ComplexF64}}}: Values of all parameters for all solutions.

• fixed_parameters::OrderedCollections.OrderedDict{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::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).

• seed::Union{Nothing, UInt32}: Seed used for the solver

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,
    func::Union{Function, String},
    name::String;
    physical
) -> Any

Creates a solution class in res using the inequality condition (parsed into Symbolics.jl input).

The new class is labelled with name and stored under res.classes[name].

By default, only physical (=real) solutions are classified, false is returned for the rest.

Example

# 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

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