Hello I’m trying to find fixed points of 2 coupled oscillator and one of the parameter is defined as sqrt (function of phase space variables)
here’s the exact system,
using DynamicalSystems, Plots, ChaosTools
function rotor_1( u, params, t)
@inbounds begin
θ1,pθ1,θ2,pθ2 = u
k,l = params
r = sqrt( 2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1) )
h = exp(-k*r) * (1+k*r) / r^3
du1 = pθ1
du2 = h * ( sin(θ1-θ2) + l*sin(θ1) )
du3 = pθ2;
du4 = -h * ( sin(θ1-θ2) + l*sin(θ2) )
end
return SVector{4}(du1,du2,du3,du4)
end
u0 = [2.4,2.7,0.,5.72]
rotor_system = CoupledODEs(rotor_1, u0, [2. , 1.8] )
v = interval(-pi, pi)
box = IntervalBox(v, v, v,v)
fp = fixedpoints(rotor_system, box)
This last fixedpoints was from ChaosTools and was giving
ERROR: DomainError with -1.3322676295501878e-15:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
...
...
...
would not be violated in the solution due to adaptivity. For example, an ODE solver or optimization
routine may check a step at `new_u` which violates the domain constraint, and if violated reject the
step and use a smaller `dt`. However, the throwing of this error will have halted the solving process.
Thus the recommended fix is to replace this function with the equivalent ones from NaNMath.jl
(https://github.com/JuliaMath/NaNMath.jl) which returns a NaN instead of an error. The solver will then
so I replaced sqrt with NaNmath.sqrt in the rotor_1 function, ie
r = NaNMath.sqrt( 2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1) )
but this was giving
julia> fp = fixedpoints(rotor_system, box)
ERROR: StackOverflowError:
Stacktrace:
[1] sqrt(x::IntervalArithmetic.Interval{Float64}) (repeats 79984 times)
@ NaNMath ~/.julia/packages/NaNMath/ceWIc/src/NaNMath.jl:18
Unsure how to proceed further too,
I have also tried with simple power of (1/2)
r = ( 2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1) )^(1/2)
last one took more than half an hour without any results.
My Motive is to find fixed points and the Eigen Value of their jacobian at that those points. From looking at the defined system, the fixed points doesn’t depend on the parameter h exp(-k*r) * (1+k*r) / r^3
since both k, r are positive in my system this factor can never go to zero so I redefined my
rotor_1 as
# r = sqrt( 2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1) )
# h = exp(-k*r) * (1+k*r) / r^3
du1 = pθ1
du2 = ( sin(θ1-θ2) + l*sin(θ1) ) # no h factor here
du3 = pθ2;
du4 = ( sin(θ1-θ2) + l*sin(θ2) ) # same no h factor her
from this I was able to find fixed points of the system but still, to find eigenvalues I need to include the h parameter and I’m not sure to which package to use to find eigenvalues of jacobian.
I know this seems to be very niche problem, but I’m feeling lost Any help or pointers are much appreciated.