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

This is a work in progress!

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

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

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

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"):

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:

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.

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

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
• FLINT
• DDEABM
• DOP853
• DVODE
• libode