# \[suggestions\] A Way to find fixed points of coupledODE with a parameter that's sqrt function
**URL:** https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858
**Category:** Modelling & Simulations
**Tags:** ode, dynamical-systems
**Created:** [November 6, 2023, 1:52pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858 "2023-11-06T13:52:46Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 6, 2023, 1:52pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858/1 "2023-11-06T13:52:46Z")
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,
```julia
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
```julia
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
```julia
r = NaNMath.sqrt( 2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1) )
```
but this was giving
```julia
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)
```julia
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
```julia
# 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.
---
### Author:  [@ocp](https://discourse.julialang.org/u/ocp)
#### Post date: [November 6, 2023, 3:06pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858/2 "2023-11-06T15:06:58Z")
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?
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 6, 2023, 4:07pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858/3 "2023-11-06T16:07:03Z")
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`
---
### Author:  [@ocp](https://discourse.julialang.org/u/ocp)
#### Post date: [November 6, 2023, 4:20pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858/4 "2023-11-06T16:20:49Z")
Sorry, I must have transposed my `θ`s.
For `θ2 = -2*atan(1/√19); θ1 = -2*atan(√19)`, r^2 = 0.
```julia
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
```
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 6, 2023, 5:18pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858/5 "2023-11-06T17:18:29Z")
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?
---
### Author:  [@ocp](https://discourse.julialang.org/u/ocp)
#### Post date: [November 6, 2023, 5:51pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858/6 "2023-11-06T17:51:14Z")
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]](https://www.wolframalpha.com/input?i=solve%5B2+%2B+1.8%5E2+-+2*1.8*%28cos%5Bx1%5D+-+cos%5Bx2%5D%29+-+2*cos%28x2-x1%29+%3D%3D+0%2C+x1%2C+x2%5D)
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 6, 2023, 7:16pm UTC](https://discourse.julialang.org/t/suggestions-a-way-to-find-fixed-points-of-coupledode-with-a-parameter-thats-sqrt-function/105858/7 "2023-11-06T19:16:28Z")
actually the solution of
```julia
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
```julia
du2 = ( sin(θ1-θ2) + l*sin(θ1) )
du4 = ( sin(θ1-θ2) + l*sin(θ2) )
```
for the case of l = 1.8
```julia
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 : )