JuliaDynamics / ChaosTools.jl

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

Problems with `basins_of_attraction` for in-place continuous systems #188

Closed Datseris closed 3 years ago

Datseris commented 3 years ago

Minimal Working Example

    ds = Systems.lorenz_iip()
    xg = yg = range(-5.,5.; length=5)
    zg = range(0,10; length=5)
    @time basin, attractors = basins_of_attraction((xg,yg,zg), ds)

I've been running this for 5 minutes now and it still hasn't terminated, even though I got a coarse grid. THe Lorenz63 model for default parameters also has only a single attractor, so this shouldn't be hard to find.

We may have a performance problem and/or a bug somewhere for continuous systems in-place methods.

awage commented 3 years ago

I will have a look next week on thin one.

Datseris commented 3 years ago

I'm thinking that we made a mistake with the PR you made for discrete in-place systems, where you modified reinit! to use copy. DiffEq doesn't do this. Instead, we should have modified our own CompleteAndReinit to do a copy. Will do so now.

EDIT: ok no you were correct, the copy in DynamicalSystemsBase.jl was definitely necessary,

Datseris commented 3 years ago

Somewhere we are doing something really wrong. I've added copy in our CompleteAndReinit just for safety. Then, with the new progress bar, I've run the above example again. It got stuck here forever: image

awage commented 3 years ago

Seems that the bugfix in #190 also solved this problem. It has nothing to do with in-place systems:

ds = Systems.lorenz_iip()
xg = yg = range(-30.,30.; length=10)
zg = range(0,50; length=10)
@time basin, attractors = basins_of_attraction((xg,yg,zg), ds)
0.069029 seconds (6.73 k allocations: 508.406 KiB)
attractors
Dict{Int16, Dataset{3, Float64}} with 1 entry:
  1 => 3-dimensional Dataset{Float64} with 87 points

The axis were not well defined. The trajectory spends too much time out of the grid and is lost for the algorithm. The bug was an incorrect initialization of the counter that was looping forever and not identifying the lost state.

With the grid that you have defined earlier we have:

ds = Systems.lorenz_iip()
xg = yg = range(-5.,5.; length=5)
zg = range(0,10; length=5)
@time basin, attractors = basins_of_attraction((xg,yg,zg), ds)
  8.511954 seconds (10.23 k allocations: 1.068 MiB)

Which is worse because the algorithm spends a lot of time outside the grid without finding the attractor (for example when the attractor is spinning around one of the focus).