[R-sig-dyn-mod] question about the limit of conversion rate

[Thomas Petzoldt]

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