Announcing rklib, a new modern Fortran library for fixed and variable-step Runge-Kutta integration. Still sort of a work in progress, but it currently includes a couple dozen methods.
Some of the algorithms exist in no other publicly-available Fortran implementation as far as I can tell.
Join me and help improve it!
I have already programed rk4 and Dormand-Prince methods. What I do not have is a git* fully operating account.
There are a number of methods in the book by Hairer et al.
I will study your codes. Many thanks for sharing them.
Amazing work! Congrats @jacobwilliams.
I don’t have time right now but should really go back to refactoring ODEPACK, eventually it would be nice to have a unified Modern Fortran package for ODEs/DAEs…
I’m glad to see the Tsitouras method is included (actually, several versions of it.) It is not very well known, even though I think it should be the “new” standard of the “Runge-Kutta-something” explicit methods, performing slightly better than Dormand-Prince.
Edit: I would recommend you add some Tsitouras methods with the FSAL property. Specifically, Tsitouras-43 and SS6-5. The former might need more integration steps in complex problems, but each step is quite cheap, time-wise. It is also very efficient for simpler problems. The latter has 9 stages, one of them free. If fixed-step methods are included in the package, I guess those two shouldn’t be missing.
You might also consider adding the 2nd and 3rd order TVD schemes of Gottlieb, Shu and Osher for time dependent problems along with some of the low storage Strong-Stability Preserving schemes (SSP). These schemes are used a lot in CFD when you need time accurate solutions to unsteady problems… There are also some classes of implicit-explicit RK schemes that have some nice properties but they are probably beyond the scope of what you want. The explicit 2nd, 3rd, and 4th order TVD and SSP schemes are usually sufficient and are economical (and easy) to implement.
See.
Gottlieb and Shu, “Total Variation Diminishing Runge_Kutta Schemes” , Mathematics of Computation," Vol 67, Number 221, Jan. 1998.
Gottlieb, S., “On High Order Strong-Stability Preserving Runge-Kutta and Multi-step Time discretizations”, Journal of Scientific Computing, Vol 25 no 112, Nov. 2005.
Spiteri and Ruuth, "A new class of optimal high-order strong-stability preserving time discretization methods, SIAM Journal of Numerical Analysis, Vol. 40, no. 2, 2002 pp 469-491.
and the cited references
Very cool. This is great to see a nice implementation in Modern Fortran.
If such kind of TVD methods a judged interesting/useful, I have some implementations here HR-WENO/src/hrweno_tvdode.f90 at main · HugoMVale/HR-WENO · GitHub. My API is most likely different from the one employed in this RK library, but I could easily make the switch.
There is also more related in some more recent papers, such as
https://dx.doi.org/10.1016/j.jcp.2009.11.006
https://doi.org/10.1016/10.1137/10080960X
https://doi.org/10.5433/1679-0375.2019v40n2p123
just four examples lying on my disc, there is so much being published in this area…
I checked again today, and there are new methods added in the library, including my favorite for such kind of problems, namely Tsitouras-54. Usually, numeric computing environments have Dormand-Prince-54 as the default initial-value problem solver, but according to my experiments Tsitouras-54 is slightly better (performance- and accuracy-wise.) As with many other RK methods, the fifth-order (7 stages, FSAL) version is probably the most “balanced” one, with a very good accuracy/execution time ratio. But I see the library also has other, more expensive methods, in case “insane” accuracy is required.
There are definitely other methods that could be added, but this is already a direct “competitor” to a similar Julia library (which I’m not going to link here.) And it’s all in modern Fortran. I would argue that such a library doesn’t really need too much of OOP, but that’s just me and it doesn’t matter for the user anyway.
Definitely recommended. Well done @jacobwilliams!
Without a doubt, the Julia one is the state of the art in this area right now. It has everything. It’s the type of library I wish people were writing in Modern Fortran. Right now, I’m just concentrating on Runge-Kutta methods, but maybe one day this library could be expanded to include other types of methods.
I’m mainly interested in the high-order methods, but I will gladly accept pull requests for any of the others mentioned above. Join me!
One of the main things I want to add is dense output. Some of these RK pairs have interpolants for dense output, so I’d like to get those in there. I’ve used code with this feature before, but I’ve never coded it up before and don’t yet fully understand how to do that.
Great job ! I woud like to propose an alternative API that could offer more flexibility. The main idea is to allow the state “vector” to be of a user-defined type rather than a real array. Typically a derived type containing several arrays, possibly of different precisions, possibly complex rather than real. I have a demo implementation here :
and an example here :
This leverages modern Fortran features such as abstract types. @jacobwilliams would you have time to have a look and comment ? Is that approach too different from yours or could it be merged into rklib ?
Agree with @dubos. We came up with a quite similar solution in our atmospheric dynamics code, and find it very convenient.