Hi all, I'm trying to use deSolve to solve a series of differential equations with rk4 mimicking a SEIR model, while including an event/function that is not solely time-dependent.
Explicitly: I want to introduce vaccination 7 days after the proportion of I2/N2 reaches 0.01. Here is the code I am using: require(deSolve) require(sfsmisc) SEIR <- function(t, x, p) { with(as.list(c(x,p)),{ dS<-b*N-d*S-beta*S*I/N-V(I, N, t)*S dE<- -d*E+beta*S*I/N -epsilon*E-V(I, N, t)*E dI<- -d*I+epsilon*E-gamma*I-mu*I-V(I, N, t)*I dR<--d*R+gamma*I+V(I, N, t)*S+V(I, N, t)*E+V(I, N, t)*I dN<-dS+dE+dI+dR list(c(dS, dE, dI, dR, dN)) }) } V <-function(I, N, t) {ifelse(t >=8 & I[t-7]/N[t-7]>0.01, 0.25, 0)} num_years <- 10.0 time_limit <-num_years*365.00 Ni <-1.0E3 b <-1/(10.0*365) d <-b beta <-0.48 epsilon <-1/4 gamma <-1/4 mu <--log(1-0.25)*gamma parms <-c(Ni=Ni, b=b, d=d, beta=beta, epsilon=epsilon, gamma=gamma, mu=mu) xstart <-c(S=999, E=0, I=1, R=0, N=1000) times <- seq(0.0, time_limit, 1.0) tol <- 1e-16 my.atol <- rep(tol,5) my.rtol <- 1e-12 out_rk4 <- as.data.frame(rk4(xstart, times, SEIR, parms)) outfilename <- paste("Basic SEIR.csv") write.csv(out_rk4,file=outfilename,row.names=FALSE, col.names=FALSE) If I remove function V and the associated parts within the differential equations, the model runs just fine. If I define V as V<-function(I, N) {ifelse(I/N >0.01, 0.25, 0) the model functions just fine. Any pointers as to how I can code a function that relies on solutions from previous time steps? thank you in advance, Aimee. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.