0

我在 R 的元种群中编写 SIR 模型，我想集成 systema，为此我使用 de deSolve 和 C 编译代码，我以前使用过这个，但在有几个参数的情况下，现在我会有 Nxm参数，其中 N 是系统的维度，所以我想要那个

/* file age3classp.c */
#include <R.h>
static double parms[3];
static double forc[1];

#define N parms[0]
#define N1 parms[1]
#define gam3 parms[2]

这个参数是向量或矩阵 NxN 可能吗？

在 C 中，我的模型将采用以下形式：

# SIR metapopulation model:
SIR <- function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    dS = c()
    dI = c()
    dR = c()
    for(i in c(1:dim)){
      dS[i] <- delta_N[i]*(S[i]+I[i]+R[i])
      dI[i] <- 10
      dR[i] <- 10
    }
    list(c(dS, dI, dR))
  }) 
}

population <- c(S <- matrix(100,ncol=N,nrow =1 ), I <- matrix(10,ncol=N,nrow =1 ),
                R <- matrix(0,ncol=N,nrow =1 ))
z <- ode(population, times, SIR, parameters)

通过这种方式，它不会将 S[i] 或其他变量识别为变量，而只是作为初始条件值。

我该怎么做才能将其识别为变量？

4

1 回答 1

0

是的，这是可能的，并且有不同的方法可以做到这一点，具体取决于您的 C 编程技能。最简单的方法是将状态和参数放在两个长向量中，然后使用编号参数和变量索引在 C 级别将其拆分。然后将方程表示为for-环。

为了提高可读性，还可以使用

  1. 索引的预处理器常量或
  2. 联合和结构（见下文）
  3. 状态向量和参数列表

在R和C级别上，状态 ( y) 始终被视为向量，但参数 ( ) 也可以作为列表向下传递，然后在C级别进行拆分。这可能很棘手，需要对R的数据结构有所了解。p

但是，我建议在R级别开始矢量化。R对矢量化模型非常快，因此加速可能无法补偿C编程工作。可以在此处找到如何实现矢量化捕食者-猎物模型的示例。

另一个想法是使用代码生成器，因此您可以查看 CRAN 包rodeo，它从公式（即 LibreOffice 或 Excel）表中创建快速 Fortran 代码。使用不需要 Fortran 知识。

可以在一篇论文 ( https://doi.org/10.1016/j.envsoft.2017.06.036 ) 和https://dkneis.github.io/上的包文档中找到有关牛仔竞技表演的更多信息

如果有人真的想在这里用C 语言编写 Lotka-Volterra-Competition 模型的小型实现，请参阅 Wikipedia with 3 states。p参数在 C 级别作为参数向量传递，而 aunion用于提高可读性：

/* file model.c */

#include <R.h>

union parvec {
  struct {
    double r[3], a[6];
  };
  double value[9];
} p;

/* initializer  */
void initmod(void (* odeparms)(int *, double *))
{
    int N = 9; /* total number of parameters */
    odeparms(&N, p.value);
}

/* Derivatives */
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip) {

  double y_sum = 0;
  for (int i = 0; i < *neq; i++) {
    y_sum = 0;
    for (int j = 0; j < *neq; j++) y_sum += p.a[i + *neq * j] * y[j];
    ydot[i] = p.r[i] * y[i] * (1 - y_sum);
  }
}

这里调用R代码：

# file call_model.R
library(deSolve)

system("R CMD SHLIB model.c")
dyn.load("model.dll")

p <- c(r = c(0.1, 0.3, 0.04), A = c(0.2, 0.3, 0.3, 0.5, 0.4, 0.2))
y <- c(X = c(2, 2, 2))
times <- seq(0, 200, by = 0.1)

out <- ode(y, times, func = "derivs", parms = p,
           dllname = "model", initfunc = "initmod")
matplot.0D(out)

dyn.unload("model.dll")

当然，更详细的解决方案是可能的。

C中具有3个方程的竞争模型

于 2021-11-29T18:26:18.407 回答