TLDR: My main challenge is just how do I write the root function that checks an arbitrary number of state variables, x, and then apply the event function such that all state variables n that have a value less than the threshold (n <= x) are acted upon by the event function?

I'm trying to use deSolve for a set of Lotka-Volterra equations, but with many state variables (i.e. not just a predator and prey but 20 interacting organisms).

I want to use a root function and event function to be constantly checking if any state variable values dip below a threshold value (e.g. 1.0) and if they do, use the event function to make that particular state variable 0. I've been messing around with a simple minimal example, but can't quite understand how to extend this to check all the state variables and just apply to the one(s) that is/are below the threshold.

The

But now I want to extend this so that it checks all the state variables (in this case just the 2), but ideally in a way that is flexible to different numbers of state variables. I have tried messing around with doing this in some sort of loop but can't seem to figure it out. My main challenge is just how do I write the root function that checks an arbitrary number of state variables, x, and then apply the event function such that all state variables n that have a value less than the threshold (n <= x) are acted upon by the event function?

Perhaps useful information is at some point I would like to implement a separate (not root-based) event function to change a parameter at some pre-set times, so ideally the solution to this problem could interface with additional event function implementation.

Help much appreciated as always!!

CodePudding user response：

One can use a vectorized version of the LV model and then write rootfun and eventfun also in vectorized style:

library(deSolve)

model <- function(t, y, parms) {
  with(parms, {
    dy <- r * y    y * (A %*% y)
    list(dy)
  })
}

## int6eraction matrix
parms <- list(
  r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
  A = matrix(c(
    0.0, 0.0, -0.2, 0.0, # prey 1
    0.0, 0.0, 0.0, -0.1, # prey 2
    0.2, 0.0, 0.0, 0.0,  # predator 1; eats prey 1
    0.0, 0.1, 0.0, 0.0), # predator 2; eats prey 2
    nrow = 4, ncol = 4, byrow = TRUE)
)

times = seq(0, 150, 1)
y0  = c(n1 = 1, n2 = 1, n3 = 2, n4 = 2)

out <- ode(y0, times, model, parms)
plot(out)

## defined as global variables for simplicity, can also be put into parms
threshold <- 0.2  # can be a vector of length(y0)
y_new     <- 1.0  # can be a vector of length(y0)

## uncomment the 'cat' lines to see what's going on
rootfun <- function (t, y, p) {
  #cat("root at t=", t, "\n")
  #cat("y old =", y, "\n")
  return(y - threshold)
}

eventfun <- function(t, y, p) {
  #cat("y old =", y, "\n")
  y <- ifelse(y <= threshold, y_new, y)
  #cat("y new =", y, "\n")
  return(y)
}

out <- ode(y0, times, model, parms, 
           events = list(func = eventfun, root = TRUE), rootfunc=rootfun)
plot(out)