Working on a problem that has lead me to the Integrator Interface of DifferentialEquations.jl and have ran into something I can’t seem to work out.
My case requires my to solve my
ODEProblem on a timescale of days, but alter the larger system yearly.
In a normal run of my system I can use
PeriodicCallback(action!, 365) over my known timespan to invoke such changes. This works without issue.
Now, I need to solve my system continuously (by that I mean in a while loop essentially), therefore I’m not aware of
tspan[end], and wish to use
step!(integrator, 365.0, true) to advance the ODE system. If I call this a couple of times however, I run into:
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable. └ @ DiffEqBase ~/.julia/packages/DiffEqBase/TjqaN/src/integrator_interface.jl:343
Initially, I thought this may be because I wasn’t defining a jacobian for my ODEs, so I’ve used ModelingToolkit.jl to do that for me, but that wasn’t the problem it seems. My system is below for you to take a look at:
using OrdinaryDiffEq using ModelingToolkit @parameters t nutrients ib ip r H₁ H₂ H₃ H₄ cb cp prmax ce mp rv cv mv @variables B(t) P(t) V(t) @derivatives D' ~ t FR = B^2 / (B^2 + H₄^2) eqs = [ D(B) ~ ib + r * (nutrients / (nutrients + H₁)) * B - cb * B^2 - prmax * FR * P, D(P) ~ ip + ce * prmax * FR * P * (V / (V + H₂)) - mp * P - cp * P^2, D(V) ~ rv * V - cv * V^2 - mv * (V * B^2 / (H₃^2 + B^2)), ] sys = ODESystem(eqs) u0 = [ B => 20.0 P => 1.8 V => 50.0 ] p = [ nutrients => 0.7 ib => 2e-5 ip => 2e-5 r => 7.5e-3 H₁ => 0.5 H₂ => 20 H₃ => 20 H₄ => 15 cb => 7.5e-5 cp => 2.75e-4 prmax => 5e-2 ce => 0.14 mp => 2.25e-3 rv => 7e-3 cv => 6e-5 mv => 7e-3 ] tspan = (0.0, 365.0) prob = ODEProblem(sys, u0, tspan, p; jac = true, sparse = true) integrator = init(prob, Tsit5())
Tsit5() handles this system without issue in my usual case with discrete, yearly callbacks.
So, this particular case will step to
t=365.0 without issue, but the second step lands at
t=796.9, whereas I’d like it to land at
I’ve tried to manually
add_tstop!(integrator, 730.0), but I suspect this does the same thing that setting the third argument of
The only workaround I can think of is let the solver step to where it needs, then use the interpolated results in my external functions after the
Are there any other options I could try?