# 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)
# \donttest{
# 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, "]")))
# }
Run the code above in your browser using DataLab