MRKISS Library
MR RK (Runge-Kutta) KISS (Keep It Super Simple)

MRKISS is a simple, tiny library with zero dependencies[*] that aims to make it easy to use and experiment with explicit Runge-Kutta methods.

From a feature standpoint this library doesn't do anything you can't find in more comprehensive packages like SUNDIALS or Julia's DifferentialEquations.jl; however, it has the rather charming feature of being super simple and self contained. This makes it ideal for sharing results and code with others others without asking them to recreate my technical computing environment or learn about new software. It also makes it easy to get things up and and running quickly when moving to a new supercomputer, cloud, or HPC cluster.

Here is a complete example in 23 lines that solves Lorenz's system and generates a CSV file of the results.

!-------------------------------------------------------------------------------------------------
program minimal

  use :: mrkiss_config,      only: rk, ik
  use :: mrkiss_solvers_nt,  only: steps_fixed_stab_nt
  use :: mrkiss_utils,       only: print_solution
  use :: mrkiss_erk_kutta_4, only: a, b, c

  implicit none

  real(kind=rk),    parameter :: y_iv(3)  = [1.0_rk, 0.0_rk, 0.0_rk]
  real(kind=rk),    parameter :: param(3) = [10.0_rk, 28.0_rk, 8.0_rk/3.0_rk]
  real(kind=rk),    parameter :: t_end    = 50.0_rk

  real(kind=rk)               :: solution(4, 10000)
  integer(kind=ik)            :: status, istats(16)

  call steps_fixed_stab_nt(status, istats, solution, eq, y_iv, param, a, b, c, t_end_o=t_end)
  call print_solution(status, solution, filename_o="minimal.csv")

contains

  subroutine eq(status, dydt, y, param)
    integer(kind=ik), intent(out) :: status
    real(kind=rk),    intent(out) :: dydt(:)
    real(kind=rk),    intent(in)  :: y(:)
    real(kind=rk),    intent(in)  :: param(:)
    dydt = [ param(1)*(y(2)-y(1)), y(1)*(param(2)-y(3))-y(2), y(1)*y(2)-param(3)*y(3) ]
    status = 0
  end subroutine eq

end program minimal

You can check out a couple larger examples:

• Three Body Problem
• Lornez Attracter

If you just want to jump in, then take a look at the Quick Start section.

The IVPs I work with are generally pretty well behavied:

• Generally non-stiff
• Time forward (\(\Delta{t} \ge 0\))
• Defined by a small (<50) set of equations expressable in closed form.

Typical examples are strange attractors and systems related to chaotic science models from celestial/classical mechanics, population dynamics, oscillating chemical reactions, and electronic circuits.

My motivation for solving IVPs generally revolve around generative art and visualization. You will actually see this in the code and feature set of the library.

Things I care about:

• Simple to use for simple problems.
• Easily create custom solvers for the, admittedly bizarre, demands of generative art.
• Graceful response to evaluation failure in derivative functions
• A good selection of predefined RK methods
• Easy to use, hardwired methods for fixed step size visualization use cases:
  • Fixed \(t\) step size solvers
  • Fixed \(\mathbf{y}\) space step size solvers
• Solutions include derivative values by default for visualization tools that perform Hermite interpolation.
• Interpolate entire solutions to new time points (Hermite & linear).
• Programmable step processing. Examples:
  • Stop the routine if the solution curve is too long in y-space
  • Stop the routine if the step delta, or some components of it, are too long in y-space
  • Stop the routine if the solution has returned to the IV
  • Stop the routine if the solution intersects itself
  • Provide an alternate y-delta and redo the step based on some condition.
  • Trigger a bisection search for a t_delta fitting some condition based on t-space and/or y-space. Examples:
    • Find t_delta so that y-delta, or some components of it, are the perfect length.
    • Find where a step crosses over a boundary in space (ex: root finding)
    • Find where a step approaches closest to a point (ex: like the problem's IV)
• Runge-Kutta Research
  • Try out new RK methods by simply feeding the solvers a Butcher tableau.
  • Directly accessible one step routines for assembling custom solvers.
  • Simple code flow to facilitate instrumentation and deep runtime analysis and reporting.
  • Individual access to each method in an embedded tableau, and control over how each is used.
  • Maple worksheets rational values, variable floating point approximations, and stability graphs for every Tableau.
  • I have included a few RK methods more for research interests than practical usefulness.
• Easy deployment & sharing
  • Easy to compile and tune for a new architecture.
  • Zero external dependencies[*] except a Fortran compiler.
  • 100% standard Fortran that works with various compilers (Intel, Cray, NAG, gfortran, clang fortran, Nvidia, etc…).
  • Simple text output that can be compressed and sent back home or shared with others.

Things I don't care about:

• Usage error checking. For example, the code makes no attempt to check that the user has supplied consistent Butcher tableau arguments, or that t_delta values are positive, etc…
• Performance. I can generally perform hundreds of thousands of RK steps in a few milliseconds for the problems I work with. This gives me a lot of performance headroom allowing me to not worry about sophisticated techniques to avoid RK steps. In fact, this library diverges from best practices in a couple significant ways:
  • I don't use interpolating polynomials for intrastep approximations. I even have a bisection routine that takes an RK step for every bisection!
  • I use generic loops to compute RK steps over the Butcher tableau instead of optimized formulas.
  • Butcher tableau arrays are not sparse. In fact, I even include the top and final row full of zeros!

Within the confines of this software, we define a system of ODEs as:

\[ \frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} = \mathbf{f}(t, \mathbf{y}) = \left[\begin{array}{c} \frac{\mathrm{d}y_1}{\mathrm{d}t} \\ \vdots \\ \frac{\mathrm{d}y_n}{\mathrm{d}t} \\ \end{array}\right] = \left[\begin{array}{c} f_1(t, \mathbf{y}) \\ \vdots \\ f_n(t, \mathbf{y}) \\ \end{array}\right] = \left[\begin{array}{c} f_1(t, [y_1, \cdots, y_n]^\mathrm{T}) \\ \vdots \\ f_n(t, [y_1, \cdots, y_n]^\mathrm{T}) \\ \end{array}\right] \]

