[R-sig-dyn-mod] how to pulse input

[Jinsong Zhao]

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."

Thanks again.

Best,
Jinsong


On 2014/4/22 3:49, Thomas Petzoldt wrote:
> Hi Jinsong Zhao,
>
> your approach seems to work (see below), but I cannot decide if it is
> "correct", because this depends on your question. In your case, the
> pulse does not occur "at time = 25", but rather during a time period
> between 25 and 26. You can add Uc as an "auxiliary variable" to the
> return statement of your model (# <---) and run it with a small time
> step to see what happens.
>
> Instead of "if" constructions in the model, pulse inputs can also be
> specified either as forcings (for a time *period*) or as events at a
> given *point of time*, see ?forcings and ?events in the manuals, or:
>
> http://desolve.r-forge.r-project.org/slides/tutorial.pdf#49
>
> Thomas
>
>
>
> library(deSolve)
>
> model <- function(time, state, parameters) {
>    with(as.list(c(state,parameters)), {
>      # rate of change
>      Dc <- Kd * Enzc
>      ### pulse input at time = 25 ###
>      if (time >= 25 & time < 26) Uc <- 10 + Dc else Uc <- DOC + Dc
>      ###
>      print(c(time, Uc))
>      EPc <- Ke * Uc
>      ELc <- Kl * Enzc
>      Re <- EPc * (1 - SUE) / SUE
>      Rm <- Km * Biomc
>      Rg <- (Uc - EPc - Re - Rm) * (1 - SUE)
>      dEnzc <- EPc - ELc
>      dBiomc <- Uc - EPc - Re - Rm - Rg
>      # return the rate of change
>      list(c(dEnzc, dBiomc), Uc = Uc) # <---
>    })
> }
>
> ### Model parameters
> parameters <- c(Kd = 1,
>                  DOC = 0.0,
>                  Ke = 0.05,
>                  Kl = 0.05,
>                  Km = 0.022,
>                  SUE = 0.5)
>
> ### State variables
> state <- c(Enzc = 0.7,
>             Biomc = 29)
>
> times <- seq(0, 30, by = .01)
> out <- ode(y = state, times = times, func = model, parms = parameters)
>
> plot(out)
>
>
>