On 4/25/2014 8:55 PM, Jinsong Zhao wrote:
Hi Thomas, Thank you very much for the reply. "auxiliary variable" help me much, so I can monitor all variable during the simulation. In my situation, events doesn't work, because Uc is not a state variable, while I consulted the help page of events that stated that "An 'event' occurs when the value of a state variable is suddenly changed, e.g. because a value is added, subtracted, or multiplied."
... ok, then use a forcing function. The following example is
a modified version taken from the ?rk help page.
SPCmod <- function(t, x, parms, input) {
with(as.list(c(parms, x)), {
import <- input(t)
dS <- import - b*S*P + g*C # substrate
dP <- c*S*P - d*C*P # producer
dC <- e*P*C - f*C # consumer
res <-
list(c(dS, dP, dC), import=import)
})
}
parms <- c(b = 0.1, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.1)
times <- seq(0, 100, length = 501)
xstart <- c(S = 1, P = 1, C = 1)
## forcing function: external signal with rectangle impulse
signal <- data.frame(times = times, import = 0)
signal$import[signal$times >= 10 & signal$times <= 11] <- 0.2
sigimp <- approxfun(signal$times, signal$import, rule = 2)
out <- ode(xstart, times, SPCmod, parms, hini = 0.1, input = sigimp)
plot(out)
Note that 'in theory' discrete events violate the assumptions of
differential equations that are continuous systems by definition.
Fortunately, most solvers of deSolve are quite robust and can handle
this in practice, but pulses with smooth transitions can be faster.
This was discussed last year, see the following thread in the archives:
https://stat.ethz.ch/pipermail/r-sig-dynamic-models/2013q2/000214.html
Hope it helps
Thomas