Instantaneous Lyapunov exponent for systems with parameter drift

src/ChaosTools.jl

@@ -26,6 +26,7 @@ include("periodicity/davidchacklai.jl")
 include("rareevents/mean_return_times/mrt_api.jl")
 
 include("chaosdetection/lyapunovs/lyapunov.jl")
+include("chaosdetection/lyapunovs/EAPD.jl")
 include("chaosdetection/lyapunovs/lyapunov_from_data.jl")
 include("chaosdetection/lyapunovs/lyapunovspectrum.jl")
 include("chaosdetection/lyapunovs/local_growth_rates.jl")

src/chaosdetection/lyapunovs/EAPD.jl

@@ -0,0 +1,116 @@
+export ensemble_averaged_pairwise_distance,lyapunov_instant
+
+"""
+    lyapunov_instant(ρ,times;interval=1:length(times)) -> λ(t)
+
+Convenience function that calculates the instantaneous Lyapunov exponent by taking the slope of 
+the ensemble-averaged pairwise distance function `ρ` wrt. to the saved time points `times` in `interval`.
+"""
+function lyapunov_instant(ρ,times;interval=1:length(times))
+    _,s = linreg(times[interval], ρ[interval]) #return estimated slope  
+    return s
+end
+
+"""
+    ensemble_averaged_pairwise_distance(ds,init_states::StateSpaceSet,T;kwargs...) -> ρ,t
+
+Calculate the ensemble-averaged pairwise distance function `ρ` for non-autonomous dynamical systems 
+with a time-dependent parameter. Time-dependence is assumed to be linear (sliding). To every member 
+of the ensemble `init_states`, a perturbed initial condition is assigned. `ρ(t)` is the natural log 
+of phase space distance between the original and perturbed states averaged over all pairs, calculated 
+for all time steps up to `T`. 
+
+This function implements the method described in 
+https://doi.org/10.1016/j.physrep.2024.09.003.
+
+## Keyword arguments
+
+* `sliding_param_rate_index = 0`: index of the parameter that gives the rate of change of the sliding parameter
+* `initial_params = deepcopy(current_parameters(ds))`: initial parameters 
+* `Ttr = 0`: transient time used to evolve initial states to reach 
+    initial autonomous attractor (without sliding)
+* `perturbation = perturbation_uniform`: if given, it should be a function `perturbation(ds,ϵ)`,
+   which outputs perturbed state vector of `ds` (preferrably `SVector`). If not given, a normally distributed 
+   random perturbation with norm `ϵ` is added.
+*  `Δt = 1`: step size 
+*  `ϵ = sqrt(dimension(ds))*1e-10`: initial distance between pairs of original and perturbed initial conditions 
+"""
+function ensemble_averaged_pairwise_distance(ds,init_states::StateSpaceSet,T;sliding_param_rate_index=0,
+    initial_params = deepcopy(current_parameters(ds)),Ttr=0,perturbation=perturbation_uniform,Δt = 1,ϵ=sqrt(dimension(ds))*1e-10)
+
+	set_parameters!(ds,initial_params)
+    N = length(init_states)
+    d = dimension(ds)
+    dimension(ds) != d && throw(AssertionError("Dimension of `ds` doesn't match dimension of states in init_states!"))
+
+    nt = length(0:Δt:T) #number of time steps
+    ρ = zeros(nt) #store ρ(t)
+    times = zeros(nt) #store t
+
+    #duplicate every state
+    #(add test particle to every ensemble member)
+    init_states_plus_copies = StateSpaceSet(vcat(init_states,init_states))
+
+    #create a pds for the ensemble
+    #pds is a ParallelDynamicalSystem
+    pds = ParallelDynamicalSystem(ds,init_states_plus_copies)
+
+	#set to non-drifting for initial ensemble
+    sliding_param_rate_index != 0 && set_parameter!(pds,sliding_param_rate_index,0.0) 
+
+    #step system pds to reach attractor(non-drifting)
+    #system starts to drift at t0=0.0
+    for _ in 0:Δt:Ttr
+        step!(pds,Δt,true)
+    end
+
+    #rescale test states 
+    #add perturbation to test states
+    for i in 1:N 
+        state_i = current_state(pds,i)
+        perturbed_state_i = state_i .+ perturbation(ds,ϵ)
+        #set_state!(pds.systems[N+i],perturbed_state_i)
+        set_state!(pds,perturbed_state_i,N+i)
+    end
+
+	#set to drifting for initial ensemble
+    set_parameters!(pds,initial_params) 
+
+    #set back time to t0 = 0
+    reinit!(pds,current_states(pds))
+
+    #calculate EAPD for each time step
+    ensemble_averaged_pairwise_distance!(ρ,times,pds,T,Δt)
+    return ρ,times
+
+end
+
+#calc distance for every time step until T
+function ensemble_averaged_pairwise_distance!(ρ,times,pds,T,Δt)
+    for (i,t) in enumerate(0:Δt:T)
+        ρ[i] = ensemble_averaged_pairwise_distance(pds)
+        times[i] = current_time(pds)
+        step!(pds,Δt,true)
+    end
+end
+
+#calc distance for current states of pds
+function ensemble_averaged_pairwise_distance(pds)
+
+    states = current_states(pds)
+    N = Int(length(states)/2)
+
+    #calculate distance averages
+    ρ = 0.0
+    for i in 1:N
+        ρ += log.(norm(states[i] - states[N+i]))
+    end
+    return ρ/N 
+
+end
+
+function perturbation_uniform(ds,ϵ)
+    D, T = dimension(ds), eltype(ds)
+    p0 = randn(SVector{D, T})
+    p0 = ϵ * p0 / norm(p0)  
+end

