Poincare section for Henon Helies system

New to Julia
, ,
1

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.

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.

Poincare_map
Poincare_map1200×900 108 KB

2

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.

Here is presented, as the last example, 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.

1 Like
3

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:

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

image
image779×261 63.1 KB

as you can see, you don’t have to manually interpolate the solution nor find the crossings yourself! :wink:

1 Like
4

Thank you. But it is necessary to solve this problem without using additional packages.

5

Thank you very much for you detailed answer and beatiful attached figure. But it is necessary to solve this problem without using additional packages.

6

Adding differential equations brings in all dependencies that dynamical systems has.

7

Sorry I don’t understand what you mean.

8

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.

1 Like
9

I get it, thank you.