deSolve: Не могу понять, как досрочно остановить решатель од с корневыми функциями

Question

Я не понимаю, как остановить решатель при выполнении определенного условия. Я подготовил фиктивную модель 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, инициировать событие, если значение равно нулю, и продолжить выполнение во всех остальных случаях. Но вместо этого корень, равный 2, этого не останавливает.

Answer 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))
}

должен работать со всеми остальными строками без изменений. Обратите внимание, что компонент триггера теперь не используется и может быть удален.