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)
Compiled code does not give different results to deSolve models implemented in R, except potential rounding errors within the limits of
atol
andrtol
.The reasons of the differences in the original post where two errors in the code. One can correct it as follows:
static double
asparms[3];
instead ofparms[4]
t
in derivs is a pointer, i.e.*t
so that the code reads as:
Here the comparison between the two simulations, somewhat adapted and generalized:
Note: set
.dll
or.so
according to your OS, or detect it with.Platform$dynlib.ext
.