[R-sig-dyn-mod] Events in deSolve

Soetaert, Karline K.Soetaert at nioo.knaw.nl
Wed Jan 12 11:34:49 CET 2011


Theo, 

The problem is that, in the current implementation of events, the "event
times" must also be present in your "output times". 
And the solvers do not check this - so in your case the events are
simply ignored.

We will change that in the next release of deSolve, i.e. probably such
that the solver will give an error message when this is the case.



In your case, if you do:

times <- 0:96
y <- ode(init_cond, times, model, params, events = list(data =
eventdat))


you will see the effect of your "events"



while 

y[times %in% c(0, 74, 78, 84, 96), ] 


will give you the output only at requested times.



Karline

-----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 Theo
Nayler
Sent: Wednesday, January 12, 2011 11:26 AM
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] Events in deSolve

Hi

I have a question regarding to the implementation of events in deSolve.
I
wish to execute a number of events at specific times but obtain an
output
from the solver at different times. However, I am having trouble doing
this
and wonder if it is possible? As an example, the following does not seem
to
implement the events:

library("deSolve")

model<-function(t, R, params)
{
# Parameters
ka<-params[1]
kel<-params[2]
k12<-params[3]
k21<-params[4]

# Equations
dR<-vector(len = 3)
dR[1]<--R[1]*ka
dR[2]<-R[1]*ka-R[2]*k12+R[3]*k21-R[2]*kel
dR[3]<-R[2]*k12-R[3]*k21

list(dR)
}

params<-c(1.2, 0.09, 1.6, 1)
init_cond<-c(v1 = 1000, v2 = 0, v3 = 0)

eventdat<-data.frame(var = c(1, 1, 1, 1),
                     time = c(4, 24, 48, 72),
                     value = c(1000, 1000, 1000, 1000),
                     method = c("add", "add", "add", "add"))
times<-c(0, 74, 78, 84, 96)

y<-lsoda(init_cond, times, model, params, events = list(data =
eventdat))


I would be grateful for any suggestions. I am new to R so apologise if I
am
missing something obvious!

Kind regards
Theo

	[[alternative HTML version deleted]]

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models



More information about the R-sig-dynamic-models mailing list