CallbackSet problem in estimating entry/exit times in a hyper-sphere

[Datseris]

Hi there,

this issue comes from the discourse post https://discourse.julialang.org/t/recording-entry-and-exit-time-of-state-space-sets-with-differentialequations-jl/46029 . There I describe the problem I am trying to solve so I won't repeat it here.

When implementing this problem, I was able to use the callback functionality to solve it for a single sphere radius ε. I thought of using CallbackSet for the case of nested spheres but this leads to problems.

I am attaching a full MWE :

using OrdinaryDiffEq: Tsit5
using DiffEqBase: ODEProblem, solve
using DiffEqBase: ContinuousCallback, CallbackSet
using StaticArrays
using Distances

εdistance(u, u0, ε::Real) = euclidean(u, u0) - ε

# Assumes εs are sorted in reverse order
function transit_time_stats(f, u0, p, εs, T;
        alg = Tsit5(), diffeq...
    )
    eT = eltype(T)
    exits = [eT[] for _ in 1:length(εs)]
    entries = [eT[] for _ in 1:length(εs)]

    # Make the magic callbacks:
    callbacks = ContinuousCallback[]
    for i in eachindex(εs)
        crossing(u, t, integ) = εdistance(u, u0, εs[i])
        negative_affect!(integ) = push!(entries[i], integ.t)
        positive_affect!(integ) = push!(exits[i], integ.t)
        cb = ContinuousCallback(crossing, positive_affect!, negative_affect!;
            save_positions = (false, false)
        )
        push!(callbacks, cb)
    end
    cb = CallbackSet(callbacks...)

    prob = ODEProblem(f, u0, (0.0, T), p)
    sol = solve(prob, alg;
        callback=cb, save_everystep = false, dense = false,
        save_start=false, save_end = false, diffeq...
    )
    return exits, entries
end

function roessler_eom(u, p, t)
    @inbounds begin
    a, b, c = p
    du1 = -u[2]-u[3]
    du2 = u[1] + a*u[2]
    du3 = b + u[3]*(u[1] - c)
    return SVector{3}(du1, du2, du3)
    end
end

u0 = SVector(
    4.705494942754781,
    -10.221120945130545,
    0.06186563933318555,
)
p = [0.2, 0.2, 5.7]
T = 100000.0 # maximum time
εs = sort!([1.0, 0.1, 0.01]; rev=true)

exits, entries = transit_time_stats(roessler_eom, u0, p, εs, T)

and I am also attaching the following figure:

Which shoes the Roessler attractor (the system I am simulating), the u0 as a scatter point, and the radii around u0. It becomes evident that at any point the trajectory must enter the ε bal and then cross it, if it cames close enough to u0.

Now, when running the code I get the following behavior: exits[1] (largest ball crossings) is okay. exits[2] is not okay:

julia> exits[2]
668935-element Array{Float64,1}:
     0.009521130925448734
   105.12850606231594
   227.804685596825
   350.51615792363293
   965.1955493652829
  1000.2144507686401
  1035.2413869702136
  1105.333586114064
  1339.5640443313348
     ⋮
 77971.27432202504
 77971.27432202504
 77971.27432202504
 77971.27432202504
 77971.27432202504

Given the system I've plotted, it is clear that there is some "wrong" behavior here, because the value 77971.27432202504 is repeated as crossing time, and thus the interpolation algorithm clearly gets stuck there.

(entries[3] is still empty, i.e. the trajectory has not yet come close enough to u0 so this part is irrelevant)