R:将参数移交给函数(deSolve-package)

R: handing over parameters to function (deSolve-package)

我最近下载了 deSolve 包来求解 ODE 模型。我根据我找到的一些示例代码编写了一个模型,但是,由于我收到有关未知参数的错误信息,因此似乎在传递参数方面存在问题。

  library(deSolve)

  model <- function(t, y, parms) {
  dY1 = -(k1 * y[1]) + (k2 * y[6]) - (k13 * y[1] * 400*sin(((2*pi)/period_ca)*t-phase_ca)+600) + (k14 * y[2]) - (k17 * y[1] * 400*sin(((2*pi)/period_dg)*t-phase_dg)+600) + (k18 *y[11]) - (k3*y[1] * AA) + (k4 * y[7])
  dY2 = (k13 * y[1] * 400*sin(((2*pi)/period_ca)*t-phase_ca)+600) - (k14 * y[2]) - (k15 * y[2] *400*sin(((2*pi)/period_dg)*t-phase_dg)+600) + (k16 * y[3]) - (k5 * y[2]) + (k6 * y[8]) - (k7 * y[2] * AA) + (k8 * y[9])
  dY3 = (k15 * y[2] * 400*sin(((2*pi)/period_dg)*t-phase_dg)+600) - (k16 * y[3]) - (k9 * y[3]) + (k10 * y[10])
  dY4 = -(k11 * y[4]) + (k12 * y[5]) + (k19 * y[11] * AA) - (k20 * y[4])
  dY5 =  (k11 * y[4]) - (k12 * y[5])
  dY6 = (k1 * y[1]) - (k2 * y[6])
  dY7 = (k3 + y[1] * AA) - (k4 * y[7])
  dY8 = (k5 * y[2]) - (k6 * y[8])
  dY9 = (k7 * y[2] * AA) - (k8 * y[9])
  dY10 = (k9 * y[3]) - (k10 * y[10])
  dY11 = (k17 * y[1] * 400*sin(((2*pi)/period_dg)*t-phase_dg)+600) - (k18 * y[11]) - (k19 * y[11] * AA) + (k20 * y[4])
  list(c(dY1, dY2, dY3,dY4, dY5, dY6,dY7, dY8, dY9, dY10, dY11))
  }

  yini <- c(y1 = 1000, y2 = 0, y3 = 0, y4 = 0, y5 = 0, y6 = 20, y7 = 0, y8 =    0, y9 = 0, y10 = 0, y11 = 0)
  times <- seq(from = 0, to = 5000, by = 0.1)
  parms <- c(AA=11000, k1=1, k2=50, k3=1.2e-7, k4=0.1, k5=1.2705, k6=3.5026, k7=1.2e-6, k8=0.1, k9=1, k10=0.1, k11=2, k12=0.2, k13=0.0006, k14=0.5, k15=7.998e-6,
           k16=8.6348, k17=6e-7, k18=0.1, k19=1.8e-5, k20=2, period_ca=100, phase_ca=0, period_dg=100,
           phase_dg=0)
  out <- ode (times = times, y = yini, func = model, parms = parms)

这里 dY1 到 dY11 代表某些系统组件的微分方程。 parms 是定义必要参数值的向量,yini 定义初始条件和乘以时间尺度。

我收到以下错误消息:

Error in func(time, state, parms, ...) : object 'k1' not found

我是 R 的新手,不明白问题的根源(我找到的所有示例代码都是以相同的方式构造的)。

您必须从 deSolve 包中查找一些示例。要在函数中使用您的参数,您需要使用 with 函数:

model <- function(t, y, parms) {
  with(as.list(c(y, parms)), {
  dY1 = ...
  .
  .
  .
  dY11 = ...

  list(c(dY1, dY2, dY3,dY4, dY5, dY6,dY7, dY8, dY9, dY10, dY11))
})
}

然而我得到了一个错误,因为集成失败,但这不是因为代码错误。也许您的 PC 上不是这种情况。

您可以尝试使用此代码求解 ODE,它会给出一些警告,但集成成功(在我的 PC 上):

out <- ode(times = times, y = yini, func = model, parms = parms, method = "bdf")