JuliaDynamics / DynamicalSystemsBase.jl

Definition of dynamical systems and integrators for DynamicalSystems.jl
Other
53 stars 27 forks source link

poincaremap direction mismatch #175

Closed rusandris closed 11 months ago

rusandris commented 1 year ago

In the documentation of PoincareMap it says:

direction = -1: Only crossings with sign(direction) are considered to belong to the surface of section. Positive direction means going from less than b to greater than b.

I think it should be the other way around, because positive direction right now means going from greater than b to less than b. This should be a small fix in the docstring, unless you really want it to work that way.

I tried this

using DynamicalSystemsBase
using Plots

using OrdinaryDiffEq: Vern9

@inbounds function lorenz_rule!(du, u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
    return nothing
end

u0 = [0, 10.0, 0]
p0 = [10, 28, 8/3]
diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)

ds = CoupledODEs(lorenz_rule!, u0, p0; diffeq)

plane= [1.0,0.0,0.0,5.0] #plane in 3D phase space with x = 5.0

pmap_minus = poincaremap(deepcopy(ds),plane;direction=-1)
pmap_plus = poincaremap(deepcopy(ds),plane;direction=+1)

psec_minus = current_state(pmap_minus)
println("time of negative direction crossing: ",current_crossing_time(pmap_minus))

psec_plus = current_state(pmap_plus)
println("time of positive direction crossing: ",current_crossing_time(pmap_plus))

traj,t = trajectory(ds,current_crossing_time(pmap_plus)+0.1;Ttr = 0,Δt = 1e-3)

plot(traj[:,1],traj[:,2],xlabel="x",ylabel="y",zlabel="z",arrow=:arrow,lc=:gray10,label="traj",lw=2)
plot!([plane[end],plane[end]],[-10,30],ls=:dash,lc=:gray80,lw=3,label="plane")
plot!([psec_minus[1]],[psec_minus[2]],st=:scatter,ms=6,label="-1",mc=:orange)
plot!([psec_plus[1]],[psec_plus[2]],st=:scatter,ms=6,mc=:green,label="+1")
Datseris commented 1 year ago

Can you please paste the output picture?

rusandris commented 1 year ago

Can you please paste the output picture?

poincare_test

Datseris commented 1 year ago

okay yeah I had fucked up what means positive. Please do a PR to change the "positive" to "negative"?

rusandris commented 1 year ago

Sure thing. While we're at it, I think it would be useful to get psection points from both directions in some cases. Do you think that PoincareMap should have an option for that? What about poincaresos returning the timepoints of all sections (also in case of datasets)? I've ran into some situations recently where these would've been useful. I can help you add these features when I find some time for it. Let me know what you think.

Datseris commented 1 year ago

Sure; I guess we should alter the direction to allow both as well somehow. Probably using 0 as the value.

rusandris commented 11 months ago

You can close this as well.