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

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),sin(u)]
# 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

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

should I define a differential equation’s callback in the main system such that it u 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  = pθ1;
du  = (h*(sin(θ1-θ2)+l*sin(θ1)))
du  = pθ2;
du  = (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 < 0
integrator.u += 2*π
else
integrator.u -= 2*π
end
end

condition(u,t,integrator) = (integrator.u <= -π  || integrator.u >= π)
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?