我有一个ODE，我想用从R的deSolve包中调用的编译C代码来解决它.讨论中的ODE是一个指数衰减模型(y'=-d*exp(g*time)*y):
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)
