Skip to content

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

3 messages · Jinsong Zhao, Thomas Petzoldt

Original
Jinsong Zhao
# 
Hi there,

I want to affect a model with pulse input. Now I do this something like 
this:

C_Only <- 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))
     })
}

Is it the correct way to set pulse input? How I can set the pulse input 
time within function argument?

Any help will be really appreciated.

Best wishes,
Jinsong
1 day later
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)
3 days later
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: