使用 disSolve 建模事件时出现未知状态变量错误答案

【问题标题】：Unkown state variable error when modelling events using deSolve使用 disSolve 建模事件时出现未知状态变量错误
【发布时间】：2020-04-02 15:00:46
【问题描述】：

我试图在广义 Lotka-Volterra 模型中对扰动事件进行建模，其中在时间 t 处，将 1 添加到变量 e。我不断遇到以下错误：

检查事件中的错误（事件、时间、Ynames、dllname）： “事件”中的未知状态变量：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))

【问题讨论】：

(1) 删除行“a
感谢@tpetzoldt 的cmets。我匆忙发布了错误的代码。该帖子现在有一个可行的示例。我同意我也可以出于相同目的使用强制功能，但想从数据框开始。在任何情况下，您能否扩展您的第 3 点和第 4 点？
第 3 点和第 4 点是相关的。如果您使用矩阵，那么命名状态变量在开始时会有点复杂。我建议要么从标量模型开始，要么使用事件函数而不是事件表。在这里可以找到一些示例：tpetzoldt.github.io/deSolve-forcing/deSolve-forcing.html

标签： r

【解决方案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)

【讨论】：