`AttractorsViaRecurrences` has trouble when systems have fast and slow dynamics

[KalelR]

Issue Description

When systems have fast and slow dynamics, AttractorsViaRecurrences can identify the slow sections as attractors. Take for concreteness the system

with parameters shown below. It's not even a so-called fast-slow system, but it has fast-slow dynamics: it has a limit cycle that passes close to two saddle points, one at the origin and another at (K,0) (It's almost a heteroclinic cycle). The trajectories slow down near the two saddle points, as they map near the stable manifolds of the saddles, and then speed up as they are ejected via the unstable manifolds of the saddles. So two sections are fast, two sections are slow.

The phase plane is like this:

Expected Behavior

The only attractor of the system is the limit cycle, but AttractorsViaRecurrences is returning the two saddle points exactly as attractors (see the MWE below) I'm actually surprised it identified the two saddle points so precisely, I don't know why that happens since they are not mapped themselves. I would have expected that it identifies some points in their neighborhood. But the idea why this is happening is that the trajectories spend a lot of time mapping near them, and it seems it's hard for the algorithm to distinguish between regions of very slow evolution and attractors.

I am using this example because it's rather simple (only one attractor) and easy to understand the dynamics. A friend of mine was also trying the method for a fast-slow system, and I think it was even messier there. I can get the details later.

MWE

MWE stands for Minimal Working Example.

using DynamicalSystems

#prey = N; predator = P 
@inbounds @inline function predator_prey(u, p, t)
    α, γ, ϵ, ν, h, K, m = p
    N, P = u
    du1 = α*N*(1 - N/K) - γ*N*P / (N+h)
    du2 = ϵ * (  ν*γ*N*P/(N+h) - m*P    )
    return SVector{2}(du1, du2)
end

γ = 2.5
h = 1 
ν = 0.5 
m = 0.4
ϵ = 1.0
α = 0.8 
K = 15 

u0 = rand(2)
ds = ContinuousDynamicalSystem(predator_prey, u0, [α, γ, ϵ, ν, h, K, m])
xg = yg = range(0, 15; length = 50)
mapper = AttractorsViaRecurrences(ds, (xg, yg))
basins, attractors = basins_of_attraction(mapper; show_progress = false)
attractors[1]
attractors[2]

#increasing a lot the parameters, takes longer
xg = yg = range(0, 15; length = 1000)
mapper = AttractorsViaRecurrences(ds, (xg, yg); mx_chk_fnd_att=100000000, mx_chk_hit_bas=10000, Δt=0.1)
basins, attractors = basins_of_attraction(mapper; show_progress = false)
attractors[1]
attractors[2]

# example trajectory 
T = 1000
Ttr= 0
Δt = 0.1
tr = trajectory(ds, T; Δt, Ttr); t = 0:Δt:T
using GLMakie
fig = Figure()
ax1 = Axis(fig[1, 1], xlabel="t", ylabel="Prey-Predator")
ax2 = Axis(fig[2, 1], xlabel="Prey", ylabel="Predator")
scatterlines!(ax2, tr[:,1], tr[:, 2], linewidth=1, color=t, marker=:circle)
lines!(ax1, t, tr[:, 1], linewidth=1, color=:black, label="Prey")
lines!(ax1, t, tr[:, 2], linewidth=1, color=:blue, label="Predator")
fig[1,2] = Legend(fig, ax1)

Is there a simple solution to this?

[Datseris]

xg = yg = range(0, 15; length = 1000), but this means that the unstable fixed points are also included in the set of initial conditions, right? Then it makes sense that the algorithm finds them. After all, these points are mapped to themselves (that's the very definition of a fixed point, unstable or not). The question that remains is why it doesn't find the limit cycle. For this, I think you can increase the steps required to stay within an attractor to believe that you belong to it. I.e., increase mx_chk_att.

Actually, I now see a scenario where the algorithm will never terminate. If the boxes are such that contain the two "found" attractors (unstable fixed points), but also contain the limit cycle, we will have a scenario where the algorithm will constantly swap between searching for a new attractor (the limit cycle) and its counter being reset when it reaches the box of the fixed points. But it won't be mapped to the fixed points if mx_ch_att is high enough, so then it goes again into the next phase of finding attractors....

Can you please check and confirm this? I guess you get it automatically with the grid you have an dby increasing mx_chk_att.

Thus, the solution is to have a grid that does not explicitly have the unstable fixed points as initial conditions and increase mx_chk_fnd_att so that the slow lingering of hte trajectory next to the saddle points doesn't consider them as attractors.

[Datseris]

Which also gives me the idea that we should have an internal counter of max number of steps that would terminate the algorithm and return an error, so that the algorithm doesn't run forever without the user realizing it. An internal counter that is always incremented, and never reset, separate from all other countrers, with really high threshold value like 1_000_000_000 or so.

[Datseris]

also notice that because of the way we scan the grid of initial conditions, the very first trajectory the algorithm initialiuzes is the one of the unstable fixed point in the bottom left. So it finds this "attractor" immediatelly :) hahaha

[Datseris]

but this is an amazing issue. thanks a lot for opening it up. It shows that we clearly need to add some more advice on potential pitfalls of the algorithm on the next paper.

