使用 deSolve 编译的 C ODE 给出与 R 不同的结果

Question

我有一个 ODE，我想使用从 R 的 deSolve 包调用的编译 C 代码来求解它。有问题的 ODE 是一个指数衰减模型 (y'=-d* exp(g* time)*y):\n但是从 R 中运行编译后的代码会为 R 的本机 deSolve 提供不同的结果。就像它们被翻转了 180\xc2\xba 一样。这是怎么回事？

\n

C代码实现

\n

/* file testODE.c */\n#include <R.h>\nstatic double parms[4];\n#define C parms[0] /* left here on purpose */\n#define d parms[1]\n#define g parms[2]\n\n/* initializer  */\nvoid initmod(void (* odeparms)(int *, double *))\n{\n  int N=3;\n  odeparms(&N, parms);\n}\n\n/* Derivatives and 1 output variable */\nvoid derivs (int *neq, double t, double *y, double *ydot,\n             double *yout, int *ip)\n{\n  // if (ip[0] <1) error("nout should be at least 1");\n  ydot[0] = -d*exp(-g*t)*y[0];\n}\n/* END file testODEod.c */\n

\n

R 实现 - Native deSolve

\n

  testODE <- function(time_space, initial_contamination, parameters){\n  with(\n    as.list(c(initial_contamination, parameters)),{\n      dContamination <- -d*exp(-g*time_space)*Contamination\n      return(list(dContamination))\n    }\n  )\n}\n\nparameters <- c(C = -8/3, d = -10, g =  28)\nY=c(y=1200)\ntimes <- seq(0, 6, by = 0.01)\ninitial_contamination=c(Contamination=1200) \nout <- ode(initial_contamination, times, testODE, parameters, method = "radau",atol = 1e-4, rtol = 1e-4)\n\nplot(out)\n

\n

R 实现 - 从 deSolve 运行编译后的代码

\n

library(deSolve)\nlibrary(scatterplot3d)\ndyn.load("Code/testODE.so")\n\nY <-c(y1=initial_contamination) ;\nout <- ode(Y, times, func = "derivs", parms = parameters,\n           dllname = "testODE", initfunc = "initmod")\n\n\nplot(out)\n

\n

Answer 1

编译的代码不会为R中实现的deSolve模型提供不同的结果，除了和限制内的潜在舍入误差之外。atolrtol

原始帖子中差异的原因是代码中有两个错误。可以按如下方式更正它：

1. 声明static double为 parms[3];而不是parms[4]
2. derivs中的时间t是一个指针，即*t

代码如下：

/* 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];
}

这是两个模拟之间的比较，经过一定程度的调整和概括：

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 !!

注意：根据您的操作系统设置.dll或，或使用 检测。.so.Platform$dynlib.ext