# Poincare section for Henon Helies system
**URL:** https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059
**Category:** New to Julia
**Tags:** question, differentialequation, dynamical-systems
**Created:** [December 21, 2024, 12:15pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059 "2024-12-21T12:15:45Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@nn-heckfy](https://discourse.julialang.org/u/nn-heckfy)
#### Post date: [December 21, 2024, 12:15pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/1 "2024-12-21T12:15:45Z")
Greetings. I’m trying to solve a problem of plotting a Poincare map for the Henon Heiles’s system. But the figure is completely painted over, but should have cavities. Maybe I’m using the Julia commands incorrectly somehow. Thank you in advance.
```julia
using DifferentialEquations, CairoMakie, LaTeXStrings, Symbolics
E=4e-2
T=1e2
tspan=(0,T)
f = Figure()
ax = Axis(f[1, 1],
title = L"H=\frac{1}{2}(p_1^2+p_2^2)+\frac{1}{2}(q_1^2+q_2^2)+q_1^2 q_2 -\frac{1}{3}q_2^3",
subtitle=latexstring("E=$E"),
xlabel = L"q_2",
ylabel = L"p_2",
)
function henon(du,u,p,t)
du[1]=u[3]
du[2]=u[4]
du[3]=-(u[1]+2*u[1]*u[2])
du[4]=-u[2]-u[1]^2+u[2]^2
end
q10=0
for q20 in -0.5:1e-2:0.5
for p20 in -0.4:1e-2:0.4
if (2*E-p20^2-(q10^2+q20^2)-2*(q10^2*q20)+(2/3)*q20^3)>0
u0=[q10,q20,sqrt(2*E-p20^2-(q10^2+q20^2)-2*(q10^2*q20)+(2/3)*q20^3),p20]
prob=ODEProblem(henon,u0,tspan)
sol=solve(prob,RK4())
# scatter!(axe,sol[1,:],(sol[3,:].^2+sol[4,:].^2)./2+(sol[1,:].^2+sol[2,:].^2)./2+(sol[1,:].^2).*sol[2,:]-(sol[2,:].^3)./3)
#Q1=sol[1,:]; Q2=sol[2,:]; P1=sol[3,:]; P2=sol[4,:]
for i in 1:(length(sol[:])-1)
if (sol[1,i]<0) && (sol[1,i+1])>0
scatter!(ax, (sol[2,i]+sol[2,i+1])/2, (sol[4,i]+sol[4,i+1])/2,color=:blue, markersize=4)
end
end
end
end
end
f
```
UPD1: For T=1e4 and q20,p20 variables range -0.5:1e-1:-0.5 it turns out to be something right. But Dataseris’s figure(see below) is much more attractive.
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [December 21, 2024, 3:50pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/2 "2024-12-21T15:50:29Z")
DynamicalSystems.jl can compute the Poincaré map, given the surface/plane of section and the initial condition:
[https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/#DynamicalSystemsBase.PoincareMap](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/#DynamicalSystemsBase.PoincareMap).
Here is presented, as the last example, [https://github.com/JuliaDynamics/DynamicalSystems.jl/blob/main/docs/src/visualizations.md](https://github.com/JuliaDynamics/DynamicalSystems.jl/blob/main/docs/src/visualizations.md) the interactive visualization of the Poincaré map orbits associated to the Hénon-Heiles dynamical system.
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [December 21, 2024, 3:53pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/3 "2024-12-21T15:53:19Z")
Hi there,
disclaimer: i didn’t look at your code to find the bug, but I am posting something that will make your life easier. Haha thankfully @empet beat me to it. I was going to suggest the `PoincareMap` structure. DynamicalSystems.jl wraps over DifferentialEquations.jl to provide dynamical systems analysis tools, such as a dedicated poincare map dynamical system. Here is the code:
```julia
using DynamicalSystems, CairoMakie
using OrdinaryDiffEqVerner: Vern9
E=4e-2
T=1e2
tspan=(0,T)
f = Figure()
ax = Axis(f[1, 1],
title = L"H=\frac{1}{2}(p_1^2+p_2^2)+\frac{1}{2}(q_1^2+q_2^2)+q_1^2 q_2 -\frac{1}{3}q_2^3",
subtitle=L"E=$($E)",
xlabel = L"q_2",
ylabel = L"p_2",
)
function henon(du,u,p,t)
du[1]=u[3]
du[2]=u[4]
du[3]=-(u[1]+2*u[1]*u[2])
du[4]=-u[2]-u[1]^2+u[2]^2
end
# Use integrator with higher tolerance (my experience much better for Hamiltonians)
ds = CoupledODEs(henon, rand(4), nothing; diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9))
pmap = PoincareMap(ds, (1, 0.0)) # hyperplane of 1st variable crossing 0.0
q10=0
for q20 in -0.5:1e-2:0.5
for p20 in -0.4:1e-2:0.4
if (2*E-p20^2-(q10^2+q20^2)-2*(q10^2*q20)+(2/3)*q20^3)>0
u0=[q10,q20,sqrt(2*E-p20^2-(q10^2+q20^2)-2*(q10^2*q20)+(2/3)*q20^3),p20]
reinit!(pmap, u0)
# get the next 10 crossings of the Poincare section
X, tvec = trajectory(pmap, 10)
scatter!(ax, X[:, [2, 4]]; markersize = 5, color = (:black, 0.1))
end
end
end
f
```
which produces
as you can see, you don’t have to manually interpolate the solution nor find the crossings yourself! 😉
---
### Author:  [@nn-heckfy](https://discourse.julialang.org/u/nn-heckfy)
#### Post date: [December 21, 2024, 6:08pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/4 "2024-12-21T18:08:15Z")
Thank you. But it is necessary to solve this problem without using additional packages.
---
### Author:  [@nn-heckfy](https://discourse.julialang.org/u/nn-heckfy)
#### Post date: [December 21, 2024, 6:09pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/5 "2024-12-21T18:09:30Z")
Thank you very much for you detailed answer and beatiful attached figure. But it is necessary to solve this problem without using additional packages.
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [December 21, 2024, 6:14pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/6 "2024-12-21T18:14:12Z")
Adding differential equations brings in all dependencies that dynamical systems has.
---
### Author:  [@nn-heckfy](https://discourse.julialang.org/u/nn-heckfy)
#### Post date: [December 22, 2024, 7:57am UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/7 "2024-12-22T07:57:20Z")
Sorry I don’t understand what you mean.
---
### Author:  [@Benny](https://discourse.julialang.org/u/Benny)
#### Post date: [December 22, 2024, 8:39am UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/8 "2024-12-22T08:39:55Z")
Datseris meant that DifferentialEquations and DynamicalSystems have an extremely large overlap in used packages, so there’s no practical reason to abstain from the latter for the former. You may have your own reasons to implement a solution from DifferentialEquations alone, so the CoupledODEs/PoincareMap example can just serve as proof that something very similar to your code gets the result (or something closer to it, the math is lost on me). DynamicalSystems documents how those are built upon DifferentialEquations, so that could be worth studying.
---
### Author:  [@nn-heckfy](https://discourse.julialang.org/u/nn-heckfy)
#### Post date: [December 22, 2024, 8:48am UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/9 "2024-12-22T08:48:17Z")
I get it, thank you.
---
### Author:  [@Ondrej\_Zelenka](https://discourse.julialang.org/u/Ondrej_Zelenka)
#### Post date: [November 10, 2025, 3:12pm UTC](https://discourse.julialang.org/t/poincare-section-for-henon-helies-system/124059/10 "2025-11-10T15:12:48Z")
Hello! While I am not a Julia programmer, I noticed this post and would like to give my two cents. Hope this is still relevant a year after the question was posted with no visible resolution.
I have observed something very similar which could potentially be the same issue when insufficient interpolation order was used. The code snippet seems to use 0th order interpolation (sort of), I recommend switching to 1st order interpolation and reducing the step size (I would not go above 3rd order interpolation).
As I said, I am not a Julia programmer, and my code may have a mistake, but please try this:
```julia-auto
scatter!(ax, (sol[1,i+1]*sol[2,i]-sol[1,i]*sol[2,i+1])/(sol[1,i+1]-sol[1,i]), (sol[1,i+1]*sol[4,i]-sol[1,i]*sol[4,i+1])/(sol[1,i+1]-sol[1,i]), color=:blue, markersize=4)
```