R: Solving ODE System and Sensitivity Equations

r2cvodes {r2sundials}

R Documentation

Solving ODE System and Sensitivity Equations

Description

r2cvodes sets up necessary structures and calls cvodes() from SUNDIALS library to solve user defined ODE system y' = f(t, y, p), y(t0) = y0, where p is a constant parameter vector. If requested, corresponding forward sensitivity equations s'[i] = df/dy s[i] + df/dp[i], s[i](t0) = dy0(p)/dp[i] (here s[i](t)=dy(t)/dp[i]) can be solved simultaneously with the original ODE system. Root finding and proceeding can be defined as well.

Usage

r2cvodes(
  yv,
  times,
  frhs,
  param = NULL,
  tstop = as.numeric(c()),
  abstol = 1e-08,
  reltol = 1e-08,
  integrator = as.integer(c()),
  maxord = 0L,
  maxsteps = 0L,
  hin = 0,
  hmax = 0,
  hmin = 0,
  constraints = as.numeric(c()),
  fjac = NULL,
  nz = 0L,
  rmumps_perm = as.integer(c()),
  nroot = 0L,
  froot = NULL,
  fevent = NULL,
  Ns = 0L,
  psens = as.numeric(c()),
  sens_init = as.numeric(c()),
  psens_bar = as.numeric(c()),
  psens_list = as.integer(c()),
  fsens = NULL,
  fsens1 = NULL,
  sens_method = as.integer(c()),
  errconS = TRUE
)

Arguments

