ODE, multiple roots and events, R

Question

I like to solve a system of coupled differential equations which involve multiple thresholds. Going through the R information this leads me to using ODE in combination with the root function and the event function.

Going through various examples, i.e. temperature model, page 14 http://cran.r-project.org/web/packages/diffEq/vignettes/ODEinR.pdf -- code pasted below --, I am able to let my model act on one threshold, i.e. finding the root/ reaching the threshold for one variable triggers an event.

library(deSolve)
yini <- c(temp = 18, heating_on=1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
 list(c(dy1, dy2))
}

rootfunc <- function(t,y,parms) c(y[1]-18, y[1]-20)

eventfunc <- function(t,y,parms) {
  y[1] <- y[1]
  y[2] <- ! y[2]  
 return(y)
}

times <- seq(from=0, to=20, by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
attributes(out)$troot

The example shows also that different roots can trigger the same event function (y[1] – 18 and y[1]-20 both trigger eventfunc). My question is, however, is it also possible to have different roots trigger different event functions? Said differently, depending on which root is found, a different eventfunc is triggered? Or alternatively, within the same eventfunct, can it do different actions depending on which root is found.

To keep it simple I first wanted to see if this can work using the same example, for instance by naming the roots and by working with an if statement? At the moment that doesn't work. Does anyone have experience with this? If you look at attributes(out) it seems like ODE does keep a record of which root is encountered $indroot (but that's after evaluation.) Thank you in advance.

# library(deSolve)
yini <- c(temp = 18, heating_on=1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
 list(c(dy1, dy2))
}

rootfunc <- function(t,y,parms) {
  yroot <- vector(len = 2)
  yroot[1] <- y[1]-18 
  yroot[2] <- y[1]-20
 return(yroot)
}

eventfunc <- function(t,y, parms) {
  y[1] <- y[1]
  ifelse(yroot[2]==2, y[2] <- y[2], y[2] <- !y[2])
 return(y)
}

times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
         rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)
attributes(out)$troot

did you ever find how to make this work? I want to trigger something based on the value of t — Herman Toothrot, Nov 22 '16 at 22:27
BTW the ifelse in the deriv funct, will create problems. It's better to move it into the event fucnt or create a smooth changes between the heat on and off state — Bakaburg, Dec 01 '21 at 13:21
An `ifelse` does not always create problems. It depends on the model structure and the strength of the impact. — tpetzoldt, Dec 02 '21 at 06:08
What is meant with `ifelse(yroot[2]==2, y[2] <- y[2], y[2] <- !y[2])`? It looks strange for two reasons: (1) assingment of a variable to itself. (2) logical negation of a numerical variable. Note also that the weblink is broken. — tpetzoldt, Dec 02 '21 at 06:22

Bakaburg · Answer 1 · 2021-12-02T08:23:35.250

1

I approached a similar problem by using a global variable set in the root or directly in the main function (useful if you trigger it based on overcoming a threshold in a particular direction.

The global flag then changes the behaviour of the event function.

Not very elegant but it works.

In your case the code would become:

trigger <- FALSE

rootfunc <- function(t,y,parms) {
  yroot <- vector(len = 2)
  yroot[1] <- y[1]-18 
  yroot[2] <- y[1]-20

  if (yroot[2] == 0) trigger <<- TRUE

  return(yroot)
}

eventfunc <- function(t,y, parms) {
  y[1] <- y[1]
  if (trigger) y[2] <- y[2] else y[2] <- !y[2]
  return(y)
}

edited Dec 02 '21 at 08:23

answered Dec 01 '21 at 13:09

Bakaburg

3,165
4
32
64

Two questions and a comment: (1) I don't think that it works this way, except if the "global flag" within `rootfunc` is assigned with the super assignment operator `<<-`, but I would not consider this as good practice here. (2) Why such a global flag? We have access to the states in both functions anyway, so one can evaluate the roots also in the `event`-function. (3) For special cases, it is possible to create a "technical" state variable, that stores additional information. – tpetzoldt Dec 02 '21 at 06:02
You're right about the <<-, I mistyped it! About the use of the global, it was a fast solution for @Linda. As you may remember I am working on a similar problem and my solution too was to add technical state variables and change it in the event fun. It was just a bit annoying to have to re-evaluate the condition in both the root and the event; my model was simple but it can get cumbersome on complex ones isn't it? – Bakaburg Dec 02 '21 at 08:30
1

You can write a function that evaluates the condition and then call it from both, the root and the event function. Super assignment is often questionable, especially for variables in the global workspace, but I would consider it as legal way if done in a the context of a `function` or `local`. – tpetzoldt Dec 02 '21 at 09:37
1

I eventually ditched the superassignment, but yes it was inside a function call – Bakaburg Dec 02 '21 at 13:20

score 1 · Answer 2 · answered Dec 02 '21 at 06:46

The state y of the system is available in both, the root function and the event function, so it can be used as a condition which event to trigger.

For more complex cases it is of course also possible to dispatch the events from a main event function to different functions for the details, the same is also possible for checking the root condition.

Thanks to @Bakaburg for spotting this interesting unanswered question.

Here a solution that also simplifies some programming construcs:

library(deSolve)
yini <- c(temp = 18, heating_on = 1)

temp <- function(t,y, parms) {
  dy1 <- ifelse(y[2] == 1, 1.0, -0.5)
  dy2 <- 0
  list(c(dy1, dy2))
}

rootfunc <- function(t, y, parms) {
  yroot <- c(y[1] - 18, y[1] - 20)
  return(yroot)
}

eventfunc <- function(t, y, parms) {
  yroot <- c(y[1] - 18, y[1] - 20)
  whichroot <- which(abs(yroot) < 1e-6) # specify tolerance
  y[2] <- if(whichroot == 2) 0 else 1
  return(y)
}

times <- seq(from=0, to=20,by=0.1)
out <- lsode(times=times, y=yini, func = temp, parms = NULL, 
             rootfun = rootfunc, events = list(func=eventfunc, root = TRUE))
plot(out, lwd=2)

root finding and alternating events

One question. From our previous discussion I understood that you cannot have conditional settings in the main fun since it would create non-differentiable changes. In my model the technical state variable stores the dT itself and is changed conditionally in the event. Does the condition in the deriv fun works here because the `if` would be triggered after the event fun is triggered which restarts the derivation, solving the discontinuity? — Bakaburg, Dec 02 '21 at 08:35

ODE, multiple roots and events, R

2 Answers2

Linked