DynamicalSystems Poincare map of sin(u[1]) and u[2] [using ProjectedDynamicalSystems]

Specific Domains Modelling & Simulations
, , , ,
1

I will take an example from a predefined dynamical system

using  PredefinedDynamicalSystems, DynamicalSystems
ds = Systems.double_pendulum( [π/2, 0, 0, π]; G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
Points, t = trajectory(ds, 50000; Ttr = 0.0, Δt =0.01)

# Points array is of this form [θ₁, ω₁, θ₂, ω₂]
# and if I need to have poincare section of plane θ₁=0

psos = poincaresos(ds, (1,0.0), Tmax = 50000.0; u0 = [π/2,0.,0.,π])
lines(psos[:,3],psos[:,4]) # plot of poincare map θ₂, vs ω₂

for a system I’m studying once θ₁ goes to negative it just keeps on going negative, so I need
to calculate poincare section of [sin(θ₁), ω₁, sin(θ₂), ω₂] ie whenever the plane sin(θ₁)=0 I need
a plot of sin(θ₂) vs ω₂

looking at the ProjectedDynamicalSystems Overarching tutorial · DynamicalSystems.jl there seem to be a way to project these and then use poincaresos but I couldn’t figure what’s complete_state and such

projection1(u) = [sin(u[1]),sin(u[3])]
# don't know what to give for complete_state
# pds = ProjectedDynamicalSystem(rotor_system_2, projection1, complete_state_1

any help or pointers are much appreciated

2

Hi, I just moved this post and a previous question of yours to the Modelling & Simulation category, as I hope you might get more relevant answers from the subscribers to that category.

1 Like
3

should I define a differential equation’s callback in the main system such that it u[1] is bracketed
in [-pi, pi] range?

so yes I tried call back function defined GitHub - awage/Basins.jl on this example

this is code for replication

function rotor_1(du, u, params, t)

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

                                          r  = ((2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1))^1/2)
                                          h  = exp(-k*r) * (1+k*r) / r^3;

                                          du[1]  = pθ1;
                                          du[2]  = (h*(sin(θ1-θ2)+l*sin(θ1)))
                                          du[3]  = pθ2;
                                          du[4]  = (h*(sin(θ2-θ1)-l*sin(θ2)))
            end

u0 = [-pi,0.,pi,0.]
params = [2.0,1.5]

# now use callback
# We have to define a callback to wrap the phase in [-π,π]
function affect!(integrator)
    if integrator.u[1] < 0
        integrator.u[1] += 2*π
    else
        integrator.u[1] -= 2*π
    end
end

condition(u,t,integrator) = (integrator.u[1] <= -π  || integrator.u[1] >= π)
cb = DiscreteCallback(condition,affect!)

diffeq = (alg = Vern7(), reltol = 1.0e-10, abstol = 1.0e-10, atol = 1.0e-10, rtol = 1.0e-10, callback=cb)

# coupled_ode
rotor_system = CoupledODEs(rotor_1, u0, params; diffeq,t0=0.0)

now let’s check the trajectory generated now

Points, t = trajectory(rotor_system, 50000.0; Ttr = 0.0, Δt =0.01)

 maximum(Points[:,1]) 
# 3.4152245610206577

 minimum(Points[:,1])
# -3.4250833268013063

 maximum(Points[:,3])
# 3.1415926536003895

 minimum(Points[:,3])
# -3.1270534882495484

so callback works but for θ1 I’m missing the bracket of limits (help needed) and now
energy checking

function hamil_cal(p) 
        θ1,pθ1,θ2,pθ2 = p
        k, l = params
        r   = (sqrt(abs(2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1))))
        H   =  ((pθ1^2+pθ2^2)/2) + (exp(-k*r)/r)
        return H
end

energy = Vector{Float64}(undef,0)
energy_abs = Vector{Float64}(undef,0)

function energy_check(Points)
       global energy = Vector{Float64}(undef,0)
       global energy_abs = Vector{Float64}(undef,0)
       @showprogress for u in Points
             push!(energy,hamil_cal(u))
             push!(energy_abs,abs(hamil_cal(u)-hamil_cal(u0)))
             end
       @show maximum(energy_abs)
       end

energy_check(Points)
# maximum(energy_abs) = 4.0548661650596785e-10

it’s seems to follow the system’s trend so reinitialize and find the poincare map

rotor_system = CoupledODEs(rotor_1, u0, params; diffeq,t0=0.0)
psos = poincaresos(rotor_system, (1,0.0), Tmax = 500000.0; u0=u0)

this gives me

julia> psos = poincaresos(rotor_system, (1,0.0), Tmax = 500000.0; u0=u0)
4-dimensional StateSpaceSet{Float64} with 1 points
 NaN  NaN  NaN  NaN

I have looked into poincaresos manpage Orbit Diagrams & PSOS · DynamicalSystems.jl
from what I see hyperplane is defined in this for plotting section of θ2,pθ2 whenever θ1 becomes 0.0

ofcourse If I look at the Points[:,1] it’s indeed crossing 0 so maybe poincaresos is not using callback defined, I don’t know how to proceed further. Is there any other way or package where I can do this?