Closed rusandris closed 1 year ago
Thank you very much! This is a great step forwards!
Can you please open a PR at DynamicalSystemsBase.jl, in which you define a function
successful_step(integ) = true # this is the default definition of the function
successful_step(integ::DEIntegrator) = ... # the code you have here
which returns true
if the latest step the integrator took was successful? For integrators that don't have this return code field the function will always return true. Then you can use that function in this code here, which will lead to cleaner code as there would be no reason to initialize has_retcode
etc.
Furthermore, instead of length(get_state(ds))
use dimension(ds)
.
lastly, you need to add a test in the test suite that indeed confirms that we get NaN
s for unstble initial conditions.
Hi! I've modified this pull request. Now it's much cleaner with the suggested successful_step
functions. And also, I hope it works for all simple integrators and tests for all retcodes, not just :Unstable
. I've tested the changes locally, but I need help with this
lastly, you need to add a test in the test suite that indeed confirms that we get NaNs for unstble initial conditions.
Thanks, I'll go through the other PR now and give some help.
For the tests of the current PR, you need to add a test that computes the lyapunovspectrum for the unstable system-parameters mentioned in the issue and ensure that the return is NaN.
Hi! I am facing this same problem, but the implemented solution on this post is not on the current version of lyapunov.jl. If you try this minimal code here
using DynamicalSystems, DifferentialEquations
function fHH!(du, u, p, t)
du[1] = (p[4] .+ p[2] * sin.(2*pi*p[1]*t) - p[8]*u[2]^(3)*u[3]*(u[1] - p[5]) - p[9]*u[4]^4*(u[1] - p[6]) - p[10]*(u[1] - p[7]))/p[3]
du[2] = 0.1 * (u[1] + 40.0)/(1.0 - exp(-(u[1] + 40.0) / 10.0))* (1.0 - u[2]) - 4.0 * exp(-(u[1] + 65.0) / 18.0) * u[2]
du[3] = 0.07 * exp(-(u[1] + 65.0) / 20.0) * (1.0 - u[3]) - 1.0 / (1 + exp(-(u[1] + 35.0) / 10.0)) * u[3]
du[4] = 0.01 * (u[1] + 55.0) / (1.0 - exp(-(u[1] + 55.0) / 10.0)) * (1.0 - u[4]) - 0.125 * exp(-(u[1] + 65.0) / 80.0) * u[4]
return nothing
end
function testlyap()
Cm = 1.0 # membrane capacitance, in uF/cm^2
g_Na = 120.0 # maximum conducances, in mS/cm^2
g_K = 36.0
g_L = 0.3
E_Na = 50.0 # Nernst reversal potentials, in mV
E_K = -77.0
E_L = -54.387
I0 = 0.0
Astim = 6660
fStim = 4.0
u0 = [-6.49963792e+01, 5.29550879e-02, 5.95994171e-01, 3.17732402e-01]
diffeq = (alg = Tsit5(), abstol = 1e-7, reltol = 1e-7)
p0 = [fStim, Astim, Cm, I0, E_Na, E_K, E_L, g_Na, g_K, g_L]
forcedHH = ContinuousDynamicalSystem(fHH!, u0, p0; diffeq)
return lyapunov(forcedHH, 10; Ttr = 0, Δt = 0.5)
end
testlyap()
With this specific abstol
and reltol
it returns a warning and simple takes too long to finish. In fact on the example here for T = 10 ms is almost instantaneous, but it returns the error After rescaling, the distance of reference and test states was not d0_lower ≤ d ≤ d0_upper as expected. Perhaps you are using a dynamical system where the algorithm doesn't work
which occurs because d
becomes zero and the program shuts down.
This is a secondary problem that might also just return a NaN
since we have coalescent orbits.
When I update my local file lyapunov.jl
with lines for testing d=0
and the successful_step(pds) || return NaN
everything works fine. Of course, for some cases the problems just occur for specific parameter, integrator (I tried a lot with Rosenbrock23
to account stiffness but still have same problem) and abstol
, reltol
values, but since I'm running over a large variety of parameters, I need to set these beforehand.
So, the question is: is this not implemented because it is provisional and we must resolve 184 or any other reason? Is there a problem with changing only my local files, thinking about reproducibility?
Yes, all functions should have the successful_step(ds) || return NaN
option when sensible. lyapunov
is one. I think I closed this because it was opened in DynamicalSystems.jl v2 and now we are in v3, so the PR would need a rework. Feel free to open a new one that does the same thing. Add your example system above as a test case in the test folder. And that should be it!
This is an attempt to fix issue 184. I've added error checking:
lyapunov
andlyapunovspectrum
now return NaN values immediately in the case of instability. Note that this version only tests for the:Unstable
retcode and only if the integrator is a usual one. In case of the simple integrator types thedt<dtmin
error is raised on instability. I also did some benchmarking onSystems.lorenz
. With error checking:Without error checking: