# Lyapunov Spectrum simple code problem

**URL:** <https://discourse.julialang.org/t/lyapunov-spectrum-simple-code-problem/70182>\
**Category:** Numerics\
**Tags:** question, diffeq, numerics, physics\
**Created:** [October 21, 2021, 10:58pm UTC](https://discourse.julialang.org/t/lyapunov-spectrum-simple-code-problem/70182 "2021-10-21T22:58:08Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gabriel\_Carita](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gabriel_carita/32/28087_2.png) [@Gabriel\_Carita](https://discourse.julialang.org/u/Gabriel_Carita)\
**Post date:** [October 21, 2021, 10:58pm UTC](https://discourse.julialang.org/t/lyapunov-spectrum-simple-code-problem/70182/1 "2021-10-21T22:58:08Z")

</div>

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](https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent)

Here is my code:

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![lbenet](https://avatars.discourse-cdn.com/v4/letter/l/35a633/32.png) [@lbenet](https://discourse.julialang.org/u/lbenet)\
**Post date:** [October 22, 2021, 12:10am UTC](https://discourse.julialang.org/t/lyapunov-spectrum-simple-code-problem/70182/2 "2021-10-22T00:10:46Z")

</div>

Hi, and welcome!

As far as I know, DynamicalSystems.jl allows you to compute the Lyapunov spectrum; see [this](https://juliadynamics.github.io/DynamicalSystems.jl/latest/chaos/lyapunovs/#Lyapunov-Spectrum). 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](https://www.cambridge.org/core/books/lyapunov-exponents/0C96D0FCA7D41027741F485A9DA629AE), Arkady Pikovsky and Antonio Politi, Cambridge U. Press (2015). We have implemented that method in [TaylorIntegration.jl](https://github.com/PerezHz/TaylorIntegration.jl), the code is [here](https://github.com/PerezHz/TaylorIntegration.jl/blob/master/src/lyapunovspectrum.jl), which exploits automatic differentiation techniques. Again, the code is open, and you can adapt your code following the [documentation](https://perezhz.github.io/TaylorIntegration.jl/latest/lorenz_lyapunov/).

---

<div class="post-metadata">

**Author:** ![Gabriel\_Carita](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gabriel_carita/32/28087_2.png) [@Gabriel\_Carita](https://discourse.julialang.org/u/Gabriel_Carita)\
**Post date:** [October 22, 2021, 12:43am UTC](https://discourse.julialang.org/t/lyapunov-spectrum-simple-code-problem/70182/3 "2021-10-22T00:43:20Z")

</div>

> [@lbenet](#):
>
> code is open, and you can adapt your code followi

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.