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

1

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. However, CVODE does not seem to automatically switch between non-stiff and stiff.

So my question is,

  1. how do you switch between non-stiff solver and stiff ODE solver?
  2. If I want to solve from time =0 to time =10. If a non-stiff solver did not return error message from t=0 to 5, however when t>5 it reports error message. Can I trust the non-stiff solver’s results from t=0 to t=5? Or do I have to use a stiff solver to solve from t=0 to 10 again?

Thanks much in advance!

2

From my experience, the LSODA family of routines are the most robust ODE solvers you will find. It is no coincidence virtually all mathematical packages just copied them. They do a great job, both switching between non-stiff or stiff methods, and performance-wise. LSODAR and DLSODAR are the ones I use the most, because I usually need the extra root functionality (in the models I was working for years, the domain of integration is not usually known in advance, so the root finding was essential for me). However, if you don’t need the extra root functionality, LSODA or DLSODA should do.

ODEPACK is written in Fortran 77 and it is spaghetti-code (a “GO TO hell”) so diving into the code to figure out how exactly it does the switching between non-stiff and stiff methods is serious trouble (I tried years ago, and gave up with a headache). You might find more information about the switching in the papers the authors published:
[1] A. C. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers -
Scientific Computing, eds. R. S. Stepleman et al., North-Holland, Amsterdam,
p. 55 (1983).
[2] L. R. Petzold, Siam J. Sci. Stat. Comput. 4, 136 (1983).
As you can see it is old software, but still used today extensively. No need to re-invent the wheel, especially if that “wheel” works remarkably well.

Now, the following is not exactly related to your question, but it might be useful if you decide to use ODEPACK: it is not the easiest software ever… That’s because it has external dependencies (which are spaghetti code as well), and because the ODE solvers themselves expect lots of mandatory arguments you need to provide. This is why years ago I wrote a Fortran 95 module, namely DLSODAR_F95, which facilitate things. The package contains everything you need to use DLSODAR in an easy way: All dependencies are included, so you won’t need to pick the specific SLATEC routines needed. The solver itself is replaced by ODE_DLSODAR, which acts as an interface with the original code, so you won’t need to deal with ODEPACK’s peculiarities. Instead, ODE_DLSODAR needs only the absolutely necessary arguments; it sets up the rest automatically, then calls the original DLSODAR to do the job. You still have full control over the original subroutine via ODE_DLSODAR’s optional arguments - but I doubt you will ever need this functionality. The package also includes helper subroutines to facilitate output of the the solver (replaces the awful error reporting of the original, and gives you more control over the results). You basically include the module DLSODAR_F95 in your code and you are ready to go (an example case showing how to do that is included in the package). You should have no issues figuring out how the software works because each subroutine is documented in the code itself, explaining everything the programmer needs to know.
Now, I realize all this is about DLSODAR, which is for double precision and includes root finding. This is the solver I use all the time, that’s why the package is about it. If you want to use the single precision variant and/or without the root finding functionality, you can basically modify the code (it should be easy), or I can do that, if needed.

2 Likes
3

Thank you very much @Pap ! Welcome to the Disclose!
I just realized that you are the author of the DLSODAR_F95 which is so cool! It is actually your words that

DLSODAR_F95 is a convenient Fortran 95 interface for ODEPACK’s most sophisticated solver, DLSODAR (by A. C. Hindmarsh and L. R. Petzold.) DLSODAR is a very powerful multi-step solver for initial-value problems, switching between non-stiff and stiff integration methods automatically. …

That somewhat convinced me to try LSODA family. It does seem to be very good solver. Its speed is close to some modern Fortran ode solver. It can automatically switch between non-stiff and stiff is particularly cool! It can also do interpolation as dvode and rkc solver, which is also cool!
I just have no idea why in cvode or dvode they seems deleted this nice non-stiff and stiff switch ability :sweat_smile:.

4

They say one picture is worth 1000 words. The pictures below show the solution of the Shampine’s ball of flame problem, for δ=0.000001. This problem looks innocently easy, but it is actually very hard for an algortithm that cannot switch to a stiff solver when needed. See this paper for details (section 3.2.1), or just look it up.

The first figure shows the solution obtained by a fairly robust adaptive-size method (specifically Runge-Kutta-Fehlberg/4th-5th order) - but without a mechanism to switch to a stiff solver; as you can see, the method performs well until the stiffness starts; then it starts to take more and more integration points, and falls the the well-know stiffness trap: The algorithm is unable to return to a normal step size, even when the stiffness is over, so it keeps taking tiny steps, even when it’s not necessary anymore. In this particular case, it solves the problem with 51200 integration steps, 99.7% of which are in the region x>107985, where the solution is basically constant, so the vast majority of the computations are totally useless.

