我喜欢求解涉及多个阈值的耦合微分方程组。通过查看 R 信息，我将 ODE 与根函数和事件函数结合使用。

通过各种示例，即温度模型，第 14 页http://cran.r-project.org/web/packages/diffEq/vignettes/ODEinR.pdf - 下面粘贴的代码 - 我能够让我的模型发挥作用在一个阈值上，即找到一个变量的根/达到阈值会触发一个事件。

library(deSolve)\nyini <- c(temp = 18, heating_on=1)\n\ntemp <- function(t,y, parms) {\n  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)\n  dy2 <- 0\n list(c(dy1, dy2))\n}\n\nrootfunc <- function(t,y,parms) c(y[1]-18, y[1]-20)\n\neventfunc <- function(t,y,parms) {\n  y[1] <- y[1]\n  y[2] <- ! y[2]  \n return(y)\n}\n\ntimes <- seq(from=0, to=20, by=0.1)\nout <- lsode(times=times, y=yini, func = temp, parms = NULL, \n         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))\nplot(out, lwd=2)\nattributes(out)$troot\n

该示例还表明，不同的根可以触发相同的事件函数（y[1] \xe2\x80\x93 …

我试图通过 deSolve 求解这个常微分方程组，dX/dt = -X*a + (YX) b + c 和 dY/dt = -Y a + (XY)*b 时间为 [0,200]，a=0.30 ，b=0.2，但时间 [50,70] 时 c 为 1，否则为 0。我一直使用的代码是，

time <- seq(0, 200, by=1)
parameters <- c(a=0.33, b=0.2, c=1)
state <- c(X = 0, Y = 0)

    two_comp <- function(time, state, parameters){
      with(as.list(c(state, parameters)), {
        dX = -X*a + (Y-X)*b + c
        dY = -Y*a + (X-Y)*b
        return(list(c(dX, dY)))
      })
    }

out <- ode(y = state, times = time, func = two_comp, parms …

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

C代码实现

/* 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 …

我正在尝试复制此处发布的湖泊食物网络模型。该模型代表两条食物链（沿海与远洋），由顶级捕食者（鱼类）连接。我已经对模型进行了编码，但是当我在 2-3 个时间步长后运行它时，模型会生成NaN. 我已经多次检查我的代码，寻找括号等问题，但找不到问题。

如果我将fish初始丰度设置为 0，模型就会运行，所以我认为问题一定出在模型的鱼组件上。

以下是方程式：

Ap = 中上层资源，Z = 中上层浮游动物，Pp = 中上层捕食者，F = 鱼类，Al = 沿岸资源，I = 无脊椎动物，Pl = 沿岸捕食者。

这是我对模型进行编码的尝试：

library(deSolve)

# define the model
vade_2005_model <- function(Time, State, Pars){
  
  with(as.list(c(State, Pars)), {

# pelagic components -----------------------------------------------

# resource
pel_res_dt <- (rPel * AP * (1 - (AP/(KT * q)))) - (aZP * ((Z * AP)/(AP + bZP)))

# zooplankton
pel_zoo_dt <-  (ef * aZP * ((Z * AP)/(AP + bZP))) …