Different result from get_steady_states in case of a 2D sweep depending on sampling rate in the same linear range

[soumyasrtk]

I get different number of solution branches from get_steady_states in the case of a 2D sweep if I change the sampling number (eg, going from LinRange(-0.25, 0.25, 80) to LinRange(-0.25, 0.25, 81) ) even when I sweep the parameter between the same values. The bug does not appear in HarmonicBalance v0.6.4, but it appears in HarmonicBalance v0.8.0

Expected behavior

The number of solution branches should have remained the same in both cases. The error does not happen when I am working with a positive [0,0.25] or a negative [-0.25,0] range, but only when I am working with a negative to positive range [-0.25,0.25].

Minimal Reproducible Example

using HarmonicBalance
@variables ω0,γ,ω,α,F₁,m,δ, t,x(t), T;
natural_equation =  d(x,t,2) +γ*d(x,t) + ω0^2*x + α*x^3;
factor = γ^2/4;
factor /= γ^2/4 + δ^2;
forces_new =  F₁*cos(ω*t) - factor*F₁*m*cos(ω*t);
diffEOM_new = DifferentialEquation(natural_equation + forces_new, x);
add_harmonic!(diffEOM_new, x, [ω]);
averagediffEOM_new = get_harmonic_equations(diffEOM_new);
fixed = (F₁=>0.02,α => 1.0,γ=>0.025,ω0=>1.,ω=>1.1)
swept = (m=> LinRange(0.0,1.0,50), δ => LinRange(-0.25, 0.25, 80)); #here, try LinRange(-0.25, 0.25, 81)
res2D = get_steady_states(averagediffEOM_new, swept, fixed);
plot_phase_diagram(res2D)

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

HarmonicBalance v0.8.0

• Output of versioninfo()

Julia Version 1.9.3

[oameye]

Hey Soumya. Should the last term in natural_equation be cubed?

[soumyasrtk]

Oh hey, sorry about that. Must've left it out while copying over here. Yes, it should be cubed.
Edited and corrected it.

[oameye]

Hey Soumya,

The problem seems to be solved if you use the total_degree solver in get_steady_states by adding the kwarg method=:total_degree. This is because by default we use method=:warmup, where we first try to find all branches with a simpler warmup exercise and then do the complete homotopy. However, if some roots slip through the cracks in the warmup, they will also not appear in the final solutions. In general, I recommend doing a :total_degree (use some threads to make it faster). We probably should make :total_degree default and recommend :warmup when performance is important for, for example, larger resolution 2d sweeps.

using HarmonicBalance

@variables ω0, γ, ω, α, F₁, m, δ, t, x(t);
natural_equation = d(x, t, 2) + γ * d(x, t) + ω0^2 * x + α * x^3;
factor = γ^2 / 4;
factor /= γ^2 / 4 + δ^2;
forces_new = F₁ * cos(ω * t) - factor * F₁ * m * cos(ω * t);
diffEOM_new = DifferentialEquation(natural_equation + forces_new, x);
add_harmonic!(diffEOM_new, x, [ω]);
averagediffEOM_new = get_harmonic_equations(diffEOM_new);
fixed = (F₁ => 0.02, α => 1.0, γ => 0.025, ω0 => 1.0, ω => 1.1)
swept = (m => range(0.0, 1.0, 50), δ => range(-0.25, 0.25, 80)); #here, try LinRange(-0.25, 0.25, 81)

res2D = get_steady_states(averagediffEOM_new, swept, fixed, method=:total_degree)
plot_phase_diagram(res2D)

[oameye]

I am confused why it did work on v1.6.4 on your side. Because, I had the same "bug" on v1.6.4. It would also surprise me, as that part of the code did note change at all.

It could be this line:
https://github.com/NonlinearOscillations/HarmonicBalance.jl/blob/25fd0ae1657f4af9f61c425fe5497192ba5569e0/src/solve_homotopy.jl#L256C3-L256C74

 complex_pert = [1E-2 * issubset(p, keys(sweep))*randn(ComplexF64) for p in problem.parameters]

I will look into it

[oameye]

I cannot seem to replicate it working on 1.6.4