4

我喜欢解决一个涉及多个阈值的耦合微分方程系统。浏览 R 信息，这使我将 ODE 与根函数和事件函数结合使用。

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

library(deSolve)
yini <- c(temp = 18, heating_on=1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
 list(c(dy1, dy2))
}

rootfunc <- function(t,y,parms) c(y[1]-18, y[1]-20)

eventfunc <- function(t,y,parms) {
  y[1] <- y[1]
  y[2] <- ! y[2]  
 return(y)
}

times <- seq(from=0, to=20, by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
attributes(out)$troot

该示例还表明，不同的根可以触发相同的事件函数（y[1] – 18 和 y[1]-20 都触发 eventfunc）。但是，我的问题是，是否也可以让不同的根触发不同的事件函数？换句话说，根据找到哪个root，触发不同的eventfunc？或者，在同一个 eventfunct 中，它是否可以根据找到的根执行不同的操作。

为了简单起见，我首先想看看这是否可以使用相同的示例，例如通过命名根和使用 if 语句？目前这行不通。这个事情谁有经验？如果您查看 attributes(out)，ODE 似乎确实记录了遇到 $indroot 的根（但这是在评估之后。）提前谢谢您。

# library(deSolve)
yini <- c(temp = 18, heating_on=1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
 list(c(dy1, dy2))
}

rootfunc <- function(t,y,parms) {
  yroot <- vector(len = 2)
  yroot[1] <- y[1]-18 
  yroot[2] <- y[1]-20
 return(yroot)
}

eventfunc <- function(t,y, parms) {
  y[1] <- y[1]
  ifelse(yroot[2]==2, y[2] <- y[2], y[2] <- !y[2])
 return(y)
}

times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
attributes(out)$troot
4

2 回答 2

1

我通过在根中或直接在主函数中使用全局变量集来解决类似的问题（如果您基于克服特定方向的阈值来触发它，这很有用。

然后全局标志会改变事件函数的行为。

不是很优雅，但它有效。

在您的情况下，代码将变为：

trigger <- FALSE

rootfunc <- function(t,y,parms) {
  yroot <- vector(len = 2)
  yroot[1] <- y[1]-18 
  yroot[2] <- y[1]-20

  if (yroot[2] == 0) trigger <<- TRUE

  return(yroot)
}

eventfunc <- function(t,y, parms) {
  y[1] <- y[1]
  if (trigger) y[2] <- y[2] else y[2] <- !y[2]
  return(y)
}
于 2021-12-01T13:09:53.573 回答
1

系统的状态y在根函数和事件函数中都是可用的，因此它可以作为触发哪个事件的条件。

对于更复杂的情况，当然也可以将事件从主事件函数分派到不同的函数以获取细节，同样可以检查根条件。

感谢@Bakaburg 发现这个有趣的未回答问题。

这是一个也简化了一些编程结构的解决方案：

library(deSolve)
yini <- c(temp = 18, heating_on = 1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
  list(c(dy1, dy2))
}

rootfunc <- function(t, y, parms) {
  yroot <- c(y[1] - 18, y[1] - 20)
  return(yroot)
}

eventfunc <- function(t, y, parms) {
  yroot <- c(y[1] - 18, y[1] - 20)
  whichroot <- which(abs(yroot) < 1e-6) # specify tolerance
  y[2] <- if(whichroot == 2) 0 else 1
  return(y)
}

times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
             rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)

根发现和交替事件

于 2021-12-02T06:46:52.953 回答