rklib: A modern Fortran library of fixed and variable-step Runge-Kutta solvers.
The focus of this library is single-step, explicit Runge-Kutta solvers for 1st order differential equations.
- The library includes a wide range of both fixed and variable-step Runge-Kutta methods, from very low to very high order.
- It is object-oriented and written in modern Fortran.
- It allows for defining a variable-step size integrator with a custom-tuned step size selection method. See
stepsize_classin the code. - The
realkind is selectable via a compiler directive (REAL32,REAL64, orREAL128). - Integration to an event is also supported. The root-finding method is also selectable (via the roots-fortran library).
- Number of fixed-step methods: 27
- Number of variable-step methods: 48
- Total number of methods: 75
| 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) |
- LS = Low storage
- SSP = Strong stability preserving
- FSAL = First same as last
- CFL = Courant-Friedrichs-Lewy
Basic use of the library is shown here (this uses the rktp86 method):
program rklib_example
use rklib_module, wp => rk_module_rk
use iso_fortran_env, only: output_unit
implicit none
integer,parameter :: n = 2 !! dimension of the system
real(wp),parameter :: tol = 1.0e-12_wp !! integration tolerance
real(wp),parameter :: t0 = 0.0_wp !! initial t value
real(wp),parameter :: dt = 1.0_wp !! initial step size
real(wp),parameter :: tf = 100.0_wp !! endpoint of integration
real(wp),dimension(n),parameter :: x0 = [0.0_wp,0.1_wp] !! initial x value
real(wp),dimension(n) :: xf !! final x value
type(rktp86_class) :: prop
character(len=:),allocatable :: message
call prop%initialize(n=n,f=fvpol,rtol=[tol],atol=[tol])
call prop%integrate(t0,x0,dt,tf,xf)
call prop%status(message=message)
write (output_unit,'(A)') message
write (output_unit,'(A,F7.2/,A,2E18.10)') &
'tf =',tf ,'xf =',xf(1),xf(2)
contains
subroutine fvpol(me,t,x,f)
!! Right-hand side of van der Pol equation
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(:),intent(in) :: x
real(wp),dimension(:),intent(out) :: f
f(1) = x(2)
f(2) = 0.2_wp*(1.0_wp-x(1)**2)*x(2) - x(1)
end subroutine fvpol
end program rklib_exampleThe result is:
Success
tf = 100.00
xf = -0.1360372426E+01 0.1325538438E+01
Running the unit tests will generate some performance plots. The following is for the variable-step methods compiled with quadruple precision (e.g, fpm test rk_test_variable_step --compiler ifort --flag "-DREAL128"): rk_test_variable_step_R16.pdf
A Fortran Package Manager manifest file is included, so that the library and test cases can be compiled with FPM. For example:
fpm build --profile release
fpm test --profile release
To use rklib within your FPM project, add the following to your fpm.toml file:
[dependencies]
rklib = { git="https://github.com/jacobwilliams/rklib.git" }By default, the library is built with double precision (real64) real values. Explicitly specifying the real kind can be done using the following processor flags:
| Preprocessor flag | Kind | Number of bytes |
|---|---|---|
REAL32 |
real(kind=real32) |
4 |
REAL64 |
real(kind=real64) |
8 |
REAL128 |
real(kind=real128) |
16 |
For example, to build a single precision version of the library, use:
fpm build --profile release --flag "-DREAL32"
To generate the documentation using FORD, run:
ford ford.md
- The library requires roots-fortran.
- The unit tests require pyplot-fortran, to generate the performance plots.
- The
coefficientsapp (not required to use the library, but used to generate some of the code) requires the mpfun2020-var1 arbitrary precision library.
All of these will be automatically downloaded by FPM.
The latest API documentation for the master branch can be found here. This was generated from the source code using FORD.
The original version of this code was split off from the Fortran Astrodynamics Toolkit in September 2022.
To add a new method to this library:
- Update the tables (either the fixed or variable one in
scripts/generate_files.py) - Run
python scripts/generate_files.pyto update all the include files. This script will generate all the boilerplate code for all the methods. It will also thisREADMEfile. - Add a step function (either in
rklib_fixed_steps.f90orrklib_variable_steps.f90). Note that you can generate a template of an RK step function using thescripts/generate_rk_code.pyscript. The two command line arguments are the number of function evaluations required and the method name (e.g.,'rk4'). Edit the template accordingly (note at the FSAL ones have a slightly different format). - Update the unit tests.
The rklib source code and related files and documentation are distributed under a permissive free software license (BSD-3).
- E. B. Shanks, "Higher Order Approximations of Runge-Kutta Type", NASA Technical Note, NASA TN D-2920, Sept. 1965.
- E. B. Shanks, "Solutions of Differential Equations by Evaluations of Functions" Math. Comp. 20 (1966).
- E. Fehlberg, "Classical Fifth-, Sixth-, Seventh-, and Eighth-Order Runge-Kutta Formulas with Stepsize Control", NASA TR R-2870, 1968.
- E. Fehlberg, "Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems", NASA Technical Report R-315, July 1, 1969.
- J. H. Verner, "Explicit Runge-Kutta Methods with Estimates of the Local Truncation Error", SIAM Journal on Numerical Analysis, 15(4), 772-790, 1978.
- T. Feagin, "High-Order Explicit Runge-Kutta Methods"
- J. C. Butcher, "A history of Runge-Kutta methods", Applied Numerical Mathematics, Volume 20, Issue 3, March 1996, Pages 247-260
- J. C. Butcher, "On Runge-Kutta Processes of High Order", Oct. 28, 1963.
- G. E. Müllges & F. Uhlig, "Numerical Algorithms with Fortran", Springer, 1996.
- K. Fox, "Numerical Integration of the Equations of Motion of Celestial Mechanics", Celestial Mechanics 33, p 127-142, 1984.
- Mathematics Source Library
- Maple worksheets on the derivation of Runge-Kutta schemes
- Index of numerical integrators
- J. Williams, Fehlberg's Runge-Kutta Methods, Feb. 10, 2018.
- C.-W. Shu, S. Osher, "Efficient implementation of essentially non-oscillatory shock-capturing schemes", Journal of Computational Physics, 77(2), 439-471, 1988.
- S. Ruuth, "Global optimization of explicit strong-stability-preserving Runge-Kutta methods." Mathematics of Computation 75.253 (2006): 183-207.
- Jiang, Guang-Shan, and Chi-Wang Shu. "Efficient implementation of weighted ENO schemes." Journal of computational physics 126.1 (1996): 202-228.
- FOODIE Fortran Object-Oriented Differential-equations Integration Environment
- FLINT Fortran Library for numerical INTegration of differential equations
- DDEABM Modern Fortran implementation of the DDEABM Adams-Bashforth algorithm
- DOP853 Modern Fortran Edition of Hairer's DOP853 ODE Solver. An explicit Runge-Kutta method of order 8(5,3) for problems y'=f(x,y); with dense output of order 7
- DVODE Modern Fortran Edition of the DVODE ODE Solver
- ODEPACK Work in Progress to refactor and modernize the ODEPACK Library
- libode Easy-to-compile, high-order ODE solvers as C++ classes