if value is in vector, perform this function
On 03-03-2013, at 00:57, Louise Stevenson <louise.stevenson at lifesci.ucsb.edu> wrote:
Hi,
I'm trying to set up R to run a simulation of two populations in which every 3.5 days, the initial value of one of the populations is reset to 1.5. I'm simulation an experiment we did in which we fed Daphnia populations twice a week with algae, so I want the initial value of the algal population to reset to 1.5 twice a week to simulate that feeding. I've use for loops and if/else loops before but I can't figure out how to syntax "if t is in this vector of possible t values, do this command, else, do this command" if that makes sense. Here's what I have (and it doesn't work):
params = c(1, 0.15, 0.164, 1)
init = c(1.5, 0.05)
t=seq(1,60, by=0.5) #all time values, experiment ran for 60 days
#feeding sequence - every "3.5 days"
feed_days = seq(1,60,by=3.5)
Daphnia <- function(t,x,params){
C_D = x[2];
C_A = 0;
for(t %in% feed_days){
if t == TRUE {
C_A = 1.5
}
else{
C_A = 0
}}
gamma = params[1]; m_D = params[2]; K_q = params[3]; q_max = params[4];
M_D = m_D * C_D
I_A = (C_D * q_max * C_A) / (K_q + C_A)
r_D = gamma * I_A
return(
list(c(
- I_A,
r_D - M_D
)))
}
library(deSolve)
results <- ode(init, t, Daphnia, params, method = "lsoda")
You have been given a correction for expression for (t %in% feed_days). But even with that correction things will not do as you seem to want. The argument "t" of function Daphnia is the integration time the ode solver is passing and almost certainly is NOT an element of the vector t defined at the start of your script. That "t" is the "the time sequence for which output is wanted" (see ode help); it is what is put into the output of ode. There is no reason to assume that the Daphnia argument t is an element of feed_days. You can easily check this by inserting a print(t) in Daphnia. So C_A will be 0 most of the time. It would certainly help if you named the elements of the init vector and the return list of Daphnia. In Daphnia x[2] is C_D. But what is x[1] (C_A?)? I think you will have to look at deSolve events but I'm not sure if that is possible or required/desired with your model. Berend