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.