The Lyapunov exponent spectrum is not smooth

[YouknowIdoknow]

Hi, sir. when I use the JiTcdde to get the Lyapunov exponent spectrum varying with the delay, the Lyapaunov exponent spectrum is not smooth, as shown in the diagram on the right. But when I input the time delays one by one, run them individually to save the final result, and then repeat this operation, I get a smooth Lyapunov exponent spectrum with changing time delays, as shown in the diagram on the left. How can I solve this problem? The following is my code

Code 1:

from jitcdde import jitcdde_lyap, y, t
import numpy as np
from symengine import sin, pi

tau = 8.0
a = 1.2
b = 0.15

def ff(t,y,tau,a,b):
   return [y(1, t - tau), y(2), -a * y(1) - a * y(2) + a * sin(2 * pi * b * y(0))]

n_lyap = 6
DDE = jitcdde_lyap(ff(t,y,tau,a,b), n_lyap = n_lyap)
DDE.constant_past( [1.0,0.0,0.0], time=0.0 )
DDE.step_on_discontinuities()
max_step = 0.05
DDE.set_integration_parameters(max_step=max_step, atol=1e-5)
pre_T = 10
DDE.integrate_blindly(pre_T, max_step)
data = []
lyaps = []
weights = []
LEs = []
for time in np.arange(0, 1000, 0.05):
# for time in np.arange(DDE.t, DDE.t+10000, 10):
    state, lyap, weight = DDE.integrate(time)
    data.append(state)
    lyaps.append(lyap)
    weights.append(weight)

lyaps = np.vstack(lyaps)
np.savetxt("LEs_t.dat", lyaps)
for i in range(n_lyap):
    Lyap = np.average(lyaps[400:, i], weights=weights[400:])
    LEs.append(Lyap)
    print(Lyap)

np.savetxt("LEs.dat", LEs)
with open('eye1.dat','ab') as f:
    np.savetxt(f, LEs, delimiter=" ")

Code 2:

from jitcdde import jitcdde_lyap, y, t
import numpy as np
from scipy.stats import sem
from symengine import sin, pi

a = 1.2
b = 0.15

def ff(t,y,tau,a,b):
   return [y(1, t - tau), y(2), -a * y(1) - a * y(2) + a * sin(2 * pi * b * y(0))]

Tau = np.linspace(0.1, 8, 159)

for tau in Tau:
    n_lyap = 6
    DDE = jitcdde_lyap(ff(t, y, tau, a, b), n_lyap=n_lyap)
    DDE.constant_past([1.0, 0.0, 0.0], time=0.0)
    DDE.step_on_discontinuities()
    max_step = 0.05
    DDE.set_integration_parameters(max_step=max_step, atol=1e-5)
    pre_T = 10
    DDE.integrate_blindly(pre_T, max_step)
    data = []
    lyaps = []
    weights = []
    LEs = []
    for time in np.arange(0, 1000, 0.05):
        state, lyap, weight = DDE.integrate(time)
        data.append(state)
        lyaps.append(lyap)
        weights.append(weight)

    lyaps = np.vstack(lyaps)
    np.savetxt("LEs_t.dat", lyaps)
    for i in range(n_lyap):
        Lyap = np.average(lyaps[400:, i], weights=weights[400:])
        LEs.append(Lyap)
        print(Lyap)

    np.savetxt("LEs.dat", LEs)
    with open('eye3.dat', 'ab') as f:
        np.savetxt(f, LEs, delimiter=" ")

[Wrzlprmft]

You clearly found something profoundly weird.

I can reproduce your results insofar as the maximum Lyapunov exponent is inexplicably small (i.e., practically 0) for certain values of τ, e.g., for τ=4.346249 but not for τ=4.34625. (I didn’t look much at the other Lyapunov exponents, but the phenomenon appears to trickle down to them.) This phenomenon persists no matter what I throw at it. My best guess so far would be a resonance phenomenon.

What I could not reproduce is the absence of this phenomenon when computing Lyapunov exponents separately. Instead I get exactly the same results as when computing them in a Python loop. Can you tell me how exactly you computed all Lyapunov exponents with Code 1? It might give me a hint at what is actually going on. Could it be that you used a different selection of τ values and thus somehow missed the resonance?

[YouknowIdoknow]

I calculate the corresponding Lyapunov exponents by adding up the time delays τ from 0.1 to 8 in intervals of 0.1, e.g. when τ = 0.1, the first set of Lyapunov exponents is obtained, when τ = 0.2, the second set of Lyapunov exponents is obtained, and so on, when τ = 8, the last set of Lyapunov exponents is obtained. Finally, all the Lyapunov indices are combined to obtain the spectrum of Lyapunov indices on the left, which varies from time to time.

