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
\[\begin{split}\begin{aligned}
y^{n+1} &=& y^n + \sum_{i=1}^s m_i k_i~,
\quad {\rm Err}^{n+1} = \sum_{i=1}^s e_i k_i\\
\nonumber
T_i &=& t^n + \alpha_i h~, \quad
Y_i =y^n + \sum_{j=1}^{i-1} a_{ij} k_j~,\\
\nonumber
A &=& \left[ \frac{1}{h \gamma} - J^T(t^n,y^n) \right]\\
\nonumber
A \cdot k_i &=& f\left( \, T_i,
\, Y_i \,\right) + \sum_{j=1}^{i-1} \frac{c_{ij}}{h} k_j + h \gamma_i
f_t\left(t^n,y^n\right)~.
\end{aligned}\end{split}\]
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
\[\begin{split}\begin{aligned}
\gamma = 1 + 1/sqrt{2} & \qquad & a_{2,1} = 1/\gamma & \qquad & c_{2,1} = -2/\gamma &\\
m_1 = 3/(2\gamma) & \qquad & m_2 = 1/(2\gamma) & \qquad & e_1 = 1/(2\gamma) &\\
e_2 = 1/(2\gamma) & \qquad & \alpha_1 = 0 & \qquad & \alpha_2 = 1 &\\
\gamma_1 = \gamma & \qquad & \gamma_2 = -\gamma &\\
\end{aligned}\end{split}\]
ROS-3
\[\begin{split}\begin{aligned}
a_{2,1} = 1 & \qquad & a_{3,1} = 1 & \qquad & a_{3,2} = 0 &\\
c_{2,1} = -1.015 & \qquad & c_{3,1} = 4.075 & \qquad & c_{3,2} = 9.207 &\\
m_1 = 1 & \qquad & m_2 = 6.169 & \qquad & m_3 = -0.427 &\\
e_1 = 0.5 & \qquad & e_2 = -2.908 & \qquad & e_3 = 0.223 &\\
alpha_1 = 0 & \qquad & \alpha_2 = 0.436 & \qquad & \alpha_3 = 0.436 &\\
\gamma_1 = 0.436 & \qquad & \gamma_2 = 0.243 & \qquad & \gamma_3 = 2.185 &\\
\end{aligned}\end{split}\]
ROS-4
\[\begin{split}\begin{aligned}
a_{2,1} = 2 & \qquad & a_{3,1} = 1.868 & \qquad & a_{3,2} = 0.234 &\\
a_{4,1} = a_{3,1} & \qquad & a_{4,2} = a_{3,2} & \qquad & a_{4,3} = 0 &\\
c_{2,1} = -7.137 & \qquad & c_{3,1} = 2.581 & \qquad & c_{3,2} = 0.652 &\\
c_{4,1} = -2.137 & \qquad & c_{4,2} = -0.321 & \qquad & c_{4,3} = -0.695 &\\
m_1 = 2.256 & \qquad & m_2 = 0.287 & \qquad & m_3 = 0.435 &\\
m_4 = 1.094 & \qquad & e_1 = -0.282 & \qquad & e_2 = -0.073 &\\
e_3 = -0.108 & \qquad & e_4 = -1.093 & \qquad & \alpha_1 = 0 &\\
\alpha_2 = 1.146 & \qquad & \alpha_3 = 0.655 & \qquad & \alpha_4 = \alpha_3 &\\
\gamma_1 = 0.573 & \qquad & \gamma_2 = -1.769 & \qquad & \gamma_3 = 0.759 &\\
\gamma_4 = -0.104
\end{aligned}\end{split}\]
RODAS-3
\[\begin{split}\begin{aligned}
a_{2,1} = 0 & \qquad & a_{3,1} = 2 & \qquad & a_{3,2} = 0 &\\
a_{4,1} = 2 & \qquad & a_{4,2} = 0 & \qquad & a_{4,3} = 1 &\\
c_{2,1} = 4 & \qquad & c_{3,1} = 1 & \qquad & c_{3,2} = -1 &\\
c_{4,1} = 1 & \qquad & c_{4,2} = -1 & \qquad & c_{4,3} = -8/3 &\\
m_1 = 2 & \qquad & m_2 = 0 & \qquad & m_3 = 1 &\\
m_4 = 1 & \qquad & e_1 = 0 & \qquad & e_2 = 0 &\\
e_3 = 0 & \qquad & e_4 = 1 & \qquad & \alpha_1 = 0 &\\
\alpha_2 = 0 & \qquad & \alpha_3 = 1 & \qquad & \alpha_4 = 1 &\\
\gamma_1 = 0.5 & \qquad & \gamma_2 = 1.5 & \qquad & \gamma_3 = 0 &\\
\gamma_4 = 0
\end{aligned}\end{split}\]
RODAS-3.1
\[\begin{split}\begin{aligned}
a_{2,1} = 0 & \qquad & a_{3,1} = 0.646601929740551 &\\
a_{3,2} = 0.409567801987914 & \qquad & a_{4,1} = 0.646601929740551 &\\
a_{4,2} = 0.409567801987914 & \qquad & a_{4,3} = 1 &\\
c_{2,1} = 4.198495621784201 & \qquad & c_{3,1} = 3.711590161613010 &\\
c_{3,2} = -1.787771994729384 & \qquad & c_{4,1} = 4.458898153216104 &\\
c_{4,2} = -2.024095448516552 & \qquad & c_{4,3} = -2.626700600119396 &\\
m_1 = 0.646601929740551 & \qquad & m_2 = 0.409567801987914 &\\
m_3 = 1 & \qquad & m_4 = 1 &\\
e_1 = 0 & \qquad & e_2 = 0 &\\
e_3 = 0 & \qquad & e_4 = 1 &\\
\alpha_1 = 0 & \qquad & \alpha_2 = 0 &\\
\alpha_3 = 1 & \qquad & \alpha_4 = 1 &\\
\gamma_1 = 0.515 & \qquad & \gamma_2 = 1.628546001287715 &\\
\gamma_3 = 0 & \qquad & \gamma_4 = 0
\end{aligned}\end{split}\]
RODAS-4
\[\begin{split}\begin{aligned}
\alpha_1 = 0 & \qquad & \alpha_2 = 0.386 & \qquad & \alpha_3 = 0.210 &\\
\alpha_4 = 0.630 & \qquad & \alpha_5 = 1 & \qquad & \alpha_6 = 1 &\\
\gamma_1 = 0.25 & \qquad & \gamma_2 = -0.104 & \qquad & \gamma_3 = 0.104 &\\
\gamma_4 = -0.036 & \qquad & \gamma_5 = 0 & \qquad & \gamma_6 = 0 &\\
a_{2,1} = 1.544 & \qquad & a_{3,1} = 0.946 & \qquad & a_{3,2} = 0.255 &\\
a_{4,1} = 3.314 & \qquad & a_{4,2} = 2.896 & \qquad & a_{4,3} = 0.998 &\\
a_{5,1} = 1.221 & \qquad & a_{5,2} = 6.019 & \qquad & a_{5,3} = 12.537 &\\
a_{5,4} = -0.687 & \qquad & a_{6,1} = a_{5,1} & \qquad & a_{6,2} = a_{5,2} &\\
a_{6,3} = a_{5,3} & \qquad & a_{6,4} = a_{5,4} & \qquad & a_{6,5} = 1 &\\
c_{2,1} = -5.668 & \qquad & c_{3,1} = -2.430 & \qquad & c_{3,2} = -0.206 &\\
c_{4,1} = -0.107 & \qquad & c_{4,2} = -9.594 & \qquad & c_{4,3} = -20.47 &\\
c_{5,1} = 7.496 & \qquad & c_{5,2} = -0.124 & \qquad & c_{5,3} = -34 &\\
c_{5,4} = 11.708 & \qquad & c_{6,1} = 8.083 & \qquad & c_{6,2} = -7.981 &\\
c_{6,3} = -31.521 & \qquad & c_{6,4} = 16.319 & \qquad & c_{6,5} = -6.058 &\\
m_1 = a_{5,1} & \qquad & m_2 = a_{5,2} & \qquad & m_3 = a_{5,3} &\\
m_4 = a_{5,4} & \qquad & m_5 = 1 & \qquad & m_6 = 1 &\\
e_1 = 0 & \qquad & e_2 = 0 & \qquad & e_3 = 0 &\\
e_4 = 0 & \qquad & e_5 = 0 & \qquad & e_6 = 1 &\\
\end{aligned}\end{split}\]
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:
\[\begin{split}\begin{aligned}
%y^{n+1} &=& y^n + \sum_{i=1}^s m_i k_i, \qquad
\delta y^{n+1} &=& \delta y^n + \sum_{i=1}^s m_i \ell_i\\
\nonumber
T_i &=& t^n + \alpha_i h~, %\quad Y_i =y^n + \sum_{j=1}^{i-1} a_{ij} k_j~,
\quad \delta Y_i = \delta y^n + \sum_{j=1}^{i-1} a_{ij} \ell_j\\
%A &=& \left[ \frac{1}{h \gamma} - J^T(t^n,y^n) \right]\\
%\nonumber
%A \cdot k_i &=&
% f\left( \, T_i,\, Y_i \,\right)
% + \sum_{j=1}^{i-1} \frac{c_{ij}}{h} k_j
% + h \gamma_i f_t\left(t^n,y^n\right)~,\\
\nonumber
A \cdot \ell_i &=&
J\left( \, T_i,\, Y_i \,\right)
\cdot \delta Y_i
+ \sum_{j=1}^{i-1} \frac{c_{ij}}{h} \ell_j\\
\nonumber
&& +
\left( H( t^n, y^n )\times k_i \right) \cdot \delta y^n
+ h \gamma_i J_t\left(t^n,y^n\right) \cdot \delta y^n\end{aligned}\end{split}\]
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
\[\begin{split}\begin{aligned}
\label{Ros_disc_adj}
%A &=& \left[ \frac{1}{h \gamma} - J^T(t^n,y^n) \right]\\
%\nonumber
A \cdot u_i
&=& m_i \lambda^{n+1} + \sum_{j=i+1}^s \left( a_{ji} v_j + \frac{c_{ji}}{h}
u_j \right)~,\\
\nonumber
v_i &=& J^T(T_i,Y_i)\cdot u_i~, \quad i = s,s-1,\cdots,1~,\\
\nonumber
\lambda^n &=& \lambda^{n+1} + \sum_{i=1}^s \left( H(t^n,y^n) \times
k_i\right)^T
\cdot u_i\\
\nonumber
&& + h J^T_t(t^n,y^n) \cdot \sum_{i=1}^s \gamma_i u_i+ \sum_{i=1}^s v_i\end{aligned}\end{split}\]
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).
Rosenbrock with mechanism auto-reduction
Integrator file: int/rosenbrock_autoreduce.f90
Mechanism auto-reduction (described in Lin et al. [2023]) expands
previous work by Santillana et al. [2010] and
Shen et al. [2020] to a computationally efficient implementation
in KPP, avoiding memory re-allocation, re-compile of the code, and
on-the-fly mechanism reduction based on dynamically determined
production and loss rate thresholds.
We define a threshold \(\delta\) which can be fixed (as in
Santillana et al. [2010]) or determined by the production and
loss rates of a “target species” scaled by a factor
\[\begin{aligned}
\delta = max(P_{target}, L_{target}) * \alpha_{target}.
\end{aligned}\]
For each species \(i\), the species is partitioned as “slow” iff.
\[\begin{aligned}
max(P_i, L_i) < \delta
\end{aligned}\]
if the species is partitioned as “slow”, it is solved explicitly
(decoupled from the rest of the mechanism) using a first-order
approximation. Otherwise, “fast” species are retained in the implicit
Rosenbrock solver.
Rosenbrock with H211b time stepping
Integrator file: int/rosenbrock_h211b_qssa.f90
H211b time stepping according to Söderlind [2003], as
implemented by Dreger et al. [2025].