Instability in computing steady states for certain parameter ranges

[oameye]

When plotting the stability diagram of the parametron with the following code,

using HarmonicBalance, ModelingToolkit, Plots

@variables t x(t) # symbolic variables
@parameters ω ω0 F0 Γ λ θ β ψ η
natural_equation = d(d(x, t), t) + Γ*d(x, t) + ω0^2*(1-λ*cos(2*ω*t + ψ))*x + β*x^3 + η*d(x, t)*x^2
forces = F0*cos(ω*t + θ)
diff_eq = DifferentialEquation(natural_equation ~ forces, x)

add_harmonic!(diff_eq, x, ω);
harmonic_eq = get_harmonic_equations(diff_eq);

varied_stability = (ω => LinRange(0.95, 1.05, 100), λ => LinRange(1e-6, 0.06, 100))
fixed_stability = Dict(ω0 => 1.0, Γ => 0.01, F0 => 0.0, β => 1.0, η => 0.3, θ => 0.0, ψ => 0.0)
result_stability = get_steady_states(harmonic_eq, varied_stability, fixed_stability);
plot_stability = plot_phase_diagram(result_stability, class= "stable", title = "η=$(fixed_stability[η])")

we find instability in computing steady states for certain parameter ranges.

This is caused by the fact random_warmup is set to true (probably) by default in get_steady_states. It seems the algorithm is not converging to all solutions, as it is doing the parameter homotopy from a system of polynomial that was adequate for a given parameter, but is not working well when you change the parameters of your system. The way that we implemented a fix is to start tracking the solutions from a trivial polynomial system which we know has the maximum number of solutions possible. This you do setting random_warmup = false in get steadystates.

It doesn't seem to be a problem we can solve, it might be more of an issue for the HomotopyContinuation.jl people.

A fix for now is setting random_warmup = false in get_steady_states function.

[jkosata]

We’ve seen this with coupled parametrons too, seems like neither false nor true is an optimal default setting. Does this noise also appear when plotting class=“all” ? Just to make sure it’s not due to misclassification …

[oameye]

It does seem to come from the Homotopy calculation

[oameye]

The next step will be to make a minimal example in the HomotopyContinuation.jl framework and open an issue at HomotopyContinuation.jl.