The goal is to find numerical values for the unknown function \(\mathbf{y}:\mathbb{R}\rightarrow\mathbb{R}^{n}\).

We define an embedded explicit Runge-Kutta method via a set of coefficients organized into a Butcher tableau:

\[ \begin{array}{l|llll} c_1 & a_{11} & a_{12} & \dots & a_{1s} \\ c_2 & a_{21} & a_{22} & \dots & a_{2s} \\ c_3 & a_{31} & a_{32} & \dots & a_{3s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \dots & a_{ss} \\ \hline \rule{0pt}{12pt} & \check{b}_1 & \check{b}_2 & \dots & \check{b}_s \\ & \hat{b}_1 & \hat{b}_2 & \dots & \hat{b}_s \\ \end{array} \]

Explicit methods, which are the focus of MRKISS, have \(c_1=0\) and \(a_{ij}=0\) for \(i\le j\).

The word embedded indicates that we actually have two explicit Runge-Kutta methods using the same \(\mathbf{a}\) matrix and \(\mathbf{c}\) vector. That is to say each \(\mathbf{b}\) vector defines a unique, explicit Runge-Kutta method. MRKISS supports both embedded and non-embedded (no \(\mathbf{\hat{b}}\) vector defined) methods.

Given \(\Delta{t}\) and initial conditions (\(t_0\) and \(\mathbf{y_0}\)), we may form an approximation of \(\mathbf{y}(t_0+\Delta{t})\) as:

\[ \mathbf{y}(t_0+\Delta{t}) \approx \mathbf{y_0}+\mathbf{\Delta\check{y}} \]

and, for embedded methods, an estimate of this approximation's error from:

\[\vert\mathbf{\Delta\check{y}} - \mathbf{\Delta\hat{y}} \vert\]

With \(\mathbf{\Delta\check{y}}\) and \(\mathbf{\Delta\hat{y}}\) (we only have \(\mathbf{\Delta\hat{y}}\) for embedded methods) computed as follows:

\[ \begin{array}{l} \mathbf{\Delta\check{y}} = \Delta{t}\sum_{i=1}^s \check{b}_i \mathbf{k}_i \\ \mathbf{\Delta\hat{y}} = \Delta{t}\sum_{i=1}^s \hat{b}_i \mathbf{k}_i \\ \end{array} \]

and the \(\mathbf{k}_i\) defined as:

\[ \mathbf{k}_i = \mathbf{f}\left(t + c_i \Delta{t},\, \mathbf{y} + \Delta{t} \sum_{j=1}^{i-1} a_{ij} \mathbf{k}_j\right) \]

In MRKISS an explicit Runge-Kutta method is specified by directly providing the Butcher tableau via arguments to subroutines.

MRKISS provides several predefined methods in modules found in the "lib/" directory. Each module defines a single tableau via parameters with names mirroring the Butcher Tableau arguments documented in the previous section. In addition, these modules also have a parameter containing the number of stages for the overall method and the number of stages for any embedded method that differs from the overall method.

• s – The number of stages for the entire method.
• s1 – The number of stages for the b1 method if it differs from s.
• s2 – The number of stages for the b2 method if it differs from s.

In some special cases an EERK may have more than two methods embedded. If so you may find variables for these additional methods following the same naming conventions. See mrkiss_eerk_cash_karp_5_4.f90 for an example.

The modules follow a simple naming conventions:

• They have one of two prefixes:

  mrkiss_eerk_

  The module contains an embedded explicit Runge Kutta method.

  mrkiss_erk_

  The module contains an explicit Runge Kutta method – i.e. it is not embedded.

• The names end with numbers indicating the orders of the b1 and b2 methods. These numbers are separated from the rest of the name by an underscore.

In addition to the parameters, the comments in these files normally include at least the following three sections:

IMO

Personal commentary about the method in question. Please note this material is simply my personal opinion.

Known Aliases

These include names used in the literature as well as names in some common ODE software.

References

I try to include the original reference if I have it. I also frequently include discussions found in other texts.

To make all this concrete, here is what one of these modules looks like (mrkiss_erk_kutta_4.f90):

!-----------------------------------------------------------------------------------------------------------------------
!> Butcher tableau for the classic 4 stage Runge-Kutta method of O(4)
!!
!! IMO: Useful for low accuracy applications; however, I find I rarely use it.
!!
!! Known Aliases: 'RK4' (OrdinaryDiffEq.jl), 'RK41' (Butcher), & 'The Runge-Kutta Method'.
!!
!! References:
!!   Kutta (1901); Beitrag Zur N\"herungsweisen Integration Totaler Differentialgleichungen; Z. Math. Phys. 46; p435-53
!!   Hairer, Norsett & Wanner (2009). Solving Ordinary Differential Equations. I: Nonstiff Problems. p138
!!   Butcher (2016); Numerical Methods for Ordinary Differential Equations. 3rd Ed; Wiley; p102
!!
module mrkiss_erk_kutta_4
  use mrkiss_config, only: rk, ik
  implicit none
  public
  integer(kind=ik), parameter :: s      = 4
  real(kind=rk),    parameter :: a(s,s) = reshape([ 0.0_rk, 0.0_rk, 0.0_rk, 0.0_rk,  &
                                                    1.0_rk, 0.0_rk, 0.0_rk, 0.0_rk,  &
                                                    0.0_rk, 1.0_rk, 0.0_rk, 0.0_rk,  &
                                                    0.0_rk, 0.0_rk, 2.0_rk, 0.0_rk], [s, s]) / 2.0_rk
  real(kind=rk),    parameter :: c(s)   = [         0.0_rk, 1.0_rk, 1.0_rk, 2.0_rk]          / 2.0_rk
  integer(kind=ik), parameter :: p      = 4
  real(kind=rk),    parameter :: b(s)   = [         1.0_rk, 2.0_rk, 2.0_rk, 1.0_rk]          / 6.0_rk
end module mrkiss_erk_kutta_4

Also note all the zeros. KISS! Seriously, it takes up a tiny bit of extra space and simplifies the code considerably…

Each embedded method defines two Runge-Kutta methods. Normally these two methods are used in conjunction to simultaneously estimate the solution and the error. In this library, the p1 & b1 method is recommended for approximating the solution while the p2 & b2 method should be used to estimate error. This is a recommendation, and is in no way enforced by the library. When the higher order method is used for the solution, we say we are using local extrapolation. Note that each of the methods in an embedded Butcher tableau may be used individually as a non-embedded method.

In addition to the module files, several maple worksheets may be found in the "rk_methods_maple/" directory. The filenames mirror the names of the modules. These worksheets contain the coefficients for the method's Butcher tableau, code to convert the coefficients into floating point values, and a plot of the method's stability region.

Throughout the code you will see subroutines, functions, and types suffixed with "_nt" or "_wt":

• _nt stands for "No T" – homogeneous problems.
• _wt stands for "With T" – non-homogeneous problems.

In the documentation below you will see "_*t" in subroutine names as shorthand to indicate both the "_nt" and "_wt" versions.

The equation to be solved is implimented in a user provided subroutine with one of the following two signatures:

For Non-Homogeneous (with t) problems:

subroutine deq_iface_wt(status, dydt, t, y, param) ! Non-Homogeneous Case (with t)
  integer(kind=ik), intent(out) :: status
  real(kind=rk),    intent(out) :: dydt(:)
  real(kind=rk),    intent(in)  :: t
  real(kind=rk),    intent(in)  :: y(:)
  real(kind=rk),    intent(in)  :: param(:)

For Homogeneous (no t) problems:

subroutine deq_iface_nt(status, dydt, y, param)    ! Homogeneous Case (no t)
  integer(kind=ik), intent(out) :: status
  real(kind=rk),    intent(out) :: dydt(:)
  real(kind=rk),    intent(in)  :: y(:)
  real(kind=rk),    intent(in)  :: param(:)

The arguments are as follows:

status ........ A status code. A positive value indicates failure.
                Do not return a value larger than 255!
dydt .......... The value of for f(t, y) is returned in this argument
t ............. The time (only for deq_iface_wt)
y ............. Values for the dependent variables
param ......... Constant parameters

This function should return the value for \( \mathbf{f}(t, \mathbf{y}) \) in dydt. The value of status should be non-positive, \((-\infty, 0]\), if everything worked, and a value between 1 and 255 inclusive, \([1, 255]\), if something went wrong. This value will be passed back via the status argument of higher level routines to indicate an error condition.

steps_adapt_etab_*t() uses traditional adaptive step size

• This solver is very similar to solvers found in other ODE packages.
• Programmable step processing
• A programmable bisection option to solve for interesting t_delta values
• Sophisticated curve length computations, and exit options when a maximum length is reached
• It can end precisely on a time value, or it can simply quit when a step goes beyond a maximum time value.
• These last two could be achieved with the programmable step processing and bisection features, but these requirements are so common that is convenient to have them directly available.

steps_fixed_stab_*t() uses fixed time steps

• Solution points separated by fixed time steps allow animations of the solution to naturally display velocity.
• This is a good place to start when writing a custom solver.
• With most modern ODE packages, this would be done with interpolation.
• This routine has the option to use Richardson extrapolation.

steps_condy_stab_*t() uses fixed (y-space) steps

• Produce solution points separated by fixed deltas in y-space, or some subset of y-space.
• This is a good place to start when writing a custom solver with a bisection step.
• A parametric plot of the first two components of a solution looks better when the points are uniformly separated.
• With most modern ODE packages, this would be done with interpolation.

Behind all of the above high level solvers are single step routines to carry out the step calculations. These are handy for creating DIY solvers.

• one_step_stab_*t() non-embedded RK methods
• one_richardson_step_stab_*t() uses Richardson extrapolation with non-embedded RK methods
• one_step_etab_*t() embedded RK methods
• one_step_rk4_*t() hardwired RK4 for unit tests
• one_step_rkf45_*t() hardwired RKF45 for unit tests

If you are interested playing around with MRKISS as quickly as possible, then this section is for you.

All of the code is in the module source files with no external dependencies at all. So you just need to call the modules from your code, and then compile/link everything together.

You can do that by just listing all the source files on the command line with most Fortran compilers. For example, we could compile the lorenz.f90 example in the examples/ directly like this:

cd examples
gfortran.exe lorenz.f90 ../src/*.f90

That said, most people will probably want to use a build system. If GNU Make is your thing, then the files in the make_include/ directory may be of help. In particular the makefile fragment include.mk provides useful targets and variables. The makefile in the examples/ directory is a good guide on how to use include.mk. In essence you do the following in your makefile:

1. Set MRKISS_PATH in your makefile to the path of the MRKISS source directory – that's the one with the include.mk file.
2. Set FC, FFLAGS, & AR if necessary – most of the time you can use the defaults.
3. Include the "include.mk" file in the MRKISS source directory.
4. Add a build rule for your program.

Your makefile will look something like this:

MRKISS_PATH = ../MRKISS

# Set FC, FFLAGS, & AR here.  The include below has the settings I use on my system.
include $(MRKISS_PATH)/tools_gfortran.mk

include $(MRKISS_PATH)/include.mk

your_program : your_program.f90 $(MRKISS_OBJ_FILES)
    $(FC) $(FFLAGS) $^ -o $@

Note the rule for your_program in the makefile above takes the lazy approach of adding every MRKISS module as a dependency regardless of if your program actually needs them all. This is how most people use the modules because it's simple. The cost might be a couple seconds of extra compile time. You can explicitly list out the modules in the makefile if you wish. Such a rule might look like the following:

your_program : your_program.f90 mrkiss_config$(OBJ_SUFFIX) mrkiss_solvers_wt(OBJ_SUFFIX) mrkiss_utils$(OBJ_SUFFIX)
    $(FC) $(FFLAGS) $^ -o $@

This section is about how I test MRKISS.

The tests/ directory contains code I primary use for testing MRKISS while the examples/ directory contains code I primarily use to demonstrate how to use MRKISS. The difference between a "test" and an "example" in MRKISS is a little bit slippery. Some of the tests, like tc1_* and tc2_*, could be considered demonstrations. In addition, I use all of the code in examples/ for tests.

The tests can be run by changing into the appropriate directory (tests/ or examples/), and building the make target tests. For example:

cd tests
make -j 16 tests

Note the -j 16 argument to make. When running all the tests, especially in the tests/ directory, I strongly recommend running in parallel.

In addition, the make files in tests/ and examples/ have numerous additional targets to run various individual tests or subsets of the test suite. These targets are documented in the subsections below.