[R-sig-dyn-mod] hmax and its role in ode

[Karline Soetaert]

Andras, 

you also have to look at the size of these negative values: none is smaller than -1e-8, which is, given your tolerances equal to zero.

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 Andras Farkas
Sent: woensdag 23 oktober 2013 13:33
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] hmax and its role in ode

Dear List

Please see this simple code:

plist <-data.frame(Ka=rlnorm(1000,0.3845992,0.09733002),c=rep(1,1000))
derivs <- function(t, state, pars) {
  with(as.list(pars), {
    dy1 <- - Ka * state[1]
    return(list(c(dy1)))
  })}

model <- function(pars, times=seq(0, 168, by = 1)) {
  state <- c(a = 600)
  return(ode(y = state, times = times, func = derivs, parms = plist[i,],method="lsoda",rtol = 1e-8,atol =1e-8,hmax=1)) }

modelout <- list()
for (i in 1:nrow(plist))
  modelout[[i]] <- model(plist)

b <-data.frame(modelout)
which(b<0)

as you will see, plenty of the values of state variable a in data frame b estimated as less than zero, which I would not expect for this first order process (unless my code is written wrong somewhere...). If we were to change times in the model function to seq(0, 168, by = 0.1) and hmax to 0.1, we will get different estimations, none of them below zero, which is what I am expecting from this simple first order process. Would you not expect to get the same results at the same time points though regardless of changing these output times? Any experience with this issue you may have?

All of your input would be greatly appreciated,

Sincerely,

Andras
	[[alternative HTML version deleted]]