Instantaneous Lyapunov exponent for systems with parameter drift (#345)

* add EAPD

* revised version and tests (still need docs)

* add some docstrings

* fix docstring part

* include EAPD test

* minor fixes, improve docstring, tests

* fix dummy parameter test

* increase tolerance for testing

* Update lyapunovs.jl: relax test

* fix comment

* add EAPD to docs

* fix citation

* add changelog and increase version number

* correct vertsion number

Updated version from v3.3.2 to v3.4 in CHANGELOG.md.

* Correct version number

---------

Co-authored-by: George Datseris <[email protected]>

Author: rusandris

CHANGELOG.md

@@ -1,5 +1,9 @@
 # main
 
+# v3.4
+
+ - Added instantaneous Lyapunov exponent for systems with parameter drift: `lyapunov_instant` uses slope of the ensemble-averaged pairwise distance function returned by `ensemble_averaged_pairwise_distance` to assess chaoticity at a certain time
+
 # v3.3.1
 
 After 7 years, we only now realized that `lyapunov` gave incorrect results for fixed points of continuous time systems. We've now fixed that. Unfortunately this decreases the computational performance of the function overall, but correctness is more important.

Project.toml

@@ -1,7 +1,7 @@
 name = "ChaosTools"
 uuid = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 repo = "https://github.com/JuliaDynamics/ChaosTools.jl.git"
-version = "3.3.1"
+version = "3.4.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

docs/src/lyapunovs.md

@@ -134,7 +134,6 @@ Colorbar(fig[1,2], obj)
 fig
 ```
 
-
 ## Lyapunov exponent from data
 
 ```@docs
@@ -192,4 +191,75 @@ fig
 
 As you can see, using `τ = 15` is not a great choice! The estimates with
 `τ = 7` though are very good (the actual value is around `λ ≈ 0.89...`).
-Notice that above a linear regression was done over the whole curves, which doesn't make sense. One should identify a linear scaling region and extract the slope of that one. The function `linear_region` from [FractalDimensions.jl](https://github.com/JuliaDynamics/FractalDimensions.jl) does this!
+Notice that above a linear regression was done over the whole curves, which doesn't make sense. One should identify a linear scaling region and extract the slope of that one. The function `linear_region` from [FractalDimensions.jl](https://github.com/JuliaDynamics/FractalDimensions.jl) does this!
+
+
+## Instantaneous Lyapunov exponent for non-autonomous systems
+```@docs
+ensemble_averaged_pairwise_distance
+```
+
+```@docs
+lyapunov_instant
+```
+
+Let's see first if the ensemble approach is equivalent to the usual time-averaging case (Benettin algorithm) in the autonomous case.
+Here is a simple example using the Henon map:
+```@example MAIN
+henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+ds = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3, 0.0])
+
+init_states = StateSpaceSet(0.2 .* rand(1000,2))
+pidx = 3 # set to dummy, not used anywhere (no drift)
+ρ,times = ensemble_averaged_pairwise_distance(ds,init_states,100,pidx;Ttr=5000)
+λ_inst = lyapunov_instant(ρ,times;interval=20:30) #fit to middle part of the curve (slope is constant until saturation)
+λ = lyapunov(ds,1000;Ttr=5000) #standard (Benettin) way
+@show λ_inst, λ   
+```
+
+Now look at the nonautonomous Duffing map with drifting ε parameter:
+```@example MAIN
+
+using CairoMakie
+using ChaosTools
+
+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
+
+duffing = duffing_drift() 
+duffing_map = StroboscopicMap(duffing,2π)
+init_states = randn(5000,2) #use an ensemble of 5000 
+pidx = 4 #ε is the fourth parameter
+ρ,times = ensemble_averaged_pairwise_distance(duffing_map,StateSpaceSet(init_states),100,pidx;Ttr=20)
+
+#measure slope of ρ at two places  
+λ_inst = lyapunov_instant(ρ,times;interval=5:10)
+λ_inst2 = lyapunov_instant(ρ,times;interval=22:25)
+
+fig,ax,obj = scatter(times, ρ, 
+    markersize = 6,
+    color = :gray10,
+    label = L"\omega = 1, \beta = 0.2, \epsilon_0 = 0.4, \alpha=0.00045",
+    axis = (xlabel = L"t", ylabel = L"\rho(t)",xlabelsize = 20,ylabelsize = 20))
+
+lines!(ax, times[5:10], ρ[5] .+ λ_inst*[0:5;], color = (:red, 0.7),linewidth = 3)
+lines!(ax, times[22:25], ρ[22] .+ λ_inst2*[0:3;],color = (:red, 0.7), linewidth = 3)
+
+text!(ax, times[5]+8, ρ[10],text = L"\lambda = %$(round(λ_inst;digits=3))",
+color = :red,align = (:left, :center),fontsize = 18)
+
+text!(ax, times[22]+8, ρ[24],text = L"\lambda = %$(round(λ_inst2;digits=3))", 
+color = :red, align = (:left, :center), fontsize = 18)
+
+axislegend(ax, position = :rb, nbanks = 2,labelsize = 18)
+fig
+```
+

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,155 @@
+export ensemble_averaged_pairwise_distance,lyapunov_instant
+
+"""
+    ensemble_averaged_pairwise_distance(ds, init_states::StateSpaceSet, T, pidx;kwargs...) -> ρ,t
+
+Calculate the ensemble-averaged pairwise distance function `ρ` for non-autonomous dynamical systems 
+with a time-dependent parameter, using the metod described by [^Jánosi,Tél2024]. Time-dependence is assumed to be a linear drift. The rate of change
+of the parameter needs to be stored in the parameter container of the system `p = current_parameters(ds)`,
+at the index `pidx`. In case of autonomous systems (with no drift), `pidx` can be set to any index as a dummy.
+To every member of the ensemble `init_states`, a perturbed initial condition is assigned.
+`ρ(t)` is the natural log of state space distance between the original and perturbed states averaged 
+over all pairs, calculated for all time steps up to `T`. 
+
+
+## Keyword arguments
+
+* `initial_params = deepcopy(current_parameters(ds))`: initial parameters 
+* `Ttr = 0`: transient time used to evolve initial states to reach 
+    initial autonomous attractor (without drift)
+* `perturbation = perturbation_normal`: 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 
+
+
+## Description
+In non-autonomous systems with parameter drift, long time averages are less useful to assess chaoticity.
+Thus, quantities using time averages are rather calculated using ensemble averages. Here, a new 
+quantity called the Ensemble-averaged pairwise distance (EAPD) is used to measure chaoticity of 
+the snapshot attractor/ state space object traced out by the ensemble [^JánosiTél2024].
+
+To any member of the original ensemble (`init_states`) a close neighbour (test) is added at an initial distance `ϵ`. Quantity `d(t)` is the 
+state space distance between a test particle and an ensemble member at time t .
+If `init_states` are randomly initialized (far from the attractor at the initial parameter), and there's no transient, 
+the first few time steps cannot be used to calculate any reliable averages.
+The function of the EAPD `ρ(t)` is defined as the average logarithmic distance between original and 
+perturbed initial conditions at every time step: `ρ(t) = ⟨ln d(t)⟩`
+
+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.
+
+## Example
+Example of a rate-dependent (linearly drifting parameter) system
+
+```julia
+#r parameter is replaced by r(n) = r0 + R*n drift term
+function drifting_logistic(x,p,n)
+    r0,R = p
+    return SVector( (r0 + R*n)*x[1]*(1-x[1]) )
+end
+
+r0 = 3.8 #inital parameter
+R = 0.001 #rate parameter
+p = [r0,R] # pidx = 2
+
+init_states = StateSpaceSet(rand(1000)) #initial states of the ensemble 
+ds = DeterministicIteratedMap(drifting_logistic, [0.1], p)
+ρ,times = ensemble_averaged_pairwise_distance(ds,init_states,100,2;Ttr=1000)
+```
+
+[^JánosiTél2024]: Dániel Jánosi, Tamás Tél, Phys. Rep. **1092**, pp 1-64 (2024)
+"""
+function ensemble_averaged_pairwise_distance(ds,init_states::StateSpaceSet,T,pidx;
+    initial_params = deepcopy(current_parameters(ds)),Ttr=0,perturbation=perturbation_normal,Δt = 1,ϵ=sqrt(dimension(ds))*1e-10)
+
+	set_parameters!(ds,initial_params)
+    original_rate = current_parameter(ds, pidx)
+    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
+    set_parameter!(pds,pidx,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_parameter!(pds,pidx,original_rate) 
+
+    #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_normal(ds,ϵ)
+    D, T = dimension(ds), eltype(ds)
+    p0 = randn(SVector{D, T})
+    p0 = ϵ * p0 / norm(p0)  
+end
+
+"""
+    lyapunov_instant(ρ,times;interval=1:length(times)) -> λ(t)
+
+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=eachindex(times))
+    _,s = linreg(times[interval], ρ[interval]) #return estimated slope  
+    return s
+end

test/chaosdetection/EAPD.jl

@@ -0,0 +1,59 @@
+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, 0.0]) # add third dummy parameter
+
+#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(1000,2))
+    pidx = 3 # set to dummy, not used anywhere (no drift)
+    ρ,times = ensemble_averaged_pairwise_distance(ds,init_states,100,pidx;Ttr=5000)
+    lyap_instant = lyapunov_instant(ρ,times;interval=20:30)
+
+    #lyapunov exponent
+    λ = lyapunov(ds,1000;Ttr=5000)
+    @test isapprox(lyap_instant,λ;atol=0.05)
+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)
+    pidx=4
+
+    ρ,times = ensemble_averaged_pairwise_distance(duffing_map,init_states_auto,100,pidx;Ttr=0)
+    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() 
+    duffing_map = StroboscopicMap(duffing,2π)
+    init_states = randn(5000,2) 
+    ρ,times = ensemble_averaged_pairwise_distance(duffing_map,StateSpaceSet(init_states),100,pidx;Ttr=20)
+    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/chaosdetection/lyapunovs.jl

@@ -64,7 +64,7 @@ end
         end
 
         ds = CoupledODEs(lorenz_rule, fill(10.0, 3), [10, 20, 8/3])
-        @test lyapunov(ds, 1000, Ttr = 100) ≈ 0 atol = 1e-3
+        @test lyapunov(ds, 1000, Ttr = 100) ≈ 0 atol = 1e-2
     end
 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")