Solving harmonic equations

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=false) # ignore the Jacobian for now
Problem(eom::HarmonicEquation; Jacobian::Matrix) # use the given matrix as the Jacobian

get_steady_states(prob::Problem, 
                    swept_parameters::ParameterRange,
                    fixed_parameters::ParameterList;
                    random_warmup=false,
                        threading=false,
                        sorting="hilbert")

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

• random_warmup: If true, 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 sweptparameters is then much faster than the initial solving. If `randomwarmup=false`, 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".

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 $ω$

# 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

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{Float64}}

  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::Any

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

classify_solutions!(res::Result, condition::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")

sort_solutions(solutions::Array; sorting) -> 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.