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.
I need to study this to see if this is something I would want to do. Mainly I donāt want to burden the user if they just want to integrate a vector of real numbers (which is my main use case).
FYI: I just tagged this library as 1.0.0. The methods included are:
| Name | Description | Properties | Order | Stages | Registers | CFL | Reference |
|---|---|---|---|---|---|---|---|
euler |
Euler | 1 | 1 | 1 | 1.0 | Euler (1768) | |
midpoint |
Midpoint | 2 | 2 | 2 | ? | ||
heun |
Heun | 2 | 2 | 2 | ? | ||
rkssp22 |
2-stage, 2nd order TVD Runge-Kutta Shu-Osher | SSP | 2 | 2 | 1 | 1.0 | Shu & Oscher (1988) |
rk3 |
3th order Runge-Kutta | 3 | 3 | 3 | ? | ||
rkssp33 |
3-stage, 3rd order TVD Runge-Kutta Shu-Osher | SSP | 3 | 3 | 1 | 1.0 | Shu & Oscher (1988) |
rkssp53 |
5-stage, 3rd order SSP Runge-Kutta Spiteri-Ruuth | SSP | 3 | 5 | 2 | 2.65 | Ruuth (2006) |
rk4 |
Classic 4th order Runge-Kutta | 4 | 4 | 4 | Kutta (1901) | ||
rks4 |
4th order Runge-Kutta Shanks | 4 | 4 | 4 | Shanks (1965) | ||
rkr4 |
4th order Runge-Kutta Ralston | 4 | 4 | 4 | Ralston (1962) | ||
rkls44 |
4-stage, 4th order low storage non-TVD Runge-Kutta Jiang-Shu | LS | 4 | 4 | 2 | Jiang and Shu (1988) | |
rkls54 |
5-stage, 4th order low storage Runge-Kutta Carpenter-Kennedy | LS | 4 | 5 | 2 | 0.32 | Carpenter & Kennedy (1994) |
rkssp54 |
5-stage, 4th order SSP Runge-Kutta Spiteri-Ruuth | SSP | 4 | 5 | 4 | 1.51 | Ruuth (2006) |
rks5 |
5th order Runge-Kutta Shanks | 5 | 5 | 5 | Shanks (1965) | ||
rk5 |
5th order Runge-Kutta | 5 | 6 | 6 | ? | ||
rkc5 |
5th order Runge-Kutta Cassity | 5 | 6 | 6 | Cassity (1966) | ||
rkl5 |
5th order Runge-Kutta Lawson | 5 | 6 | 6 | Lawson (1966) | ||
rklk5a |
5th order Runge-Kutta Luther-Konen 1 | 5 | 6 | 6 | Luther & Konen (1965) | ||
rklk5b |
5th order Runge-Kutta Luther-Konen 2 | 5 | 6 | 6 | Luther & Konen (1965) | ||
rkb6 |
6th order Runge-Kutta Butcher | 6 | 7 | 7 | Butcher (1963) | ||
rk7 |
7th order Runge-Kutta Shanks | 7 | 9 | 9 | Shanks (1965) | ||
rk8_10 |
10-stage, 8th order Runge-Kutta Shanks | 8 | 10 | 10 | Shanks (1965) | ||
rkcv8 |
11-stage, 8th order Runge-Kutta Cooper-Verner | 8 | 11 | 11 | Cooper & Verner (1972) | ||
rk8_12 |
12-stage, 8th order Runge-Kutta Shanks | 8 | 12 | 12 | Shanks (1965) | ||
rkz10 |
10th order Runge-Kutta Zhang | 10 | 16 | 16 | Zhang (2019) | ||
rko10 |
10th order Runge-Kutta Ono | 10 | 17 | 17 | Ono (2003) | ||
rkh10 |
10th order Runge-Kutta Hairer | 10 | 17 | 17 | Hairer (1978) |
| Name | Description | Properties | Order | Stages | Registers | CFL | Reference |
|---|---|---|---|---|---|---|---|
rkbs32 |
Bogacki & Shampine 3(2) | FSAL | 3 | 4 | 4 | Bogacki & Shampine (1989) | |
rkssp43 |
4-stage, 3rd order SSP | SSP, LS | 3 | 4 | 2 | 2.0 | Kraaijevanger (1991), Conde et al. (2018) |
rkf45 |
Fehlberg 4(5) | 4 | 6 | 6 | Fehlberg (1969) | ||
rkck54 |
Cash & Karp 5(4) | 5 | 6 | 6 | Cash & Karp (1990) | ||
rkdp54 |
Dormand-Prince 5(4) | FSAL | 5 | 7 | 7 | Dormand & Prince (1980) | |
rkt54 |
Tsitouras 5(4) | FSAL | 5 | 7 | 7 | Tsitouras (2011) | |
rks54 |
Stepanov 5(4) | FSAL | 5 | 7 | 7 | Stepanov (2022) | |
rkpp54 |
Papakostas-PapaGeorgiou 5(4) | FSAL | 5 | 7 | 7 | Papakostas & Papageorgiou (1996) | |
rkpp54b |
Papakostas-PapaGeorgiou 5(4) b | FSAL | 5 | 7 | 7 | Papakostas & Papageorgiou (1996) | |
rkbs54 |
Bogacki & Shampine 5(4) | 5 | 8 | 8 | Bogacki & Shampine (1996) | ||
rkss54 |
Sharp & Smart 5(4) | 5 | 7 | 7 | Sharp & Smart (1993) | ||
rkdp65 |
Dormand-Prince 6(5) | 6 | 8 | 8 | Dormand & Prince (1981) | ||
rkc65 |
Calvo 6(5) | 6 | 9 | 9 | Calvo (1990) | ||
rktp64 |
Tsitouras & Papakostas NEW6(4) | 6 | 7 | 7 | Tsitouras & Papakostas (1999) | ||
rkv65e |
Verner efficient (9,6(5)) | FSAL | 6 | 9 | 9 | Verner (1994) | |
rkv65r |
Verner robust (9,6(5)) | FSAL | 6 | 9 | 9 | Verner (1994) | |
rkv65 |
Verner 6(5) | 6 | 8 | 8 | Verner (2006) | ||
dverk65 |
Verner 6(5) āDVERKā | 6 | 8 | 8 | Verner (?) | ||
rktf65 |
Tsitouras & Famelis 6(5) | FSAL | 6 | 9 | 9 | Tsitouras & Famelis (2006) | |
rktp75 |
Tsitouras & Papakostas NEW7(5) | 7 | 9 | 9 | Tsitouras & Papakostas (1999) | ||
rktmy7 |
7th order Tanaka-Muramatsu-Yamashita | 7 | 10 | 10 | Tanaka, Muramatsu & Yamashita (1992) | ||
rktmy7s |
7th order Stable Tanaka-Muramatsu-Yamashita | 7 | 10 | 10 | Tanaka, Muramatsu & Yamashita (1992) | ||
rkv76e |
Verner efficient (10:7(6)) | 7 | 10 | 10 | Verner (1978) | ||
rkv76r |
Verner robust (10:7(6)) | 7 | 10 | 10 | Verner (1978) | ||
rkss76 |
Sharp & Smart 7(6) | 7 | 11 | 11 | Sharp & Smart (1993) | ||
rkf78 |
Fehlberg 7(8) | 7 | 13 | 13 | Fehlberg (1968) | ||
rkv78 |
Verner 7(8) | 7 | 13 | 13 | Verner (1978) | ||
dverk78 |
Verner āMapleā 7(8) | 7 | 13 | 13 | Verner (?) | ||
rkdp85 |
Dormand-Prince 8(5) | 8 | 12 | 12 | Hairer (1993) | ||
rktp86 |
Tsitouras & Papakostas NEW8(6) | 8 | 12 | 12 | Tsitouras & Papakostas (1999) | ||
rkdp87 |
Dormand & Prince RK8(7)13M | 8 | 13 | 13 | Prince & Dormand (1981) | ||
rkv87e |
Verner efficient (8)7 | 8 | 13 | 13 | Verner (1978) | ||
rkv87r |
Verner robust (8)7 | 8 | 13 | 13 | Verner (1978) | ||
rkev87 |
Enright-Verner (8)7 | 8 | 13 | 13 | Enright (1993) | ||
rkk87 |
Kovalnogov-Fedorov-Karpukhina-Simos-Tsitouras 8(7) | 8 | 13 | 13 | Kovalnogov, Fedorov, Karpukhina, Simos, Tsitouras (2022) | ||
rkf89 |
Fehlberg 8(9) | 8 | 17 | 17 | Fehlberg (1968) | ||
rkv89 |
Verner 8(9) | 8 | 16 | 16 | Verner (1978) | ||
rkt98a |
Tsitouras 9(8) A | 9 | 16 | 16 | Tsitouras (2001) | ||
rkv98e |
Verner efficient (16:9(8)) | 9 | 16 | 16 | Verner (1978) | ||
rkv98r |
Verner robust (16:9(8)) | 9 | 16 | 16 | Verner (1978) | ||
rks98 |
Sharp 9(8) | 9 | 16 | 16 | Sharp (2000) | ||
rkf108 |
Feagin 8(10) | 10 | 17 | 17 | Feagin (2006) | ||
rkc108 |
Curtis 10(8) | 10 | 21 | 21 | Curtis (1975) | ||
rkb109 |
Baker 10(9) | 10 | 21 | 21 | Baker (?) | ||
rks1110a |
Stone 11(10) | 11 | 26 | 26 | Stone (2015) | ||
rkf1210 |
Feagin 12(10) | 12 | 25 | 25 | Feagin (2006) | ||
rko129 |
Ono 12(9) | 12 | 29 | 29 | Ono (2006) | ||
rkf1412 |
Feagin 14(12) | 14 | 35 | 35 | Feagin (2006) |
Sure. It should be doable to have both a simple interface and a more extensible interface. Like predefined types that can be used when working with a plain real array. I can come back with something in this vein.
@jacobwilliams: This is a commendable body of work, and I hope that the code package will be adopted widely. Being applied by many users will help with giving the package a good workout and correcting errors.
An instance where your integrator package can be used is the CAMB cosmology package, where they still use, in essence the old Netlib DVERK code, with about 50 GOTO statements in about 500 lines of code. Since I had recently looked at the CAMB software, after a user had posted a problem to the Intel Fortran site, I tried out the DVERK in your package, and found an error!
Line 1561 of file rklib_variable_steps.f90 now contains
call me%f(t+h, x,f1)
and should instead be
call me%f(t, x, f1)
Unfortunately, this sort of error can go unnoticed because the only consequence may be more integration steps used and small changes in the results ā changes that do not shout āerror here!ā unless one has at hand independently obtained reference results. Moreover, this error will have no effect if the ODEs being solved are autonomous.
For the test case that I used, which involved three coupled ODEs, with the error left uncorrected, I obtained:
tf = 12.906594
xf = 1.5631650764E-01 8.1839336759E-01 6.2787801657E-01
used 1712 function evaluations
and with the error corrected, I obtained:
tf = 12.906594
xf = 1.5656410324E-01 8.1904586941E-01 6.2743142363E-01
used 568 function evaluations
Searching for the same error pattern in the file revealed more instances. I do not know enough about the numerous R-K schemes to declare these as errors without looking up published references in each case, so users should take the following list as a hint for further checking:
771: call me%f(t+h, x,f1)
856: call me%f(t+h, x,f1)
951: call me%f(t+h, x,f1)
1474: call me%f(t+h, x,f1)
1561: call me%f(t+h, x,f1)
1783: call me%f(t+h, x,f1)
1907: call me%f(t+h, x,f1)
2345: call me%f(t+h, x,f1)
2702: call me%f(t+h, x,f1)
3711: call me%f(t+h, x,f1)
3839: call me%f(t+h, x,f1)
4826: call me%f(t+h, x,f1)
5631: call me%f(t+h, x,f1)
8221: call me%f(t+h, x,f1)
Woah! I will definitely look into this. Thanks!
Perhaps the folks in the Runge-Kutta Club (https://www.math.auckland.ac.nz/~butcher/RKclub/) would appreciate your work.
All these issues should be fixed now (in 1.0.1 release). Thanks for pointing this out to me!
Do you have benchmarks of these solvers compared to other diffeq libraries? Iād be interested to see comparisons vs DifferentialEquations.jl and Sundials.