/rklib

Fixed and variable-step Runge-Kutta solvers in Modern Fortran

Primary LanguageFortranBSD 3-Clause "New" or "Revised" LicenseBSD-3-Clause

rklib

rklib: A modern Fortran library of fixed and variable-step Runge-Kutta solvers.

Language GitHub release CI Status codecov last-commit

Description

The focus of this library is single-step, explicit Runge-Kutta solvers for 1st order differential equations.

Novel features:

  • 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_class in the code.
  • The real kind is selectable via a compiler directive (REAL32, REAL64, or REAL128).
  • Integration to an event is also supported. The root-finding method is also selectable (via the roots-fortran library).

Available Runge-Kutta methods:

  • Number of fixed-step methods: 27
  • Number of variable-step methods: 48
  • Total number of methods: 75

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)

Properties key:

  • LS = Low storage
  • SSP = Strong stability preserving
  • FSAL = First same as last
  • CFL = Courant-Friedrichs-Lewy

Example use case

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_example

The result is:

Success
tf = 100.00
xf = -0.1360372426E+01  0.1325538438E+01

Example performance comparison

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

Compiling

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

3rd Party Dependencies

  • The library requires roots-fortran.
  • The unit tests require pyplot-fortran, to generate the performance plots.
  • The coefficients app (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.

Documentation

The latest API documentation for the master branch can be found here. This was generated from the source code using FORD.

Notes

The original version of this code was split off from the Fortran Astrodynamics Toolkit in September 2022.

For developers

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.py to update all the include files. This script will generate all the boilerplate code for all the methods. It will also this README file.
  • Add a step function (either in rklib_fixed_steps.f90 or rklib_variable_steps.f90). Note that you can generate a template of an RK step function using the scripts/generate_rk_code.py script. 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.

License

The rklib source code and related files and documentation are distributed under a permissive free software license (BSD-3).

References

See also

  • 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