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

Question

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

3

我有一个捕食-猎物模型，参数和初始值已指定。在这里，我使用两种方法解决微分方程：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)

enter image description here

我试图使用deSove软件包来解决这个ODE系统，并期望得到类似的输出。代码可以运行，但提供了不同的答案。（参考《生态学入门》中的例子）

##################
#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时出了错误，但是我无法找到我的错误。感谢任何花时间查看此问题的人。

- Wes McClintick

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

你的问题是从梯度函数返回项的顺序错了 - 如果你的初始条件是{n0，p0}，你需要 list(c(dndt，dpdt))。

我进行了一些额外的化妆调整。

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")