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

Question

我有一个 ODE，我想使用从 R 的 deSolve 包调用的编译 C 代码来求解。有问题的 ODE 是一个指数衰减模型 (y'=-d* exp(g* time)*y)：但是运行 R 中的编译代码给出了与 R 的原生 deSolve 不同的结果。就好像它们被翻转了 180º。怎么回事？

C代码实现

/* 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 实现 - 本机 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 实现 - 运行来自 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)

Answer 1

编译代码不会给 与 deSolve 模型在 R 中实现的不同结果，除了atol 和 rtol.

范围内的潜在舍入误差

与原文不同的原因post其中代码中有两个错误。可以更正如下：

将 static double 声明为 parms[3]; 而不是 parms[4]
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 !!

注意：根据你的OS设置.dll或.so，或者用.Platform$dynlib.ext检测。

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

Compiled C ODE gives different results to R's using deSolve

c

r

desolve

C代码实现

R 实现 - 本机 deSolve

R 实现 - 运行来自 deSolve 的编译代码

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

Compiled C ODE gives different results to R's using deSolve

c

r

desolve

C代码实现

R 实现 - 本机 deSolve

R 实现 - 运行 来自 deSolve 的编译代码

R 实现 - 运行来自 deSolve 的编译代码