# Handling Dynamical system with trigonometric function
**URL:** https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429
**Category:** Modelling & Simulations
**Tags:** trigonometry, differentialequation, dynamical-systems
**Created:** [November 19, 2023, 5:44am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429 "2023-11-19T05:44:16Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 5:44am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/1 "2023-11-19T05:44:16Z")
```julia
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
julia> sin(pi-0) + sin(0)
1.2246467991473532e-16
```
so this was causing [Fixed points & Periodicity · ChaosTools.jl](https://juliadynamics.github.io/ChaosTools.jl/dev/periodicity/#Fixed-points) 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
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 19, 2023, 7:49am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/2 "2023-11-19T07:49:33Z")
> [@hisacro](#):
>
> `function rotor_1( u, params, t)`
Maybe this rewrite of `rotor_1` helps (untested, just playing around):
```julia
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.
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 19, 2023, 8:07am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/3 "2023-11-19T08:07:51Z")
> [@nsajko](#):
>
> ```julia-auto
> y = fma(c1, c2, fma(s1, s2, 1))
> r = (sqrt ∘ abs ∘ fma)(l, l - 2*(c1 - c2), -2*y)
>
> ```
I have a hunch this might be an improvement:
```julia
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
```
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 8:13am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/4 "2023-11-19T08:13:11Z")
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
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 19, 2023, 8:16am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/5 "2023-11-19T08:16:10Z")
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).
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 8:32am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/6 "2023-11-19T08:32:18Z")
okay this is what I meant, I found Symbolics.ji can do this
```julia
> 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!
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 19, 2023, 9:26am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/7 "2023-11-19T09:26:08Z")
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.
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 9:33am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/8 "2023-11-19T09:33:37Z")
it was +2 in the `r^2` rather than -2 so I added d+4.0 and this is my new rotor\_2 now
```julia
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!~~
```julia
> 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
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 19, 2023, 9:35am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/9 "2023-11-19T09:35:47Z")
OK then, you’re welcome!
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 9:40am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/10 "2023-11-19T09:40:10Z")
> [@nsajko](#):
>
> of `k`, `l`, `θ1`, `θ2` be drawn from?
just saw this reply sorry,
`θ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
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 9:44am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/11 "2023-11-19T09:44:22Z")
> [@nsajko](#):
>
> To be specific, I want to know what values can `x` in my code above take.
`x = k*r` and `r` can take values from `0 to (l+2)` so `maximum(x)=(l+2)*k`
sorry for the correction
---
### Author:  [@nsajko](https://discourse.julialang.org/u/nsajko)
#### Post date: [November 19, 2023, 9:47am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/12 "2023-11-19T09:47:46Z")
> [@hisacro](#):
>
> `exp(-x) * (1 + x)`
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.
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 9:51am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/13 "2023-11-19T09:51:16Z")
oh thank you! will look into this further
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [November 19, 2023, 9:56am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/14 "2023-11-19T09:56:24Z")
@nsajko thanks a lot for your support with this question!
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 10:45am UTC](https://discourse.julialang.org/t/handling-dynamical-system-with-trigonometric-function/106429/15 "2023-11-19T10:45:30Z")
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!
```julia
# 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
# newer
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
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