# Handling Dynamical system with trigonometric function

using DynamicalSystems, GLMakie, OrdinaryDiffEq
function rotor_1( u, params, t)
#                      @inbounds begin

θ1,pθ1,θ2,pθ2  = u
k,l = params

r  = (sqrt(abs(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(θ2-θ1)-l*sin(θ2)))
#                      end

return SVector{4}(du1,du2,du3,du4)
end

k = 2.  #exponential factor
l = 1.5 #seperation distance

params = [k, l]
u0=[pi,0.0,-pi,0.] #indeed a fixed point of the system

diffeq = (alg = Vern7(), reltol = 1e-10, abstol = 1e-10 , atol=1e-10,rtol=1e-10)
rotor_system =  CoupledODEs(rotor_1, u0, params; diffeq, t0 = 0.0)
Points, t = trajectory(rotor_system, 8000.0; Ttr = 0.0, Δt =0.1)

lines(t,columns(Points)[1]) # contradicts the results after time steps 100


this is my system as it depends on sin(θ2-θ1) and one obvious thing is

julia> sin(pi-0) + sin(0)
1.2246467991473532e-16


so this was causing Fixed points & Periodicity · ChaosTools.jl to not find the exact ones and misinterpreting the results.

I believe this has to do with floating error not the julia or dynamicalsystems.ji itself but how to tackle it?

I have looked into ~/.julia/packages/PredefinedDynamicalSystems/Af98t/src/continuous_famous_systems.jl but couldn’t find one that has 2 angle parameters inside trignometric function.

Any help is much appreciated

Maybe this rewrite of rotor_1 helps (untested, just playing around):

function rotor_1(u, params, t)
θ1, pθ1, θ2, pθ2 = u
k, l = params

(s1, c1) = sincos(θ1)
(s2, c2) = sincos(θ2)

y = fma(c1, c2, fma(s1, s2, 1))
r = (sqrt ∘ abs ∘ fma)(l, l - 2*(c1 - c2),  -2*y)
x = k * r
h = exp(-x) * (1 + x) / r^3

s = s1*c2 - c1*s2

du1 = pθ1
du2 =  h * fma(s1, l, s)
du3 = pθ2
du4 = -h * fma(s2, l, s)

SVector{4}(du1, du2, du3, du4)
end


My goal was to use trig identities to eliminate expressions like sin(a - b) and cos(a - b), to delay catastrophic cancellation. Then it was possible to use fma (fused multiply-add) in many more places, which could very well help with accuracy, too. I think this should also be faster, given the significant degree of deduplication.

1 Like

I have a hunch this might be an improvement:

r = let
a = fma(s1, s2, 1)
b = fma(c1, c2, a)
c = fma(l, c1 - c2, b)
d = fma(l, l, -2*c)
(sqrt ∘ abs)(d)
end

1 Like

is there a way I can see that simplification of fma with keeping variables as it is to check it matches with original equation of motion (something like maxima can do)

idea seems promising

I’m not sure what you mean.

If you’re asking about the semantics of FMA, fma(a, b, c) is defined like this: fl({a \cdot b} + c), where fl denotes a function that rounds a real number into a floating-point number. The point being that there’s only a single rounding, as opposed to the two that would happen if you just did a*b + c instead: fl(fl(a \cdot b) + c).

1 Like

okay this is what I meant, I found Symbolics.ji can do this

> using Symbolics
> @variables θ1,θ2

> fma_check(a,b,c)=a*b+c
fma_check (generic function with 1 method)

> (s1, c1) = sincos(θ1)
(sin(θ1), cos(θ1))

> (s2,c2) = sincos(θ2)
(sin(θ2), cos(θ2))

> fma_check(c1,s2,1)
1 + sin(θ2)*cos(θ1)



will let you know after testing a bit! thanks!

1 Like

What range may each of k, l, θ1, θ2 be drawn from? There’s another idea I’m thinking about. To be specific, I want to know what values can x in my code above take.

it was +2 in the r^2 rather than -2 so I added d+4.0 and this is my new rotor_2 now

function rotor_2(u, params, t)

θ1, pθ1, θ2, pθ2 = u
k, l = params

(s1, c1) = sincos(θ1)
(s2, c2) = sincos(θ2)

# fma(a,b,c) = a*b + c

r = let
a = fma(s1, s2, 1)
b = fma(c1, c2, a)
c = fma(l, c1 - c2, b)
d = fma(l, l, -2*c)
(sqrt ∘ abs)(d+4.0)
end

x = k * r
h = exp(-x) * (1 + x) / r^3

s = s1*c2 - c1*s2

du1 = pθ1
du2 =  h * fma(s1, l, s)
du3 = pθ2
du4 = -h * fma(s2, l, s)

return SVector{4}(du1, du2, du3, du4)
success!
end


success!

> u0
4-element Vector{Float64}:
-3.141592653589793
0.0
3.141592653589793
0.0

> sampling_time = 0.1
0.1

> total_time = 500000.0
500000.0

> transient_time = 0.0
0.0

> Points_2, t = trajectory(rotor_system_2, total_time; Ttr = transient_time, Δt =sampling_time)

4-dimensional StateSpaceSet{Float64} with 5000001 points
-3.14159   0.0           3.14159   0.0
-3.14159  -7.90889e-116  3.14159  -8.27251e-116
-3.14159  -1.58178e-115  3.14159  -1.6545e-115
-3.14159  -2.37267e-115  3.14159  -2.48175e-115
-3.14159  -3.16355e-115  3.14159  -3.30901e-115

> Points_1 # previous system
4-dimensional StateSpaceSet{Float64} with 5000001 points
-3.14159   0.0          3.14159   0.0
-3.14159   2.02673e-19  3.14159  -1.41871e-18
-3.14159   4.05346e-19  3.14159  -2.83743e-18
-3.14159   6.0802e-19   3.14159  -4.25614e-18



and I have checked with energy conservation previously it was evolving with an error of
1.9974651099818885e-10 now it’s 0.0

I will mark it has solved, hearty thanks ; ) and have a great day