test/chaosdetection/EAPD.jl

@@ -0,0 +1,57 @@
+using ChaosTools, Test
+using LinearAlgebra
+
+henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+henon() = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3])
+
+#test if ensemble averaging gives the same as 
+#the usual lyapunov exponent for autonomous system
+@testset "time averaged and ensemble averaged lyapunov exponent" begin
+    ds = henon()
+
+    #eapd slope
+    init_states = StateSpaceSet(0.2 .* rand(5000,2))
+    ρ,times = ensemble_averaged_pairwise_distance(ds,init_states,100;Ttr=1000,sliding_param_rate_index=0)
+    lyap_instant = lyapunov_instant(ρ,times;interval=10:20)
+
+    #lyapunov exponent
+    λ = lyapunov(ds,1000;Ttr=1000)
+    @test isapprox(lyap_instant,λ;atol=0.01)
+end
+
+#test sliding Duffing map 
+#-------------------------duffing stuff-----------------------
+#https://doi.org/10.1016/j.physrep.2024.09.003
+
+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
+
+@testset "Duffing map" begin
+    #----------------------------------hamiltonian case--------------------------------------
+    duffing = duffing_drift(;β = 0.0,α=0.0,ε0=0.08) #no dissipation -> Hamiltonian case
+    duffing_map = StroboscopicMap(duffing,2π)
+    init_states_auto,_ = trajectory(duffing_map,5000,[-0.85,0.0];Ttr=0) #initial condition for a snapshot torus
+    #set system to sliding
+    set_parameter!(duffing_map,4,0.0005)
+
+    ρ,times = ensemble_averaged_pairwise_distance(duffing_map,init_states_auto,100;Ttr=0,sliding_param_rate_index=4)
+    lyap_instant = lyapunov_instant(ρ,times;interval=50:60)
+    @test isapprox(lyap_instant,0.87;atol=0.01) #0.87 approximate value from article
+
+    #-----------------------------------dissipative case------------------------------------
+    duffing = duffing_drift() #no dissipation -> Hamiltonian case
+    duffing_map = StroboscopicMap(duffing,2π)
+    init_states = randn(5000,2) 
+    ρ,times = ensemble_averaged_pairwise_distance(duffing_map,StateSpaceSet(init_states),100;Ttr=20,sliding_param_rate_index=4)
+    lyap_instant = lyapunov_instant(ρ,times;interval=2:20)
+    @test isapprox(lyap_instant,0.61;atol=0.01) #0.61 approximate value from article
+
+end

test/runtests.jl

@@ -25,6 +25,7 @@ testfile("chaosdetection/gali.jl")
 testfile("chaosdetection/partially_predictable.jl")
 testfile("chaosdetection/01test.jl")
 testfile("chaosdetection/expansionentropy.jl")
+testfile("chaosdetection/EAPD.jl")
 
 testfile("stability/fixedpoints.jl")
 testfile("periodicity/periodicorbits.jl")