Alternative method for chaotic saddles in high dimensions

[Datseris]

I just came across this paper:

https://doi.org/10.1063/1.4973235

Very interesting. It appears to be an alternative method to the one @reykboerner is using in PR #61 . We should implement it here at some point in the future.

Also relevant for @raphael-roemer .

[awage]

Hi! At some point I was about to push a PR to Chaostools with the stagger and step method [Sweet et. al. (2001)].

I have a code to reproduce one figure of the original paper:

[awage]

I leave the code here. It is undocumented and unoptimized. But it's a starting point!

using DynamicalSystemsBase
# using ChaosTools
using LinearAlgebra:norm
using Random
using Plots

function F!(du, u ,p, n)
    x,y,u,v = u
    A = 3; B = 0.3; C = 5.; D = 0.3; k = 0.4; 
    du[1] = A - x^2 + B*y + k*(x-u)
    du[2] = x
    du[3] = C - u^2 + D*v + k*(u-x)
    du[4] = u
    return 
end

function stagger_trajectory(x0 ,ds, δ, Tm, R_min, R_max)
    integ = integrator(ds)
    T = escape_time!(x0, integ, R_min, R_max)
    xi = deepcopy(x0) 
    while T < Tm 
        Tp = 0; xp = zeros(length(x0))
        while Tp ≤ T
            r = rand_u(δ,length(x0))
            xp = xi .+ r 
            while _is_in_grid(xp, R_min, R_max) == false
                r = rand_u(δ,length(x0))
                xp = xi .+ r 
            end
            # @show xp
            Tp = escape_time!(xp, integ, R_min, R_max)
        end
        xi = xp; T = Tp
    end
    return xi
end

function rand_u(δ, n)
    a = -log10(δ)
    s = (15-a)*rand() + a
    u = (rand(n).- 0.5)
    u = u/norm(u)
    return 10^-s*u
end

function stagger_and_step(x0 ,ds, δ, Tm, N, R_min, R_max)
    integ = integrator(ds)
    xi = stagger_trajectory(x0, ds, δ, Tm, R_min, R_max) 
    δ = 1e-10
    v = Vector{Vector{Float64}}(undef,N)
    v[1] = xi
    @show xi
    for n in 1:N
        if escape_time!(xi, integ, R_min, R_max) > Tm
            set_state!(integ,deepcopy(xi))
            step!(integ)
        else
            Tp = 0; xp = zeros(length(xi))
            while Tp ≤ Tm
                r = δ*(rand(length(xi)).- 0.5)
                xp = xi .+ r 
                while _is_in_grid(xp, R_min, R_max) == false
                    r = δ*(rand(length(xi)).- 0.5)
                    xp = xi .+ r 
                end
                Tp = escape_time!(xp, integ, R_min, R_max)
            end
            set_state!(integ,deepcopy(xp))
            step!(integ)
        end 
        @show xi = get_state(integ)
        v[n] = xi
    end
    return v
end

function escape_time!(x0, integ, R_min, R_max) 
    set_state!(integ,deepcopy(x0))
    integ.t = 1
    x = get_state(integ)
    k = 1; max_step = 10000;
    while _is_in_grid(x, R_min, R_max) 
        k > max_step && break
        step!(integ)
        x = get_state(integ)
        k += 1
    end
    return integ.t
end

function _is_in_grid(u, R_min, R_max) 
    iswithingrid = true
    @inbounds for i in eachindex(R_min)
        if !(R_min[i] ≤ u[i] ≤ R_max[i])
            iswithingrid = false
            break
        end
    end
    return iswithingrid
end

# The region should not contain any attractors. 
R_min = [-4; -4.; -4.; -4.] 
R_max = [4.; 4.; 4.; 4.]

# Initial condition.
x0 = [(R_max[k] - R_min[k])*rand() + R_min[k] for k in 1:length(R_min)]

@show x0
df = DiscreteDynamicalSystem(F!, x0) 
xi = stagger_trajectory(x0, df, 2., 30, R_min, R_max) 
@show Tp = escape_time!(xi, integrator(df), R_min, R_max)

v = stagger_and_step(xi, df, 2., 30, 10000, R_min, R_max) 
v = hcat(v...)'
plot(v[:,1], v[:,3], seriestype = :scatter, markersize = 1)

[Datseris]

awesome thanks!!! @reykboerner is the stagger and step method the same as the one you are trying to implement?

[awage]

These are different method. Stagger and step is more general I think and does not use the same idea. It is much slower because it involves some random search in a domain.

[Datseris]

Oh okay. So then there are at least three methods that can find saddle/edge states:

1. What @reykboerner is implementing in #61
2. The "stagger and step", for which we have the draft code @awage left in this comment: https://github.com/JuliaDynamics/Attractors.jl/issues/118#issuecomment-1852129291
3. The method I outlined in the start of this issue, from the paper by Sala, Leitao and Altmann.

[awage]

I will work on the code of the stagger and step and submit a PR. I don't know how well it works for ODEs.

There is another method for invariant set that is easy to implement but the results looks difficult to analyze:

https://www.ams.org/journals/notices/202206/noti2489/noti2489.html?adat=June/July%202022&trk=2489&galt=none&cat=feature&pdfissue=202206&pdffile=rnoti-p936.pdf

This is the Lagrangian descriptors method. You need to integrate forward and backward the dynamical system up to some time tau and do a simple transformation using a what they call a descriptor. You obtain a picture where the stable and unstable manifolds appear. But you have to extract the information from this picture with a processing.

[Datseris]

cc @reykboerner have you seen this discussion?

[reykboerner]

Sounds interesting! I am not that familiar with these alternative algorithms (stagger-and-step method, algorithms by Sala et al.) but can add another potential candidate to the list: non-invasive feedback control

This is a way to stabilize unstable equilibria and thus explore them in simulation, as explained in Wuyts and Sieber 2023 (see Eq. 19 and Fig. 6) and Sieber, Omel’chenko and Wolfrum 2014. I may try applying this method in the future.