1 Like

OK then, you’re welcome!

θ1, θ2 I’m bracketing to -pi to pi can be 0 to 2*pi too

k is a positive definte parameter right now I have fixed it to 2.0

l is also positive definite parameter (separation distance), in my analysis I’m making for 1.5 but the evident chaotic part is l<2.0 and above that it’s showing quasi periodic

1 Like

x = k*r and r can take values from 0 to (l+2) so maximum(x)=(l+2)*k

sorry for the correction

1 Like

I guess you already solved this then, but in case you’d need even better accuracy or performance, a small improvement could maybe be achieved by approximating this function (\exp(-x) \cdot (1 + x)) with a polynomial.

1 Like

oh thank you! will look into this further

1 Like

@nsajko thanks a lot for your support with this question!

2 Likes

actually this still doesn’t solve! the parameter values here were [k=2.5, l = 89] only for this
the results were like this for parameter [k=2.5, l=1.5] both fms defined system and only cos sin are matching!

sorry! this is without my proper testing I came to an conclusion. Damn me!

# earlier
function rotor_1( u, params, t)
#                      @inbounds begin

θ1,pθ1,θ2,pθ2  = u
k,l = params

r  = (sqrt(abs(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(θ2-θ1)-l*sin(θ2)))
#                      end

return SVector{4}(du1,du2,du3,du4)
end

function rotor_2(u, params, t)

θ1, pθ1, θ2, pθ2 = u
k, l = params

(s1, c1) = sincos(θ1)
(s2, c2) = sincos(θ2)

# fma(a,b,c) = a*b + c

r = let
a = fma(s1, s2, 1)
b = fma(c1, c2, a)
c = fma(l, c1 - c2, b)
d = fma(l, l, -2*c)
(sqrt)(d+4.0)
end

x = k * r
h = exp(-x) * (1 + x) / r^3

s = s1*c2 - c1*s2

du1 = pθ1
du2 =  h * fma(s1, l, s)
du3 = pθ2
du4 = -h * fma(s2, l, s)

return SVector{4}(du1, du2, du3, du4)

end

###############################
k = 2.5  #exponential factor
l = 1.5 #seperation distance

total_time = 500000.0
sampling_time = 0.1
transient_time = 0.0
################################

diffeq = (alg = Vern7(), reltol = 1e-10, abstol = 1e-10 , atol=1e-10,rtol=1e-10)

rotor_system_1 =  CoupledODEs(rotor_1, u0, params; diffeq, t0 = 0.0)
rotor_system_2 =  CoupledODEs(rotor_2, u0, params; diffeq, t0 = 0.0)

Points_1, t = trajectory(rotor_system_1, total_time; Ttr = transient_time, Δt =sampling_time)
Points_2, t = trajectory(rotor_system_2, total_time; Ttr = transient_time, Δt =sampling_time)



this is whole comparison code

julia> Points_1
4-dimensional StateSpaceSet{Float64} with 5000001 points
-3.14159   0.0           3.14159    0.0
-3.14159   3.61313e-19   3.14159   -2.52919e-18
-3.14159   7.22626e-19   3.14159   -5.05838e-18
-3.14159   1.08394e-18   3.14159   -7.58757e-18
-3.14159   1.44525e-18   3.14159   -1.01168e-17

julia> Points_2
4-dimensional StateSpaceSet{Float64} with 5000001 points
-3.14159  0.0          3.14159   0.0
-3.14159  3.61313e-19  3.14159  -2.52919e-18
-3.14159  7.22626e-19  3.14159  -5.05838e-18
-3.14159  1.08394e-18  3.14159  -7.58757e-18
-3.14159  1.44525e-18  3.14159  -1.01168e-17

julia> Points_2[end]
4-element SVector{4, Float64} with indices SOneTo(4):
-2.1002090499357218
0.2249507918476913
1.5381130463951052
0.11349819349001435

julia> Points_1[end]
4-element SVector{4, Float64} with indices SOneTo(4):
-4.091961018961626
-0.09646451216587304
-0.5122524725797137
-0.23686290047212855


in the end they are evolving differently, there’s very slight variation but this should be a fixed point of the system, I believe the 2nd and 4th variables are making these to blow up