# [suggestions] A Way to find fixed points of coupledODE with a parameter that's sqrt function

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.

I don’t have a direct answer, but it does appears that there \theta_1 and \theta_2 for which r=0, for example θ1 = -2*atan(1/√19) and θ2 = -2*atan(√19). At those points, the floating point is probably getting in the way and giving really small negative numbers.

However, it seems that that since h is Inf when r==0 there might be another underlying problem with your equations of motion?

1 Like

but r is defined as 2 + (l^2) + 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1) )

so there’s non vanishing term l^2 with the points θ1 = -2*atan(1/√19) , θ2 = -2*atan(√19) and l is the parameter I’m varying to analyze stability.

Actually r in the system is magnitude of relative distance between 2 rotors, it’s positive definite and is zero only in the case of point particles ie l=0
but in my case l>0.5

Sorry, I must have transposed my θs.

For θ2 = -2*atan(1/√19); θ1 = -2*atan(√19), r^2 = 0.

julia> l = 1.8
julia> r2(θ1, θ2) = 2 + (l^2) + 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1)
julia> θ2 = -2*atan(1/√19); θ1 = -2*atan(√19)
julia> r2(θ1, θ2)
-2.220446049250313e-16

julia> r2(θ2, θ1) # == (2*l)^2 = 12.96
12.96

1 Like

totally missed it, I can see now! this is causing r to be complex number.
mind asking what transcendental solver do you use to find the that value?

No problem. I used WolframAlpha to find the roots: solve[2 + 1.8^2 - 2 * 1.8*(cos[x1] - cos[x2]) - 2*cos(x2-x1) == 0, x1, x2]

1 Like

actually the solution of

r  = ( 2 + l^2 -  2*l*(cos(θ1)- cos(θ2))  - 2*cos(θ2-θ1) )^(1/2)


for l < = 2.0 is also fixed points! of differential equations

du2  =  ( sin(θ1-θ2) + l*sin(θ1) )
du4  =  ( sin(θ1-θ2) + l*sin(θ2) )


for the case of l = 1.8


julia> l = 1.8
1.8

julia> # solving for fixed points
function trans(z)
θ1,θ2 = z

SVector(sin(θ1-θ2)+l*sin(θ1), sin(θ2-θ1)-l*sin(θ2));

end
trans (generic function with 1 method)

julia> rts = roots(trans, IntervalBox(-pi..pi, 2))
11-element Vector{Root{IntervalBox{2, Float64}}}:
Root([-0.451027, -0.451026] × [-2.69057, -2.69056], :unique)
Root([-3.1416, -3.14159] × [3.14159, 3.1416], :unknown)
Root([3.14159, 3.1416] × [-3.1416, -3.14159], :unknown)
Root([-3.1416, -3.14159] × [-5.12118e-08, 3.42732e-08], :unknown)
Root([-2.19864e-09, 3.58359e-09] × [-4.85786e-09, 8.04318e-09], :unique)
Root([0.451026, 0.451027] × [2.69056, 2.69057], :unique)
Root([-3.38977e-08, 8.74439e-09] × [3.14159, 3.1416], :unknown)
Root([3.14159, 3.1416] × [-5.83401e-08, 2.12105e-08], :unknown)
Root([-3.1416, -3.14159] × [-3.1416, -3.14159], :unknown)
Root([3.14159, 3.1416] × [3.14159, 3.1416], :unknown)
Root([-3.85559e-08, 5.74371e-08] × [-3.1416, -3.14159], :unknown)

julia> rts[1]
Root([-0.451027, -0.451026] × [-2.69057, -2.69056], :unique)

julia>  -2*atan(1/√19), -2*atan(√19) # solution to r = 0
(-0.4510268117962624, -2.6905658417935308)

julia> rts[6]
Root([0.451026, 0.451027] × [2.69056, 2.69057], :unique)

julia> +2*atan(1/√19), +2*atan(√19)  # solution to r = 0
(0.4510268117962624, 2.6905658417935308)


so this was why the eigen values were blowing up even with NaNMath.sqrt

Thanks for pointing it out, have a great day : )