0
180
ViewsC ODE compilado da resultados diferentes a R usando deSolve
Tengo un ODE que me gustaría resolver usando un código C compilado llamado desde el paquete deSolve de R. La ODE en cuestión es un modelo de decaimiento exponencial (y'=-d* exp(g* time)*y): pero ejecutar el código compilado desde dentro de R da resultados diferentes a la deSolve nativa de R. Es como está ahí están volteadas 180º. ¿Qué está sucediendo?
implementación de código C
/* file testODE.c */ #include <Rh> 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 */Implementación de R - deSolve nativo
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)Implementación de R: ejecute código compilado desde 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
El código compilado no da resultados diferentes a los modelos deSolve implementados en R , excepto posibles errores de redondeo dentro de los límites de atol y rtol .
Las razones de las diferencias en la publicación original fueron dos errores en el código. Se puede corregir de la siguiente manera:
- Declare
static doublecomoparms[3];en lugar deparms[4] - El tiempo
ten derivadas es un puntero, es decir,*t
para que el código se lea como:
/* file testODE.c */ #include <Rh> #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]; }Aquí la comparación entre las dos simulaciones, algo adaptada y generalizada:
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 !! 
Nota: configure .dll o .so según su sistema operativo, o detéctelo con .Platform$dynlib.ext .