PSA: How to help yourself debug differential equation solving issues

[ChrisRackauckas]

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:

https://diffeq.sciml.ai/dev/basics/faq/

How do I debug why the differential equation solver is diverging?

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 TRBDF2(), KenCarp4(), or QNDF().

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,Tsit5()) to 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 CVODE_Adams, CVODE_BDF, IDA, and ARKODE .
• ODEInterfaceDiffEq.jl, a wrapper for the classic Hairer Fortran codes like dorpi5, dop853, radau, rodas, etc.
• LSODA.jl, a wrapper for the classic lsoda algorithm.
• MATLABDiffEq.jl, a wrapper for the MATLAB ODE solvers ode45, ode15s, etc.
• 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 u goes to zero, that means that u can 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 f function 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 f, then f no 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 dtmin trying to reduce the randomness to zero (if you do need randomness, use an SDE or RODE solver).
  • If your f function is modifying u, then calling f with 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 f function 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 + dt1 before trying t + dt2 and 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 “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.

What about performance issues?

For understanding how to solve differential equations with the top-notch performance, follow the tutorials which mention how to get high performance:

https://diffeq.sciml.ai/dev/tutorials/advanced_ode_example/

https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html

https://mitmath.github.io/18337/lecture2/optimizing

Optimizing Serial Code in Julia 1: Memory Models, Mutation, and Vectorization

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!

Benevolent onlookers: how to help

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.