# DynamicalSystems Poincare map of sin(u\[1\]) and u\[2\] \[using ProjectedDynamicalSystems\]
**URL:** https://discourse.julialang.org/t/dynamicalsystems-poincare-map-of-sin-u-1-and-u-2-using-projecteddynamicalsystems/106443
**Category:** Modelling & Simulations
**Tags:** question, plotting, makie, differentialequation, dynamical-systems
**Created:** [November 19, 2023, 9:00pm UTC](https://discourse.julialang.org/t/dynamicalsystems-poincare-map-of-sin-u-1-and-u-2-using-projecteddynamicalsystems/106443 "2023-11-19T21:00:20Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 19, 2023, 9:00pm UTC](https://discourse.julialang.org/t/dynamicalsystems-poincare-map-of-sin-u-1-and-u-2-using-projecteddynamicalsystems/106443/1 "2023-11-19T21:00:20Z")
I will take an example from a predefined dynamical system
```julia
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](https://juliadynamics.github.io/DynamicalSystems.jl/dev/tutorial/#DynamicalSystemsBase.ProjectedDynamicalSystem) there seem to be a way to project these and then use poincaresos but I couldn’t figure what’s complete\_state and such
```julia
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
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [November 19, 2023, 11:02pm UTC](https://discourse.julialang.org/t/dynamicalsystems-poincare-map-of-sin-u-1-and-u-2-using-projecteddynamicalsystems/106443/2 "2023-11-19T23:02:45Z")
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.
---
### Author:  [@hisacro](https://discourse.julialang.org/u/hisacro)
#### Post date: [November 20, 2023, 2:40am UTC](https://discourse.julialang.org/t/dynamicalsystems-poincare-map-of-sin-u-1-and-u-2-using-projecteddynamicalsystems/106443/3 "2023-11-20T02:40:38Z")
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](https://github.com/awage/Basins.jl/blob/master/#user-content-5---computation-of-the-saddle-embedded-in-the-boundary-experimental-feature) on this example
this is code for replication
```julia
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
```julia
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
```julia
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
```julia
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
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](https://juliadynamics.github.io/DynamicalSystems.jl/v2.3/chaos/orbitdiagram/#Poincar%C3%A9-Surface-of-Sectionand)
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?