Handling unstable true solution issues with the Integrator Interface

Specific Domains Numerics
1

Hi All,

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 t=730.0.

I’ve tried to manually add_tstop!(integrator, 730.0), but I suspect this does the same thing that setting the third argument of step! to true.

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 step! i.e. integrator.sol(730.0).

Are there any other options I could try?

2

Note:

step!(integrator,dt[,stop_at_tdt=false])

https://diffeq.sciml.ai/latest/basics/integrator/#Initialization-and-Stepping-1

So what you’re looking for is:

step!(integrator, 730.0, stop_at_tstop=true)
3

Hi Chris,

Yep, this is exactly what I’m using when I come up against these problems. Note that stop_at_tstop is a variable, not a keyword so I’ve just used step!(integrator, 365.0, true) in my examples above.

1 Like
4

OK, so I think I’ve identified the problem, which could ultimately be a bug in the sense that it may not be tested for.

Take the same setup as above and run it like so:

...
tspan = (0.0, 365.0)
prob = ODEProblem(sys, u0, tspan, p; jac = true, sparse = true)
integrator = init(prob, Tsit5())
step!(integrator, 365.0, true)
#integrator.t = 365.0
step!(integrator, 365.0, true)
#Warning: dt <= dtmin. Aborting.

On the other hand:

...
tspan = (0.0, 1.0)
prob = ODEProblem(sys, u0, tspan, p; jac = true, sparse = true)
integrator = init(prob, Tsit5())
step!(integrator, 365.0, true)
#integrator.t = 365.0
step!(integrator, 365.0, true)
#integrator.t = 730.0

So essentially I shouldn’t try to step when t=tspan[end]. Did I miss this point in the docs somewhere or should I file an issue for this?

5

That just looks like an untested edge case. solve! goes to tspan[end], and then step! could continue from there. Or step! could do all of it, but I think if you step! to the end then continue, it might hit a weird floating point bug since it’ll be 1e-16 away from the true end or something like that. Worth an issue.

FWIW, you should probably just be doing tspan = (0.0,Inf) here.

1 Like
6

Great. That clears it up. Thanks!

7

Hey @ChrisRackauckas, I have a quick follow-up question to this.

In regards to mutating the value of u in the integrator, the documentation suggests using callbacks. I’m using that method normally, but in this particular instance, when I’m only using the step! interface, I never call solve! and therefore cannot provide a callback (as far as I can see).

Modifications within a callback affect! surrounded by saves provides an error-free handling of the discontinuity.

(my emphasis). Since I am saving at specific, discrete time points - is it OK if modify u directly at these points before step!-ing further?

8

Yeah, that’s exactly what DiscreteCallback is doing.

1 Like