Message-ID: <51700BD0.1060900@tu-dresden.de>
Date: 2013-04-18T15:05:52Z
From: Thomas Petzoldt
Subject: [R-sig-dyn-mod] metapopulation SIR model
In-Reply-To: <1366295911.3783.YahooMailNeo@web171506.mail.ir2.yahoo.com>
Dear Vincent,
R is a vectorized language. Loops are often not necessary, and named
vectors work differently. Try the example below. In addition, it would
be also possible to write the model completely in matrix formulation.
Thomas
n <- 3
parameters <- c(beta=1000,gamma=365/13)
S <- rep(0.99, times=n)
I <- rep(0.001, times=n)
R <- rep(0.009, times=n)
state <- c(S, I, R)
sirmodel <- function(t, state, parameters) {
S <- state[1:n]
I <- state[(n+1):(2*n)]
R <- state[(2*n+1):(3*n)]
with (as.list(parameters),{
dS <- -beta*S*I
dI <- beta*S*I-gamma*I
dR <- gamma*I
# return the rate of change
list(c(dS, dI, dR))
})
}
times <-seq(from=0,to=60/365,by=1/365/4)
print(system.time(
out <- ode(y = state, times = times, func = sirmodel, parms = parameters)
))
head(out)
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-", mfrow=c(3,3))