1

当满足特定条件时，我对如何停止求解器感到困惑。我准备了一个虚拟 SIR 模型，一旦 I 隔室达到某个值，它就会停止。但在我的代码中，求解器只是继续：

library(deSolve)
library(dplyr)

pars <- c(beta = .1, gamma = .04)

init <- c(S = 100, I = .01, R = 0, trig = 0)

rootFun <- function(t, y, pars) {
    r <- 1
    if (y['I'] > 10 & y['trig'] == 0) r <- 0
    if (y['I'] > 80) r <- 2
    
    if (r == 2) print('should finish')
    
    return(r)
}

eventFun <- function(t, y, pars) {
    message('First threshold passed!')
    
    y['trig'] <- 1
    
    y
}

derFun <- function(t, y, pars) {
    with(as.list(c(y, pars)), {
        dS = -S * I * beta
        dI = S * I * beta - I * gamma
        dR = I * gamma
        
        list(c(dS, dI, dR, 0))
    }) 
}

ode(y = init, func = derFun, parms = pars, times = 1:100, events = list(func = eventFun, root = TRUE, terminalroot = 2),
    rootfun = rootFun) %>% invisible()

如果根计算为 2，求解器应该停止，如果计算为 0，则触发一个事件，并在所有其他情况下继续。但是根为 2 并不能阻止它。

4

1 回答 1

2

在 event(root=>action) 机制中，事件位于状态的连续函数的根。在您的情况下，根函数是y['I']-10and y['I']-80，rootfun是这些函数的列表（或返回其值列表的函数）。

该函数在解曲线的所有部分上频繁评估以检测符号变化（对于某些步进器，如果组件是分段常数并且根函数恰好为零，它也可能起作用）。然后用包围方法细化符号变化的间隔。除了提供这些值之外，根函数中不应该也不能发生任何处理。

对状态的动作被编码在 中eventfun，它在事件之后返回新的状态。在内部，集成在根处停止，并以返回的状态作为初始值重新启动。

terminalroot终止是用变量编码的。它是一个索引并确定哪个根函数提供终止事件。

所以

rootFun <- function(t, y, pars) {
return(c(y['I']-10, y['I']-80))
}

应该与所有其他行保持不变。请注意，触发器组件现在未使用，可以删除。

于 2021-12-16T10:35:45.150 回答