# Handling unstable true solution issues with the Integrator Interface

**URL:** <https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [June 11, 2020, 3:58pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208 "2020-06-11T15:58:35Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [June 11, 2020, 3:58pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/1 "2020-06-11T15:58:36Z")

</div>

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:

```julia
┌ 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:

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 11, 2020, 5:24pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/2 "2020-06-11T17:24:16Z")

</div>

Note:

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

```

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

So what you’re looking for is:

```julia
step!(integrator, 730.0, stop_at_tstop=true)

```

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [June 11, 2020, 5:49pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/3 "2020-06-11T17:49:00Z")

</div>

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](https://github.com/SciML/DiffEqBase.jl/blob/803298e174e05c3ab5a139a4738670a3c73fb176/src/integrator_interface.jl#L565), not a keyword so I’ve just used `step!(integrator, 365.0, true)` in my examples above.

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [June 11, 2020, 6:02pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/4 "2020-06-11T18:02:35Z")

</div>

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:

```julia
...
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:

```julia
...
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 11, 2020, 6:05pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/5 "2020-06-11T18:05:22Z")

</div>

> [@Libbum](#):
>
> 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?

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.

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [June 11, 2020, 6:07pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/6 "2020-06-11T18:07:34Z")

</div>

Great. That clears it up. Thanks!

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [June 15, 2020, 9:27am UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/7 "2020-06-15T09:27:15Z")

</div>

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

In regards to mutating the value of `u` in the integrator, [the documentation](https://docs.sciml.ai/stable/basics/integrator/#Note-about-mutating-1) 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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 15, 2020, 12:06pm UTC](https://discourse.julialang.org/t/handling-unstable-true-solution-issues-with-the-integrator-interface/41208/8 "2020-06-15T12:06:19Z")

</div>

> [@Libbum](#):
>
> Since I am saving at specific, discrete time points - is it OK if modify `u` directly at these points before `step!` -ing further?

Yeah, that’s exactly what DiscreteCallback is doing.