`yv`	const numeric vector, initial values of state vector (`y0`). Its length defines the number of equations in the ODE system and is refered hereafter as 'Neq'.
`times`	const numeric vector, time point values at which the solution is stored
`frhs`	R function or XPtr pointer to Rcpp function, defines rhs `f(t, y, p)`, i.e. returns (or fills in place) first derivative vector of length `Neq` (cf. details for argument list)
`param`	any R object (default `NULL`), parameter object (can be a vector, list, environment, ...) passed through to user defined functions namely to `frhs`. It can be useful to get the code more readable as named parameters are authorized. (Starting from this one, the following parameters are optional)
`tstop`	const numeric scalar (default `as.numeric(c())`) if non empty, defines a time point beyond which the calculation is not made (e.g. out of definition domain). If a vector is given, only the first value is used.
`abstol`	const numeric scalar (default `1.e-8`), absolute tolerance in ODE solving (used for values which are close to 0)
`reltol`	const double (default `1.e-8`), relative tolerance in ODE solving (used for values which are not so close to 0)
`integrator`	integer scalar (default `as.integer(c())`), defines which integration scheme should be used by cvodes: implicit (`CV_BDF`) or explicit (`CV_ADAMS`). Default empty vector is equivalent to `CV_BDF`. Constants `CV_BDF` and `CV_ADAMS` (as other mentioned constants) are defined by the package and are available for user.
`maxord`	const integer scalar (default `0`), defines maximal order of time scheme. Default 0 values is equivalent to 5 for `integrator=CV_BDF` or to 12 for `integrator=CV_ADAMS`
`maxsteps`	const integer scalar (default `0`), maximum of internal steps before reaching next time point from 'times'. Default 0 is equivalent to 500.
`hin`	const numeric scalar (default `0`), value of the initial step size to be attempted. Default 0.0 corresponds to `cvodes`' default value which is internally estimated at `t0`.
`hmax`	const numeric scalar (default `0`), maximum absolute value of the time step size (>= 0.0). The default 0 value is equivalent to Inf.
`hmin`	const numeric scalar (default `0`), minimum absolute step size (>= 0.0)
`constraints`	const numeric vector (default `as.numeric(c())`), if non empty, defines constraint flags. If constraints[i] is 0 then no constraint is imposed on `y[i]`; 1.0 then `y[i]` will be constrained to be `y[i]` >= 0.0; -1.0 then `y[i]` will be constrained to be `y[i]` =< 0.0.; 2.0 then `y[i]` will be constrained to be `y[i]` > 0.0.; -2.0 then `y[i]` will be constrained to be `y[i]` < 0.0. If a time step respecting hmin and avoiding constraint violation cannot be found, an error is triggered.
`fjac`	R function of XPtr pointer to Rcpp function (default `NULL`), users supplied function that calculates Jacobian matrix df/dy which can be dense or sparse (cf. details for parameter list). If `fjac` is not supplied, a SUNDIALS internal approximation to Jacobian is used when needed.
`nz`	const integer scalar (default `0`), number of non zero elements in Jacobian. If > 0, this parameter triggers a sparse Jacobian usage. In this case, the previous parameter `fjac`, must be a function (or a pointer to a function) with appropriate parameter list (cf. details) filling a sparse matrix. It is advised to have `nz` accounting not only for non zeros in the Jacobian itself but also for diagonal terms even if they are zeros. In this case, a supplementary memory allocation can be avoided. If a sparse Jacobian is used, corresponding sparse linear system is solved with the help of rmumps package.
`rmumps_perm`	integer scalar (default `as.integer(c())`), defines permutation method that will be used by rmumps during symbolic analysis before solving sparse linear systems. Default value is equivalent to `RMUMPS_PERM_AUTO`. Possible values (defined in rmumps package) are `RMUMPS_PERM_AMF`, `RMUMPS_PERM_AMD`, `RMUMPS_PERM_AUTO`, `RMUMPS_PERM_QAMD`, `RMUMPS_PERM_PORD`, `RMUMPS_PERM_METIS`, `RMUMPS_PERM_SCOTCH`. An appropriate choice of permutation type can reduce memory requirements as well as calculation speed for sparse system solving. If a vector is supplied, only the first value is read. If `nz=0`, this parameter is not used.
`nroot`	const integer scalar (default `0`), defines number of roots that user wishes to track.
`froot`	R function of XPtr pointer to Rcpp function (default `NULL`), user defined function calculating root vector of length `nroot`. When at least one component of this vector crosses 0, a root finding event is triggered and corresponding root handling function is called (cf. `fevent`). This defines a very flexible root tracking mechanism. User can define a root based both on time values and state vector (y).
`fevent`	R function of XPtr pointer to Rcpp function (default `NULL`), defines function for root proceeding. When a root is encountered, this function is called. Its return value defines what to do with this root. User's reaction (i.e. the return value of `fevent`) can be one of three types: `R2SUNDIALS_EVENT_IGNORE` - do nothing (can be helpful if, for example, 0 is crossed not in pertinent sens); `R2SUNDIALS_EVENT_HOLD` - the time point at which the root happened as well as the corresponding `rootsfound` vector are added to the root matrix delivered with output attributes (cf. details). This time point is also added to the ODE solution which can lead to a new time point, originally absent in 'times' vector; `R2SUNDIALS_EVENT_STOP` - stops the ODE solving. If it happens before reaching the last value in 'times', the ODE solution can have less time points than initially defined in 'times' parameter. `fevent` is allowed to modify the vector state y in an arbitrary way. It can be helpful for modeling some controlling events for example adding some compound to chemical mix or modifying speed vector if an obstacle is hitted. If such modification takes place, the ODE system is restarted to handle in appropriate way induced discontinuities.
`Ns`	const integer scalar (default `0`), number of parameters for which forward sensitivity system has to be solved.
`psens`	numeric vector (default `as.numeric(c())`), if not empty, defines a vector of parameters for which (or for some components of which) forward sensitivity system has to be solved. If its length > `Ns`, then vector `psens_list` of length `Ns` must define which are the `Ns` components of `psens` with respect to which sensitivities are to be computed. Note that parameters used in sensitivity calculations must be put in `psens` vector (and not somewhere in 'param' object) only if user wish to rely on SUNDIALS internal procedure for estimation of sensitivity rhs. If user provides its own function for sensitivity rhs, he is free to use either of 'param' or `psens` for passing and accessing sensitivity parameters.
`sens_init`	numeric matrix (default `as.numeric(c())`), matrix with `Neq` rows and `Ns` columns defining initial condition for sensititivy system. Default value is equivalent to `Neq` x `Ns` zero matrix which means that `y0` is not dependent on `p`.
`psens_bar`	numeric vector (default `as.numeric(c())`), a vector of `Ns` positive scaling factors. The values provided in this vector, describe orders of magnitude for `psens` components. Having such estimations can help to increase the accuracy for sensitivity vector calculations or for the right-hand side calculation for the sensitivity systems based on the internal difference-quotient function. If not defined, it is equivalent to set all `psens_bar` components to 1
`psens_list`	const integer vector (default `as.integer(c())`), if `length(psens) > Ns`, this index vector can be used to indicate which are components of `psens` that are concerned with sensitivity calculations. Default value is valid only when `length(psens)==Ns`, i.e. `psens_list` is not used. Values in `psens_list` are 1-based (as regular integer indexes in R)
`fsens`	R function of XPtr pointer to Rcpp function (default `NULL`), user defined sensitivity rhs for all `psens` at once.
`fsens1`	R function of XPtr pointer to Rcpp function (default `NULL`), user defined sensitivity rhs for `is`-th `psens` component (`is` is passed as input parameter to `fsens1`, cf. details). Only one of `fsens` or `fsens1` can be supplied. If none of them is supplied, SUNDIALS' internal difference-quotient function is used. In this case, sensitivity parameters must be passed only through `psens` argument.
`sens_method`	integer scalar (default `as.integer(c())`), constant defining the method used for solving sensitivity system. Allowed values are `CV_SIMULTANEOUS` or `CV_STAGGERED`. Default value is equivalent to `CV_SIMULTANEOUS`. `CV_SIMULTANEOUS` means that the state and sensitivity variables are corrected at the same time. `CV_STAGGERED` means that the correction step for the sensitivity variables takes place at the same time for all sensitivity equations, but only after the correction of the state variables has converged and the state variables have passed the local error test
`errconS`	constant logical scalar (default `TRUE`), specifies whether sensitivity variables are to be included (`TRUE`) or not (`FALSE`) in the error control mechanism

