Adaptive SDESolver yields "Warning: dt <= dtmin. Aborting", unsure why, large non-adaptive dt works fine

Specific Domains Modelling & Simulations
, ,
1

I have a simple model, which I am trying to simulate, however, I am getting

┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase /home/SLCU/torkel.loman/.julia/packages/SciMLBase/GW7GW/src/integrator_interface.jl:345

errors. I am quite sure I run the same model like half a year/a year ago sometime, and it was all fine (but I have been unable to dig up a working example from then). Either way, if I try a large fixed dt it seems fine. Generally, I haven’t had problems with this combination of model/parameters, so am unsure what is going on.

MWE:

using StochasticDiffEq

function f(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1] = v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1) - σ
    du[2] = (1/τ)*(σ-A)
end
function g(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1,1] = η*sqrt(v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1))
    du[1,2] = -η*sqrt(σ)
    du[1,3] = 0
    du[1,4] = 0
    du[2,1] = 0
    du[2,2] = 0
    du[2,3] = η*sqrt((1/τ)*σ)
    du[2,4] = -η*sqrt((1/τ)*A)
end;

sprob = SDEProblem(f,g,[1.,1.],(0.,2000.),[5.0, 5.0, 0.1, 0.1, 3.0, 0.1],noise_rate_prototype=zeros(2,4))
sol = solve(sprob,ImplicitEM());

the solver fails pretty much independent of the seed. If I plot it:

using Plots
plot(sol)

image

the output looks normal, and nothing special happens when it crashes. I’m not sure if it is useful, but I can plot the dt over time:

plot(sol.t[2:end],sol.t[2:end]-sol.t[1:end-1])

image

Finally, if I try fixing the stepsize, it works well, even for rather large dts:

sprob = SDEProblem(f,g,[1.,1.],(0.,2000.),[5.0, 5.0, 0.1, 0.1, 3.0, 0.1],noise_rate_prototype=zeros(2,4))
sol = solve(sprob,ImplicitEM();adaptive=false,dt=0.25);

Finally, this output looks like:

using Plots
plot(sol)

image

The output of

using Pkg
Pkg.status()

is

      Status `~/Desktop/temporary_julia_environment/Project.toml`
  [91a5bcdd] Plots v1.27.4
  [789caeaf] StochasticDiffEq v6.46.0

and the Julia version 1.7.0.

Just clarifying the actual question: Why does the adaptive time-stepper failing, and is there some way to make it able to succesfully simulate this model?

2

I am attempting to circumvent this by using force_dtmin=true, but am unsure whenever I am doing something potentially unwise…

3

While it can solve with a higher dt, is the error below the adaptive tolerance in a strong sense? Most likely it’s just due to not being able to hit the required tolerance.

1 Like
4

You are right, it seems to have to do with the error. If I add abstol=1e-1 to the adaptive time-stepper it works without errors/warnings.

Just checking so that I know what is (probably going on): Basically, the solution is noisy. Even if the adaptive stepper tries to reduce the dt to keep the errors within the limits, the fluctuations are too much and it eventually fails. However, if I increase the error tolerance then the threshold is passed.

In the non-adaptive time-stepper, there’s a fixed dt, so it never actually checks the error, and this is never a problem. If I set my fixed dt too high, I get a

┌ Warning: Newton steps could not converge and algorithm is not adaptive. Use a lower dt.
└ @ SciMLBase /home/SLCU/torkel.loman/.julia/packages/SciMLBase/GW7GW/src/integrator_interface.jl:357

error, but that is a different phenomenon(?)

This is useful, I have never really paid much attention to the error (tolerance), but maybe this illustrates how it becomes important?

Thansk a lot!

5

that’s different, that’s the solution diverging.

The error is measured in strong error as well, so it’s not necessarily the error of the average.

1 Like
6

Thanks, I think I get what is going on, this has been really useful :+1: