[R] Difficulty coding time-forced functions in deSolve
Aimee Kopolow
alj27 at georgetown.edu
Wed Apr 2 05:51:34 CEST 2014
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.
More information about the R-help
mailing list