How to get the first 30 Lyapunov exponents of this modle?

[YouknowIdoknow]

Hi, Sir. Here I am studying a laser system model, I want to calculate the Kolmogorov–Sinai entropy, which is the sum of all positive Lyapunov exponents of the chaotic system. I have no idea how to calculate the 30 Lyapunov exponents since the model with delay time. Can you help me get the first 30 Lyapunov exponents? The initial value is [0.01 0.01 0.01]

[Wrzlprmft]

What do you have so far? I have no problem to help you with a specific problem, but I won’t do the legwork for you. Also, I cannot imagine that you need need thirty Lyapunov exponents. With such a model, you usually have a few positive Lyapunov exponents, one zero, a few moderately negative ones. All further Lyapunov exponents are extremely negative, very difficult to compute as a consequence, and not informative.

[YouknowIdoknow]

This is my script:

from jitcdde import jitcdde_lyap, y, t, anchors, past_y
import numpy as np
from scipy.stats import sem
from symengine import cos, sin, pi, Symbol
from ttictoc import tic, toc

# --------------------------the system parameters----------------------------
rc = 5.36e11
rn = 7.53e9
rs = 5.96e9
J = 1.222
rp = 1.91e10
b = 3.2
si = 0.09
m = 0.56
fm = 36e9
fi = 41.8e9
kf = 0.026
# τ = 2.4e-9
τ = Symbol("tau")

f = [
    0.5 * (rc * rn / (rs * J) * y(2) - rp * (y(0) * y(0) + y(1) * y(1) - 1)) * (y(0) + b * y(1)) + si * (
                1 + 2 * m * cos(2 * pi * fm * t)) * rc * cos(2 * pi * fi * t) + kf * rc * y(0, t - τ),
    0.5 * (rc * rn / (rs * J) * y(2) - rp * (y(0) * y(0) + y(1) * y(1) - 1)) * (y(1) - b * y(0)) - si * (
                1 + 2 * m * cos(2 * pi * fm * t)) * rc * sin(2 * pi * fi * t) + kf * rc * y(1, t - τ),
    -(rs + rn * (y(0) * y(0) + y(1) * y(1))) * y(2) - rs * J * (y(0) * y(0) + y(1) * y(1) - 1) + rs * rp / rc * J * (
                y(0) * y(0) + y(1) * y(1)) * (y(0) * y(0) + y(1) * y(1) - 1),
]

n_lyap = 10
DDE = jitcdde_lyap(
    f,
    n_lyap=n_lyap,
    control_pars=[τ],
    max_delay=10,
    simplify=True,
    verbose=False,
)
print(DDE.delays)
DDE.compile_C(simplify=True, do_cse=True)
taus = np.linspace(1e-9, 2e-9, 5)

dt = 1e-12
times = np.arange(dt, 400e-9, dt)
pre = int(40e-9 // dt)

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

    DDE.constant_past([0.01, 0.01, 0.01], 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('LEs_tau.dat', 'ab') as f:
        np.savetxt(f, [[tau, *LEs]], delimiter="\t")

    DDE.purge_past()

print(toc())

[Wrzlprmft]

Okay, and what is your problem with it? Why did you not apply the answer to Issue #44?

[YouknowIdoknow]

You mean that I need rescale parameters such that the dynamical variables and timescale of the dynamics are roughly in the order of magnitude of 1?

[Wrzlprmft]

Yes. For Lyapunov exponent there is no way around that because the tangent vectors will be scaled in the order of magnitude of 1.

[YouknowIdoknow]

OK，I try to do it. Thank you so much.