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

[Thomas Petzoldt]

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)