Strange behavior with DifferentialEquations.jl

General Usage
1

@ChrisRackauckas, I know you are busy today, but for afterwards, I have what “seems” to be a bug in DifferentialEquations.jl. Here is a MWE:

using DifferentialEquations

## Data generation
function lotka!(du, u, p, t)
    α, β, γ, δ = p
    du[1] = α*u[1] - β*u[2]*u[1]
    du[2] = γ*u[1]*u[2]  - δ*u[2]
    return du
end

# Define the experimental parameter
tspan = (0.0f0,8.9f0)

u0 = Float32[0.44249296,4.6280594]
p_ = Float32[1.3, 0.9, 0.8, 1.8]
prob = ODEProblem(lotka!, u0,tspan, p_)
solution = solve(prob, Tsit5(), abstol=1e-12, reltol=1e-12, saveat = 0.1)

# Ideal data
X = Array(solution)
t = solution.t

When I execute these lines, I get that t.solution is a 68-element Vector{Float32}. But I expected a 90-element vector with times [0., 0.1, 0.2, ....., 6.9].
I am on a Mac with Ventura and an M1 chip. Here Project.toml:

[deps]
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"

Finally, the results of versioninfo():

Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.3.0)
  CPU: 8 × Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 2 on 4 virtual cores
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 2

Any insights are appreciated. I am perfectly happy to assume there is something I do not understand.

2

I believe this has to do with passing both abstol/reltol and saveat. If you fix the timesteps it can take, it doesn’t make sense to give it a tolerance to hit.

3

You are correct. That makes perfect sense. Thanks!

While this solved the problem, I think there is more to this than meets the eye. I thought that the integration proceeds with adaptive time steps to maintain the required tolerances, and saves at the specified values. In other words, lower tolerances would lead to saved values that are less precise. I feel that keep the tolerances in the call should not affect the number of times the solution is saved.
So I feel there is a need for a deeper explanation. Thanks!

4

saveat has nothing to do with time steps.

When you run the code you get a warning:

┌ Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).
└ @ SciMLBase C:\Users\accou\.julia\packages\SciMLBase\cxGMJ\src\integrator_interface.jl:504

so if you check solution.retcode it says:

julia> solution.retcode
ReturnCode.MaxIters = 4

so it exited early. The reason is because it couldn’t solve a Float32 problem to 12 digits: Float32 only has 8 digits of accuracy so it’s not really possible to have reltol that low (@Julius_Martensen make note :sweat_smile: )

The following is more reasonable and works without early exit:

solution = solve(prob, Tsit5(), abstol=1e-7, reltol=1e-7, saveat = 0.1)

# Ideal data
X = Array(solution)
t = solution.t
1 Like