Instantaneous Lyapunov exponent for systems with parameter drift

[rusandris]

This is a non-autonomous version of Lyapunov exponent for systems with parameter drift, based on Dániel Jánosi, Tamás Tél, 2024 .

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 ρ(t) : the local slope can be considered an instantaneous Lyapunov exponent.

This was previously discussed in #337 .

Example (low-precision version of Fig.11 a) from the article above):

using DynamicalSystems
using LinearAlgebra
using StatsBase
using Plots,LaTeXStrings


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(1000,2) #initalize ensemble
ρ,times = EAPD(duffing_map,init_states,100;Ttr=20,sliding_param_rate_index=4); #calc EAPD
lyapunov_instant(ρ,times;interval=1:20) #slope of first 20 time steps

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

[Datseris]

can you also post the pictures as the artivcle is not open access?

[Datseris]

@rusandris is this PR still a draft or ready to review? Furthermore, are the test failures a result of this PR or unrelated?

[rusandris]

can you also post the pictures as the artivcle is not open access?

I'm not sure if I can just post the screenshots from the original article here, but anyways here are my version of Fig.11.

The slope values came out to be only approximately equal to the ones in the article (for an ensemble of N=10000), but trend-wise it seems to be okay.

[rusandris]

@rusandris is this PR still a draft or ready to review? Furthermore, are the test failures a result of this PR or unrelated?

I've added a test for the regular Henon map (autonomous) just to check if the instantaneous Lyapunov exponent coincides with the regular Lyap exp from the Benettin method (this is basically a sanity check that the time- and ensemble averages are the same).

Oops, I forgot to include my EAPD.jl tests in runtests.jl so the failures are definitely coming from somewhere else. I might need to look into those as well.

[rusandris]

Tests seem to be ok, the error that keeps popping up comes from elsewhere:

Lorenz stable: Test Failed at /home/runner/work/ChaosTools.jl/ChaosTools.jl/test/chaosdetection/lyapunovs.jl:67
  Expression: ≈(lyapunov(ds, 1000, Ttr = 100), 0, atol = 0.001)
   Evaluated: -0.0015200001782929603 ≈ 0 (atol=0.001)

I suggest this needs to be dealt with in a separate PR.

[codecov-commenter]

Codecov Report

:white_check_mark: All modified and coverable lines are covered by tests. :white_check_mark: Project coverage is 62.21%. Comparing base (928bdb0) to head (82fdf05). :warning: Report is 17 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #345      +/-   ##
==========================================
+ Coverage   54.86%   62.21%   +7.35%     
==========================================
  Files          22       26       +4     
  Lines         904     1064     +160     
==========================================
+ Hits          496      662     +166     
+ Misses        408      402       -6     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[rusandris]

Thanks @rusandris ; Let's finish this PR then! You need to:

1. Add the new method to the documentation.

I've added to lyapunovs.md, and it looks fine for a first version (I managed to build it locally), but I couldn't figure out how to render paper references properly yet

[Datseris]

Sorry for being late, I had some personal life problems. I am thinking whether it is possible to combine this with the RateSystem of https://github.com/JuliaDynamics/CriticalTransitions.jl/pull/117 which creates a RateSystem that could be used directly here. But to do this, I must first ask @rusandris : does this functionality here only work / make sense if the forcing is linear? Or the methodology would work with any kind of forcing?

[rusandris]

does this functionality here only work / make sense if the forcing is linear? Or the methodology would work with any kind of forcing?

In the original review paper the parameter drift was chosen to be linear for simplicity. The method itself - calculating the ensemble averaged distance between trajectory pairs - doesn't depend on this linearity. You just have to be careful when taking the slope of ρ(t), as it only can be considered as Lyapunov exponent if the distance between pairs is typically small. And that depends on your step size, not the functional form of the drift. I might be wrong but that's how I understand it right now.

Looking at the code, the function ensemble_averaged_pairwise_distance does use the assumption that the drift is linear. Particularly, pidx parameter is for the index of the drifting rate. But that can be generalized if needed, especially if we want to make this compatible with the RateSystem.

[Datseris]

I'll merge this in for now, as the rate system will take a while more. Thanks @rusandris !