Details

The package r2sundials was designed to avoid as much as possible memory reallocation in callback functions (frhs and others). C++ variants of these functions are fully compliant with this design principle. While R counterparts are not (as per R design). Here, we define callback function interfaces that user has to abide to. Pointers to C++ variants to be passed to r2cvodes() can be obtained with the help of RcppXPtrUtils. See examples for illustrations of such use.
Right hand side function frhs provided by user calculates derivative vector y'. This function can be defined as classical R function or a Rcpp/RcppArmadillo function. In the first case, it must have the following list of input arguments
frhs(t, y, param, psens)
and return a derivative vector of length Neq. Here t is time point (numeric scalar), y current state vector (numeric vector of length Neq), param and psens are passed through from r2cvodes() arguments.
In the C++ case, it is defined as
int (*frhs)(double t, const vec &y, vec &ydot, RObject &param, NumericVector &psens)
and return an integer status flag, e.g. CV_SUCCESS. For other possible status flags see the original SUNDIALS documentation (search for cvs_guide.pdf). The derivatives are stored in-place in ydot vector. See examples section for a usage sample.
fjac is a function calculating Jacobian matrix. Its definition varies depending on 1) kind of used Jacobian: dense or sparse and 2) on programming language used: R or C++ (i.e. Rcpp/RcppArmadillo).

For dense Jacobian calculated in R, the arguments are:
fjac(t, y, ydot, param, psens)
and the expected return value is Neq x Neq dense Jacobian matrix df/dy.
For dense Jacobian calculated in C++ the definition is following:
int (*fjac)(double t, const vec &y, const vec &ydot, mat &J, RObject &param, NumericVector &psens, vec &tmp1, vec &tmp2, vec &tmp3)
It must return a status flag. The resulting Jacobian is stored in-place in the the matrix J. Auxiliary vectors tmp1 to tmp3 are of length Neq and are made available for intermediate storage thus avoiding memory reallocation at each call to fjac().
For sparse Jacobian calculated in R, the arguments are:
fjac(t, yv, ydotv, param, psens)
The return value is a list with fields i (row indices), p (column pointers) and v (matrix values) defining the content of sparse Jacobian in CSC (condensed sparse column) format. The values stored in i and p vectors are supposed to be 1-based, as it is common in R language.
For sparse Jacobian calculated in C++ the definition is following:
int (*fjac)(double t, const vec &y, const vec &ydot, uvec &i, uvec &p, vec &v, int n, int nz, RObject &param, NumericVector &psens, vec &tmp1, vec &tmp2, vec &tmp3)
here n=Neq, nz is passed through from r2cvodes() arguments. The resulting sparse Jacobian is stored in-place in vectors i, p, v corresponding to the CSC (Compressed Sparse Column) format. Their respective dimensions are nz, n+1 and nz. The values stored in i and p must be 0 based as per usage in C++. The return value is a status flag.

