+
Skip to content

jacobwilliams/rklib

Repository files navigation

rklib

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

Language GitHub release CI Status codecov last-commit

Description

This is a work in progress!

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.

Available methods:

Method name Type Order Number of Stages Reference
euler Fixed-step 1 1 Euler (1768)
midpoint Fixed-step 2 2 ?
heun Fixed-step 2 2 ?
rkssp22 Fixed-step, SSP 2 2 Shu & Oscher (1988)
rk3 Fixed-step 3 3 ?
rkssp33 Fixed-step, SSP 3 3 Shu & Oscher (1988)
rk4 Fixed-step 4 4 Kutta (1901)
rks4 Fixed-step 4 4 Shanks (1965)
rkls44 Fixed-step, LS 4 4 Jiang and Shu (1988)
rkssp54 Fixed-step, SSP 4 5 Ruuth (2006)
rks5 Fixed-step 5 5 Shanks (1965)
rkb6 Fixed-step 6 7 Butcher (1963)
rk7 Fixed-step 7 9 Shanks (1965)
rk8-10 Fixed-step 8 10 Shanks (1965)
rkcv8 Fixed-step 8 11 Cooper & Verner (1972)
rk8-12 Fixed-step 8 12 Shanks (1965)
rkz10 Fixed-step 10 16 Zhang (2019)
rko10 Fixed-step 10 17 Ono (2003)
rkh10 Fixed-step 10 17 Hairer (2003)
rkbs32 Variable-step 3 4 (FSAL) Bogacki & Shampine (1989)
rkf45 Variable-step 4 6 Fehlberg (1969)
rkck54 Variable-step 5 6 Cash & Karp (1990)
rkdp54 Variable-step 5 7 (FSAL) Dormand & Prince (1980)
rkt54 Variable-step 5 7 (FSAL) Tsitouras (2011)
rks54 Variable-step 5 7 (FSAL) Stepanov (2022)
rkdp65 Variable-step 6 8 Dormand & Prince (1981)
rkc65 Variable-step 6 9 Calvo (1990)
rktp64 Variable-step 6 7 Tsitouras & Papakostas (1999)
rkv65e Variable-step 6 9 (FSAL) Verner (1994)
rkv65r Variable-step 6 9 (FSAL) Verner (1994)
rktf65 Variable-step 6 9 (FSAL) Tsitouras & Famelis (2006)
rktp75 Variable-step 7 9 Tsitouras & Papakostas (1999)
rktmy7 Variable-step 7 10 Tanaka, Muramatsu & Yamashita (1992)
rkv76e Variable-step 7 10 Verner (1978)
rkv76r Variable-step 7 10 Verner (1978)
rkf78 Variable-step 7 13 Fehlberg (1968)
rkv78 Variable-step 7 13 Verner (1978)
rktp86 Variable-step 8 12 Tsitouras & Papakostas (1999)
rkdp87 Variable-step 8 13 Prince & Dormand (1981)
rkv87e Variable-step 8 13 Verner (1978)
rkv87r Variable-step 8 13 Verner (1978)
rkk87 Variable-step 8 13 Kovalnogov, Fedorov, Karpukhina, Simos, Tsitouras (2022)
rkf89 Variable-step 8 17 Fehlberg (1968)
rkv89 Variable-step 8 16 Verner (1978)
rkt98a Variable-step 9 16 Tsitouras (2001)
rkv98e Variable-step 9 16 Verner (1978)
rkv98r Variable-step 9 16 Verner (1978)
rkf108 Variable-step 10 17 Feagin (2006)
rkc108 Variable-step 10 21 Curtis (1975)
rkb109 Variable-step 10 21 Baker (?)
rks1110a Variable-step 11 26 Stone (2015)
rkf1210 Variable-step 12 25 Feagin (2006)
rko129 Variable-step 12 29 Ono (2006)
rkf1412 Variable-step 14 35 Feagin (2006)

Notes:

  • LS = Low storage
  • SSP = Strong stability preserving
  • FSAL = First same as last

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 :: x0 = 0.0_wp !! initial x value
  real(wp),parameter :: dx = 1.0_wp !! initial step size
  real(wp),parameter :: xf = 100.0_wp !! endpoint of integration
  real(wp),dimension(n),parameter :: y0 = [0.0_wp,0.1_wp] !! initial y value

  type(rktp86_class) :: prop
  real(wp),dimension(n) :: yf
  character(len=:),allocatable :: message

  call prop%initialize(n=n,f=fvpol,rtol=[tol],atol=[tol])
  call prop%integrate(x0,y0,dx,xf,yf)
  call prop%status(message=message)

  write (output_unit,'(A)') message
  write (output_unit,'(A,F7.2/,A,2E18.10)') &
              'xf =',xf ,'yf =',yf(1),yf(2)

contains

  subroutine fvpol(me,x,y,f)
  !! Right-hand side of van der Pol equation

  implicit none

  class(rk_class),intent(inout)     :: me
  real(wp),intent(in)               :: x
  real(wp),dimension(:),intent(in)  :: y
  real(wp),dimension(:),intent(out) :: f

  real(wp),parameter :: mu  = 0.2_wp

  f(1) = y(2)
  f(2) = mu*(1.0_wp-y(1)**2)*y(2) - y(1)

  end subroutine fvpol

end program rklib_example

The result is:

xf = 100.00
yf = -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

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

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

License

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

References

See also

点击 这是indexloc提供的php浏览器服务,不要输入任何密码和下载