From 9771c55e6c2484208368d6c4f75cde80399693fd Mon Sep 17 00:00:00 2001 From: rusandris Date: Sat, 2 Jul 2022 20:18:13 +0300 Subject: [PATCH 1/4] Fix error checking in lyapunov and lyapunovspectrum --- src/chaosdetection/lyapunovs/lyapunov.jl | 7 +++++++ src/chaosdetection/lyapunovs/lyapunovspectrum.jl | 9 +++++++++ 2 files changed, 16 insertions(+) diff --git a/src/chaosdetection/lyapunovs/lyapunov.jl b/src/chaosdetection/lyapunovs/lyapunov.jl index f8a04c42..c0b17cb5 100644 --- a/src/chaosdetection/lyapunovs/lyapunov.jl +++ b/src/chaosdetection/lyapunovs/lyapunov.jl @@ -95,10 +95,14 @@ end inittest_default(D) = (state1, d0) -> state1 .+ d0/sqrt(D) function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt) + has_retcode = hasfield(typeof(pinteg),:sol) # transient t0 = pinteg.t while pinteg.t < t0 + Ttr step!(pinteg, Δt) + if has_retcode + pinteg.sol.retcode ==:Unstable && return NaN + end d = λdist(pinteg) lt ≤ d ≤ ut || rescale!(pinteg, d/d0) end @@ -113,6 +117,9 @@ function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt) # evolve until rescaling while lt ≤ d ≤ ut step!(pinteg, Δt) + if has_retcode + pinteg.sol.retcode ==:Unstable && return NaN + end d = λdist(pinteg) pinteg.t ≥ t0 + T && break end diff --git a/src/chaosdetection/lyapunovs/lyapunovspectrum.jl b/src/chaosdetection/lyapunovs/lyapunovspectrum.jl index 387139e0..7e3fe07b 100644 --- a/src/chaosdetection/lyapunovs/lyapunovspectrum.jl +++ b/src/chaosdetection/lyapunovs/lyapunovspectrum.jl @@ -81,6 +81,7 @@ function lyapunovspectrum(ds::DS{IIP, S, D}, N, Q0::AbstractMatrix; end function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = false) + has_retcode = hasfield(typeof(integ),:sol) if show_progress progress = ProgressMeter.Progress(N; desc = "Lyapunov Spectrum: ", dt = 1.0) end @@ -89,6 +90,10 @@ function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = t0 = integ.t while integ.t < t0 + Ttr step!(integ, Δt) + if has_retcode + integ.sol.retcode ==:Unstable && return fill(NaN,length(get_state(integ))) + end + Q, R = _buffered_qr(B, get_deviations(integ)) set_deviations!(integ, Q) end @@ -99,6 +104,10 @@ function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = t0 = integ.t for i in 1:N step!(integ, Δt) + if has_retcode + integ.sol.retcode ==:Unstable && return fill(NaN,length(get_state(integ))) + end + Q, R = _buffered_qr(B, get_deviations(integ)) for j in 1:k @inbounds λ[j] += log(abs(R[j,j])) From 53e2e1ac86bddfa4e72348ddef99578607e05288 Mon Sep 17 00:00:00 2001 From: rusandris Date: Wed, 13 Jul 2022 11:41:06 +0300 Subject: [PATCH 2/4] use successful_step functions for error checking --- src/chaosdetection/lyapunovs/lyapunov.jl | 9 ++------- src/chaosdetection/lyapunovs/lyapunovspectrum.jl | 10 +++------- 2 files changed, 5 insertions(+), 14 deletions(-) diff --git a/src/chaosdetection/lyapunovs/lyapunov.jl b/src/chaosdetection/lyapunovs/lyapunov.jl index c0b17cb5..f8f6e3b5 100644 --- a/src/chaosdetection/lyapunovs/lyapunov.jl +++ b/src/chaosdetection/lyapunovs/lyapunov.jl @@ -95,14 +95,11 @@ end inittest_default(D) = (state1, d0) -> state1 .+ d0/sqrt(D) function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt) - has_retcode = hasfield(typeof(pinteg),:sol) # transient t0 = pinteg.t while pinteg.t < t0 + Ttr step!(pinteg, Δt) - if has_retcode - pinteg.sol.retcode ==:Unstable && return NaN - end + successful_step(pinteg) || return NaN d = λdist(pinteg) lt ≤ d ≤ ut || rescale!(pinteg, d/d0) end @@ -117,9 +114,7 @@ function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt) # evolve until rescaling while lt ≤ d ≤ ut step!(pinteg, Δt) - if has_retcode - pinteg.sol.retcode ==:Unstable && return NaN - end + successful_step(pinteg) || return NaN d = λdist(pinteg) pinteg.t ≥ t0 + T && break end diff --git a/src/chaosdetection/lyapunovs/lyapunovspectrum.jl b/src/chaosdetection/lyapunovs/lyapunovspectrum.jl index 7e3fe07b..920ccc3c 100644 --- a/src/chaosdetection/lyapunovs/lyapunovspectrum.jl +++ b/src/chaosdetection/lyapunovs/lyapunovspectrum.jl @@ -81,7 +81,7 @@ function lyapunovspectrum(ds::DS{IIP, S, D}, N, Q0::AbstractMatrix; end function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = false) - has_retcode = hasfield(typeof(integ),:sol) + if show_progress progress = ProgressMeter.Progress(N; desc = "Lyapunov Spectrum: ", dt = 1.0) end @@ -90,9 +90,7 @@ function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = t0 = integ.t while integ.t < t0 + Ttr step!(integ, Δt) - if has_retcode - integ.sol.retcode ==:Unstable && return fill(NaN,length(get_state(integ))) - end + successful_step(integ) || return fill(NaN,length(get_state(integ))) Q, R = _buffered_qr(B, get_deviations(integ)) set_deviations!(integ, Q) @@ -104,9 +102,7 @@ function lyapunovspectrum(integ, N, Δt::Real, Ttr::Real = 0.0, show_progress = t0 = integ.t for i in 1:N step!(integ, Δt) - if has_retcode - integ.sol.retcode ==:Unstable && return fill(NaN,length(get_state(integ))) - end + successful_step(integ) || return fill(NaN,length(get_state(integ))) Q, R = _buffered_qr(B, get_deviations(integ)) for j in 1:k From 59354435141fb2efdcc61ca451b0e4e847e77ba6 Mon Sep 17 00:00:00 2001 From: rusandris Date: Thu, 11 Aug 2022 12:47:24 +0300 Subject: [PATCH 3/4] fix silent inf loop in lyapunov and add tests --- test/chaosdetection/lyapunov_exponents.jl | 56 +++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/test/chaosdetection/lyapunov_exponents.jl b/test/chaosdetection/lyapunov_exponents.jl index a70eb506..5874d6e4 100644 --- a/test/chaosdetection/lyapunov_exponents.jl +++ b/test/chaosdetection/lyapunov_exponents.jl @@ -150,4 +150,60 @@ end end +#test lyapunov and lyapunovspectrum for unstable and stable ds + +@testset "lyapunov tests for unstable systems" begin + + + ############## define unstable dynamical systems ################## + function divergent_chaotictest!(du,u, p, t) + k1, k2 = p + x, y, z = u + du[1] = y + du[2] = z + du[3] = k1-7*y+x^2+k2*z + end + + u0 = rand(3) + para = [2, -0.9] #k2=k6=0 + + ds_cont_unstable = ContinuousDynamicalSystem(divergent_chaotictest!, u0, para) + ds_disc_unstable = DiscreteDynamicalSystem(divergent_chaotictest!, u0, para) + + diffeq = (alg=Tsit5(),) + + T = 50 + + @test isnan(lyapunov(ds_disc_unstable,T)) + @test all(isnan,lyapunovspectrum(ds_disc_unstable,T)) + + @test_logs (:warn, "Instability detected. Aborting") isnan(lyapunov(ds_cont_unstable,T;diffeq)) + @test_logs (:warn, "Instability detected. Aborting") all(isnan,lyapunovspectrum(ds_cont_unstable,T;diffeq)) + +end + + +@testset "lyapunov tests for stable systems" begin + + ################### define stable dynamical systems #################### + + ds_cont_stable = Systems.lorenz() + ds_disc_stable = Systems.henon() + + + T = 50 + diffeq = (alg=Tsit5(),) + + @test lyapunov(ds_cont_stable,T) != NaN + @test lyapunov(ds_disc_stable,T) != NaN + @test all(isfinite,lyapunovspectrum(ds_cont_stable,T)) + @test all(isfinite,lyapunovspectrum(ds_disc_stable,T)) + + @test lyapunov(ds_cont_stable,T;diffeq) != NaN + @test lyapunov(ds_disc_stable,T;diffeq) != NaN + @test all(isfinite,lyapunovspectrum(ds_cont_stable,T;diffeq)) + @test all(isfinite,lyapunovspectrum(ds_disc_stable,T;diffeq)) + +end + end From 84abed260a1a4cc191ced65ebd398eac798f512a Mon Sep 17 00:00:00 2001 From: rusandris Date: Thu, 11 Aug 2022 12:50:38 +0300 Subject: [PATCH 4/4] fix silent inf loop in lyapunov and add tests --- src/chaosdetection/lyapunovs/lyapunov.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/chaosdetection/lyapunovs/lyapunov.jl b/src/chaosdetection/lyapunovs/lyapunov.jl index f8f6e3b5..cd8e89d3 100644 --- a/src/chaosdetection/lyapunovs/lyapunov.jl +++ b/src/chaosdetection/lyapunovs/lyapunov.jl @@ -119,6 +119,7 @@ function lyapunov(pinteg, T, Ttr, Δt, d0, ut, lt) pinteg.t ≥ t0 + T && break end # local lyapunov exponent is the relative distance of the trajectories + isnan(d) && return NaN a = d/d0 λ += log(a) rescale!(pinteg, a)