# Chaostools reinit! function

**URL:** <https://discourse.julialang.org/t/chaostools-reinit-function/26046>\
**Category:** Modelling & Simulations\
**Created:** [July 5, 2019, 4:27pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046 "2019-07-05T16:27:00Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 5, 2019, 4:27pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/1 "2019-07-05T16:27:00Z")

</div>

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

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

```

```julia
lyapunovs(integ, N, 1.0,1000.0)

```

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

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

```

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

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 5, 2019, 5:20pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/2 "2019-07-05T17:20:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 5, 2019, 5:21pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/3 "2019-07-05T17:21:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 5, 2019, 5:24pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/4 "2019-07-05T17:24:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 5, 2019, 5:27pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/5 "2019-07-05T17:27:55Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 5, 2019, 5:35pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/6 "2019-07-05T17:35:09Z")

</div>

`lyapunovs` might have more tolerances.

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 5, 2019, 7:02pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/7 "2019-07-05T19:02:05Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/b/1/b14cbc40521f8a250849b90312e4c09ff0dd9c38.png)

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

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

```

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

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

```

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

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

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [July 5, 2019, 8:32pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/8 "2019-07-05T20:32:14Z")

</div>

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](https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/src/continuous.jl#L217) )

```julia

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.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [July 5, 2019, 8:34pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/9 "2019-07-05T20:34:14Z")

</div>

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**.

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 7, 2019, 6:09pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/10 "2019-07-07T18:09:03Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/f/b/fbf1dd6089769d21a2466724182fe2545336b8e2.png)

They finally “rejoin”, converge to the Lyapunov exponent, if we consider very long times:  
 ![3-chain_ts_1000000_reinit](https://global.discourse-cdn.com/julialang/original/3X/7/3/73b6fd15270cb0898b0d5880d9720f2d207fc389.png)

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](https://global.discourse-cdn.com/julialang/original/3X/e/8/e823b69e089f7cb77f789cef9949c6c60180881b.png)

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

> [@](#):
>
> 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?

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

> [@](#):
>
> 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.

That is correct. Oseledets theorem gives theoretical answer to this ([Oseledets theorem - Wikipedia](https://en.wikipedia.org/wiki/Oseledets_theorem)). 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!

> [@](#):
>
> But I honestly would not expect to find the exponent with accuracy within more than 1% of the true value.

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.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [July 7, 2019, 9:26pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/11 "2019-07-07T21:26:36Z")

</div>

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

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/](https://juliadynamics.github.io/DynamicalSystems.jl/dev/chaos/choosing/)

---

<div class="post-metadata">

**Author:** ![jamblejoe](https://avatars.discourse-cdn.com/v4/letter/j/ee7513/32.png) [@jamblejoe](https://discourse.julialang.org/u/jamblejoe)\
**Post date:** [July 7, 2019, 10:50pm UTC](https://discourse.julialang.org/t/chaostools-reinit-function/26046/12 "2019-07-07T22:50:38Z")

</div>

Hey @Datseris,

> [@](#):
>
> thanks for this very detailed reply, I am happy to see someone so well versed in the topic that they know Oseledet’s theorem 🙂

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

> [@](#):
>
> 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/](https://juliadynamics.github.io/DynamicalSystems.jl/dev/chaos/choosing/)

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 🙂