DifferentialEquations.jl - Solving ODE with "alg_hints=[:stiff]" is 20x slower, any alternatives?

General Usage
, , ,
1

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?

2

You can find recommended solvers here in the documentation

It is common that alg_hints is slower, as it might try different approaches at first. Which is especially relevant if the ODE solution itself is fast.

3 Likes
3

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.

4

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).

Warning: Instability detected. Aborting

5

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?

2 Likes
6

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

1 Like
7

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.

8

Speaking from my own experience, the most likely reason could be a small implementation mistake. Especially since it’s running in R.

If you could share the code or a MWE, people here can probably help you finding the discrepancy faster.

1 Like
9

It’s the exact same solver. If the steps aren’t exactly the same, then there’s an implementation error. Did you try deSolveDiffEq.jl to double check?