saving the trajectory from Lyapunov spectrum calculation

[rseydam]

I am trying to get the trajectory for which the spectrum is calculated. I tried it with the following:

_diffeq = (alg = Tsit5(), abstol = 1e-7, reltol = 1e-7, saveat = (1.0:0.1:10.0) )
ds = ContinuousDynamicalSystem(rhs!, u0, p; diffeq)
λ1 = @time Lyapunov spectrum(ds, 1000, p.k; Δt=1.0, Ttr= 1000.0, show_progress = true);_

Then when looking in ds.integ.sol the solution is not retained. How do I keep it?

[pclus]

I agree this option should be included. If you want the trajectory and the Lyapunov exponents right now the only option I found is to make separate computations, which basically means computing the trajectory twice...

[Datseris]

Hi there, apologies for the late reply,

DynamicalSystems.jl operates directly on integrators and never saves intermediate states unless required by the nonlinear dynamics algorithm. That is why currently the saveat option is over-written.

However, it is easy to change this. You need to make a pull requests that edits this code:

https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/main/src/core_systems/continuous_time_ode.jl#L93-L95

The Pull Request should change the order and expand remaining..., special_kwargs..., after it defines save_start = false, save_end = false, save_everystep = false,, so that these options are not overwritten.

Notice however that accessing ds.integ.sol is not public API and will not be officially supported. Use at your own regard. Also, you need to first make tds = TangentDynamicalSystem(...) and then access tds.integ.sol.

I agree this option should be included. If you want the trajectory and the Lyapunov exponents right now the only option I found is to make separate computations, which basically means computing the trajectory twice...

Sure, this is useful, but not part of the algorithm, and not always useful. We have to be careful here. If at every implemented algorithm we allowed every option that someone found useful, then the majority of the source code would be auxiliary options instead of the core code for the nonlinear dynamics algorithm. And source code clarity is one of the cornerstones of DynamicalSystems.jl. The solution above, of simply allowing saveat to be passed to the integrator, is a simple solution for people to be happy enough I hope.

To alter the particular function to achieve what you want, is also super easy though:

1. You copy paste this function: https://github.com/JuliaDynamics/ChaosTools.jl/blob/main/src/chaosdetection/lyapunovs/lyapunovspectrum.jl#L67-L99 and rename it to mylyapunovspectrum,
2. You add the state_history = [current_state(tds)] at the start of the function.
3. In the core step loop, right after https://github.com/JuliaDynamics/ChaosTools.jl/blob/main/src/chaosdetection/lyapunovs/lyapunovspectrum.jl#L89 you add the line push!(state_history, current_state(tds)).
4. You change the return statement to return λ, state_history so that you obtain the exponent and the trajectory.

[rseydam]

Hi, Thank you for the detailed response! Thank you also for explaining why the code the way it is and how to potentially adapt it. In the meantime, I found a way that also works for me. I have included a SavingCallback in _diffeq. In this way, one can have a look at the final bit of the trajectory or potentially other things.

[Datseris]

Oh, that's also another solution, albeit less performant than my suggestion of doing a pull requesti n DynamicalSystemsBase.jl. I think i'll add a documentation section soon explaining how the package works internally so people aren't surprised as you were here. But fixing the (incorrect) over-write of keywords such as saveat is also a good idea in general.