froot calculates a root vector, i.e. a vector whose components are tracked for 0 crossing during the time course in ODE solving. If written in R, its call follows the following pattern:
froot(t, y, param, psens)
and it must return a numeric vector of length 'nroot'. If written in C++, it is defined as
int (*froot)(double t, const vec &y, vec &vroot, RObject &param, NumericVector &psens)
The tracked values are stored in-place in vroot. The returned value is a status flag.
fevent handles the event of root finding. If written in R, the calling pattern is
fevent(t, yvec, Ns, ySm, rootsfound, param, psens)
and the return value is a list with named component "ynew", "flag" and "sens_init". The last item is required only when Ns > 0. Current value of sensitivity matrix dy/dp can be found in parameter ySm. Integer vector 'rootsfound' of length 'nroot' provides information on vroot components that triggered the root event. If rootsfound[i] != 0, it means that vroot[i] is a root otherwise it is not. Moreover, the sign of rootsfound[i] is meaningful. If rootsfound[i] > 0 then vroot[i] is increasing at 0 crossing. Respectively, if rootsfound[i] < 0 then vroot[i] is decreasing. The vector 'ynew' in the output list can define a new state vector after event handling (for example, an abrupt change in velocity direction and/or magnitude after an obstacle hit). The field 'flag' in the output list is authorized to take only three values: R2SUNDIALS_EVENT_IGNORE, R2SUNDIALS_EVENT_HOLD and R2SUNDIALS_EVENT_STOP described here-before. The matrix sens_init is used for a possible restarting of sensitivity calculation. It must contain a derivative matrix dynew/dp of size Neq x Ns. If this field is absent (NULL) then a zero matrix is assumed.
If written in C++, this function is defined as
int (*fevent)(double t, const vec &y, vec &ynew, int Ns, std::vector<vec> &ySv, const ivec &rootsfound, RObject &param, NumericVector &psens)
The new state vector can be stored in-place in ynew and the status flag indicating what to do with this event is the return value. If sensitivity calculation is going on (i.e. Ns > 0), the current value of sensitivity vectors can be found in ySv and their new values can be stored in-place in the same parameter. Note that if ynew is different from the vale of y when the root was found, the ODE is restarted from this time point to handle correctly the discontinuity. As a result, there will be two columns corresponding to the same time point: one with the state vector at root finding and one with ynew values. The same is true for sensitivity output, if it is part of the problem.
fsens calculates rhs for sensitivity system. If written in R, it must be defined as
fsens(Ns, t, y, ydot, ySm, param, psens)
and return a dataframe in which i-th column correspond to s'[i] sensitivity derivative vector. Among other parameters, it receives ySm which is a Neq x Ns matrix having the current values of sensitivity vector (i-th vector is in i-th column).
If written in C++, it has to be defined as
int (*fsens)(int Ns, double t, const vec &y, const vec &ydot, const std::vector<vec> &ySv, std::vector<vec> &ySdotv, RObject &param, NumericVector &psens, const vec &tmp1v, const vec &tmp2v)
Note a slight difference in the input parameters compared with the R counterpart. Here ySv plays the role of ySm and is not a matrix but a vector of Armadillo vectors. To access m-th component of s[i], one can simply do ySv[i][m] and the whole s[i] is selected as ySv[i]. Such data structure was retained to keep as low as possible new memory reallocation. The resulting sensitivity derivatives are to be stored in-place in ySdotv according to the same data organization scheme as in ySv. This function returns a status flag.
fsens1 does the same as fsens but provides derivatives of sensitivity vectors on one-by-one basis. This second form is provided for user's convenience as in some cases the code can become more readable if it calculates only one vector s'[i] at a time. If written in R, this function has to be defined as
fsens1(Ns, t, y, iS, ydot, yS, param, psens)
here iS is the index of calculated vector s'[iS] and yS contains the current value of s[iS]. If written in C++, this functions has to be defined as
int (*fsens1)(int Ns, double t, const vec &yv, const vec &ydotv, int iS, const vec &ySv, vec &ySdotv, RObject &param, NumericVector &psens, const vec &tmp1v, const vec &tmp2v)
The result, i.e. s'[iS] is to be stored in-place in ySdotv vector. This function returns a status flag.

