Hi! I was trying out the lyapunovspectrum function in Julia, for a few standard dynamical systems like Rossler.
I compared the results obtained by using lyapunovspectrum to the spectrum obtained by manually initiating a tangent_integrator and evolving it and rescaling it successively, However the results obtained are largely different, could anyone help me understand the problem.
Thank You.
Here’s my code ~~
## rossler generation
""" generate rossler system
u0~ initial state
"""
function generate_rossler_system(u0; p=[0.2, 0.2, 5.7])
##eom~~
@inline @inbounds function rossler_eom( u, p, t)
a = p[1]; b = p[2]; c = p[3]
du1 = -u[2] - u[3]
du2 = u[1] + a*u[2]
du3 = b + (u[1] - c)*u[3]
return SVector{3}(du1, du2, du3)
end
##jacobian~~
@inline @inbounds function rossler_jac( u, p, t)
a ,b ,c = p
jac = @SMatrix [0 -1 -1;
1 a 0;
(0) (0) (u[1]-c)]
return jac
end
rossler_o = ContinuousDynamicalSystem(rossler_eom, u0, p, rossler_jac)
#rossler_p = ContinuousDynamicalSystem(rossler_eom, u0, perturbed_parameters, rossler_jac)
return rossler_o
end
## tangent integration for deviation vectors ~
function deviation(system, t_rescale, t_steps;Ttr= 100.0, ϵ = 0.0001 )
dim = length(system.u0)
t_int = tangent_integrator(system, dim)
step!(t_int, Ttr, true)
dev = get_deviations(t_int)
set_deviations!(t_int, ϵ*Matrix(qr(dev).Q))
lambda_spectrum = zeros((dim, t_steps))
for i in 1:t_steps
step!(t_int, t_rescale, true)
dev = get_deviations(t_int)
for dv in 1:dim
lambda_spectrum[dv, i] = (1/t_rescale)*log(norm(dev[:,dv])/ϵ )
end
set_deviations!(t_int, ϵ*Matrix(qr(dev).Q))
end
return lambda_spectrum
end
But the results came out as follows
Damn, I wish there was a way to track keywords in discourse.com…
Hi!
Before anything else, notice that your Lyapunov spectrum is incorrect even when using DynamicalSystems.jl. You get 0.09, 0.09, -5 but for a continuous system, at least one exponent must be zero. The most likely error is that you have written the Jacobian wrong. Indeed, the entry J[3,1] must be u[3], not zero.
Try again after fixing the Jacobian and let me know.
2 Likes
Hey! thanks for your reply. I reimplemented the code by modifying the Jacobian but there still is significant differnce between the spectrum obtained from the deviation constructed by me and the lyapunovs function available from your package …
as shown below ~ || ** the rest of the code is the same except the jacobian definition
p.s : I have also found this issue in case where the jacobian is obtained automatically using ForwardDiff ! Do you think it could be the function deviation that I have defined that has error, but I really couldn’t figure the bug by myself!
Also I would like ask this question, for any arbitrary ContinousDynamicalSystem we don’t have a prior estimate of the LLE, so what I mostly resort to in my code is to have an estimate of the LLE using the lyapunov function and then use the reciprocal of the LLE as the Transient Time ; Ttr = 1/ LLE to obtain the spectrum. Given the reciprocal of LLEs is a charactersitic time-scale of chaos in the system, evolving the time for for this muh of time would be enough to settle the trajectory on the attractor (** if any). What do you think of this ?
You can compare with our implementation and see where you deviatie from it: https://github.com/JuliaDynamics/ChaosTools.jl/blob/master/src/chaosdetection/lyapunovs.jl#L93-L133 . The algorithm for computing the exponents relies in making a QR decomposition of the deviation vectors, not getting their norm. When our book on nonlinear dynamics is out, have a look at Section 3.2 and this will explain to you why.
No, this is not valid, at least not in the general case. The time to converge to an attractor and the exponential divergence timescale on the attractor are two different things. A better way to quantify convergence time to attractor would be to calculate the mean dissipation rate of hte dynamical system over all the state space.
1 Like
Hey! I was informed recently that this book authored by you is soon to be released and congratulations on your work. I was wondering if there was any way to have access to it either through any open source medium or through my institutional access, as it would be very much essential for my projects.
Thank you
The printed copy is not out yet that’s why I haven’t done an announcement yet
1 Like
There seems to be an issue importing the fixedpoints function from ChaosTools in my code
What exactly is wrong with this ?
