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 options include the error norm (max, 2-norm, RMS, or user-defined), the scaling of the relative tolerance, the step size limits and safety factor, the PI (Lund) stabilized step size control, and how the two error estimates of rkdp85 are combined. With the appropriate options, the step size control of Hairer's DOPRI5 and DOP853 codes is reproduced (see the documentation of stepsize_class).
• 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).
• Dense output (interpolation within a step) for some of the variable-step methods. See Dense output.

Available Runge-Kutta methods:

• Number of fixed-step methods: 28
• Number of variable-step methods: 48
• Total number of methods: 76

Fixed-step methods:

Variable-step methods:

Properties key:

• LS = Low storage
• SSP = Strong stability preserving
• FSAL = First same as last
• DENSE = Has a dense output (continuous extension). See Dense output.
• 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

Dense output

Variable-step methods with the DENSE property have an interpolating polynomial over each step, so the solution can be computed at any time within a step without taking extra integration steps. It is off by default (so there is no overhead), and is enabled when initializing the integrator:

• dense_output=.true. : enables the interpolate and dense_interval methods, and uses the dense output for root finding in integrate_to_event (which saves function evaluations).
• dense_output_step=dt : the report function is called at a fixed interval (t0, t0+dt, t0+2*dt, ...), plus the final time (or the event), regardless of the steps the integrator takes. This also enables the dense output.

For example, to report every 0.1 time units:

type(rkdp85_class) :: prop
call prop%initialize(n=n,f=fvpol,rtol=[tol],atol=[tol],report=report,dense_output_step=0.1_wp)
call prop%integrate(t0,x0,dt,tf,xf)

Or, to compute points yourself within the report function (which is called after each step):

subroutine report(me,t,x)
  class(rk_class),intent(inout)    :: me
  real(wp),intent(in)              :: t
  real(wp),dimension(:),intent(in) :: x
  logical :: available
  real(wp) :: t_prev
  real(wp),dimension(size(x)) :: x_mid
  call me%dense_interval(available,t0=t_prev) ! the step that was just taken: [t_prev, t]
  if (available) call me%interpolate((t_prev+t)/2.0_wp,x_mid) ! e.g., the midpoint
end subroutine report

When to use it:

• When the output is needed at specific times (e.g., evenly-spaced points for plotting, or to compare with measurements). This is much more efficient than limiting the step size (or using fixed_step_mode), since the integrator still takes the largest steps that meet the error tolerance.
• For event finding with an expensive derivative function.
• Not needed if only the final state (or the points at the integration steps) are required.

Notes: Only the most recent step can be interpolated, and interpolate must not be called from the derivative function. Some methods (e.g., rkdp85) need additional function evaluations for the dense output, but these are only done for steps where it is used.

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

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).
• If the method has a dense output (continuous extension), add DENSE to its properties in the table, and add a <name>_dense_coefficients routine to rklib_variable_dense.f90. This returns the coefficients as they are published, in either the polynomial form $b_i(\theta)$ (dense_form_polynomial) or the form used in DOPRI5/DOP853 (dense_form_dopri). If the dense output requires additional function evaluations, also add the method to the dense_extra_stage_methods list, and add a <name>_dense_extra_stages routine that computes them (these are only computed when the dense output is used). When enabled (with dense_output=.true. in initialize), the interpolated state can then be computed with the interpolate method (e.g., from the report function, using dense_interval to get the interval of the step that was just taken), and it is used for event location in integrate_to_event.
• The dense output routines for the Verner methods (rkv65e, rkv65r, rkv76e, rkv76r, rkv87e, rkv87r, rkv98e, rkv98r) and for rkdp87, rkev87 and dverk78 are generated by python scripts/generate_dense.py, from Verner's coefficient files (downloaded from the Internet Archive, since Verner's site is no longer online) and from Peter Stone's Maple worksheets. The script checks the coefficients against the step coefficients in rklib_variable_steps.f90, converts the interpolants to the form used in DOPRI5/DOP853 (which has a much smaller roundoff error for these coefficients), and updates the marked block in rklib_variable_dense.f90. Use --check to also verify all the order conditions of the interpolants. See the docstring of that script for details.
• 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