[R] Difficulty coding time-forced functions in deSolve

[Thomas Petzoldt]

Hi,

does the following help you?

Thomas


require(deSolve)

SEIR <- function(t, x, p) {

   if (t < 7)
     xlag <- x
   else
     xlag <- lagvalue(t - 7)

   V <- ifelse(xlag[3]/sum(xlag) > 0.01, 0.25, 0)

   ## uncomment the following for printing some internal information
   #cat("t=", t, " -- ")
   #cat(xlag, " -- ", V, "\n")

   N <- sum(x)
   with(as.list(c(x,p)),{
     dS<-  b*N - d*S - beta*S*I/N - V*S
     dE<- -d*E+beta*S*I/N - epsilon*E - V*E
     dI<- -d*I + epsilon*E - gamma*I - mu*I - V*I
     dR<- -d*R + gamma*I + V*S + V*E + V*I
     list(c(dS, dE, dI, dR), N = unname(N), V = unname(V))
   })
}

num_years  <- 1
time_limit <- num_years*365.00

b       <- 1/(10.0*365)
d       <- b
beta    <- 0.48
epsilon <- 1/4
gamma   <- 1/4
mu      <- -log(1-0.25)*gamma
parms   <- c(b=b, d=d, beta=beta, epsilon=epsilon, gamma=gamma, mu=mu)

xstart <-c(S=999, E=0, I=1, R=0)
times   <- seq(0.0, time_limit, 1.0)

out <-  dede(xstart, times, SEIR, parms, rtol=1e-8, atol=1e-8)

plot(out)