# Solving ODE System and Sensitivity Equations `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. ```r 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 ## Returns 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 in stats vector corresponds to SUNDIALS function pattern CVodeGet, e.g. "NumRhsEvals" was obtained with CVodeGetNumRhsEvals() call. For detailed meaning of each statistics, user is invited to refer to [SUNDIALS documentation](https://computing.llnl.gov/) (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` ## 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). In this package, 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 ¶m,NumericVector &psens)` and return an integer status flag, e.g. `CV_SUCCESS`. For other possible status flags see the original [SUNDIALS documentation](https://computing.llnl.gov/) (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¶m, 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 ¶m, 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 ¶m,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 &ySv, const ivec &rootsfound, RObject ¶m,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, conststd::vector &ySv, std::vector &ySdotv, RObject ¶m,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, intiS, const vec &ySv, vec &ySdotv, RObject ¶m, 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. ## Examples ```r # 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 ¶m, 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 "), verbose=FALSE) # use 'cacheDir=""' 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 ¶m, 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 "), verbose=FALSE) # use 'cacheDir=""' 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 ¶m, 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 "), verbose=FALSE) # use 'cacheDir=""' 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 &ySv, const ivec &rootsfound, RObject ¶m, 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 "), verbose=FALSE) # use 'cacheDir=""' 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 ¶m, 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 "), 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 ¶m, 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 "), verbose=FALSE) # use 'cacheDir=""' 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 ¶m, 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 "), verbose=FALSE) # use 'cacheDir=""' 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, "]"))) ```

r2cvodes function