Hello!
I’ve encountered a weird behavior of lyapunovspectrum when I was investigating my nonlinear system. A MWE can be found here:
using ChaosTools, DifferentialEquations, CairoMakie
function couple(x, y, vThr, k, hOffset)
return (vThr - x)/(1 + exp(k*(hOffset-y)))
end
@inbounds function odesys(u::T, p, t=0.0) where T
u1, v1, w1, u2, v2, w2 = u
α, β, d, vThr, k, hOffset = p
du1 = u1*(1 - u1 - α*v1 - β*w1) + d * couple(u1, u2, vThr, k, hOffset)
dv1 = v1*(1 - v1 - α*w1 - β*u1) + d * couple(v1, v2, vThr, k, hOffset)
dw1 = w1*(1 - w1 - α*u1 - β*v1) + d * couple(w1, w2, vThr, k, hOffset)
du2 = u2*(1 - u2 - α*w2 - β*v2) + d * couple(u2, u1, vThr, k, hOffset)
dv2 = v2*(1 - v2 - α*u2 - β*w2) + d * couple(v2, v1, vThr, k, hOffset)
dw2 = w2*(1 - w2 - α*v2 - β*u2) + d * couple(w2, w1, vThr, k, hOffset)
return T([du1, dv1, dw1, du2, dv2, dw2])
end
# initial point
u1 = 0.08020870490181783
v1 = 0.042614670405671906
w1 = 0.7652877851548388
u2 = 0.7652877851548388
v2 = 0.042614670405671906
w2 = 0.08020870490181783
ic = [u1, v1, w1, u2, v2, w2]
# param values
alpha = 0.2451
beta = 2.324
d = 0.3606476945244957
W = 0.1
k = 50.0
h = 0.5
pars = [alpha, beta, d, W, k, h]
# integrator parameters
alg = RK4()
# alg = DP8()
abstol = 1e-9
reltol = 1e-9
dtmax = 5e-1
dtmin = 1e-20
diffeq = (dtmax = dtmax,
dtmin = dtmin,
alg = alg,
abstol = abstol,
reltol = reltol,
)
# doing DS stuff and looking for Lyapunov exponents
dsys = CoupledODEs(odesys, ic, pars; diffeq)
tanDsys = TangentDynamicalSystem(dsys)
lyaps = lyapunovspectrum(tanDsys, 100_000; Ttr=3000.)
@show lyaps
lle = lyapunov(dsys, 100_000.)
@show lle
# plot projections
totalTime = 6_000
X, t = trajectory(dsys, totalTime)
u1, v1, w1, u2 , v2 , w2 = columns(X)
fig1 = lines(u1, w1)
save("fig1.png", fig1)
While the projection of a found attractor looks like a simple limit cycle, there is a weird discrepancy between the calculations of the Largest Lyapunov Exponent.
The output on my computer is the following:
lyaps = [0.001347582572478865, 0.0013624061655748648, 8.583893939022226e-6, -0.1916855026931509, -0.8523570829071934, -1.276863571175826]
lle = 3.3207502234233297e-6
The lle variable (LLE computed by Benettin algorithm from lyapunov function) is close to zero and corresponds to a limit cycle, while the top Lyapunov exponent of spectrum vector lyaps is noticeably larger than zero. Increasing the time of LE computations does not seem to change that behaviour. Is there any caveat that I might be overlooking?