Non-autonomous version of Lyapunov exponent (EADP)

[rusandris]

In non-autonomous systems with parameter drift, long time averages are less useful to assess chaoticity.

There are several different ways to address this, one is related to the ensemble approach. Here, a new quantity called the Ensemble-averaged pairwise distance (EAPD) is used. An analog of the classical largest Lyapunov exponent can be extracted from the EAPD) function: the local slope can be considered an instantaneous Lyapunov exponent.

See for example here , section 3.3 :

The function of the EAPD in continuous time is defined as
$\rho(t) = \langle ln r(t) \rangle$
where the average is taken over an ensemble of trajectories. To any member of the ensemble a test particle is ordered at an initial distance $r_0$ . Quantity $r(t)$ is the dimensionless phase space distance between a test particle and an ensemble member at time $t$.

I was wondering if we could use ParallelDynamicalSystem for handling large ensembles of trajectories of the drifting system.
The test particles would just be a copy of the ensemble (can be even inside the same ParallelDynamicalSystem or separately). For estimating the slope, we could use the already existing tools from FractalDimensions.jl FractalDimensions.slopefit.

As far as I know, there isn't a function for this right now. Let me know if this is a good idea, and whether ChaosTools is the right package. (I've seen CriticalTransitions has some functionality for rate-dependent phenomena but seems more niche).

[Datseris]

Isn't this what https://juliadynamics.github.io/DynamicalSystemsDocs.jl/chaostools/stable/lyapunovs/#ChaosTools.local_growth_rates does? If you pass in a drifting system, does it work?

(if yes, we can add it to the docs as an example and also mention in the function docstring it can be used for EADP)

[rusandris]

If you pass in a drifting system, does it work?

Unfortunately not, or at least not directly. To calculate EAPD at each period (time step), one needs to make sure that all ensemble members are stepped at the same time and that states as well as time are updated between distance calculations.
If I'm not mistaken, local_growth_rates re-initializes the ensemble members to t=0 (more precisely, if one uses a StroboscopicMap, the underlying CoupledODEs are reinitialized) at each call, which is fine for its intended purpose, but I'm not sure how to use it in this case.

For this drifting system at least

using DynamicalSystems

function duffing_drift(u0 = [0.1, 0.25]; ω = 1.0, β = 0.2, ε0 = 0.4, α=0.00045)
    return CoupledODEs(duffing_drift_rule, u0, [ω, β, ε0, α])
end

@inbounds function duffing_drift_rule(x, p, t)
    ω, β, ε0, α = p
    dx1 = x[2]
    dx2 = (ε0+α*t)*cos(ω*t) + x[1] - x[1]^3 - 2β * x[2]
    return SVector(dx1, dx2)
end

my rough implementation for EAPD is this

function EAPD(ds,init_states::Matrix,T;Ttr,ϵ=sqrt(2)*1e-10)
    N = size(init_states)[1]
    ρ = Float64[]
    times = Int64[]

    #save drifting rate value
    α = current_parameters(ds)[4]

    #add test particle to every ensemble member
    init_states_plus_pert = StateSpaceSet(vcat(init_states,init_states))

    #create a pds for the ensemble
    #pds is a ParallelDynamicalSystem
    pds = ParallelDynamicalSystem(ds,init_states_plus_pert)

    set_parameter!.(pds.systems,4,0.0) #set to non-drifting for initial ensemble
    #setting parameter for all systems in one go doesn't work?

    #step system pds to reach attractor(non-drifting)
    #system starts to drift at t0=0.0
    for _ in 1:Ttr
        step!(pds)
    end
    
    #add perturbation to test particles
    for i in 1:N 
        state_i = current_state(pds.systems[i])
        perturbed_state_i = state_i .+ perturbation(ds,ϵ)
        set_state!(pds.systems[N+i],perturbed_state_i)
    end

    set_parameter!.(pds.systems,4,α) #set to drifting for initial ensemble

    #set back time to t0 = 0
    reinit!(pds,current_states(pds))

    #calculate EAPD at each time step
    for _ in 1:T
        ρ_t = EAPD(pds)
        push!(ρ,ρ_t)
        push!(times,current_time(pds))
        step!(pds)
    end

    return ρ,times

end

function EAPD(pds)

    states = current_states(pds)                        
    N = Int(length(states)/2)

    #calculate distance averages
    dists = [norm(states[i] - states[N+i]) for i in 1:N]
    ρ = mean(log.(dists))
    return ρ 

end

function perturbation(ds,ϵ)
    D, T = dimension(ds), eltype(ds)
    Q0 = randn(SVector{D, T})
    Q0 = ϵ * Q0 / norm(Q0)  
end

which successfully reproduces Fig.11 a) from the article above :

duffing = duffing_drift()
duffing_map = StroboscopicMap(duffing,2π)

init_states = randn(1000,2)
@time ρ,times = EAPD(duffing_map,init_states,100;Ttr=20)

pl = plot(times,ρ,ylabel=L"\rho(t)",lc=:gray10,lw=2,xlabel=L"t",legend=false,guidefontsize=28)

I'm not sure if ParallelDynamicalSystem was meant to be hold thousands of systems. Maybe there is a better way to store all the data of the ds objects without redundancy (vectors of parameters).

[Datseris]

okay, so 1) you can definitely contribute the EADP method as one more function (and in teh docstring highlight the delicate difference with the local growth rates).

For 2) , yeah this is a mistake. At the moment there is specialized parallel integrator for CoupledODEs and a generic fallback that copies systems for anything else. But Stroboscopicmap is actually just an ODE so we can reuse the same parallel code with small tweak in the source code. Are you willing to contribute this via a PR?

[rusandris]

But Stroboscopicmap is actually just an ODE so we can reuse the same parallel code with small tweak in the source code.

See the initial PR here .