Неизвестная ошибка переменной состояния при моделировании событий с использованием deSolve

Question

Я пытаюсь смоделировать тревожное событие в обобщенной модели Лотки-Вольтерры, где в момент времени t 1 добавляется к переменной e . Я продолжаю сталкиваться со следующей ошибкой:

Ошибка в контрольных событиях (события, времена, Ynames, dllname): неизвестная переменная состояния в 'event': e

Моя модель следующая:

lvg<-function(t, N, e, param){
    e <- 0
    dNdt <- N * r + N * (a %*% N) - N * e
    list(c(dNdt))    
}

, где N - численность популяции вида i , r - скорость роста, a - матрица взаимодействия, и e - это событие. r и a указываются в качестве предыдущих параметров, событие указывается в кадре данных. Упрощенная версия выглядит следующим образом:

#set parameters
S<- 10 #  number of species
r <- rep(1.1, S) # growth rates
a <- matrix (nrow = S, ncol = S) #interaction matrix
a[lower.tri(a)] <- -0.001
a[upper.tri(a)] <- -0.001
diag(a) <- -0.01

parms <- list (r, a) #put parameters in a list
N0 <- rep(100, S) #initial values for species abundances
ts<-seq(0, 100, 1) # time steps for solver

#create data frame for event
eventdat <- data.frame(var = c("e", "e"), time = c(10, 20), value = c(1, 1), method = c("add"))

lvout<-lsoda(N0, ts, lvg, parms, events = list(data = eventdat)) 

Answer 1

Здесь подход с функцией события вместо таблицы событий, которая в общем случае более гибкая, в данном случае проще. Также обратите внимание, что число состояний и значений параметров были изменены, чтобы получить более типичную модель L & V:

library(deSolve)

## multi-species Lotka-Volterra
lvg <- function(t, N, param) {
  with(param, {
    dNdt <- r * N + N * (a %*% N)
    list(c(dNdt))
  })
}

## simplified to 4 species, you can add more
S <- 4
N0 <- c(1,1,1,1)

## parameter list
parms <- list(
  r = c(r1 = 0.5, r2 = 0.5, r3 = -0.5, r4 = -0.5),
  a = matrix(c(
    0.0, 0.0, -0.5, 0.0, # prey 1
    0.0, 0.0, 0.0, -0.2, # prey 2
    0.5, 0.0, 0.0, 0.0,  # predator 1; eats prey 1
    0.0, 0.2, 0.0, 0.0), # predator 2; eats prey 2
    nrow = 4, ncol = 4, byrow = TRUE),
  e = rep(0.5, S)
)

ts <- seq(0, 100, 1) # time steps for solver
te <- c(20, 40)     # event times

## event function is more flexible than an event table
eventfun <- function(t, N, param){
  with (as.list(param), {
    N <- N - N * e
    return(c(N))
  })
}

## simulation without events
lvout<-lsoda(N0, ts, lvg, parms)
plot(lvout)

## simulation with events
lvout<-lsoda(N0, ts, lvg, parms, events = list(func = eventfun, time = te))
plot(lvout)