Value

numeric matrix, ODE solution where each column corresponds to a state vector at a given time point. The columns (their number is refered to as Nt) are named by time points while the rows heritates the names from yv. If no names are found in yv, the rows are simply named 'V1', 'V2' and so on. After a normal execution and without root handling, column number is equal to the length of times. However, if root handling is used, it can add or remove some time points from times. So the user must not assume that column number of output matrix is equal to length(times). Instead, actual number of time points for which the solution was calculated can be retrieved from an attribute named "times". Moreover, several attributes are defined in the returned matrix. We have mentioned "times", the others are:

stats: Some statistics about various events that could happen during cvodes run like the number of rhs calls, Jacobian calls, number of internal time steps, failure number and so on. Any component <name> in stats vector corresponds to SUNDIALS function pattern CVodeGet<name>, e.g. "NumRhsEvals" was obtained with CVodeGetNumRhsEvals() call. For detailed meaning of each statistics, user is invited to refer to SUNDIALS documentation (search for cvs_guide.pdf);
roots: matrix with row number nroot+1 and column number equal to number of roots found by the cvodes() and retained by the user. Each column is a composite vector made of time point and rootsfound vector described here-before.
sens: sensitivity 3D array with dimensions Neq x Nt x Ns

Examples

# Ex.1. Solve a scalar ODE describing exponential transition form 0 to 1
# y'=-a*(y-1), y(0)=0, a is a parameter that we arbitrary choose to be 2.
# define rhs function (here in R).
frhs_exp=function(t, y, p, psens) -p$a*(y-1)
# define parameter list
p=list(a=2)
# define time grid from 0 to 5 (arbitrary units)
ti=seq(0, 5, length.out=101)
# define initial state vector
y0=0
# we are set for a very simple r2cvodes() call
res_exp=r2sundials::r2cvodes(y0, ti, frhs_exp, param=p)
# compare the result to theoretical values: 1-exp(-a*t)
stopifnot(diff(range(1-exp(-p$a*ti) - res_exp)) < 1.e-6)

