1

I want to create a Dynamic model of butterfly ecology using deSolve. the simulation runs over several simulation years and some events are triggered by the day of the year (so I added one state variable of days ). in order to trigger those events I want to use an ifelse statement and it works fine, until I try to put in the ifelse statement an operation involving another state variable: D.egg.sus=(ifelse(days<270,(400 * adult.sus),0)). When I do so, the simulation runs, but it seems to ignore the ifelse statement. can anyone help me please? here is my full code:

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)

2 Answers2

1

I don't know about five years ago, but at the time of writing ifelse works just fine with deSolve. Your issue seems to be that the returned value of your condition did not return as you wanted. Instead you might want to use a flag variable or save the return from your ifelse to a variable that you can then use in your model.

Here is a small example demonstrating how you can use a flag in your model parameters

library(deSolve)

# Our model function, first-order
# One parameter is a flag that is used by the ifelse to set Ka to zero if TRUE.

onecomp <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    Ka = ifelse(flag == TRUE, 0, Ka) # Use ifelse to check for negative values
    
    dX <- - X*Ka
    dY <- X*Ka - Y*Ke
    list(c(dX, dY))
  })
}

times <- seq(0, 24, by = 0.01)
parameters <- c(Ka = 0.8 , Ke = 0.2, flag = FALSE)
state <- c(X = 100 , Y = 0)

# Test for TRUE
out <- ode(y = state, times = times, func = onecomp, parms = parameters)
plot(out)

# Test for FALSE, where we expect no transfer.
parameters <- c(Ka = 0.8 , Ke = 0.2, flag = TRUE)
out <- ode(y = state, times = times, func = onecomp, parms = parameters)
plot(out)

Created on 2021-01-13 by the reprex package (v0.3.0)

mhovd
  • 3,724
  • 2
  • 21
  • 47
1

The model in the question has several issues:

  1. You can use the simulation time directly instead of a state variable days, because simulation time in the function is given as t. Then just use the modulo operator %% and you don't need events anymore.
  2. the parameters are all hard-coded, so use parms=NULLin the ode function.
  3. line breaks are wrong. R continues lines if (and only if) they are not yet syntactically complete. Therefore, remove obsolete parentheses and, for example, put the - operator at the end of the line.
  4. Use of a random number e.g. rgamma within an ODE function is a very bad idea, especially for solvers with automatic time steps. ODEs are deterministic by definition. One may consider a fixed time-step solver instead, e.g. method="euler"with a very small time step or (much better) to provide the random values as an external input (forcing).
  5. If you use an external input, you can avoid the ifelse anyway.
tpetzoldt
  • 5,338
  • 2
  • 12
  • 29
  • ... and looking once again, I am not sure if this is really an ODE. Can it be possible that you have a difference equation in mind? – tpetzoldt Jan 16 '21 at 19:11