在 deSolve 中的时间步更改参数值

Question

我正在尝试使用 deSolve 求解 R 中的 ODE。使用以下代码，我希望参数 gamma0 在时间步长 0、1、2、3、4、5、6、7、8、9 和 10 取值 5，否则取 0。但是，print(gamma0) 显示 gamma0 保持在 0。

这是我的 ODE：

library(deSolve) 
param <- c(a = 0.1, b = 1) 
yini <- c(alpha0 = 6, beta0 = 2) 

mod <- function(times, yini, param) { 

  with(as.list(c(yini, param)), { 

    gamma0 <- ifelse(times %in% seq(0,10,1), 5, 0) 

    ## print(gamma0) 

    dalpha0 <- - a*alpha0 + gamma0 
    dbeta0 <- a*alpha0 - b*beta0 
    return(list(c(dalpha0, dbeta0))) 

  })} 

times <- seq(from = 0, to = 10, by = 1/24) 
out <- ode(func = mod, times = times, y = yini, parms = param) 
plot(out, lwd = 2, xlab = "day")

我做错了什么？

Answer 1

我得到的结果与您略有不同。如果我取消注释 print(gamma0) 它打印出 5 两次然后打印出 513 个零。不难从表面上追溯原因，尽管您想要的可能比我在这里提供的要多。

在您有（注释掉的）语句 print(gamma0) 的地方，输入以下行：

cat("g:", gamma0, "  t:", times, "\n")

和运行代码。您会看到它显示的前两次是 0。因为那些在您的列表中 seq(0,10,1) gamma0 是 5。之后，显示的时间值会发生变化。请注意，打印的 none 来自您的原始时间列表 seq(from = 0, to = 10, by = 1/24) 并且其中 none 是整数，因此 none满足将 gamma0 设置为 5 的条件。ode 正在对时间做一些事情（插值？），但它不仅仅是使用您提供的值。事实上，它并没有打印出 241 个 gamma0 和 times 的值。它打印出 515 个这样的值。我注意到结果 out 确实有 241 个值。

根据你的问题，我认为你假设 ode 实际上会在你的 times 处计算函数。它不是。它将 times 视为连续变量。但是你的情况

gamma0 <- ifelse(times %in% seq(0,10,1), 5, 0)

仅测试 11 个特定值 - 不是值范围。连续变量不太可能恰好达到这些值。

Answer 2

这是对您的函数的一个非常简单的修改。如果您有兴趣知道自己做错了什么，可以看下面。

mod <- function(times, yini, param) { 

  dt = times[2] - times[1]
  with(as.list(c(yini, param)), { 

    gamma0 <- ifelse(times <= 10*dt, 5, 0) 

    ## print(gamma0) 

    dalpha0 <- - a*alpha0 + gamma0 
    dbeta0 <- a*alpha0 - b*beta0 
    return(list(c(dalpha0, dbeta0))) 

  })}

编辑

与 G5W 的回答相同，问题在于您将 times 与什么进行比较。

写作时

times %in% seq(0,10,1)

你指的不是时间步长。您只需参考 times.

的值

因此，如果您想在前 10 个时间步使用它，您只需要使用我的代码或任何考虑 dt.

的代码即可

但是这里有一个问题要问你：

如果您不需要 gamma0 根据 times 进行更改并希望它在前 11 (10) 个时间步长为 5，那么为什么要将它与 times？为什么不简单地将这些时间步长设置为 5？

Answer 3

M-- 的答案对我不起作用，如果你试试这个呢？

mod <- function(times, yini, param) {
  with(as.list(c(yini, param)), {

    if (times < 10) {
      gamma0 = 5
    } else if (times >= 10) {
      gamma0 = 0
    }

    dalpha0 <- -a * alpha0 + gamma0
    dbeta0 <- a * alpha0 - b * beta0
    return(list(c(dalpha0, dbeta0)))

  })
}

这是一个 Headvise 类型的函数

在 deSolve 中的时间步更改参数值

Change parameter values at time step in deSolve

parameters

r

ode

differential-equations