Chaostools reinit! function

Specific Domains Modelling & Simulations
1

Hi,
so I am calculating Lyapunov exponents for a lot of initial values, parameters for specific systems. I want to make heavy use of the reinit! function for integrator objects, especially the tangent_integrator from the DynamicalSystems package. A thing I noticed is, that the reinit! function does not seem to reinit every property of the integrator object. See the following example

function lorenz_iip(du,u,p,t)
    sigma = p[1]
    rho = p[2]
    beta = p[3]

    # f
    du[1] = sigma*(u[2]-u[1])
    du[2] = u[1]*(rho-u[3])-u[2]
    du[3] = u[1]*u[2]-beta*u[3]

end

function lorenz_jac_iip(J,u,p,t)

    sigma = p[1]
    rho = p[2]
    beta = p[3]

    J[1,1] = -sigma
    J[1,2] = sigma
    J[1,3] = 0
    J[2,1] = rho-u[3]
    J[2,2] = -1
    J[2,3] = -u[1]
    J[3,1] = u[2]
    J[3,2] = u[1]
    J[3,3] = -beta
end

Random.seed!(42)
N=10000
tol = 1e-12
w0 = Matrix(LinearAlgebra.qr(Random.rand(3,3)).Q)
param = [16.,45.92,4.]
x_start = [19.,20.,50.]
cds = ContinuousDynamicalSystem(lorenz_iip,x_start,param,lorenz_jac_iip)
integ = tangent_integrator(cds; abstol=tol, reltol=tol)

lyapunovs(integ, N, 1.0,1000.0)

3-element Array{Float64,1}:
1.499473772711736
0.0001330755279767297
-22.499144329887265

reinit!(integ, x_start, w0)
lyapunovs(integ, N, 1.0,1000.0)

3-element Array{Float64,1}:
1.5036130016404863
2.5060873222979467e-5
-22.50343335858572

integ = tangent_integrator(cds; abstol=tol, reltol=tol)
lyapunovs(integ, N, 1.0,1000.0)

3-element Array{Float64,1}:
1.499473772711736
0.0001330755279767297
-22.499144329887265

For more chaotic systems the difference can of course be much bigger. What is the reason here?

1 Like
2

IIRC it doesn’t reinit the PI-adaptive control variables.

3

I never heard of these. I guess they control the stepsize and some inital conditions/params are chosen at random?

4

They control the dts, and so if your tolerance is high it can be the jiggle in the result

5

I chose a tolerance of 1e-12 - is this still high? I find the jiggle quite significant…

6

lyapunovs might have more tolerances.

7

I just checked different timeseries of Lyapunov exponents for the Lorenz system. I always initialized a new tangent integrator object for each calculation and printed the timeseries for the highest Lyapunov exponent, skipping the first 5000 steps:

lorenz_ts_different_pre_iters

Different results after reinitializing the integrator object seem to lie in the tolerance:

reinit!(integ, x_start, w0)
lyapunovs(integ, N, 1.0,1000.0)

3-element Array{Float64,1}:
1.5022180529877711
4.0746343093052365e-6
-22.502156732106904

reinit!(integ, x_start, w0)
lyapunovs(integ, N, 1.0,1000.0)

3-element Array{Float64,1}:
1.5012094743096116
-1.1035548815286304e-5
-22.501447250951443

reinit!(integ, x_start, w0)
lyapunovs(integ, N, 1.0,1000.0)

3-element Array{Float64,1}:
1.5029591604277805
3.2801567874854977e-6
-22.503013094270816

I probably was a bit over-concerned about this, because my original models had higher relative deviations for reinitialization of the integrator :smiley:

1 Like
8

Hey there @jamblejoe , thanks for opening this. I don’t know exactly where these differences are coming from, but I don’t think they are specifically due to DynamicalSystems or the Lyapunov exponent algorithm. The function reinit! that you quote does the following: (from https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/continuous.jl#L217 )


function DiffEqBase.reinit!(integ::AbstractODEIntegrator, u0::AbstractVector,
    Q0::AbstractMatrix; kwargs...)

    set_state!(integ, u0)
    set_deviations!(integ, Q0)
    reinit!(integ, integ.u; kwargs...)
end

I.e. it is just a wrapper of the true reinit! of DiffEqBase. Therefore anything that the true reinit! does, this function does as well. My bet would be just what Chris said, i.e. the adaptivity parameters. Do you get the same different results if you make non adaptive integration, with e.g. adaptive = false as a keyword or just using a non-adaptive solver all together?

The second comment I’d like to make is that, theoretically, the convergence of the Lyapunov exponents is defined only for infinite time, i.e. there is no guarantee that you will indeed find convergence in finite time.

1 Like
9

Ah, a third comment to make: from my experience while developing DynamicalSystems, I’ve computed λ for various systems. In general, I would expect an error bar of around 1 to 10% of the true exponent value, depending for how long you integrate for. But I honestly would not expect to find the exponent with accuracy within more than 1% of the true value.

But that’s just from my experience, and I’ve never run the algorithm for truly long times.

10

Hi @Datseris,
thanks for pointing the way to the definition of the reinit! function.

After crashing my computer by running out of memory I finished some tests now.

What I realized when using the standard solvers for my models is the “breakout” of some paths. This can be seen in the graphic below where 5 times the same initial conditions were solved with subsequent reinit! calls:
3-chain_ts_100000_reinit

They finally “rejoin”, converge to the Lyapunov exponent, if we consider very long times:
3-chain_ts_1000000_reinit

After some experiments with different solvers I switched to vern9() and got rid of the breakouts of paths for shorter calculation times:
3-chain_ts_100000_reinit_vern9

This speeds up the calculation time by roughly 25-30% as well.

I have not tried with adaptive = false but I assume vern9() uses non-adaptive stepsizes.

That is correct. Oseledets theorem gives theoretical answer to this (Oseledets theorem - Wikipedia). Of course, in practical examples, systems are rarely truely ergodic, but if one is in a non-chaotic region, Lyapunov exponents should go to zero anyways.
In practice, the time, when one would “see” a convergence, can be way beyond experimental times!

This depends of course on the “speed” of the system, so relative errors samller than 1% should be found only for very long calculation times.

1 Like
11

thanks for this very detailed reply, I am happy to see someone so well versed in the topic that they know Oseledet’s theorem :slight_smile:

A quick comment: Vern9 is indeed an adaptive solver.

EDIT: It is expected that solving the system is faster with Vern9 for the Lorenz system, see the discussion here: https://juliadynamics.github.io/DynamicalSystems.jl/dev/chaos/choosing/

1 Like
12

Hey @Datseris,

Thank you very much :slight_smile: I have the feeling that already so much has been done on chaos in dynamical systems and I am just scratching the surface :smiley:

Thanks that you pointed that out for me! I will probably run some more tests on my models.

So far I just had a look at Lyapunov exponents. I did not got into detail about other chaos indicators, especially not in your package, but I hope to do so soon :slight_smile:

1 Like