The second figure shows the result obtained with DLSODAR. Notice how the algorithm is able to return to a normal step size when the stiffness is over. In this case, DLSODAR needs only 224 integration points. Why? because it switches to a stiff solver when needed but it is able to switch back to a normal solver when the stiff part is over.

The pdf link above contains several other examples, more complex than this one. You will probably conclude you can’t go wrong with LSODA, LSODAR, DLSODA, DLSODAR. They are the state-of-the-art in the topic. :slight_smile:

Stiff_1b
Example_7b

3 Likes
5

Damn! That is decent! Thank you @Pap !

I see, in this example you mentioned, I know that with ITASK = 1 and I set the initial x = 0, final x = 20000. So LSODA solve the whole range and automatically switch between non-stiff and stiff during the process.

Now, just one question.
If I use the interpolation function DINTDY in LSODA to solve for the whole x range you mentioned in your above example. Can LSODA still automatically switch between non-stiff and stiff during the process?

I mean like, I set the initial x = 0, final x = 20000. I set ITASK = 2 in LSODA, so it solve one-step at a time, and I can use DINTDY to do interpolation at the points I want from 0 to 20000. Say I interpolate at x = 0, 100, 200, 300, 400, …, 19900, 20000. So you know what I mean. I am just saying with interpolation enabled, can LSODA automatically switch between non-stiff and stiff during the process? Thanks you very much!

PS.
By the way, about interpolation function, @ivanpribec has a nice example about dvode and perhaps RKC solver too,

The mechanism of interpolation for different solvers are more less the same.

6

DINTDY is used internally by the algorithm, but it is also available to the programmer for interpolation. In the latter case, which is the one you are interested, it is used after each call of LSODA, and it is not part of the integration. Therefore it is not involved in the stiff/non-stiff automation.

You can call DINTDY after each integration step to get an interpolated solution for any desired points found in this step interval (that’s what @ivanpribec does in his example code). However you can also just save the solution for each point selected by LSODA, then use this information to call DINTDY separately, or call another interpolation subroutine of your choice. Both methods are equally correct.

For a generic package like my DLSODAR_F95, which is not designed for a specific problem but just for any problem, the user can save the integration points and the solution there or not; no interpolation is done in the main integration loop (although it could be added). In most cases, I need the solution at the integration points, so I save it in a matrix anyway. Having the solution at the points the algorithm selected guarantees a smooth plot while having the solution only at arbitrary user-selected points does not.

Once I have the solution at the integration points DINTDY picked, I can do whatever I want. If I need an estimated interpolated value at a given point t, I can call DINTDY accordingly, or I can use another interpolation method of my choice. However if you don’t care what points LSODA picked, and only care about the solution at specific points of your interest, @ivanpribec’s method is definitely more practical.

1 Like
7

I’m not sure if you can call DVINDY separately. For dense output (interpolation) the routine DVINDY will typically use a form of Hermite interpolation that requires also derivative information (obtained by evaluating the right-hand side callback). This will typically be stored in the work array (e.g. RWORK) for the last step or so. That’s the reason why rwork appears in the call:

call dvindy (sgrid(n), 0, rwork(21), neq, ygrid(:,n), iflag)

I guess it’s possible to evaluate the derivatives later, after the solution has already been obtained, and pack them into an array which matches the expectations of DVINDY, but it kind of defeats the purpose of dense output. If you know you want the solution at a certain integration point, you should just communicate it directly to the integration routine beforehand, instead of calling the interpolation routine yourself.

A convenient way to obtain a reasonably-good interpolant only from the integration points points is to use cubic spline interpolation. However, this doesn’t come with the same accuracy guarantees as the “built-in” interpolation mechanisms for dense output.

1 Like
8

@ivanpribec I see your point. By “saving the solution for later use” I meant saving the derivatives as well - all the derivatives that are involved in the differential equations, that is (essentially the part of RWORK that contains this information in each step). That’s what I always do - even if I only need the function(s) themselves. I am so used to it that when I say “solution” I assume everybody will interpret is as “the derivatives are part of it”. Having the derivatives at hand is always a good thing - for example it helps interpolating with cubic splines or whatever interpolating method needs them. You don’t always know where you want to interpolate beforehand, that’s why I need a way to interpolate after integration is done.

Perhaps I should add another optional argument in my ODE_DLSODAR to let the user decide if he/she wants to get another output containing the solution at a user-defined grid. However I considered this a dangerous thing, because the user may choose to get just the solution at the grid but not at the points the algorithm chooses - and the grid the user picks may be dense or not. If the grid doesn’t contain at least one point in each integration interval the algorithm picked, the resulting plots may not show what the solution really is. And people love to pick the wrong grid. :laughing:.

1 Like
9

