Thank you Dr Soetaert, that puts these in better prospective... these values though are later used in calculating further state variables in a model with 9 state variables. Eventually, the solver will "crash on me" rightfully I am sure:-) with the following error message:

DLSODA-  At T (=R1) and step size H (=R2), the
      corrector convergence failed repeatedly
      or with ABS(H) = HMIN
In above message, R =
[1] 9.600744e+01 4.166607e-08
DLSODA-  At T (=R1) and step size H (=R2), the
      corrector convergence failed repeatedly
      or with ABS(H) = HMIN
In above message, R =
[1] 1.200240e+02 7.012878e-10
DLSODA-  At T (=R1) and step size H (=R2), the
      corrector convergence failed repeatedly
      or with ABS(H) = HMIN
In above message, R =
[1] 1.200848e+02 5.380739e-08
DLSODA-  At T (=R1) and step size H (=R2), the
      corrector convergence failed repeatedly
      or with ABS(H) = HMIN
In above message, R =
[1] 9.602714e+01 3.257526e-08
DLSODA-  At T (=R1) and step size H (=R2), the
      corrector convergence failed repeatedly
      or with ABS(H) = HMIN
In above message, R =
[1] 7.208421e+01 6.005561e-08
DLSODA-  At T (=R1) and step size H (=R2), the
      corrector convergence failed repeatedly
      or with ABS(H) = HMIN...

 Would you suggest a method to avoid those negative or equal to zero (they should not be zero or bellow zero in reality) values from being generated? 

thank you for your patience and input

Andras



On Wednesday, October 23, 2013 7:44 AM, Karline Soetaert <Karline.Soetaert@nioz.nl> wrote:
 
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@r-project.org [mailto:r-sig-dynamic-models-bounces@r-project.org] On Behalf Of Andras Farkas
Sent: woensdag 23 oktober 2013 13:33
To: r-sig-dynamic-models@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]]
	[[alternative HTML version deleted]]

