Lyapunov exponent of stroboscopic map perhaps incorrect?

[Datseris]

Alright, with new design of DynamicalSystems.jl v3 it is rather straightforward to compute the maximum lyapunov exponent of arbitrary dynamical systems. Below, I am attaching a code that computes the Lyapunov exponent of the duffing oscillator first as a normal system and then as a stroboscopic map:

using ChaosTools

# %% Stroboscopic
function duffing_rule(x, p, t)
    ω, f, d, β = p
    dx1 = x[2]
    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
    return SVector(dx1, dx2)
end

u0 = [0.1, 0.25]
p0 = [1.0, 0.3, 0.2, -1]
T = 2π/1.0

duffing_raw = CoupledODEs(duffing_rule, u0, p0)
duffing_map = StroboscopicMap(duffing_raw, T)

λspec = lyapunovspectrum(duffing_raw, 10000)
λmax = lyapunov(duffing_raw, 10000)
λmax_smap = lyapunov(duffing_map, 10000)

@show λspec
@show λmax, λmax*T, λmax_smap

The output is:

λspec = [0.15894265996257015, -0.35894237397571577]
(λmax, λmax * T, λmax_smap) = (0.15884881282285493, 0.9980765267914823, 2.3432028585153932)

Now, as you see above, the exponent we get from the stroboscopic map does not match the period times the exponent of the normal system. Theory says it "should". I am attaching here two pages form the book of Politi and Pikovsky

Lyapunov_exponents.pdf

where they say that for a stroboscopic map one expects that the Lyapunov exponent is the period times the exponent of the normal time system.

However, I've checked the code extensively and I am not sure there is a mistake. The code does what its supposed to do, treating the stroboscopic map as a deterministic iterated map. I hope there is some easy to find bug somewhere so that we don't have to dig deep into theory to check if the statement of the book is actually correct..

[Datseris]

This is guaranteed wrong; the stroboscopic map exponent depends strongly on d0 which doesn't make sense if d0 is small enoough.

[Datseris]

It is probably the internal integrator. If I use Vern9 with tolerance 1e-15 the result is correct. In this case the integrator takes tiny tiny steps and hence does not miss the stroboscopic poincare plane by any meaningful amount. In small tolerances the integrator exceeds the plane by a lot.

I believe that either I don't understand what step!(integ, T, true) does, or that it is not do what I think it does: step, and stop, the integrator at exactly T.

[awage]

I also had some problems with step! when the default time step is too big. What does the integrator is iterating with small steps until a stopping condition is met (in add_tstop!):

function step!(integ::DEIntegrator, dt, stop_at_tdt = false)
    (dt * integ.tdir) < 0 * oneunit(dt) && error("Cannot step backward.")
    t = integ.t
    next_t = t + dt
    stop_at_tdt && add_tstop!(integ, next_t)
    while integ.t * integ.tdir < next_t * integ.tdir
        step!(integ)
        integ.sol.retcode in (ReturnCode.Default, ReturnCode.Success) || break
    end
end

I do not remember the system, but the add_tstop! condition was systematically ignored. Maybe it depends on the solver. In any case, it seems that OrdinaryDiffEq does the job of checking that it stops exactly at T in the function handle_tstop!(integrator) .

It is worth having a look at what is going on with the precision. This way we can also give some guideline on the parameters that should be set for some systems.

[rusandris]

If I use Vern9 with tolerance 1e-15 the result is correct.

I tried this out too, with lower tolerances:

using ChaosTools,OrdinaryDiffEq

# %% Stroboscopic
function duffing_rule(x, p, t)
    ω, f, d, β = p
    dx1 = x[2]
    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
    return SVector(dx1, dx2)
end

u0 = [0.1, 0.25]
p0 = [1.0, 0.3, 0.2, -1]
T = 2π/1.0

diffeq = (alg=Vern9(),abstol = 1e-15,reltol = 1e-10)
duffing_raw = CoupledODEs(duffing_rule, u0, p0;diffeq)
duffing_map = StroboscopicMap(duffing_raw, T)

λspec = lyapunovspectrum(duffing_raw, 10000)
λmax = lyapunov(duffing_raw, 10000)
λmax_smap = lyapunov(duffing_map, 10000)

@show λspec
@show λmax, λmax*T, λmax_smap

Here are the results:

λspec = [0.15721696267070875, -0.3572169626725329]
(λmax, λmax * T, λmax_smap) = (0.15788920209828142, 0.9920471147862301, 0.9835419055970506)

So the result is indeed correct with these diffeq settings.

In small tolerances the integrator exceeds the plane by a lot.

This is what I don't see so far. Using default tolerances, the exact stepping seems to be pretty accurate in terms of time.

duffing_raw = CoupledODEs(duffing_rule, u0, p0)

for i in 1:10000
  step!(duffing_raw,T,true)
end

current_time(duffing_raw)/T ≈ 10000 #true

I think I'm missing something. Any new ideas on this issue?

[Datseris]

I havent seen this issue in a long time... Is there an issue in the end? It appears that if you use very precise integrator the results are correct. Now, why aren't the results correct for tolerances of e.g., 10-8? it is still pretty low.

[rusandris]

Side note. It seems like the lyapunov exponent computed from StroboscopicMap is very sensitive to solver tolerances in general. Here's an example for the Lorenz-system:

using ChaosTools, PredefinedDynamicalSystems

ds_cont = Systems.lorenz()
diffeq = (alg=Tsit5(),abstol = 1e-15,reltol = 1e-10)
ds_cont2 = CoupledODEs(Systems.lorenz_rule,[0.0, 10.0, 0.0], [10.0,28.0,8/3];diffeq)

T = 2π/1.0
strob_map = StroboscopicMap(ds_cont, T)
strob_map2 = StroboscopicMap(ds_cont2, T)

The difference between the stroboscopic exponents is quite large:

λ_strob = lyapunov(strob_map, 100)
λ_strob2 = lyapunov(strob_map2, 100)
@show λ_strob,λ_strob2 # => (11.40971149484336, 5.449189231407157)

The difference in the continuous exponents is small, even if the total time and Δt is approximately the same as in the stroboscopic case.

λ_cont = lyapunov(ds_cont, 100*T;Δt=T)
λ_cont2 = lyapunov(ds_cont2, 100*T;Δt=T)
@show λ_cont,λ_cont2 # => (0.8645614697777381, 0.8656378311723896)

[Datseris]

Okay then this is clearly a bug. Either the way we communicate current_state for a stroboscopic map is incorrect, or the way lyapunov handles current_state is incorrect, or the when lyapunov does the re-scaling of the satellite trajectory is incorrect. It could be something else entirely but one of these three things would be my candidates to check. I don't have the capacity to work on this but if someone does it would a well-appreciated fix.