JuliaDynamics / InteractiveDynamics.jl

Fast, general-purpose interactive applications for complex systems
https://juliadynamics.github.io/InteractiveDynamics.jl/dev/
MIT License
171 stars 25 forks source link

Problem when using Rodas5 as the differential equation solver #77

Open Lucalino opened 2 years ago

Lucalino commented 2 years ago

Hi, I encountered a MethodError, when using Rodas5 to compute the timeseries of a 3D dynamical system. It works fine with other solvers such as Vern9. The example code reads

using GLMakie
using OrdinaryDiffEq
using DynamicalSystems
using InteractiveDynamics

function my_system(x, p, t)
    @unpack a, b, c = p
    dx = exp(x[2] - 0.5) * (1.0 - x[1]) - a * tanh(10.0 * (x[1] - 0.5)) - a
    dy = a * tanh(10.0 * (x[1] - 0.5)) + a - exp(x[2] - 0.5) + (x[3] - x[2]) * 2.0
    dz = 3.0 - exp(x[3] - 0.5) - (x[3] - x[2]) * 0.5
    return SVector(dx, dy, dz)
end

p = Dict(
       :a => 2.0,
       :b => 70.0,
       :c => 4.0e8,
)

ps = Dict(
    :a => 1.9:0.01:2.1
)

x0 = [0.3, 0.5, 0.4]
ds = ContinuousDynamicalSystem(my_system, x0, p)

diffeq = (alg = Rodas5(), adaptive = false, dt = 0.001, reltol = 1e-8, abstol = 1e-8)
tr = trajectory(ds, 100; diffeq...)

u0s =  [ds.u0]

fig, obs = interactive_evolution_timeseries(
    ds, u0s, ps; tail = 10000, diffeq, idxs = (1, 2, 3), 
    lims =((0.0,1.0), (0.0, 2.0), (0.0, 2.0))
)

The error message reads:

ERROR: MethodError: Cannot `convert` an object of type Bool to an object of type SVector{3, Float64}
Closest candidates are:
  convert(::Type{T}, ::Intervals.AnchoredInterval{P, T, L, R} where {L<:Intervals.Bounded, R<:Intervals.Bounded}) where {P, T} at /Users/luca/.julia/packages/Intervals/ua9cq/src/anchoredinterval.jl:181
  convert(::Type{T}, ::Intervals.Interval{T, L, R} where {L<:Intervals.Bound, R<:Intervals.Bound}) where T at /Users/luca/.julia/packages/Intervals/ua9cq/src/interval.jl:253
  convert(::Type{SVector{N, T}}, ::CartesianIndex{N}) where {N, T} at /Users/luca/.julia/packages/StaticArrays/OWJK7/src/SVector.jl:46

I don't think, DynamicalSystems.jl is the problem because the trajectory (tr) is computed correctly. The programme only fails when interactive_evolution_timeseries() is used.

I can provide the stacktrace on request.

Datseris commented 2 years ago

quick comment on MWE: When making a MWE try to minimize everything else that does not have a connection with the code. E.g., do not use unecessary packages like DataFrames and use the minimal packages that make the code run, e.g., replace DifferentialEquations with OrdinaryDiffEq and replace GLMakie with Makie.