Description Usage Arguments Details Value Examples
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31  r2cvodes(
yv,
times,
frhs,
param = NULL,
tstop = as.numeric(c()),
abstol = 1e08,
reltol = 1e08,
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
)

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 
param 
any R object (default 
tstop 
const numeric scalar (default 
abstol 
const numeric scalar (default 
reltol 
const double (default 
integrator 
integer scalar (default 
maxord 
const integer scalar (default 
maxsteps 
const integer scalar (default 
hin 
const numeric scalar (default 
hmax 
const numeric scalar (default 
hmin 
const numeric scalar (default 
constraints 
const numeric vector (default 
fjac 
R function of XPtr pointer to Rcpp function (default 
nz 
const integer scalar (default 
rmumps_perm 
integer scalar (default 
nroot 
const integer scalar (default 
froot 
R function of XPtr pointer to Rcpp function (default 
fevent 
R function of XPtr pointer to Rcpp function (default 
Ns 
const integer scalar (default 
psens 
numeric vector (default 
sens_init 
numeric matrix (default 
psens_bar 
numeric vector (default 
psens_list 
const integer vector (default 
fsens 
R function of XPtr pointer to Rcpp function (default 
fsens1 
R function of XPtr pointer to Rcpp function (default 
sens_method 
integer scalar (default

errconS 
constant logical scalar (default 
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 ¶m, 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 inplace 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 inplace 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 1based, 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 inplace 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 inplace 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 herebefore. 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 ¶m, NumericVector &psens)
The new state vector can be stored inplace 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 inplace 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 ith 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 (ith vector is in ith 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 ¶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 mth 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 inplace 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 onebyone 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 ¶m, NumericVector &psens, const vec &tmp1v, const vec &tmp2v)
The result, i.e. s'[iS] is to be stored inplace in ySdotv
vector. This function returns a status flag.
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:
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);
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 herebefore.
sensitivity 3D array with dimensions Neq
x Nt
x Ns
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211  # Ex.1. Solve a scalar ODE describing exponential transition form 0 to 1
# y'=a*(y1), 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*(y1)
# 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: 1exp(a*t)
stopifnot(diff(range(1exp(p$a*ti)  res_exp)) < 1.e6)
# 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="using namespace arma;\n#include <r2sundials.h>", cacheDir="lib", verbose=FALSE)
# 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.e14)
# 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 5th 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="using namespace arma;\n#include <r2sundials.h>", cacheDir="lib", verbose=FALSE)
# 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] = tp["tmax"]; // t==p["tmax"]
return(0);
}
', depends=c("RcppArmadillo","r2sundials","rmumps"),
includes="using namespace arma;\n#include <r2sundials.h>", cacheDir="lib", verbose=FALSE)
# 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 ¶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="using namespace arma;\n#include <r2sundials.h>", cacheDir="lib", verbose=FALSE)
# 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="using namespace arma;\n#include <r2sundials.h>", cacheDir="lib", 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="using namespace arma;\n#include <r2sundials.h>", cacheDir="lib", verbose=FALSE)
# 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="using namespace arma;\n#include <r2sundials.h>", cacheDir="lib", verbose=FALSE)
# 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, "]")))

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.