[awage]

I found a fix for this. The problem is not trivial since you have very long transient near the saddles. So the trick is to use a very fine grid so the algorithm does not detect recurrences near unstable points. For this, we can use: A sparse matrix!

Anyway we have to think more about this problem. Maybe we can use adaptive grids for this kind of problems.

xg = yg = range(0, 15; length = 10000)
mapper = AttractorsViaRecurrences(ds, (xg, yg); sparse = true, mx_chk_fnd_att=1000, mx_chk_hit_bas=1000, Δt=0.1)
mapper(rand(2))
u = mapper.bsn_nfo.attractors[1]
plot(u[:,1], u[:,2], seriestype = :scatter)

[KalelR]

xg = yg = range(0, 15; length = 1000), but this means that the unstable fixed points are also included in the set of initial conditions, right? Then it makes sense that the algorithm finds them. After all, these points are mapped to themselves (that's the very definition of a fixed point, unstable or not). The question that remains is why it doesn't find the limit cycle. For this, I think you can increase the steps required to stay within an attractor to believe that you belong to it. I.e., increase mx_chk_att.

Makes sense!

Actually, I now see a scenario where the algorithm will never terminate. If the boxes are such that contain the two "found" attractors (unstable fixed points), but also contain the limit cycle, we will have a scenario where the algorithm will constantly swap between searching for a new attractor (the limit cycle) and its counter being reset when it reaches the box of the fixed points. But it won't be mapped to the fixed points if mx_ch_att is high enough, so then it goes again into the next phase of finding attractors....

Can you please check and confirm this? I guess you get it automatically with the grid you have an dby increasing mx_chk_att.

I think that's what's happening. I ran

xg = range(0, 15; length = 1000); yg = range(0.0, 10; length=1000); mapper = AttractorsViaRecurrences(ds, (xg, yg),  mx_chk_att=1e6, mx_chk_fnd_att=1e6, mx_chk_hit_bas=1e6, mx_chk_loc_att = 1e6) #doesnt find any

and the algorithm seems frozen after it finds one of the saddles:

This occurs when I set mx_chk_fnd_bas = 1e6, but not when only mx_chk_att=1e6, mx_chk_fnd_att=1e6.

Thus, the solution is to have a grid that does not explicitly have the unstable fixed points as initial conditions and increase mx_chk_fnd_att so that the slow lingering of hte trajectory next to the saddle points doesn't consider them as attractors.

So, removing the saddle points from the grid makes it so that the algorithm returns pretty quickly but with no attractors at all being found. I've tried the following parameters:

xg = yg = range(1, 14; length = 500); mapper = AttractorsViaRecurrences(ds, (xg, yg), mx_chk_att=10000000, mx_chk_fnd_att=10000, mx_chk_hit_bas=10000) #doesnt find any
xg = yg = range(1, 14; length = 1000); mapper = AttractorsViaRecurrences(ds, (xg, yg), mx_chk_att=10000000, mx_chk_fnd_att=10000, mx_chk_hit_bas=10000) #doesnt find any
xg = yg = range(1, 14; length = 1000); mapper = AttractorsViaRecurrences(ds, (xg, yg), mx_chk_att=1e6, mx_chk_fnd_att=1e6, mx_chk_hit_bas=1e6, mx_chk_loc_att = 1e6) #doesnt find any
xg = range(1, 14; length = 1000); yg = range(1, 10; length=1000); mapper = AttractorsViaRecurrences(ds, (xg, yg), Δt=1.0, mx_chk_att=1e6, mx_chk_fnd_att=1e6, mx_chk_hit_bas=1e6, mx_chk_loc_att = 1e6) #doesnt find any
xg = range(0.1, 14.9; length = 1000); yg = range(0.1, 10; length=1000); mapper = AttractorsViaRecurrences(ds, (xg, yg), Δt=1.0, mx_chk_att=1e6, mx_chk_fnd_att=1e6, mx_chk_hit_bas=1e6, mx_chk_loc_att = 1e6) #doesnt find any

I don't get why, though... By fixing Δt=1.0, and knowing from the typical trajectories that the system spends about 30 to 50 times near the saddle, we know that mx_chk_fnd_att and mx_chk_att should be bigger than 50 at least. Is this thinking right? Then setting to 1e6 should be enough to pass the saddle points.

[KalelR]

I found a fix for this. The problem is not trivial since you have very long transient near the saddles. So the trick is to use a very fine grid so the algorithm does not detect recurrences near unstable points. For this, we can use: A sparse matrix!

Great!

Anyway we have to think more about this problem. Maybe we can use adaptive grids for this kind of problems.

Yeah, that would be a more elegant solution than always having to refine the grid until we get the attractor! haha Maybe adaptively refining the grid, or increasing the amount of checks a lot, when the speed dx/dt of the trajectory becomes very small?

[Datseris]

So, removing the saddle points from the grid makes it so that the algorithm returns pretty quickly but with no attractors at all being found.

I am confused by this. Why...? I can't comprehend it. let's discuss it in our call tomorrow.

Adaptive grids sound good but also seem really hard to code. Adaptive integration time seems much easier doable!