[Wrzlprmft]

I calculate the corresponding Lyapunov exponents by adding up the time delays τ from 0.1 to 8 in intervals of 0.1, e.g. when τ = 0.1, the first set of Lyapunov exponents is obtained, when τ = 0.2, the second set of Lyapunov exponents is obtained, and so on, when τ = 8, the last set of Lyapunov exponents is obtained. Finally, all the Lyapunov indices are combined to obtain the spectrum of Lyapunov indices on the left, which varies from time to time.

And how do you enter the τs and run the scripts? I presume that you did not manually alter the value and run the script 80 times?

[YouknowIdoknow]

No, I do manually alter the value without changing any other paprameters and run the script 80 times. After runing one time, I will copy the Lyapunov exponents into a txt file. Furthermore, I have plot the bifurcation diagram varying with τ. Comparing with the Lyapunov exponent spectrum, it can be consistent with each other.

[Wrzlprmft]

For τ ∈ { 0.1, 0.2, …, 8.0}, I do indeed also obtain a smooth spectrum. So, the issue is not due to manually running the script, but due to using these particular values for τ (or rather not using them). However, I am still clueless why on earth this is so strongly connected to the values of τ. I can change the integration steps, integration parameters, control parameters; I can even change the dynamics, including simply speeding it up, and the issue remains.

While I dig deeper into this, does this model have a name or are there any known peculiar dynamical properties that may be of relevance to this?

Finally, instead of manually entering numbers, you can loop over something like:

taus =  [ float(f"{a}.{b}") for a in range(8) for b in range(10) ]

[Wrzlprmft]

I somewhat narrowed the problem down: If the delay only has 14 relevant digits, the problem appears not to occur. Consequentially, a makeshift fix for you is to do something like:

for tau in taus:
    tau = round(tau,14)
    …

Apart from that, consider my mind boggled and my curiosity piqued. I am investigating.

[YouknowIdoknow]

Thank you. Sir, This moldel is our team developed. I have copied the code into my sc ript and run it. But it does not work.

[Wrzlprmft]

You have to skip the 0.0 in the beginning. Anyway, the rounding solution I posted later should be more convenient in any respect.

[YouknowIdoknow]

Ok. Thank you, sir. Moreover, how can I speed up this script? Beacause I have many number which are needed to operate by this script.

[Wrzlprmft]

Moreover, how can I speed up this script?

Here is how I would do it:

Lots of code; click to expand.

``` Python 3 from jitcdde import jitcdde_lyap, y, t import numpy as np from scipy.stats import sem from symengine import sin, pi, Symbol a = 1.2 b = 0.15 τ = Symbol("tau") f = [ y(1,t-τ), y(2), a* (-y(1) - y(2) + sin( 2*pi*b*y(0) ) ), ] n_lyap = 3 DDE = jitcdde_lyap( f, n_lyap=n_lyap, control_pars=[τ], max_delay=8, verbose=False ) taus = np.linspace(1, 5, 1723) dt = 5 times = np.arange(dt,1000,dt) pre = 40//dt for tau in taus: tau = round(tau,14) DDE.constant_past([1.0, 0.0, 0.0], time=0.0) DDE.set_parameters(tau) DDE.max_delay = tau DDE.delays = [tau] DDE.adjust_diff() lyaps = [] weights = [] for time in times: _, lyap, weight = DDE.integrate(time) lyaps.append(lyap) weights.append(weight) lyaps = np.vstack(lyaps) LEs = [] for i in range(n_lyap): Lyap = np.average(lyaps[pre:,i], weights=weights[pre:]) LEs.append(Lyap) print(tau,*LEs) with open('eye2.dat', 'ab') as f: np.savetxt(f, [[tau,*LEs]], delimiter="\t") DDE.purge_past() ```

Mind that the number of Lyapunov exponents can be a major factor. So, ensure that you really need six of them.

[YouknowIdoknow]

Good job! Thank you so much

[Wrzlprmft]

I found the culprit and it all boils down to this issue in SymEngine. Due to this, one of the Jacobians needed for computing the Lyapunov exponents gets miscomputed.

For now, rounding all your delays to 14 digits is the best makeshift solution I can offer.

[YouknowIdoknow]

Ok, I get it.

[Wrzlprmft]

The underlying issue in Symengine has been fixed, so getting their latest version should solve this. (See this on how to do that.)