What is best ODE solver for stiff problem?

1

Dear all,

I know @ivanpribec has a stiff3 ODE solver for stiff problem,

Is there some other stiff ODE solver too?

@prince_mahajan point to Hairer solver too,

http://www.unige.ch/~hairer/software.html

How is stiff3 algorithm compare with Hairer method? In terms of accuracy and speed?

Thank you very much in advance!

1 Like
2

DIRK solvers are rather thoroughly investigated in the recent paper

Without any special needs, I found the suggested default solver ā€œESDIRK4(3)6L[2]SAā€ pretty good.
Speed greatly depends on how fast the underlying implicit equation can be solved (by Newton?) and how good the adaptive time step control can be tuned.

5 Likes
3

Thank you @martin !
Welcome to the community!

2 Likes
4

What is best ODE solver for stiff problem?

It depends on the problem. DifferentialEquations.jl has a great summary of recommended methods https://diffeq.sciml.ai/stable/solvers/ode_solve/#ode_solve

I do atmospheric chemical kinetics involving 10,000s ODEs, which are very stiff. For this problem, one of the best methods is CVODE BDF with the banded matrix solver. CVODE is written in C but has an excellent Fortran interface.

2 Likes
7

Have you tried SEULEX from Hairer? For my problem it seems to be 2-3 times faster than DVODE when its stiff and on par with it when its mildly or not stiff.

8

Iā€™ve had good results with ROCK2 from Hairer and RKC from Netlib. When they fail I tend to find MEBDF next best, followed by VODEPK. ā€˜Best codeā€™ is probably ā€˜best code for my problemā€™, so you may need to try a variety.

9

@bendelBoy Would it be possible for you to provide a simple driver for MEBDF? I downloaded one from here but cant seem to get it working properly.

10

Iā€™ve been using an older version, when the main routine was called STIFF.

Mine was taken from here ā€“ algorithm 703. [Collected Algorithms of the ACM]. That includes a driver.

My own routine isnā€™t suitable ā€“ built around my need, with extensions to the original. I have lots of parameters that I pass through for use with an unsteady-state process simulator.

case (8) ! MBEDF with full Jacobian (8)
mf = 22
LOUT = 10 ! For error messages
select case (m1(2))
case (1, -1); index = 1
case (2, -2); index = 2
end select

TempOutput = OutputTime
call DRIVER(ny, t0, m2(3), y, OutputTime, TempOutput,

  • m2(1), MF, INDEX, LOUT, iw1,
  • w(1), iw2 , w(iw1+1), MBND,
  • np, ns, ni, nx, np1, np2, npi,
  • o, z1, z2, s, p1, p2, pi, x, m1(5), m2(7), m3,
  • w(iw3), w(iw4), StartTime, StopTime)

m1(4) = index ! Preserve original error code

12

Thanx! The netlib version works.