Has anyone used DVODE and is it good?

[Pap]

I had in mind the usual numerical methods that explicitly require Jacobians one way or another (analytically or numerically).

[oscardssmith]

Right, what I’m saying is that the fastest solvers now can numerically compute jacobians much more accurately using dual numbers (or reverse mode automatic differentiation), and that both of these approaches produce much higher quality results than the traditional method (finite difference).

[CRquantum]

Thanks @oscardssmith !
I see that FLINT and RKC does not explicitly require Jacobian, but they are not particularly suitable for stiff problem.
Uhm in your experience, how do you switch between non-stiff and stiff ode solver?

How to switch between non-stiff ODE solver and stiff ODE solver?

Dear all, I know many of you have very decent experience in Fortran ODE solvers. But most ODE solver are either only for non-stiff ODEs or just for stiff ODEs. I only see LSODA in ODEPACK actually can automatically switch between non-stiff and stiff methods. In general, I found that the performance of LSODA is actually not bad and it seems faster than DVODE. Also can see below, Choosing a Solver — Odes 2.6.2 documentation (scikits-odes.readthedocs.io) But perhaps CVODE may be even better. H…

Because in real case, even if we use the same model of ODE, if some parameters of the ODE changes, the problem may become a stiff problem. Therefore being able to efficiently switch between non-stiff and stiff can be very useful.

[ivanpribec]

@CRquantum

One more very naïve question @Pap , hmmm, it seems there are some ode solver does not require jacobian, does those solver reliable?
Like RKC, RKC: An explicit solver for parabolic PDEs - ScienceDirect
which can be found as rkc.f in ode

Also the classic textbook fourth-order Runge-Kutta discretization doesn’t require Jacobians. This doesn’t mean it’s unreliable per se. However if your problem is stiff, it can become unstable unless you use a very small time-step. If you take many steps, it means you will also accumulate more round-off error over time. Choosing the right solver is essentially a balancing act between accuracy, stability, and performance. Which solver to choose depends highly on your problem and resources (in terms of time, hardware, available libraries).

Perhaps you could give us a bit more detail, how large is your system of ODEs, how strongly are the equations coupled, where does the problem originate from, etc.?

The RKC code is an explicit method, which automatically adapts the number of stages and/or time-step to reach a certain level of accuracy/stability. As described by the authors it is designed for non-stiff and mildly stiff problems resulting from discretization of PDE’s using the method of lines. An example would be a 2D reaction-diffusion equation discretized using a 5-point stencil in space.

The Jacobian matrix is needed for implicit time-stepping methods. For non-linear problems, the implicit time discretization results in nonlinear systems of equations. The linearization of the nonlinear system around the current point in time is where the Jacobian appears from. The benefit of implicit methods is their stability for stiff problems. This comes at the expense of 1) having to calculate the Jacobian either analytically, via finite-differences, or auto-diff, and 2) solve a linear system of equations for each stage/time-step.

Depending on your problem, you could also take a look at implicit-explicit methods (IMEX). These use operator splitting and combine properties of both explicit and implicit stepping for different terms in the system of ODEs. An example is the IRKC code.

[CRquantum]

@ivanpribec

Perhaps you could give us a bit more detail, how large is your system of ODEs, how strongly are the equations coupled, where does the problem originate from, etc.?

Thank you very much @ivanpribec !

Sure, overall, the ode system I am dealing with is not big at all, like 3 to 20 ODEs. The problem is that it is a Monte Carlo Expectation Maximization algorithm. Just say in one iteration, need to solve such ODE system like 10^5 times. 1 time means, I need to solve, say, the 3-ODE system, at 24 different point, like at t = 0, 2, 4, 6, …,46, 48. So 24 observations (Of course I do interpolation at these points, instead of solve them piece wisely). Need to solve such 24-observation 3-ODE system 10^5 times in one iteration. So the speed of ODE solver is very importance.

The 3-ode model is like,

#Prim with range
Ka, 0, 15 # so the range of Ka is (0,15), etc.
Vmax0, 0, 60
Km, 0, 60
Vc0, 0, 5
FA1, 0.1, 1
KCP, 0, 15
KPC, 0, 15

#Cov
wt
age

#Sec
VM = Vmax0 * wt**0.75
V = Vc0 * wt

#Diffeq
XP(1) = -Ka*X(1)
XP(2) = Ka*X(1)-VM/(Km*V+X(2))*X(2)-KCP*X(2)+KPC*X(3)+90
XP(3) = KCP*X(2) - KPC*X(3)

Such ODE is interesting because if any of the primary parameter Ka, Vmax0, … KPC is negative number, the ODE will highly likely just fail.
I currently use FLINT + LSODA. Using FLINT first, if it fails I call LSODA. I sometimes use RKC + LSODA too. The speed is roughly like, FLINT + LSODA take 1.7s, RKC + LSODA take 1.9s, LOSDA take 2s. DVODE if in stiff mode, took 4s or so.