0

I'm having an issue with solving a periodic diff eq model.

The trajectory of the model is cyclic and I want to stop it when a defined of cycles have been run.

For this purpose I run the following root function in deSolve::ode:

rootfun <- function (t, y, pars) {
        #message("rootfun")

        change <- detectChange(t, y)

        roots <- c(1, 1)

        if (any(change %in% c('turn_on', 'turn_off', 'recording_off'))) roots[1] <- 0
        if ('exit' %in% change) {
            message('exit')
            roots[2] <- 0
        }

        roots
    }

called in:

ode(y = stateVec, times = seq(tStart, tStop, by = tRes), func = dTempFun, parms = NULL,
                                                     events = list(func = eventfun, root = T, terminalroot = 2), verbose = T,
                                                     rootfunc = rootfun)

with detectChange being a function that given the state of the system when the root is triggered defines what eventfun should do. The root has two components, when the first is zero eventfun is run, while is the second is zero ode should stop running.

In the output, I would expect to see the "exit" message appear only once and then the result of the ode to be returned.

Instead, I have the following output:

eventfun

-------------------- 
Time settings 
--------------------
 
  Normal computation of output values of y(t) at t = TOUT 

-------------------- 
Integration settings 
--------------------
 
  Model function an R-function:  
  Jacobian not specified 

DLSODAR- A switch to the BDF (stiff) method has occurred    
at T (R1), the new step size is (R2): 315.184, 0.0208815                
eventfun
turn_on
root found at time 320.5
DLSODAR- A switch to the BDF (stiff) method has occurred    
at T (R1), the new step size is (R2): 321.846, 0.00643464                
eventfun
turn_off
root found at time 326.768
DLSODAR- A switch to the BDF (stiff) method has occurred    
at T (R1), the new step size is (R2): 327.483, 0.00632845                
eventfun
recording_off
root found at time 327.638
DLSODAR- A switch to the BDF (stiff) method has occurred    
at T (R1), the new step size is (R2): 327.996, 0.00735142                
eventfun
turn_on
root found at time 335.069
exit
exit
DLSODAR-  One or more components of g has a root  
      too near to the initial point.

("eventfun" is a flag that warns me when the event function was run, "turn_on", "turn_off", "recording_off" are the type of event triggered by the first component of the root, while "exit" is printed only if the second component is zero.)

In the end, we can see two "exit" messages followed by an error complaining of two consecutive roots. If I print the root at every cycle it's possible to see that the second root component is zero both times.

Why is this happening? should ode not stop the first time the root is zero at the second element? I do not have the same problem with simpler models like in deSolve: Can't understand how to early stop ode solver with root functions.

An hack I found to make ode stop properly was to force the two "exit" statuses be separated by a non-exit root event.

Bakaburg
  • 3,165
  • 4
  • 32
  • 64
  • Can you please provide a minimal fully reproducible example? I have an idea, why two exits occur -- and they are essentially the same. – tpetzoldt Feb 20 '22 at 19:18
  • unfortunately I cannot reproduce it using a simpler model, it works as expected in that case. But if you tell me your idea I can test it out! – Bakaburg Feb 20 '22 at 21:52
  • 1
    Then give us at least a simplified example with the root and `detectChange` functions as above, even if it does not reproduce the double exit. – tpetzoldt Feb 21 '22 at 06:18

0 Answers0