A new Runge-Kutta Fortran library

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!

42 Likes

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.

2 Likes

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ā€¦

2 Likes

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.

3 Likes

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

4 Likes

Very cool. This is great to see a nice implementation in Modern Fortran.

1 Like

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.

1 Like

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ā€¦

2 Likes

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!

1 Like

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.

8 Likes

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 ?

3 Likes

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).

1 Like

FYI: I just tagged this library as 1.0.0. The methods included are:

Fixed-step methods:

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)

Variable-step methods:

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)
15 Likes

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.

1 Like

@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)
2 Likes

Woah! I will definitely look into this. Thanks!

1 Like

Perhaps the folks in the Runge-Kutta Club (https://www.math.auckland.ac.nz/~butcher/RKclub/) would appreciate your work.

1 Like

All these issues should be fixed now (in 1.0.1 release). Thanks for pointing this out to me!

1 Like

Do you have benchmarks of these solvers compared to other diffeq libraries? Iā€™d be interested to see comparisons vs DifferentialEquations.jl and Sundials.

1 Like