# Ex. 2. Same problem but frhs is written in C++
library(RcppXPtrUtils)
ptr_exp=cppXPtr(code='
int rhs_exp(double t, const vec &y, vec &ydot, RObject &param, NumericVector &psens) {
  NumericVector p(param);
  ydot[0] = -p["a"]*(y[0]-1);
  return(CV_SUCCESS);
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
 includes=c("// [[Rcpp::plugins(cpp14)]]", "using namespace arma;", "#include <r2sundials.h>"), 
 verbose=FALSE) # use 'cacheDir="<yourdir>"' to keep compilation results between R sessions.
# For ease of use in C++, we convert param to a numeric vector instead of a list.
pv=c(a=p$a)
# new call to r2cvodes() with XPtr pointer ptr_exp.
res_exp2=r2sundials::r2cvodes(y0, ti, ptr_exp, param=pv)
stopifnot(diff(range(res_exp2 - res_exp)) < 1.e-14)

# Ex.3. Bouncing ball simulation.
# A ball falls from a height y=5 m with initial vertical speed vy=0 m/s
# and horizontal speed vx=1 m/s. The forces acting on the ball are 1) the gravity
# (g=9.81 m/s^2) and 2) air resistance f=-k_r*v (k_r=0.1 N*s/m).
# When the ball hits the ground, it bounces instantly retaining k=0.9 part
# of its vertical and horizontal speed. At the bounce, the vertical speed
# changes its sign to the opposite while horizontal speed keeps the original sign.
# Simulation should stop after the 5-th bounce or at tmax=10 s which ever comes first.
# This example illustrates usage of root finding and handling.
# We give here an example of callback functions for root handling in C++.
yv=c(x=0, y=5, vx=1, vy=0) # initial state vector
pv=c(g=9.81, k_r=0.1, k=0.9, nbounce=5, tmax=10) # parameter vector
ti=seq(0, 20, length.out=201L) # time grid

# rhs
ptr_ball=cppXPtr(code='
int rhs_ball(double t, const vec &y, vec &ydot, RObject &param, NumericVector &psens) {
  NumericVector p(param);
  ydot[0] = y[2]; // dx/dt=vx
  ydot[1] = y[3]; // dy/dt=vy
  ydot[2] = -p["k_r"]*y[2]; // dvx/dt= -k_r*vx
  ydot[3] = -p["g"] - p["k_r"]*y[3]; // dvy/dt=-g -k_r*vy
  return(CV_SUCCESS);
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
 includes=c("// [[Rcpp::plugins(cpp14)]]", "using namespace arma;", "#include <r2sundials.h>"), 
 verbose=FALSE) # use 'cacheDir="<yourdir>"' to keep compilation results between R sessions.

# root function
ptr_ball_root=cppXPtr(code='
int root_ball(double t, const vec &y, vec &vroot, RObject &param, NumericVector &psens) {
  NumericVector p(param);
  vroot[0] = y[1]; // y==0
  vroot[1] = t-p["tmax"]; // t==p["tmax"]
  return(0);
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
 includes=c("// [[Rcpp::plugins(cpp14)]]", "using namespace arma;", "#include <r2sundials.h>"), 
 verbose=FALSE) # use 'cacheDir="<yourdir>"' to keep compilation results between R sessions.

# event handler function
ptr_ball_event=cppXPtr(code='
int event_ball(double t, const vec &y, vec &ynew, int Ns, std::vector<vec> &ySv,
          const ivec &rootsfound, RObject &param, NumericVector &psens) {
  NumericVector p(param);
  static int nbounce=0;
  if (rootsfound[1] != 0) // time is out
     return(R2SUNDIALS_EVENT_STOP);
  if (rootsfound[0] > 0)
    // cross 0 in ascending sens, can happen when y < 0 in limits of abstol
    return(R2SUNDIALS_EVENT_IGNORE);
  ynew=y;
  if (++nbounce < p["nbounce"]) {
    // here nbounce=1:4
    ynew[2] *= p["k"]; // horizontal speed is lowered
    ynew[3] *= -p["k"]; // vertical speed is lowered and reflected
    return(R2SUNDIALS_EVENT_HOLD);
  } else {
    // here nbounce=5
    nbounce=0; // reinit counter for possible next calls to r2cvodes
    return(R2SUNDIALS_EVENT_STOP);
  }
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
 includes=c("// [[Rcpp::plugins(cpp14)]]", "using namespace arma;", "#include <r2sundials.h>"), 
 verbose=FALSE) # use 'cacheDir="<yourdir>"' to keep compilation results between R sessions.

# ODE solving and plotting
res_ball <- r2sundials::r2cvodes(yv, ti, ptr_ball, param=pv, nroot=2L,
  froot=ptr_ball_root, fevent=ptr_ball_event)
plot(res_ball["x",], res_ball["y",], xlab="X [m]", ylab="Y [m]",
  t="l", main="Bouncing ball simulation")

# Ex.4. Robertson chemical reactions
# This example is often used as an illustration and a benchmark for stiff ODE.
# We will demonstrate here:
# • how to use sparse Jacobian (not really meaningful for 3x3 sytem but just to give a hint);
# • how to make sensitivity calculations with user provided rhs.
#
# Let simulate the following chemical system of 3 compounds y1, y2 and y3
#  y1' = -k1*y1 + k3*y2*y3
#  y2' =  k1*y1 - k2*y2*y2 - k3*y2*y3
#  y3' =  k2*y2*y2
# Jacobian df/dy is
# 
# | -k1 |      k3*y3       |  k3*y2 |
# |-----+------------------+--------|
# |  k1 | -2*k2*y2 - k3*y3 | -k3*y2 |
# |-----+------------------+--------|
# |  0  |     2*k2*y2      |   0    |


yv <- c(y1=1, y2=0, y3=0) # initial values
pv <- c(k1 = 0.04, k2 = 3e7, k3 = 1e4) # parameter vector
ti=10^(seq(from = -5, to = 11, by = 0.1)) # exponential time grid

# pointer to rhs function
ptr_rob=cppXPtr(code='
int rhs_rob(double t, const vec &y, vec &ydot, RObject &param, NumericVector &psens) {
  NumericVector p(param);
  ydot[0] = -p["k1"]*y[0] + p["k3"]*y[1]*y[2];
  ydot[2] = p["k2"]*y[1]*y[1];
  ydot[1] = -ydot[0] - ydot[2];
  return(CV_SUCCESS);
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
 includes=c("// [[Rcpp::plugins(cpp14)]]", "using namespace arma;", "#include <r2sundials.h>"), 
 verbose=FALSE) 
# pointer to sparse jacobian function
ptr_rob_jacsp=cppXPtr(code='
int spjac_rob(double t, const vec &y, const vec &ydot, uvec &ir, uvec &pj, vec &v, int n, int nz,
              RObject &param, NumericVector &psens, vec &tmp1, vec &tmp2, vec &tmp3) {
  if (nz < 8)
    stop("spjac_rob: not enough room for non zeros, must have at least 8, instead got %d", nz);
  NumericVector p(param);
  int i=0;
  pj[0] = 0; // init pj
  // first column
  ir[i] = 0;
  v[i++] = -p["k1"];
  ir[i] = 1;
  v[i++] = p["k1"]; 
  pj[1] = i;
  // second column
  ir[i] = 0;
  v[i++] = p["k3"]*y[2];
  ir[i] = 1;
  v[i++] = -p["k3"]*y[2]-2*p["k2"]*y[1];
  ir[i] = 2;
  v[i++] = 2*p["k2"]*y[1];
  pj[2] = i;
  // third column
  ir[i] = 0;
  v[i++] = p["k3"]*y[1];
  ir[i] = 1;
  v[i++] = -p["k3"]*y[1];
  ir[i] = 2;
  v[i++] = 0; // just to have the main diagonal fully in Jacobian
  pj[3] = i;
  return(0);
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
 includes=c("// [[Rcpp::plugins(cpp14)]]", "using namespace arma;", "#include <r2sundials.h>"), 
 verbose=FALSE) # use 'cacheDir="<yourdir>"' to keep compilation results between R sessions.
# pointer to sensitivity rhs function
ptr_rob_sens1=cppXPtr(code='
int sens_rob1(int Ns, double t, const vec &y, const vec &ydot, int iS, const vec &yS, vec &ySdot,
              RObject &param, NumericVector &psens, vec &tmp1, vec &tmp2) {
  // calculate (df /dy)s_i(t) + (df /dp_i) for i = iS
  NumericVector p(param);
  // (df/dy)s_i(t)
  ySdot[0] = -p["k1"]*yS[0] + p["k3"]*y[2]*yS[1] + p["k3"]*y[1]*yS[2];
  ySdot[1] = p["k1"]*yS[0] - (p["k3"]*y[2]+2*p["k2"]*y[1])*yS[1] - p["k3"]*y[1]*yS[2]; 
  ySdot[2] = 2*p["k2"]*y[1]*yS[1];
  // + (df/dp_i)
  switch(iS) {
    case 0:
      ySdot[0] -= y[0];
      ySdot[1] += y[0];
      break;
    case 1:
      ySdot[1] -= y[1]*y[1];
      ySdot[2] += y[1]*y[1];
      break;
    case 2:
      ySdot[0] += y[1]*y[2];
      ySdot[1] -= y[1]*y[2];
  }
  return(CV_SUCCESS);
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
 includes=c("// [[Rcpp::plugins(cpp14)]]", "using namespace arma;", "#include <r2sundials.h>"), 
 verbose=FALSE) # use 'cacheDir="<yourdir>"' to keep compilation results between R sessions.
# Note that we don't use psens param for sensitivity calculations as we provide our own fsens1.
res_rob <- r2sundials::r2cvodes(yv, ti, ptr_rob, param=pv, nz=8, fjac=ptr_rob_jacsp, Ns=3,
                              fsens1=ptr_rob_sens1)
# plot ODE solution
layout(t(1:3)) # three sublots in a row
for (i in 1:3)
   plot(ti, res_rob[i,], log="x", t="l", xlab="Time", ylab=names(yv)[i])
# plot sensitivities
layout(matrix(1:9, nrow=3)) # 9 subplots in a square
for (j in 1:3) # run through pv
   for (i in 1:3) # run through y
      plot(ti, attr(res_rob, "sens")[i,,j], log="x", t="l", xlab="Time",
           ylab=parse(text=paste0("partialdiff*y[", i, "]/partialdiff*k[", j, "]")))

[Package r2sundials version 6.5.0-5 Index]