Handling Dynamical system with trigonometric function

Specific Domains Modelling & Simulations
, ,
1 
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

2

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
3

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
4

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

5

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
6

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
7

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.

8

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
9

OK then, you’re welcome!

10

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

1 Like
11

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
12

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
13

oh thank you! will look into this further

1 Like
14

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

2 Likes
15

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

# 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> 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