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

[Daniel Reed]

Hi Andras:

One possible kludge to ensure positive numbers is to transform the state variable using natural logs and solve for that. The idea is built around the fact that 

d(ln x) / dt = (1/x) * (dx/dt)

In your case, the model function might look something like this:

Model<-function(t, logx, pars){
	x<-exp(logx)
	dxdt<- (-Ka * x)
	return(list(dxdt/x))
}

When solving the model, you initially pass the natural log of the state variable to ode() instead of the state variable itself. Also, you'd need to transform the solution afterwards (i.e., exp[solution]). As I say, it's a kludge, but hopefully it helps...

Cheers,
Daniel

_______________________________
Daniel C. Reed, PhD.
Postdoctoral Fellow,
Dept. of Earth & Environmental Sciences,
University of Michigan,
Ann Arbor, MI, USA.
email: reeddc at umich.edu
web: www.danielreed.org



On Oct 23, 2013, at 8:15 AM, Karline Soetaert <Karline.Soetaert at nioz.nl> wrote:

> Sometimes (but exceptionall), what I do, is in the model function itself only calculate with the max of the state variable and zero.
> This way I avoid using the tiny negative numbers, although they are generated by the solver
> 
> Karline
> 
> From: Andras Farkas [mailto:motyocska at yahoo.com]
> Sent: woensdag 23 oktober 2013 13:56
> To: Karline Soetaert; Special Interest Group for Dynamic Simulation Models in R
> Subject: Re: [R-sig-dyn-mod] hmax and its role in ode
> 
> 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 at nioz.nl<mailto:Karline.Soetaert at 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 at r-project.org<mailto:r-sig-dynamic-models-bounces at r-project.org> [mailto: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<mailto: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]]
> 
> 
> 	[[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