deSolve：无法理解如何使用根函数提前停止 ode 求解器

Posted 2023-02-18

【中文标题】deSolve：无法理解如何使用根函数提前停止 ode 求解器

【英文标题】：deSolve: Can't understand how to early stop ode solver with root functions

【发布时间】：2022-01-19 09:55:24

【问题描述】：

我对如何在满足特定条件时停止求解器感到困惑。我准备了一个虚拟 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 并没有阻止它。

【参考方案1】：

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

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

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

终止使用terminalroot 变量编码。它是一个索引，决定了哪个根函数提供了终止事件。

所以

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

应该在所有其他行不变的情况下工作。请注意，触发器组件现在未使用，可以删除。

【讨论】：

谢谢，我完全误解了它应该如何工作！