Dear all,
Has anyone used DVODE to solve ODEs?
https://www.radford.edu/~thompson/vodef90web/
What is the drawback of DVODE? Is it slow or inaccurate or something?
Thank you very much in advance!
2 Likes
Hey @CRQuantum, I find it slightly demeaning to tacitly assume the code has some drawbacks. After all, the solver was developed by accomplished mathematicians including Alan Hindmarsh, which authored a dozen other ODE solvers. Perhaps you can search for the paper which introduced the method to learn what the authors had to say about it, and then check the papers which cited it.
If the code was developed a few decades ago it likely doesn’t follow the latest programming models or benefit from the latest compiler performance tweaks. Without specifying what are the alternatives, we are left with an incomplete comparison.
Asking others about their experience with using a certain ODE solver code is of course more than welcome.
4 Likes
Besides what @ivanpribec said, I believe dvode is also at the heart of the SciPy’s ODE integrator: scipy.integrate.ode — SciPy v1.7.1 Manual. So I would assume it is rather very good.
We should create an fpm package out of it, add examples how to use it and maintain under fortran-lang.
2 Likes
DVODE, DLSODE, and related are good. They need to be modernized. Unfortunately, it was done, but only by converting it to C ( CVODE | Computing (llnl.gov)).
DDEABM is also similar: i have modernized that one: jacobwilliams/ddeabm: Modern Fortran implementation of the DDEABM Adams-Bashforth algorithm (github.com)
2 Likes
I’m sure DVODE is great. But, CVODE is probably the best maintained flavor of VODE, with a bunch of different matrix sparsity options. It’s written in C but has a good Fortran interface.
2 Likes
Thank you all! @ivanpribec @certik @jacobwilliams @nicholaswogan
Yes @jacobwilliams,
I actually have both @prince_mahajan 's FLINT and your DOP853 and DDEABM at hand. The structure of both of your code are similar, both are modernized. I learned a lot from your code!
Perhaps according to FLINT’s benchmark, FLINT’s DOP853 is slightly faster. But I have not done through test myself.
My plan is to put FLINT, your DOP853 and DDEABM inside a package I am developing (which needs solving ODE to make prediction), and also add at least an ODE solver for stiff problem. Currently I may just add the f90 version DVODE.
@jacobwilliams ,
May I ask, what do you mean by modernized, Uhm, does that mean making Dvode object-oriented?
Will that make DVODE faster?
Thank you very much!
@nicholaswogan
Thank you very much for your suggestion of CVODE, which is a part of Sundials. It should be great.
Uhm, I currently may just try to use a Fortran ODE instead. It is relatively easy for me to use Sundials for my own research. But I am developing a package, adding Sundials inside my package may be a little bit not not too easy right now.
By the way, has anyone by any chance tried both DVODE and CVODE, if so, is CVODE faster than DVODE for the same ODEs?
What we have done and agreed upon for other packages such as FFTPACK, modernize means: compile with fpm, use free form, possibly add into modules, ensure all interfaces are present, remove implicit typing, add tests, documentation, examples and so on.
1 Like
Time difference between CVODE and DVODE will depend on the problem but I’m guessing it will be small, like within a factor of 2. The exception is if you want to do some really big problem. Then CVODE will win/allow a solution to be calculated because it has stuff like GMRES.
2 Likes
@certik @jacobwilliams
Alright. Uhm, do you think the current DVODE by Byrne and Thompson is modernized? link is below,
File Listing (radford.edu)
I have downloaded the DVODE.f90 distribution file there,
uhm, it seem it is somehow modernized a little bit.
It uses module and f90 format, and implicit none.
In windows it looks like below, and it can compile and run without problem,
Seem not too bad.
@CRquantum im not able to access the code via the link your originally posted. Do you know an alternative link?
Thank you very much for the message @nicholaswogan ,
true, strangely, the original link does not work any more. Fortunately, I have downloaded their whole website (the latest version) of dvode before.
Now I pushed it to a new gitlab repo, you could download it from the link below,
In the folder VODE_F90 Support Page , double click index.html,
then you should see the whole original website, and you can download all the corresponding files locally.
Hope that helps 
2 Likes
I have not test the speed of dvode, but I guess it should not be too bad.
dvode can handle both stiff and non-stiff problem.
The only issue seems to be that dvode does not natively support interpolation.
Thanks!! This looks great. I think the code does have interpolation. Section 5 of the documentation, which is all in dvode_f90_m.f90, talked about interpolation.
1 Like
Great! Thank you so much for letting me know! I did not read the file carefully enough LOL.
Right, from line 1554 there are instructions about how to do interpolation using subroutine DVINDY.
Then the dvode seems a quite versatile ODEs solver!
Hi @nicholaswogan ,
Have you by any chance tried the interpolation function in DVODE?
It is the subroutine DVINDY.
It is a little strange that, I solve the ODE ( just say Y(t), I want Y at given t ) within the time t range, say from T0 to Tf.
Then I can call DVINDY to interpolate Y for t>Tf. (But this is strange, it is more like extrapolation than interpolation)
However , if I just want to interpolate Y at t within T0 to Tf, then DVINDY give me some error, something like t is not within the range like below,
Error in DVINDY, T(=R1) is illegal. T is not
In the above message, R1 = 0.6000000000000D+01
in interval TCUR - HU(= R1) to TCUR(=R2)
In the above message, R1 = 0.2364435159699D+02 R2 = 0.7499440289459D+02
Just curious, have you meet similar issues?
I think DVINDY can correctly do interpolation, but I just do not know which part did I do is wrong… LOL.
The interpolation routine DVINDY has the following requirement, which is stated in the source code:
For valid results, T must lie between TCUR - HU and TCUR.
In other words, if interpolation is desired, it should be performed after each successful integration step, which covers Tcur - Hu to Tcur (both evaluated with the current value of Tcur).
If, instead, you perform a number of successive integration steps and only then try to interpolate over the entire range of the independent variable, that will not be possible. The capability for doing that could have been provided, but that would have required a large amount of supporting information to be stored. At the time the software was published (1989), how much RAM did our PCs have?
You can verify for yourself the relation between T and Tcur by adding the following print statement at the beginning of subroutine DVINDY_CORE:
print '(A,3ES12.4)','T, TP, TN = ',T,TP,TN
and running the two example problems in your Gitlab repository.
1 Like
Thank you very much @mecej4 !
What I want is, call DVODE between time [0, 24], then interpolate points at [0,2,4,6,8,10,…22,24].
So, I mean if I call
CALL DVODE_F90(..., T=0,TOUT=24,...)
Where T=0 is the initial time, TOUT=24 is the final time.
So the integration is in the range of [0,24].
I would expect then I should be able to do things like
CALL DVINDY(Time,...)
where Time should be within [0,24], am I right?
But DVINDY does not seem to interpolate as I said above.
It seems after I do
CALL DVODE_F90(..., T,TOUT,...)
T becomes TOUT, then DVINDY can only interpolate at time between like [24-xxx, 24+yyy], such as [23, 75]. But Definitely not the original range [0,24] which is what I really want.
Overall, it seems DVINDY is not particularly useful. It can only interpolate at its own particular range TCUR - HU and TCUR. Instead of the original range [0,24].
I may just CALL DVODE_F90 between [0,2], then [2,4], then [4,6], … [22,24].
Are you saying adding the print you mentioned at the position below in dvode_f90_m.f90?
Line 6941 is where I added the print you mentioned.
But here TP seems has not been calculated yet?
Or, do I add the print at line 6947 below?
Do you know if there is quick way that I can add something in DVODE so that all those points are stored? Then I guess DVINDY will be more useful.
Edit: the follow-up post from @mecej4 addresses the actual problem of output at regular intervals, without the need to call DVINDY.
To use interpolation you should consult the instructions under itask in dvode (original or _f90):
! ITASK = An index specifying the task to be performed.
! Input only. ITASK has the following values and meanings.
! 1 means normal computation of output values of y(t) at
! t = TOUT (by overshooting and interpolating).
! 2 means take one step only and return.
! 3 means stop at the first internal mesh point at or
! beyond t = TOUT and return.
! 4 means normal computation of output values of y(t) at
! t = TOUT but without overshooting t = TCRIT.
! TCRIT must be input as RWORK(1). TCRIT may be equal to
! or beyond TOUT, but not behind it in the direction of
! integration. This option is useful if the problem
! has a singularity at or beyond t = TCRIT.
! 5 means take one step, without passing TCRIT, and return.
! TCRIT must be input as RWORK(1).
I assume you are currently using itask = 1. Probably you’d want to use 2 (or 5 in case you have a known singularity) and integrate in “one-step” mode. This is a fairly common pattern in classic Fortran ODE solvers.
A similar mechanism is used also in RKC (Runge-Kutta-Chebyshev) with which I’m more familiar. The rkc subroutine has the interface:
subroutine rkc(neqn,f,y,t,tend,rtol,atol,info,work,idid)
The array of integers info is used to control the procedure. Specifically info(1) fulfill the same role as itask:
c INFO(1): RKC integrates the initial value problem from T to TEND.
c This is done by computing approximate solutions at points
c chosen automatically throughout [T, TEND]. Ordinarily RKC
c returns at each step with an approximate solution. These
c approximations show how y behaves throughout the interval.
c The subroutine RKCINT can be used to obtain answers anywhere
c in the span of a step very inexpensively. This makes it
c possible to obtain answers at specific points in [T, TEND]
c and to obtain many answers very cheaply when attempting to
c locating where some function of the solution has a zero
c (event location). Sometimes you will be interested only in
c a solution at TEND, so you can suppress the returns at each
c step along the way if you wish.
c
c INFO(1) = 0 Return after each step on the way to TEND with a
c solution Y(*) at the output value of T.
c
c ` = 1 Compute a solution Y(*) at TEND only.
c
An example of integration with interpolation to produce output at regularly spaced time intervals can be found in rkc-exa.f.
1 Like
No, your expectations do not match the design of the VODE integrator package.
Look at the code in example1.f90 and see how it calls the integrator inside a loop. For your desired calculation you would have to set up a loop that is run 12 times, with Tout = 2, 4,…,24.
Note also that the driver program in example1.f90 does not itself call DVINDY. Such calls could be needed if, for example, you wanted to obtain many more values for plotting the solutions (which is sometimes labelled “dense output”).
1 Like
I think the purpose of DVINDY is actually to get derivatives of y as stated in this comment:
! CALL DVINDY(,,,,,) Provide derivatives of y, of various
! (See below.) orders, at a specified point T, if
! desired. It may be called only after
! a successful return from DVODE.
and not interpolation per se.
While the solution with itask=2 might work, I agree with what @mecej4 just suggested.
1 Like