我想用deSolve创建一个动态的蝶类生态模型。模拟运行了几年,一些事件是在一年中的某一天触发的(因此我添加了一个状态变量days )。为了触发这些事件,我想使用ifelse语句,它工作得很好,直到我尝试在ifelse语句中加入一个涉及另一个状态变量的操作:D.egg.sus=(ifelse(days<270,(400 * adult.sus),0))。当我这样做时,模拟运行,但它似乎忽略了ifelse语句。有谁能帮帮我吗?下面是我的完整代码:
days = 1
egg.sus = 0
larvae.sus = 0
pupae.sus = 0
adult.sus = 1000
state = c(days = days, egg.sus=egg.sus, larvae.sus=larvae.sus,
pupae.sus=pupae.sus, adult.sus=adult.sus)
model = function(t, state, parameters)
{
with(as.list(c(state, parameters)),
{
D.Days = 1
D.egg.sus =
( ifelse(days<270, (400*adult.sus) ,0)) ## This is the line causing trouble
(- egg.sus/5)
(- egg.sus * rbeta(1, 6.038892/5,1.4612593)*.95)
D.larvae.sus =
(+ egg.sus/5)
(- larvae.sus * rbeta(1, 0.248531/14,0.2094379)*0.95)
(- larvae.sus/14)
D.pupae.sus =
(+ larvae.sus/14)
(- pupae.sus * rbeta(1, 0.022011/15, 1.43503))
(- pupae.sus/15)
D.adult.sus =
(+ pupae.sus/15)
(- adult.sus/30)
list(c( D.Days, D.egg.sus, D.larvae.sus,D.pupae.sus, D.adult.sus))
}
)}
events <- data.frame(var = c('days'),
time = seq(364,73000,by=365) ,
value = 0,
method = "rep")
require(deSolve)
times = seq(1,900, by = 1)
out = ode(y=state, times = times, func = model, parms = parameters, events = list(data=events))
dev.cur()
plot(out, col = 2)发布于 2021-01-17 02:48:35
问题中的模型有几个问题:
days,因为函数中的模拟时间给定为t。然后使用模运算符%%,你就不再需要事件了。parms=NULL。-运算符放在行尾。rgamma )是一个非常糟糕的想法,特别是对于具有自动时间步长的求解器。根据定义,ODE是确定性的。可以考虑使用固定时间步长的求解器,例如,使用非常小的时间步长的method="euler",或者(更好的)提供随机值作为外部输入。如果使用外部输入,则无论如何都可以避免使用ifelse。https://stackoverflow.com/questions/33533400
复制相似问题