rklib

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

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:

Variable-step methods:

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:

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

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

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