Debugging differential equation solver issues almost always boils down to doing the same thing, so this is a summary to help you out. For more information on specific issues, check out the FAQ section of the DifferentialEquations.jl documentation which highlights common issues and questions:
dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
Instability detected. Aborting
There, now this post will show up in a search. What do you do when you see these kinds of problems?
If you find yourself in differential equation solver trouble, the first thing to do is lower the tolerance. In many cases, decreasing the tolerance improves the stability, and so trying something like
abstol=1e-10,reltol=1e-10 can tell you if it’s a tolerance issue. If that’s not working, try switching to a more stable solver, something for stiff equations like
If none of those things are working, the problem is probably your model. If you want to blame Julia solvers, then there’s a good way to prove it: try a non-Julia solver. It’s as simple as changing
solve(prob,CVODE_BDF()) and now you’re using the classic Sundials C/C++ library. That’s not all, there’s a whole list:
- Sundials.jl, a wrapper for the C/C++ SUNDIALS library though
- ODEInterfaceDiffEq.jl, a wrapper for the classic Hairer Fortran codes like
- LSODA.jl, a wrapper for the classic
- MATLABDiffEq.jl, a wrapper for the MATLAB ODE solvers
- SciPyDiffEq.jl, a wrapper for SciPy’s
odeint(LSODA) and other methods (LSODE, etc.).
- deSolveDiffEq.jl, a wrapper for the commonly used R library.
Note: these solver packages are not included by default due to their extra dependencies, so you’lld need to do things like
]add Sundials; using Sundials before using the Sundials solvers
If you do this at lower tolerances and every solver fails… I’ll just put it in meme form.
So yes, if you’re model is failing with every major solver, including all of the major ones called into unmodified from C/C++ and Fortran, then the problem isn’t the solver it’s your model. I know you checked it, I know you think you debugged it, but the likelihood that every single solver ever created in all languages is incorrect instead of the code you wrote a few hours ago is very slim.
So what can it be? Here’s a list of common issues:
- Double check for sign issues. And then check again. If you have some solution piece that’s blowing up to infinity, why? Which term in its derivative gets really large? Is it supposed to be large? Why did it get large? Almost all issues boil down to just checking term by term in this fashion. Double check your model.
- Double check your assumptions on the model. If the derivative goes to zero as
ugoes to zero, that means that
ucan never go negative, right? Wrong!
u' = -sqrt(u)hits zero in finite time. Just because you believe that the system you’re modeling has a property (like being positive) doesn’t mean that the true solution to your model actually has this property. Try to prove the properties or, try and look at the solution as it begins to violate the property, look at the values in the derivative, and understand which terms are causing it.
- Double check if you’re doing anything that violates assumptions of being an ODE. As an ODE, the right-hand side
ffunction should always give you the same result:
u' = f(u,p,t)needs to be uniquely defined. These issues are just fundamental to the mathematics: if you do these things in
fno longer defines an ODE so of course it cannot be solved!
- If you put randomness into
f, then the adaptivity will think your ODE is being solved at a high error because the derivative keeps changing, and therefore it will fail and hit
dtmintrying to reduce the randomness to zero (if you do need randomness, use an SDE or RODE solver).
- If your
ffunction is modifying
u, then calling
fwith different stepsizes is not deterministic and the solver is likely to fail. If you need to do this, you should be using a callback.
- If your
ffunction is caching values from a previous step, just think about what this means. If you change
dt, then you’re changing
f. In that sense,
u'is no longer defined since it’s now dependent on how it’s being solved! Even worse, adaptive ODE solvers do not always move forwards in time: sometimes they try
t + dt1before trying
t + dt2and choosing what to do. So if you’re assuming that the cached values are from the last step, that’s not actually the case: those values can be coming from a fake (too incorrect according to error estimates) future!
- If you put randomness into
If you “know” your model can’t have an issue because it’s translated from MATLAB/Python/R, then try using MATLABDiffEq.jl/SciPyDiffEq.jl/deSolveDiffEq.jl’s wrapper solvers, which call the solvers of those respective languages/packages. If your differential equation is the same, then it should give exactly the same results, down to the step. Almost certainly if the Julia solvers are diverging these will diverge as well, and you can use this to help narrow down translation problems.
For understanding how to solve differential equations with the top-notch performance, follow the tutorials which mention how to get high performance:
If you’ve gone through these tutorials but still have performance issues, or just have clarifying questions to ask, please feel free to ask! But do check out these tutorials since they cover most of what people need!
There’s a large volume of differential equation questions pouring in because there’s a lot of adoption in undergraduate and graduate courses. This is a good thing! But it’s also time consuming and >99% of questions come down to the issues listed above and have a formulaic solution. If you’ve read this post, you now have enough information to lead someone through debugging model issues so please feel free to share this and walk someone through the step by step debugging process. If it turns out to be one of the few cases that need extra attention, then point people towards opening an issue.