Numerical methods¶
The KPP numerical library contains a set of numerical integrators selected to be very efficient in the low to medium accuracy regime (relative errors \(\sim 10^{-2} \dots 10^{-5}\)). In addition, the KPP numerical integrators preserve the linear invariants (i.e., mass) of the chemical system.
KPP implements several Rosenbrock methods: ROS–2 [Verwer et al., 1999], ROS–3 [Sandu et al., 1997], RODAS–3 [Sandu et al., 1997], ROS–4 [Hairer and Wanner, 1991], and RODAS–4 [Hairer and Wanner, 1991]. For each of them KPP implements the tangent linear model (direct decoupled ensitivity) and the adjoint models. The implementations distinguish between sensitivities with respect to initial values and sensitivities with respect to parameters for efficiency.
Note that KPP produces the building blocks for the simulation and also for the sensitivity calculations. It also provides application programming templates. Some minimal programming may be required from the users in order to construct their own application from the KPP building blocks.
In the following sections we introduce the numerical methods implemented in KPP. The symbols used in the formulas are explained in Table 20. Symbols used in numerical methods.
Symbol |
Description |
---|---|
\(s\) |
Number of stages |
\(t^n\) |
Discrete time moment |
\(h\) |
Time step \(h=t^{n+1}-t^n\) |
\(y^n\) |
Numerical solution (concentration) at \(t^n\) |
\(\delta y^n\) |
tangent linear solution at \(t^n\) |
\(\lambda^n\) |
Adjoint numerical solution at \(t^n\) |
\(f(\cdot,\cdot)\) |
The ODE derivative function: \(y'=f(t,y)\) |
\(f_t(\cdot,\cdot)\) |
Partial time derivative \(f_t( t,y)=\partial f(t,y)/\partial t\) |
\(J(\cdot,\cdot)\) |
The Jacobian \(J( t,y)=\partial f(t,y)/\partial y\) |
\(J_t(\cdot,\cdot)\) |
Partial time derivative of Jacobian \(J_t( t,y)=\partial J(t,y)/\partial t\) |
\(A\) |
The system matrix |
\(H(\cdot,\cdot)\) |
The Hessian \(H(t,y) =\partial^2 f(t,y)/\partial y^2\) |
\(T_i\) |
Internal stage time moment for Runge-Kutta and Rosenbrock methods |
\(Y_i\) |
Internal stage solution for Runge-Kutta and Rosenbrock methods |
\(k_i\), \(\ell_i\), \(u_i\), \(v_i\) |
Internal stage vectors for Runge-Kutta and Rosenbrock methods, their tangent linear and adjoint models |
\(\alpha_i\), \(\alpha_{ij}\), \(a_{ij}\), \(b_i\), \(c_i\), \(c_{ij}\), \(e_i\), \(m_i\) |
Method coefficients |
Rosenbrock methods¶
Integrator file: int/rosenbrock.f90
An \(s\)-stage Rosenbrock method (cf. Section IV.7 in [Hairer and Wanner, 1991]) computes the next-step solution by the formulas
where \(s\) is the number of stages, \(\alpha_i = \sum_j \alpha_{ij}\) and \(\gamma_i = \sum_j \gamma_{ij}\). The formula coefficients (\(a_{ij}\) and \(\gamma_{ij}\)) give the order of consistency and the stability properties. \(A\) is the system matrix (in the linear systems to be solved during implicit integration, or in the Newton’s method used to solve the nonlinear systems). It is the scaled identity matrix minus the Jacobian.
The coefficients of the methods implemented in KPP are shown below:
ROS-2¶
Stages (\(s\)): 2
Funcion calls: 2
Order: 2(1)
Stability properties: L-stable
Method Coefficients:
ROS-3¶
Stages (\(s\)): 3
Funcion calls: 2
Order: 3(2)
Stability properties: L-stable
Method Coefficients:
ROS-4¶
Stages (\(s\)): 4
Funcion calls: 3
Order: 4(3)
Stability properties: L-stable
Method Coefficients:
RODAS-3¶
Stages (\(s\)): 4
Funcion calls: 3
Order: 3(2)
Stability properties: Stiffly-accurate
Method Coefficients:
RODAS-4¶
Stages (\(s\)): 6
Funcion calls: 5
Order: 4(3)
Stability properties: Stiffly-accurate
Method Coefficients:
Rosenbrock tangent linear model¶
Integrator file: int/rosenbrock_tlm.f90
The Tangent Linear method is combined with the sensitivity equations. One step of the method reads:
The method requires a single n times n LU decomposition per step to obtain both the concentrations and the sensitivities.
KPP contains tangent linear models (for direct decoupled sensitivity analysis) for each of the Rosenbrock methods (ROS–2, ROS–3, ROS–4, RODAS–3, and RODAS–4). The implementations distinguish between sensitivities with respect to initial values and sensitivities with respect to parameters for efficiency.
Rosenbrock discrete adjoint model¶
Integrator file: int/rosenbrock_adj.f90
To obtain the adjoint we first differentiate the method with respect to \(y_n\). Here \(J\) denotes the Jacobian and \(H\) the Hessian of the derivative function \(f\). The discrete adjoint of the (non-autonomous) Rosenbrock method is
KPP contains adjoint models (for direct decoupled sensitivity analysis) for each of the Rosenbrock methods (ROS-2, ROS-3, ROS-4, RODAS-3, RODAS-4).
Runge-Kutta (aka RK) methods¶
A general \(s\)-stage Runge-Kutta method is defined as (see Section II.1 of [Hairer et al., 1993])
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
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
For \(b_i \ne 0\) the Runge-Kutta adjoint can be rewritten as another Runge-Kutta method:
Backward differentiation formulas¶
Backward differentiation formulas (BDF) are linear multistep methods with excellent stability properties for the integration of chemical systems (cf. [Hairer and Wanner, 1991], Section V.1). The \(k\)-step BDF method reads
where the coefficients \(\alpha_i\) and \(\beta\) are chosen such that the method has order of consistency \(k\).
The KPP library contains two off-the-shelf, highly popular implementations of BDF methods, described in the following sections:
LSODE¶
Integrator file: int/lsode.f90
LSODE, the Livermore ODE solver [Radhakrishnan and Hindmarsh, 1993], implements backward differentiation formula (BDF) methods for stiff problems. LSODE has been translated to Fortran90 for the incorporation into the KPP library.
VODE¶
Integrator file: int/dvode.f90
VODE [Brown et al., 1989] uses another formulation of backward differentiation formulas. The version of VODE present in the KPP library uses directly the KPP sparse linear algebra routines.