I'm trying to add an if then statement to condition the initial value of one of my state variables, and am using deSolve. Essentially, I want to introduce the 3rd ODE (in this case a 3rd species into a population) after the start of the simulation.
Here is what the code looks like without the condition:
Antia_3sp_Model <- function(t,y,p1){
# Parms
ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5]
# State vars
Pi <- y[1]; Pj <- y[2]; I <- y[3]
# ODEs
dPi = ri*Pi - k*Pi*I
dPj = rj*Pj - k*Pj*I
dI = p*I*(Pi/(Pi + o) + Pj/(Pj + o))
list(c(dPi,dPj,dI))
}
# Parm vals
ri <- 0.3; rj <- 0.2; k <- 0.001; p <- 1; o <- 1000 # Note that r can range btw 0.1 and 10 in this model
parms <- c(ri,rj,k,p,o)
# Inits
Pi0 <- 1; Pj0 <- 1; I0 <- 1
N0 <- c(Pi0,Pj0,I0)
# Time pt. sol'ns
TT <- seq(0.1,200,0.1)
# Sim
results <- lsoda(N0,TT,Antia_3sp_Model,parms,verbose = TRUE)
Here's what I have so far, after trying to add in an if then statement, saying that before time = 50, the initial value of the 3rd state variable will be 0, and that at or above time = 50, the initial value of the 3rd state variable will be 1.
Antia_3sp_Model <- function(t,y,p1){
# Parms
ri <- p1[1]; rj <- p1[2]; k <- p1[3]; p <- p1[4]; o <- p1[5]
# State vars
Pi <- y[1]; Pj <- y[2]; I <- y[3]
if (t[i] < t[50]){
Pj0 = 0
}
else if (t[i] >= t[50]){
Pj0 = 1
}
# ODEs
dPi = ri*Pi - k*Pi*I
dPj = rj*Pj - k*Pj*I
dI = p*I*(Pi/(Pi + o) + Pj/(Pj + o))
list(c(dPi,dPj,dI))
}
# Parm vals
ri <- 0.3; rj <- 0.2; k <- 0.001; p <- 1; o <- 1000 # Note that r can range btw 0.1 and 10 in this model
parms <- c(ri,rj,k,p,o)
# Inits
Pi0 <- 1; Pj0 <- 1; I0 <- 1
N0 <- c(Pi0,Pj0,I0)
# Time pt. sol'ns
TT <- seq(0.1,200,0.1)
# Sim
results <- lsoda(N0,TT,Antia_3sp_Model,parms,verbose = TRUE)
Any suggestions?
Please let me know if I should add any additional information, and thank you so much for reading! :)
Aucun commentaire:
Enregistrer un commentaire