Lyapunov spectrum, HyperChaos, DynamicalSystems.jl

Specific Domains Modelling & Simulations
1

I calculate lyapunov spectrum in hyperchaos. Why i have such spectrum:

[0.015581182326147337,
 0.0028325903213687266,
 1.0762649529273773e-6,
 -0.007004756972893843,
 -5.7866575928839605,
 -11.883439707053556]

For hyperchaos need two positive exponent lyapunov, one zero exponent and other negative. Why third exponent don’t zero and not infinitesimal value?
Code:

function sigma(x)
    return @fastmath 1.0 / ( 1.0 + exp( -10.0 * ( x  - ( - 0.25 ) ) ) )
end

function HR(u, p, t)
        
    a, b, c, d, s, xr, r,  I, vs, k1, k2, el_link  = p
    x1, y1, z1, x2, y2, z2 = u
    
    du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + el_link * ( x2 - x1 )
    du2 = c - d * x1 ^2 - y1
    du3 = r * ( s * ( x1 - xr ) - z1 )
    
    du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + el_link * ( x1 - x2 )
    du5 = c - d * x2 ^2 - y2
    du6 = r * ( s * ( x2 - xr ) - z2 )
    return SVector(du1, du2, du3,
                    du4, du5, du6)
end

condition = SA[-0.6925302979542225,
 -2.2083565485726155,
  3.4251812063410423,
 -0.49710798909310244,
 -0.24112901660181446,
  4.585643148476261]

tspan = (0.0, 500000.0)

a = 1.
b = 3.
c = 1.
d = 5.
xr = -1.6
r = 0.01 # 0.01
s = 5.
I = 4.
xv = 2.

k1= -0.17
k2 = k1
k = 0.155

p = SA[a, b, c, d, s, xr, r, I, xv, k1, k2, k]

k_space = range(0.155, 0.161, step = 0.001)

spectrum_array = zeros(6, length(k_space))
condition_array = zeros(6, length(k_space))

int(x) = floor(Int, x)
length(k_space)

for (i, k) in enumerate(k_space)
    if i == 1
        global initialcondition =  condition
    end
    println("Initial condition: $initialcondition"); flush(stdout)
    println("k: $k"); flush(stdout)
    
    condition_array[:, i] = initialcondition
    
    p = SA[a, b, c, d,
        s, xr, r, I, xv, k1, k2, k]
    
    prob = ODEProblem(HR, initialcondition, tspan, p)
    sol = solve(prob, Vern9(), abstol = 1e-11, reltol = 1e-11, dense=false, maxiters = 30000000)
    println("solve complete"); flush(stdout)
    ds_HR = ContinuousDynamicalSystem(HR, initialcondition, p )
    spectrum = lyapunovspectrum(ds_HR, tspan[2]; diffeq = (alg = Vern9(),
                                                            abstol = 1e-11, reltol = 1e-11,
                                                            maxiters = 30000000
                                                            ))
    spectrum_array[1:6, i] = spectrum[1:6]
    
    println("Spectrum: ", spectrum_array[1:6, i]); flush(stdout)
    
    initialcondition = sol[length(sol.u)]
    println(">>>>>>>>>>>>>")
    println("")
end

image for spectrum:

image
image3251×833 41.1 KB

1 Like
2

Sounds like a @Datseris question

3

tl;dr: you don’t get exact answers from floating-point arithmetic.

There are numerous approximations in this calculation that will result in not-exactly-zero values for things that are exactly zero in principle. There’s finite-timestep integration of ODEs, approximating an infinite-time average with a long finite-time sum, maybe using finite-differencing to approximate a Jacobian, doing linear algebra on moderately poorly-conditioned linear algebra in finite-precision arithmetic, and plain old floating-point round-off. For this calculation it looks like you’re getting about seven digits of precision (1e-06 compared to 11.88), which is as good as you can get if the Jacobians are 1st-order finite-difference (square root of machine epsilon, 1e-16).

5 Likes
4

This is of course absolutely off topic but I have to admit that sometimes, somehow I like this forum; some of the questions and answers are just simply outstanding. :- )

1 Like
5

that’s exactly what you have though :smiley: @John_Gibson 's answer is what you need to read in detail why you get this.

1 Like
6

is it possible to get more precise values, for example if I set the Jacobian?

7

Yes, providing the analytic Jacobian is probably the most important way to improve the precision of the results. The other thing I’d do is increase the finite integration time to get a better approximation of the infinite-time average.

I’m not sure of the details of DynamicalSystem.jl’s Jacobian approximation, but there might be a way to increase the finite-difference approximation of the Jacobian from the first-order df/dx = (f(x+\epsilon) - f(x))/\epsilon + O(\epsilon) to the second-order df/dx = (f(x+\epsilon/2) - f(x-\epsilon/2))/\epsilon + O(\epsilon^2). With the right epsilon that could get you from seven digits to ten.

3 Likes
8

Thank you, I will try to calculate the Jacobian analytically and set it

1 Like
9

We use ForwardDiff.jl

2 Likes