JuliaDynamics / ChaosTools.jl

Tools for the exploration of chaos and nonlinear dynamics
https://juliadynamics.github.io/DynamicalSystemsDocs.jl/chaostools/stable/
MIT License
186 stars 35 forks source link

Improve UPO storage in `periodicorbits.jl` #325

Closed JonasKoziorek closed 4 months ago

JonasKoziorek commented 4 months ago

I noticed that while using periodicorbits algorithm with parameter roundtol that is high (eg. 8 ), it detects many duplicate UPOs. MVE:

using ChaosTools
using PredefinedDynamicalSystems: henon
begin
    ds = henon()
    xs = LinRange(-3.0, 3.0, 60)
    ys = LinRange(-10.0, 10.0, 60)
    seeds = [SVector{2}(x,y) for x in xs for y in ys]
    o = 5
    indss, signss = lambdaperms(dimension(ds))
    λ = 0.005
    results1 = sort!(periodicorbits(ds, o, seeds, λ, indss, signss; roundtol = 8))
end

Output:

10-element Vector{SVector{2, Float64}}:
 [-1.13135448, -0.33940634]
 [-1.13135448, -0.33940633]
 [-1.13135447, -0.33940636]
 [-1.13135447, -0.33940635]
 [0.63135447, 0.18940633]
 [0.63135447, 0.18940634]
 [0.63135448, 0.18940634]
 [0.63135448, 0.18940635]
 [0.63135449, 0.18940636]
 [0.63135449, 0.18940636]

My proposed change solves this:

using ChaosTools
using PredefinedDynamicalSystems: henon
begin
    ds = henon()
    xs = LinRange(-3.0, 3.0, 60)
    ys = LinRange(-10.0, 10.0, 60)
    seeds = [SVector{2}(x,y) for x in xs for y in ys]
    o = 5
    indss, signss = lambdaperms(dimension(ds))
    λ = 0.005
    results2 = sort!(periodicorbits(ds, o, seeds, λ, indss, signss; spacetol = 1e-7))
end

Output:

julia> results2
2-element Vector{SVector{2, Float64}}:
 [-1.1313544826122608, -0.33940632522716674]
 [0.6313544679502884, 0.18940632554538597]

I have implemented FP as a binary tree. This way, new fixed points are inserted into FP only if they are spacetol away from the remaining fixed points in FP. Thus duplicates are avoided for a good choice of spacetol.

Additionally, binary tree removes the need for linear search through FP even though this doesn't seem to be a bottleneck.

Let me know if there is something to fix. This is my first PR :).

JonasKoziorek commented 4 months ago

I renamed it as VectorWithEpsRadius since it's used with SVectors.

Should the changelog include info about both PRs? I will make another PR this week adding the Davidchack-Lai algorithm. It will internally also make use of VectorWithEpsRadius. Perhaps all three could be released together.

Datseris commented 4 months ago

we can add the changelog in your 3rd PR, it's fine.