0
184
ViewsCompiled C ODE gives different results to R's using deSolve
I have an ODE which I would like to solve using compiled C code called from R's deSolve package. The ODE in question is I an exponential decay model (y'=-d* exp(g* time)*y): But running the compiled code from within R gives different results to R's native deSolve. It's as is there they are flipped 180º. What's going on?
C code implementation
/* file testODE.c */
#include <R.h>
static double parms[4];
#define C parms[0] /* left here on purpose */
#define d parms[1]
#define g parms[2]
/* initializer */
void initmod(void (* odeparms)(int *, double *))
{
int N=3;
odeparms(&N, parms);
}
/* Derivatives and 1 output variable */
void derivs (int *neq, double t, double *y, double *ydot,
double *yout, int *ip)
{
// if (ip[0] <1) error("nout should be at least 1");
ydot[0] = -d*exp(-g*t)*y[0];
}
/* END file testODEod.c */
R implementation - Native deSolve
testODE <- function(time_space, initial_contamination, parameters){
with(
as.list(c(initial_contamination, parameters)),{
dContamination <- -d*exp(-g*time_space)*Contamination
return(list(dContamination))
}
)
}
parameters <- c(C = -8/3, d = -10, g = 28)
Y=c(y=1200)
times <- seq(0, 6, by = 0.01)
initial_contamination=c(Contamination=1200)
out <- ode(initial_contamination, times, testODE, parameters, method = "radau",atol = 1e-4, rtol = 1e-4)
plot(out)
R implementation - Run compiled code from deSolve
library(deSolve)
library(scatterplot3d)
dyn.load("Code/testODE.so")
Y <-c(y1=initial_contamination) ;
out <- ode(Y, times, func = "derivs", parms = parameters,
dllname = "testODE", initfunc = "initmod")
plot(out)
1 answers
Answer question0
Compiled code does not give different results to deSolve models implemented in R, except potential rounding errors within the limits of atoland rtol.
The reasons of the differences in the original post where two errors in the code. One can correct it as follows:
- Declare
static doubleasparms[3];instead ofparms[4] - Time
tin derivs is a pointer, i.e.*t
so that the code reads as:
/* file testODE.c */
#include <R.h>
#include <math.h>
static double parms[3];
#define C parms[0] /* left here on purpose */
#define d parms[1]
#define g parms[2]
/* initializer */
void initmod(void (* odeparms)(int *, double *)) {
int N=3;
odeparms(&N, parms);
}
/* Derivatives and 1 output variable */
void derivs (int *neq, double *t, double *y, double *ydot,
double *yout, int *ip) {
ydot[0] = -d * exp(-g * *t) * y[0];
}
Here the comparison between the two simulations, somewhat adapted and generalized:
library(deSolve)
testODE <- function(t, y, parameters){
with(
as.list(c(y, parameters)),{
dContamination <- -d * exp(-g * t) * contamination
return(list(dContamination))
}
)
}
system("R CMD SHLIB testODE.c")
dyn.load("testODE.dll")
parameters <- c(c = -8/3, d = -10, g = 28)
Y <- c(contamination = 1200)
times <- seq(0, 6, by = 0.01)
out1 <- ode(Y, times, testODE,
parms = parameters, method = "radau", atol = 1e-4, rtol = 1e-4)
out2 <- ode(Y, times, func = "derivs", dllname = "testODE", initfunc = "initmod",
parms = parameters, method = "radau", atol = 1e-4, rtol = 1e-4)
plot(out1, out2) # no visible difference
summary(out1 - out2) # differences should be (close to) zero
dyn.unload("testODE.dll") # always unload before editing .c file !!
Note: set.dll or .so according to your OS, or detect it with .Platform$dynlib.ext.