Error for harmonics with decimal prefactors

[curio-sitas]

Hi,

Some frequency driven systems exhibits a "Period-Doubling" behavior, where a response exists at half the excitation frequency. Would it be possible to search for half frequency reponses ?

I tried :

add_harmonic!(system_eqs, I, [0.5ω,ω])

But it ends up with the following error :

ERROR: MethodError: no method matching isless(::ComplexF64, ::Int64)
Closest candidates are:
  isless(::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at C:\Users\brsinquin\.julia\packages\StatsBase\XgjIN\src\statmodels.jl:90
  isless(::ArfLike, ::Integer) at C:\Users\brsinquin\.julia\packages\Arblib\PdCpA\src\predicates.jl:131
  isless(::Static.StaticInteger{X}, ::Real) where X at C:\Users\brsinquin\.julia\packages\Static\Ldb7F\src\Static.jl:455
  ...

when executing : harmonic_eqs = get_harmonic_equations(system_eqs)

[jkosata]

Hello! We've definitely calculated subharmonics before (from the harmonic balance point of view, there is no "subharmonic", simply two harmonics). It works if you use ω/2 instead of 0.5ω . Better to use fractions than decimal prefactors, as you are relying on the harmonics being multiples of each other! That said, this needs a better error message :)

[curio-sitas]

Hello! We've definitely calculated subharmonics before (from the harmonic balance point of view, there is no "subharmonic", simply two harmonics). It works if you use ω/2 instead of 0.5ω . Better to use fractions than decimal prefactors, as you are relying on the harmonics being multiples of each other! That said, this needs a better error message :)

Thank you for your response, it worked. Since i want to study the harmonics of the half-system excitation frequency i need to add non integer fractionnal frequency Anzats such as : add_harmonic!(laser_eqs, I, [ω/2,ω,ω*(3//2)])

The parsing works with ω/2 and also with ω*(3//2) ! ( Just noting that it is not straight forward to find the way to write the fractions for the parser to work.)

But i end up with this error when calculating the steady states :

Computing mixed cells... 1835    Time: 0:00:02
  mixed_volume:  11169
ERROR: InexactError: trunc(Int32, -45513548967005834)
Stacktrace:
  [1] throw_inexacterror(f::Symbol, #unused#::Type{Int32}, val::Int64)
    @ Core .\boot.jl:614
  [2] checked_trunc_sint
    @ .\boot.jl:636 [inlined]
  [3] toInt32
    @ .\boot.jl:673 [inlined]
  [4] Int32
    @ .\boot.jl:763 [inlined]
  [5] convert
    @ .\number.jl:7 [inlined]
  [6] cconvert
    @ .\essentials.jl:412 [inlined]
  [7] set_si!
    @ .\gmp.jl:203 [inlined]
  [8] _broadcast_getindex_evalf
    @ .\broadcast.jl:670 [inlined]
  [9] _broadcast_getindex
    @ .\broadcast.jl:643 [inlined]
 [10] getindex
    @ .\broadcast.jl:597 [inlined]
 [11] macro expansion
    @ .\broadcast.jl:961 [inlined]
 [12] macro expansion
    @ .\simdloop.jl:77 [inlined]
 [13] copyto!
    @ .\broadcast.jl:960 [inlined]
 [14] copyto!
    @ .\broadcast.jl:913 [inlined]
 [15] copy
    @ .\broadcast.jl:885 [inlined]
 [16] materialize
    @ .\broadcast.jl:860 [inlined]
 [17] solve!(BSS::HomotopyContinuation.BinomialSystemSolver)
    @ HomotopyContinuation C:\Users\brsinquin\.julia\packages\HomotopyContinuation\I1faM\src\binomial_system.jl:179
 [18] solve
    @ C:\Users\brsinquin\.julia\packages\HomotopyContinuation\I1faM\src\binomial_system.jl:163 [inlined]
 [19] iterate(iter::HomotopyContinuation.PolyhedralStartSolutionsIterator, ::Tuple{Int64, Int64})
    @ HomotopyContinuation C:\Users\brsinquin\.julia\packages\HomotopyContinuation\I1faM\src\polyhedral.jl:81
 [20] _collect(cont::UnitRange{Int64}, itr::HomotopyContinuation.PolyhedralStartSolutionsIterator, #unused#::Base.HasEltype, isz::Base.SizeUnknown)
    @ Base .\array.jl:725
 [21] collect
    @ .\array.jl:712 [inlined]
 [22] #solve#279
    @ C:\Users\brsinquin\.julia\packages\HomotopyContinuation\I1faM\src\solve.jl:501 [inlined]
 [23] solve(args::HomotopyContinuation.ModelKit.System; show_progress::Bool, threading::Bool, catch_interrupt::Bool, target_parameters::Vector{ComplexF64}, stop_early_cb::Function, transform_result::Nothing, transform_parameters::typeof(identity), flatten::Nothing, target_subspaces::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ HomotopyContinuation C:\Users\brsinquin\.julia\packages\HomotopyContinuation\I1faM\src\solve.jl:487
 [24] _get_raw_solution(problem::Problem, parameter_values::Vector{Vector{Float64}}; sweep::OrderedCollections.OrderedDict{Symbolics.Num, Vector{Float64}}, random_warmup::Bool, threading::Bool, show_progress::Bool)
    @ HarmonicBalance C:\Users\brsinquin\.julia\packages\HarmonicBalance\iRw05\src\solve_homotopy.jl:249
 [25] get_steady_states(prob::Problem, swept_parameters::OrderedCollections.OrderedDict{Symbolics.Num, Vector{Float64}}, fixed_parameters::OrderedCollections.OrderedDict{Symbolics.Num, Float64}; random_warmup::Bool, threading::Bool, show_progress::Bool, sorting::String, classify_default::Bool)
    @ HarmonicBalance C:\Users\brsinquin\.julia\packages\HarmonicBalance\iRw05\src\solve_homotopy.jl:101
 [26] get_steady_states
    @ C:\Users\brsinquin\.julia\packages\HarmonicBalance\iRw05\src\solve_homotopy.jl:89 [inlined]
 [27] #get_steady_states#88
    @ C:\Users\brsinquin\.julia\packages\HarmonicBalance\iRw05\src\solve_homotopy.jl:130 [inlined]
 [28] get_steady_states
    @ C:\Users\brsinquin\.julia\packages\HarmonicBalance\iRw05\src\solve_homotopy.jl:130 [inlined]
 [29] #get_steady_states#89
    @ C:\Users\brsinquin\.julia\packages\HarmonicBalance\iRw05\src\solve_homotopy.jl:131 [inlined]
 [30] get_steady_states(eom::HarmonicEquation, swept::Pair{Symbolics.Num, LinRange{Float64, Int64}}, fixed::Tuple{Pair{Symbolics.Num, Float64}, Pair{Symbolics.Num, Int64}, Pair{Symbolics.Num, Float64}})
    @ HarmonicBalance C:\Users\brsinquin\.julia\packages\HarmonicBalance\iRw05\src\solve_homotopy.jl:131
 [31] top-level scope
    @ d:\Home\brsinquin\Desktop\test\test.jl:16

Here is my whole code :

using HarmonicBalance, Plots, DifferentialEquations
@variables r,m, ε, t, ω, I(t), n(t) # declare constant variables and a function x(t)

laser_eqs = DifferentialEquation([d(I,t,1) - n*I, d(n,t,1) - ε*(r*(1+m*cos(ω*t))-(n+1)*(I+1))] .~ [0,0], [I,n])
add_harmonic!(laser_eqs, I, [ω/2,ω,ω*(3//2), 2ω]) # specify the ansatz x = u(T) cos(ωt) + v(T) sin(ωt)
add_harmonic!(laser_eqs, n, [ω/2,ω,ω*(3//2), 2ω]) # specify the ansatz x = u(T) cos(ωt) + v(T) sin(ωt)
harmonic_eqs = get_harmonic_equations(laser_eqs)

fixed = (ε => 0.027, r=>2, m=>0.5)   # fixed parameters
varied = ω => LinRange(0.1, 0.9, 100)           # range of parameter values
result = get_steady_states(harmonic_eqs, varied, fixed)

plot(result, "sqrt(u2^2 + v2^2)")

[jkosata]

This looks like a problem with the solver from HomotopyContinuation.jl that occurs when solving huge systems. I have only seen it when the number of paths ran above 10⁸ or so.

Here you're looking at 16 equations of order <= 2 , which should give a max of 2^16 = 65536 solutions. By internal symmetries this reduces to 11169 and runs ok for me:

Funnily enough, your output seems like the solver is also trying to run 11169 paths ... Do other, similar-sized problems (e.g. from https://github.com/NonlinearOscillations/HarmonicBalance-notebooks ) run for you? Do you by any chance have a 32-bit OS?