Lyapunov Spectrum simple code problem

Specific Domains Numerics
, , ,
1

Hello there,

I’m trying to figure out how to solve my lyapunov spectrum handwrite code using the QR method for a simple example (lorenz system), there is something wrong with my code ? I’m trying to reproduce the values obtained in this topic : https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent

Here is my code:

using DynamicalSystems
using Plots, DataFrames, DifferentialEquations, LinearAlgebra
# Lorenz system
function loop(u, p, t)
    σ = p[1]; ρ = p[2]; β = p[3]
    du1 = σ*(u[2]-u[1])
    du2 = u[1]*(ρ-u[3]) - u[2]
    du3 = u[1]*u[2] - β*u[3]
    return SVector{3}(du1, du2, du3)
end
# Jacobian:
function loop_jac(u, p)
    σ, ρ, β = p
    J = @SMatrix [-σ  σ  0;
    ρ - u[3]  (-1)  (-u[1]);
    u[2]   u[1]  -β]
    return J
end

#initial parameters
u0 = [1,1,1]
p0 = [10,28,8/3]
t_span = (0,100.0)
prob = ODEProblem(loop, u0, t_span,p0)
solution = solve(prob, RK4(), reltol=1e-11, abstol=1e-11, dense=true)

global U = [1  0  0; 0  1  0;0  0  1]

df = DataFrame(solution)

lyapuA=[]
lyapuB=[]
lyapuC=[]

append!(lyapuA,1)
append!(lyapuB,1)
append!(lyapuC,1)


AA = 0
BB = 0
CC = 0
for i in 2:length(df.timestamp)
    v0 =solution.u[i]
    dt = df.timestamp[i] - df.timestamp[i-1]
    U_n = (U + loop_jac(v0,p0)*dt)*U
    Q, R = qr(U_n)
    U = Q

    AA = log(abs(diag(R)[1]))/(df.timestamp[i])
    BB = log(abs(diag(R)[2]))/(df.timestamp[i])
    CC = log(abs(diag(R)[3]))/(df.timestamp[i])

    append!(lyapuA,AA)
    append!(lyapuB,BB)
    append!(lyapuC,CC)
end

p2 = plot(df.timestamp,lyapuA)
p3 = plot(df.timestamp,lyapuB)
p4 = plot(df.timestamp,lyapuC)
plot(p2,p3,p4)
2

Hi, and welcome!

As far as I know, DynamicalSystems.jl allows you to compute the Lyapunov spectrum; see this. I do not know the details of the implementation, but the code is open!

Skimming throw the article you linked above, I think it uses some linear approximation for small enough time steps. You can read about a more generic algorithm in Lyapunov Exponents: A Tool to Explore Complex Dynamics, Arkady Pikovsky and Antonio Politi, Cambridge U. Press (2015). We have implemented that method in TaylorIntegration.jl, the code is here, which exploits automatic differentiation techniques. Again, the code is open, and you can adapt your code following the documentation.

3

Thank you Ibenet. I would like to implement it my self to understand the details better, actually i need to understand/implement another chaotic indicator (MEGNO) for more complex problems based on LCN. I’ll check the details you have shared thank you.