I am solving several thousand systems of ODEs, with an explicit Jacobian, using DifferentialEquations.jl. These systems are similar to the competitive Lotka-Volterra model. In some cases (depending on the parameter values), Tsit5() method works well, but in some others the command solve() throws many “stiffness” and “Instability” warning messages.

sol = solve(odeprob, Tsit5())

Increasing maxiters and tolerance does not help.

I found that using “alg_hints=[:stiff]” produces a stable implementation, but it becomes 20 times slower.

sol = solve(odeprob; alg_hints=[:stiff])

Is there a faster alternative to using “alg_hints=[:stiff]” while still considering stiffness?

what happens if you just solve without passing an algorithm or stiffness hint? you probably want to be using one of the methods that automatically switches between stiff and non stiff solvers.

It continues throwing instability warning messages (and sometimes it enters into some sort of infinite loop, as it does not return a solution after a while).

A stiff ODE solver on non-stiff ODEs is slower. That’s why non-stiff ODE solvers exist. You need to use then when you can. Is the auto-switch method not working well here?

If there are several thousands of unknowns, stiff methods will probably need to invert the jacobian (or some matrix related to it), and for this size, it can be slow. You can use GMRES for the linear solver and see if it helps

Thanks, Chris. I agree. In my case, the non-stiff solver just enters in some sort of infinite loop (once it starts throwing warning messages, it never stops). The non-stiff solver works, although it is 10x slower (naturally). The auto-switch method does not work either, it behaves the same as a non-stiff solver (throwing infinite warning messages).

Intriguingly, LSODA.jl has the same problem (infinite warning messages), while LSODA works well in R (yet, 50-100 times slower).

I am still investigating what is going on. Stiffness is a potential problem, or some extreme parameter values that lead to overflows.