deSolve: Can't understand how to early stop ode solver with root functions-CodePudding

I am confused about how to stop the solver when a certain condition is met. I prepared a dummy SIR model that should stop once the I compartment reaches a certain value. But in my code the solver simply continues on:

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()

The solver should stop if the root evaluates to 2, trigger an event if evaluates to zero and continue in all the other cases. But instead the root being 2 does not stop it.

CodePudding user response：

In the event(root=>action) mechanism, the event is located at a root of a continuous function of the state. In your case, the root functions would be y['I']-10 and y['I']-80, rootfun is the list of these functions (or the function returning the list of their values).

This function gets evaluated frequently on all parts of the solution curve to detect a sign change (it might also work, for some steppers, if a component is piecewise constant and the root function hits exactly zero). The interval of a sign change is then refined with a bracketing method. Apart from providing these values, no processing should and can happen in the root function.

The action on the state is encoded in eventfun, it returns the new state after the event. Internally, the integration is stopped at the root and restarted anew with the returned state as initial value.

Termination is encoded with the terminalroot variable. It is an index and determines which root function provides the termination event.

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

should work with all the other lines unchanged. Note that the trigger component is now unused and could be removed.