Hello!
I am trying to calculate the spectrum of Lyapunov exponents and get a message about instability. When I calculate a solution with the same settings method of integration, there is no instability message. How can this be fixed?
Thank you for help!
Code:
function HR_mem(u, p, t)
function sigma(x)
return 1.0 / ( 1.0 + exp( -10.0 * ( x - ( - 0.25 ) ) ) )
end
memristor(z, k1_me, k2_me) = k1_me + k2_me * z^2
a, b, c, d, s, xr, r, I, vs, k1, k2, k1_me, k2_me = p
x1, y1, z1, x2, y2, z2, z = u
du1 = y1 + b * x1 ^ 2 - a * x1 ^3 - z1 + I - k1 * ( x1 - vs ) * sigma(x2) + memristor(z, k1_me, k2_me)*(x2 - x1)
du2 = c - d * x1 ^2 - y1
du3 = r * ( s * ( x1 - xr ) - z1 )
du4 = y2 + b * x2 ^ 2 - a * x2 ^3 - z2 + I - k2 * ( x2 - vs ) * sigma(x1) + memristor(z, k1_me, k2_me)*(x1 - x2)
du5 = c - d * x2 ^2 - y2
du6 = r * ( s * ( x2 - xr ) - z2 )
du7 = x1 - x2
return SVector(du1, du2, du3, du4, du5, du6, du7)
end
a = 1.0; b = 3.0; c = 1.0; d = 5.0;
xr = -1.6; r = 0.01; s = 5.0; I = 4.0; xv = 2.0;
k1= -0.17; k2 = -0.17;
k1_me= 0.2
k2_me = 0.42474916387959866 # 0.4180602006688963
u0 = SA[-1.1913014173660204, -6.354990997836504, 3.5694446819312655, -1.1913014173660204, -6.354990997836504, 3.5694446819312655, 57.58493480354987]# [-1.5, 0.0, 0.0, -2.5, 0.0, 0.0, 0.0]
p = [a, b, c, d, s, xr, r, I, xv, k1, k2, k1_me, k2_me];
tstart = 0.0
tend = 5000.0
tspan = (tstart, tend)
tstep = 0.001
prob = ODEProblem(HR_mem, u0, tspan, p)
sol = solve(prob, RK4(), abstol = 1e-14, reltol = 1e-14, maxiters = 5000000);
integ_set = (alg = RK4(), adaptive = true, abstol = 1e-14, reltol = 1e-14)
ds = CoupledODEs(HR_mem, u0, p, diffeq = integ_set)
Λs = lyapunovspectrum(ds, 100)