Runge-Kutta (aka RK) methods

A general \(s\)-stage Runge-Kutta method is defined as (see Section II.1 of Hairer et al. [1987])

\[\begin{split}\begin{aligned} y^{n+1} &=& y^n + h \sum_{i=1}^s b_i k_i~,\\ \nonumber T_i &=& t^n + c_i h~, \quad Y_i = y^n + h \sum_{j=1}^{s} a_{ij} k_j~,\\ \nonumber k_i &=& f\left( \, T_i, \, Y_i \,\right)~,\end{aligned}\end{split}\]

where the coefficients \(a_{ij}\), \(b_i\) and \(c_i\) are prescribed for the desired accuracy and stability properties. The stage derivative values \(k_i\) are defined implicitly, and require solving a (set of) nonlinear system(s). Newton-type methods solve coupled linear systems of dimension (at most) \(n \times s\).

The Runge-Kutta methods implemented in KPP are summarized below:

3-stage Runge-Kutta

Integrator file: int/runge_kutta.f90

Fully implicit 3-stage Runge-Kutta methods. Several variants are available:

RADAU-2A: order 5
RADAU-1A: order 5
Lobatto-3C: order 4
Gauss: order 6

RADAU5

Integrator file: int/radau5.f90

This Runge-Kutta method of order 5 based on RADAU-IIA quadrature is stiffly accurate. The KPP implementation follows the original implementation of Hairer and Wanner [1991], Section IV.10. While RADAU5 is relatively expensive (when compared to the Rosenbrock methods), it is more robust and is useful to obtain accurate reference solutions.

SDIRK

Integrator file: int/sdirk.f90,

SDIRK is an L-stable, singly-diagonally-implicit Runge-Kutta method. The implementation is based on Hairer and Wanner [1991]. Several variants are available:

Sdirk 2a, 2b: 2 stages, order 2

Sdirk 3a: 3 stages, order 2

Sdirk 4a, 4b: 5 stages, order 4

SDIRK4

Integrator file: int/sdirk4.f90

SDIRK4 is an L-stable, singly-diagonally-implicit Runge-Kutta method of order 4. The implementation is based on Hairer and Wanner [1991].

SEULEX

Integrator file: int/seulex.f90

SEULEX is a variable order stiff extrapolation code able to produce highly accurate solutions. The KPP implementation is based on the implementation of Hairer and Wanner [1991].

RK tangent linear model

The tangent linear method associated with the Runge-Kutta method is

\[\begin{split}\begin{aligned} %y^{n+1} &=& y^n + h \sum_{i=1}^s b_i k_i~,\\ \delta y^{n+1} &=& \delta y^n + h \sum_{i=1}^s b_i \ell_i~,\\ \nonumber %Y_i &=& y^n + h \sum_{j=1}^{s} a_{ij} k_j~,\\ \delta Y_i& =& \delta y^n + h \sum_{j=1}^{s} a_{ij} \ell_j~,\\ \nonumber %k_i &=& f\left( \, T_i, \, Y_i \,\right)~,\\ \ell_i &=& J\left(T_i, \, Y_i \right) \cdot \delta Y_i ~.\end{aligned}\end{split}\]

The system is linear and does not require an iterative procedure. However, even for a SDIRK method (\(a_{ij}=0\) for \(i>j\) and \(a_{ii}=\gamma\)) each stage requires the LU factorization of a different matrix.

RK discrete adjoint model

The first order Runge-Kutta adjoint is

\[\begin{split}\begin{aligned} u_i &=& h \, J^T(T_i,Y_i)\cdot \left( b_i \lambda^{n+1} + \sum_{j=1}^s a_{ji} u_j \right)\\ %\quad i = 1 \cdots s\\ \nonumber \lambda^{n} &=& \lambda^{n+1} +\sum_{j=1}^s u_j~.\end{aligned}\end{split}\]

For \(b_i \ne 0\) the Runge-Kutta adjoint can be rewritten as another Runge-Kutta method:

\[\begin{split}\begin{aligned} u_i &=& h \, J^T(T_i,Y_i)\cdot \left( \lambda^{n+1} + \sum_{j=1}^s \frac{b_j \, a_{ji}}{b_i} u_j \right)\\ %~, \quad i = 1 \cdots s\\ \nonumber \lambda^{n} &=& \lambda^{n+1} +\sum_{j=1}^s b_j \, u_j~.\end{aligned}\end{split}\]