[R-sig-dyn-mod] Modeling dampened oscillator with random resets

Hi Felix,

With simecol, which builds on deSolve you could do the following (based on an example Thomas Petzoldt send me a couple of years ago):

library(simecol)

osci <- odeModel( main=function(time, init, parms,...){
  x <- init
  p <- parms
  ## same as above, but now we have and input, or external stimulus
  S <- approxTime1(inputs, time, rule=2)["s.in"]
  dx <- x[2]
  dy <- parms["eta"] * x[1] + parms["zeta"] * x[2] + S
  list(c(dx, dy))
}, 
                  times=c(from=0, to=200, by=0.2),
                  parms=c(eta=-0.05, zeta=-0.02, x=1, y=0),
                  initfunc = function(obj) {
                    init(obj)[c("x", "y")] <- parms(obj)[c("x", "y")]
                    return(obj)
                  }
)


## Let's give the system a kick at 70 and 120:
 
kick <- as.matrix(data.frame(time=1:200, 
                             s.in=c(rep(0,69), 
                                    1, 
                                    rep(0,49), 1, 
                                    rep(0,80))))
inputs(osci) <- kick

mySim <-  out(sim(osci))
plot(y ~ time, data=mySim, t="l")


I hope this helps.

Best regards

Markus

-----Original Message-----
From: r-sig-dynamic-models-bounces at r-project.org [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Felix Sch?nbrodt
Sent: 09 January 2013 14:33
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] Modeling dampened oscillator with random resets

Dear list,

using deSolve (or another package): is it possible to model a dampened oscillator that is reset to new states at specific points?

Background: I want to model a system that tries to approach equilibrium, but is disturbed at single events. It would roughly look like following screenshot:
https://dl.dropbox.com/u/4472780/WP/oscReset.PNG

In this screen shot, events at t=70 and t=120 disturb the oscillator. In this screen shot I simply pasted together 3 independent plots; but if the disturbance is directly modeled it should somehow be reflected in the dx etc. (e.g., at the second event I would expect that the curve first continues going up before it is regulated down again).

Here's my function for the dampened oscillator:

library("deSolve")

	osci <- function (t, states, parms) {
	 with(as.list(c(states,  parms)), {
	   dx <- y
	   dy <- eta * x + zeta * y
	   list(c(dx, dy))
	 })
	}

	states <- c(x=1, y=0)
	out <- data.frame(ode(states, osci, times = 1:100, parms = c(eta = -0.05, zeta = -0.02)))
	plot(out)


Best,
Felix
_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
**********************************************************************
The information in this E-Mail and in any attachments is...{{dropped:27}}

Dear list,

using deSolve (or another package): is it possible to model a
dampened oscillator that is reset to new states at specific points?

Yes, this is possible, using the event mechanism:


library(deSolve)
osci <- function (t, states, parms) {
   with(as.list(c(states,  parms)), {
     dx <- y
     dy <- eta * x + zeta * y
     list(c(dx, dy))
   })
}

states <- c(x = 1, y = 0)
out <- ode(states, osci, times = 1:200,
   parms = c(eta = -0.05, zeta = -0.02))
plot(out)


###  see ?events for details
eventdat <- data.frame(var = "y", time = c(70, 120),
                        value = c(0.7, 1), method = "add")

out <- ode(states, osci, times = 1:200,
   parms = c(eta = -0.05, zeta = -0.02),
   events = list(data = eventdat))
plot(out, which = "y")
abline(v = eventdat$time, col = "red", lty = "dashed")




## Thomas P.

Thanks! Works like a charm.
Felix


Am 09.01.2013 um 16:07 schrieb Thomas Petzoldt <Thomas.Petzoldt at tu-dresden.de>:

On 09.01.2013 15:33, Felix Sch?nbrodt wrote:

Dear list,

using deSolve (or another package): is it possible to model a
dampened oscillator that is reset to new states at specific points?

Yes, this is possible, using the event mechanism:


library(deSolve)
osci <- function (t, states, parms) {
 with(as.list(c(states,  parms)), {
   dx <- y
   dy <- eta * x + zeta * y
   list(c(dx, dy))
 })
}

states <- c(x = 1, y = 0)
out <- ode(states, osci, times = 1:200,
 parms = c(eta = -0.05, zeta = -0.02))
plot(out)


###  see ?events for details
eventdat <- data.frame(var = "y", time = c(70, 120),
                      value = c(0.7, 1), method = "add")

out <- ode(states, osci, times = 1:200,
 parms = c(eta = -0.05, zeta = -0.02),
 events = list(data = eventdat))
plot(out, which = "y")
abline(v = eventdat$time, col = "red", lty = "dashed")




## Thomas P.

Thanks! Works like a charm.
Felix

Thanks to Karline Soetaert who implemented this mechanism and also to 
Woodrow Setzer. Instead of a data frame it is also possible (and even 
more powerful) to use a function.

Thomas