r - 在捕食者猎物系统的生态建模中正确使用 deSolve

Question

我有一个具有指定参数和初始值的捕食者-猎物模型。我在这里以两种方式求解微分方程 1. 使用 for 循环和 2. 使用 deSolve 包。我相信 for 循环是正确的，应该会给出如下图所示的输出。

##################
#For loop attempt
##################
r=0.1; k=100; v=40; s=.1; tbeg=0; tend=1200; nints=1200
c=.06; a=.12; predator0=c(15); prey0=c(50)

dt=(tend-tbeg)/nints
prey=matrix(0,nints+1,length(prey0))
predator=matrix(0, nints+1, length(predator0))
predator[1, ]=predator0
time=numeric(nints+1)
prey[1, ]=prey0
for(i in 1:nints) {
    dprey=r*prey[i, ]-(r*prey[i, ]*prey[i, ]/k)-(s*prey[i, ]*predator[i, ])/(v+prey[i, ])*dt
    prey[i+1, ]=prey[i, ]+dprey
    dpredator=(a*prey[i, ]*predator[i, ])/(v+prey[i, ])-(c*predator[i, ])*dt
    predator[i+1, ]=predator[i, ]+dpredator
    time[i+1]=time[i]+dt}
matplot(time, prey, type="l", lty=1, main="Case 1:  Predator and Prey Populations over Time", ylab="Population", xlab="Time")
points(time, predator, type="l", lty=2)

我正在尝试使用包 deSove 来解决这个 DE 系统并期望得到类似的输出。代码运行，但提供了不同的答案。 (以下“生态学入门”示例)

##################
#deSolve attempt 
##################
library(deSolve)
#the case of a prey with a logistic growth and predator functional response
predprey_FuncResp  <- function(t, y, parms) {
 n0 <- y[[1]]
 p0 <- y[[2]]
 with(as.list(parms), {
    dpdt <- (a*n0*p0)/(v+n0) - c*p0
#    dndt <- r*n0 - (r*n0^2)/k - (s*n0*p0)/(v+n0)
    dndt <- r*n0*(1-n0/k) - (s*n0*p0)/(v+n0)
    return(list(c(dpdt, dndt)))
    })
 }
parms <- c(a=.12, c=.06 , r=.1,s=.1,k=100, v=40)
 Tmax = 1000 # time horizon  
 TimeStep = 1 # integration time step
 Time <- seq(0, Tmax, by = TimeStep)  # the corresponding vector
LV.out <- lsoda(c(n0 = 50, p0 = 15), Time, predprey_FuncResp, parms)
matplot(LV.out[,1],LV.out[,-1], type='l')

我假设我没有正确使用 deSolve 但看不到我的错误。感谢任何花时间看这个的人。

Answer 1

您的问题是，当您从梯度函数返回项时，您需要切换项的顺序——如果您的初始条件是 ，则需要 list(c(dndt, dpdt)) >{n0,p0}.

我做了一些额外的外观调整。

library(deSolve)
predprey_FuncResp  <- function(t, y, parms) {
    with(as.list(c(y,parms)), {
         dpdt <- (a*n0*p0)/(v+n0) - c*p0
         dndt <- r*n0*(1-n0/k) - (s*n0*p0)/(v+n0)
         return(list(c(dndt, dpdt)))
    })
}
parms <- c(a=.12, c=.06 , r=.1,s=.1,k=100, v=40)
Tmax = 1200 # time horizon  
TimeStep = 1 # integration time step
Time <- seq(0, Tmax, by = TimeStep)
LV.out <- ode(c(n0 = 50, p0 = 15), Time, predprey_FuncResp,parms)
par(las=1,bty="l")
matplot(LV.out[,1],LV.out[,-1], type='l', xlab="time", ylab="density")