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¶
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¶
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¶
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, rosenbrock-rodas3, 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 files: int/kpp_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/kpp_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/kpp_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/kpp_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/kpp_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.
Integrator inputs and outputs¶
In order to offer more control over the integrator, the KPP-generated subroutine provides the Optional integrator input parameters. Each of them is an array of 20 elements that allow the fine-tuning of the integrator.
Similarly, to obtain more information about the integration, the subroutine provides the Optional integrator output parameters, which are also also arrays of 20 elements.
Optional integrator input parameters¶
Optional integer (ICNTRL_U
) and real
(RCNTRL_U
) input parameters subroutine INTEGRATE
.
Setting any elements to zero will activate their default values. Array
elements not listed here are either not used or integrator-specific options.
Details can be found in the comment lines of the individual integrator files
$KPP_HOME/int/*.f90
.
-
ICNTRL_U(1)
¶
= 1
: \(F = F(y)\), i.e. independent of t (autonomous)= 0
: \(F = F(t,y)\), i.e. depends on t (non-autonomous)This option is only available for some of the integrators.
-
ICNTRL_U(2)
¶
The absolute (
ATOL
) and relative (RTOL
) tolerances can be expressed by either a scalar or individually for each species in a vector:= 0
:NVAR
-dimensional vector= 1
: scalar
-
ICNTRL_U(3)
¶
Selection of a specific method (only available for some of the integrators).
-
ICNTRL_U(4)
¶
Maximum number of integration steps.
-
ICNTRL_U(5)
¶
Maximum number of Newton iterations (only available for some of the integrators).
-
ICNTRL_U(6)
¶
Starting values of Newton iterations (only avaialble for some of the integrators).
= 0
: Interpolated= 1
: Zero
-
ICNTRL_U(15)
¶
This determines which
Update_*
subroutines are called within the integrator.= -1
: Do not call anyUpdate_*
subroutines= 0
: Use the integrator-specific default values> 1
: A number between 1 and 7, derived by adding up bits with values 4, 2, and 1. The first digit (4) activatesUpdate_SUN
. The second digit (2) activatesUpdate_PHOTO
. The third digit (1) activatesUpdate_RCONST
. |For example
ICNTRL(15)=6)
(4+2) will activate the calls toUpdate_SUN
andUpdate_PHOTO
, but not toUpdate_RCONST
.
-
RCNTRL_U(1)
¶
Hmin
, the lower bound of the integration step size. It is not recommended to change the default value of zero.
-
RCNTRL_U(2)
¶
Hmax
, the upper bound of the integration step size.
-
RCNTRL_U(3)
¶
Hstart
, the starting value of the integration step size.
-
RCNTRL_U(4)
¶
FacMin
, lower bound on step decrease factor.
-
RCNTRL_U(5)
¶
FacMax
, upper bound on step increase factor.
-
RCNTRL_U(6)
¶
FacRej
, step decrease factor after multiple rejections.
-
RCNTRL_U(7)
¶
FacSafe
, the factor by which the new step is slightly smaller than the predicted value.
-
RCNTRL_U(8)
¶
ThetaMin
. If the Newton convergence rate is smaller than ThetaMin, the Jacobian is not recomputed (only available for some of the integrators).
-
RCNTRL_U(9)
¶
NewtonTol
, the stopping criterion for Newton’s method (only available for some of the integrators).
-
RCNTRL_U(10)
¶
Qmin
(only available for some of the integrators).
-
RCNTRL_U(11)
¶
Qmax
. IfQmin < Hnew/Hold < Qmax
, then the step size is kept constant and the LU factorization is reused (only available for some of the integrators).
Optional integrator output parameters¶
Optional integer (ISTATUS_U
) and real (RSTATUS_U
)
output parameters of subroutine INTEGRATE
. Array elements not
listed here are either not used or are integrator-specific options.
Details can be found in the comment lines of the individual integrator files
$KPP_HOME/int/*.f90
.
-
ISTATUS_U(1)
¶
Number of function calls.
-
ISTATUS_U(2)`
¶
Number of Jacobian calls.
-
ISTATUS_U(3)
¶
Number of steps.
-
ISTATUS_U(4)
¶
Number of accepted steps.
-
ISTATUS_U(5)
¶
Number of rejected steps (except at very beginning).
-
ISTATUS_U(6)
¶
Number of LU decompositions.
-
ISTATUS_U(7)
¶
Number of forward/backward substitutions.
-
ISTATUS_U(8)
¶
Number of singular matrix decompositions.
-
RSTATUS_U(1)
¶
Texit
, the time corresponding to the computed \(Y\) upon return.
-
RSTATUS_U(2)
¶
Hexit
: the last accepted step before exit.
-
RSTATUS_U(3)
¶
Hnew
: The last predicted step (not yet taken. For multiple restarts, useHnew
asHstart
in the subsequent run.