Thank you very much @Pap , and also @ivanpribec !
Uhm, let me confirm. Let us say LSODA and your auto non-stiff/stiff switch example.

  1. Just say we will solve the total range of x which is [0, 20000]. I only want to know the solution at x= 20000. There are two choices of ITASK. ITASK = 1 and ITASK = 2. The former (ITASK = 1) solves the whole range automatically and give the solution at x=20000. The latter (ITASK=2) is basically solve one adaptive step at a time.
    Now let us forget about interpolation.
    So in the case of ITASK =2, I can just let LSODA give an output at each adaptive step (and do nothing else),
    therefore finally at x = 20000, both ITASK = 1 and ITASK = 2 should give the same result, is it right? I mean so the only difference between ITASK = 1 and 2 is simply that ITASK = 2 clearly show the solution at each internal adaptive step.

  2. Or ITASK = 2 may give very different result from ITASK =1 or even stuck by the stiffness trap?
    Because if ITASK = 2 only solve one adaptive step at a time, it does not ‘see’ the whole range of x, therefore LSODA in this case cannot auto switch between non-stiff and stiff?

In real case, for example, I want to know the solution at x=0, 100, 200, 300, …19900, 20000.
One way to do it, is to use ITASK = 1, so I call LSODA from 0 to 100, 100 to 200, 200 to 300, …, 19900 to 20000. You know, piece-wisely. So I have the solution at the points I want.
However it may be slower than doing interpolation. So I actually wanted to do interpolation. I set the whole range as 0 to 20000, and I set ITASK = 2. Therefore I can also get the (interpolated) solution at x=0, 100, 200, 300, …19900, 20000. But I just want to make sure that even if I use ITASK = 2, LSODA can still switch between non-stiff and stiff.
In my program, what I did is, I use ITASK = 2 and do interpolation. However if ISTATE return negative value, then I solve the problem piece-wisely simply using ITASK =1 for each piece.

But anyway, it seems the best way, perhaps as you have already pointed out. I just use ITASK = 1, however, we need LSODA somehow be able to save its solution at each internal adaptive step. In the end, we should be able to use those saved solution to do interpolation ourselves or use some interpolation subroutine.

10

Thanks @ivanpribec !
Uhm, may I ask, in your experience, how do you switch between non-stiff and stiff ODE solver?

I mean for example, what I do is, say I want to solve the range x from 0 to 10. I begin with non-stiff solver, however if x=5 it gives some error information, then I use stiff solver to solve x from 0 to 10.

Now the question is, if the non-stiff solver give error information at x=5, can I trust its result from x=0 to 5?
Or should I totally discard its results, and use stiff solver to solve the whole range which is 0 to 10 again?

11

It is a very interesting paper of yours. It gaves a way to decide id a set of ODEs is stiff (by comparing eigenvalues of Jacobian matrix). I wonder how can one judge in case of a single equation. The algorithm in LSODA apparently knows it somehow.

12

This is was a quite reasonable question - and yes, of course ITASK=2 will do what you expect it to do. It is common practice for pretty much every repetitive method to take care of that. After one step is taken all necessary information needed for the next step are there. This is why LSODA includes arrays like RWORK and IWORK.

You could integrate piece-wise as you described, with ITASK=1 for each segment. But I don’t see a real reason to do that. This is usually something we do when we have to (for example when there is an expected discontinuity inside the integration interval). Also, it’s not like you will avoid interpolation that way. Don’t forget that LSODA does interpolate anyway; INTDY is not just a subroutine for user’s convenience. It is also used internally in LSODA. So whatever you do, interpolation is there.

If you want the solution at a set of points like x=0, 100, 200, … and you know that set of points beforehand, go with ITASK=2 and interpolate to get the solution at the points you want, as in @ivanpribec’s example (it is for DVODE, but the principle is the same). In my case, I almost never have that information beforehand. So I save the solution in each step (again, with ITASK=2 but quite often I use ITASK=5 as well). After integration is done I create one or more interpolation functions, so I can have an interpolated value whenever I need it. I used such interpolating functions as an input for other numerical methods as well - for example, to compute eigenvalues.

You mentioned ISTATE errors, so I must emphasize you should not just bypass the problematic point blindly. You should investigate what happened first, and why. Plot the solution before ISTATE=-1. Is the solution correct at least there? If not, check your program and make sure the right-hand side of the differential equation(s) (and the Jacobian) was typed correctly. Even the slightest mistake there means you are solving a problem different than the one you wanted to solve. And it is very easy to make such a mistake, if the expressions are complex. In the vast majority of the cases I had issues, this was the reason.
If the solution is correct before ISTATE=-1, investigate what is the behavior of the function(s) near the point the problem occurred. Any discontinuity there, maybe? If so, try LSODAR with a constraint function to stop integration before the problematic point, or do as you said, go with ITASK=1 until a point close (but before) the problem. Then continue after that point, again with ITASK=1 or whatever.

What exactly needs to be done to treat the problematic point depends on the case. One thing is sure, you should not bypass the problem before making sure you understand what caused it.

1 Like