Fix error checking in lyapunov and lyapunovspectrum

[rusandris]

This is an attempt to fix issue 184. I've added error checking: lyapunov and lyapunovspectrum 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 the dt<dtmin error is raised on instability.
I also did some benchmarking on Systems.lorenz.
With error checking:

@btime lyapunov(ds,100, Ttr = 2000 ;diffeq)
  23.199 ms (666 allocations: 50.45 KiB)
0.9056821285687465

@btime lyapunovspectrum(ds,100, Ttr = 2000 ;diffeq)
  27.496 ms (57 allocations: 6.91 KiB)
3-element Vector{Float64}:
   0.9147419224955904
  -0.005725267193229689
 -14.575596505931658

Without error checking:

 @btime lyapunov(ds,100, Ttr = 2000 ;diffeq)
  23.201 ms (666 allocations: 50.45 KiB)
0.9056821285687465

@btime lyapunovspectrum(ds,100, Ttr = 2000 ;diffeq)
  27.469 ms (57 allocations: 6.91 KiB)
3-element Vector{Float64}:
   0.9147419224955904
  -0.005725267193229689
 -14.575596505931658

[Datseris]

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 NaNs for unstble initial conditions.

[rusandris]

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.

[Datseris]

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.

[Jaderpolli]

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?

[Datseris]

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!