Hello everyone, I am writing a code for calculating the Lyapunov indicators of the system:
ks = 0.0:0.0001:1.0; λs = zeros(length(ks), 6)
for (i, k) in enumerate(ks)
set_parameter!(ds, 14, k)
λs[i, :] .= lyapunovspectrum(ds, 300000; Δt = 0.01, show_progress = true)
end
As far as I know, the lyapunovspectrum function integrates the system each time for each new value of the parameter k. So, I want to make sure that with each new integration of the system, the state of the system at the last step of the previous integration is taken as the initial conditions. In python, I wrote this:
initial_state = np.array([-1.5, 0.0, 0.0, -2.5, 0.0, 0.0])
for j, k in enumerate(ks):
twoHRs = [
y(1) - a * y(0) ** 3 + b * y(0) ** 2 - y(2) + I - k1 * (y(0) - v_s) * (1 / (1 + sy.exp(-lambd * (y(3) - Theta)))) + k * (y(3) - y(0)),
c - d * y(0) ** 2 - y(1),
r * (s * (y(0) - x_r) - y(2)),
y(4) - a * y(3) ** 3 + b * y(3) ** 2 - y(5) + I - k2 * (y(5) - v_s) * (1 / (1 + sy.exp(-lambd * (y(0) - Theta)))) + k * (y(0) - y(3)),
c - d * y(3) ** 2 - y(4),
r * (s * (y(3) - x_r) - y(5))
]
print('k = ', k)
n = len(twoHRs)
ODE = jitcode_lyap(twoHRs, n_lyap=n)
ODE.set_integrator("vode")
ODE.set_initial_value(initial_state,0.0)
times = range(1,300000,1)
lyaps = []
for time in times:
lyaps.append(ODE.integrate(time)[1])
# converting to NumPy array for easier handling
lyaps = np.vstack(lyaps)
#we take values from the last point for new initial conditions
final_state = np.array([ODE.y_dict.get(y(0)), ODE.y_dict.get(y(1)), ODE.y_dict.get(y(2)), ODE.y_dict.get(y(3)),
ODE.y_dict.get(y(4)), ODE.y_dict.get(y(5))])
initial_state = final_state
temp = np.zeros(6)
for i in range(n):
lyap = np.average(lyaps[1000:,i])
temp[i] = lyap
stderr = sem(lyaps[1000:,i]) # Note that this only an estimate
print("%i. Lyapunov exponent: % .4f ± %.4f" % (i+1,lyap,stderr))
lce1[j] = temp[0]
lce2[j] = temp[1]
lce3[j] = temp[2]
lce4[j] = temp[3]
lce5[j] = temp[4]
lce6[j] = temp[5]
Tell me, please, is it possible to do the same on Julia?