[R] lsoda with arbitrary zero thresholds (with psuedo-solution)

[Martin Henry H. Stevens]

I have  a new and less distressing, but potentially more interesting, 
problem.
I realized the major flaw my old "solution" and now have a solution 
that kind of works but is rather inelegant and I think may be 
problematic in difficult systems.
Borrowing from the lsoda example again I once again highlight the code 
that I have changed to put in place arbitrary thresholds:

parms <- c(k1=0.04, k2=1e4, k3=3e7)
my.atol <- c(1e-6,  1e-10,  1e-6)
times <- seq(0,1000)
lsexamp <- function(t, y, p)
   {
     if(y[1] < .4) yd1 <- -y[1] ### These if, else statements are new
     else yd1 <- -p["k1"] * y[1] + p["k2"] * y[2]*y[3]
     if(y[3] < .4) yd3 <- -y[3]   ### These if,else statements are new
     else yd3 <- p["k3"] * y[2]^2
     list(c(yd1,-yd1-yd3,yd3),c(massbalance=sum(y)))
   }
out <- lsoda(c(.5,0,.5),times,lsexamp, parms, rtol=1e-4, atol= my.atol, 
hmax=.1)
matplot(out[,1],out[,2:5], type="l")
out[dim(out)[1],] # The intent of my could was to cause population 1 to 
fall to zero as soon as it reached < 0.4. However, the populations 1 
and 2 reach approximations of 0 (4e-281 and 5e-11).

So, I have two questions:
Can I set thresholds in a more elegant and simpler way?
Are the approximate zero values close enough?

Thank you kindly, as ever.
Sincerely,
Hank


On Jun 9, 2004, at 12:45 PM, Martin Henry H. Stevens wrote:

> using R 2.0.0
> I am trying to do some population modeling with lsoda, where I set 
> arbitrary zero population sizes when values get close to zero, but am 
> having no luck.
> As an example of what I have tried, I use code below from the help 
> page on lsoda in which I include my modification bordered by ###
>
> parms <- c(k1=0.04, k2=1e4, k3=3e7)
> my.atol <- c(1e-6,  1e-10,  1e-6)
> times <- seq(0,)
>
> lsexamp <- function(t, y, p)
>   { ### The next line is where I try to insert the threshold
> ifelse(y < 0.4,  0, y)
> ###### all else is unchanged
>     yd1 <- -p["k1"] * y[1] + p["k2"] * y[2]*y[3]
>     yd3 <- p["k3"] * y[2]^2
>     list(c(yd1,-yd1-yd3,yd3),c(massbalance=sum(y)))
>   }
> out <- lsoda(c(.5,0,.5),times,lsexamp, parms, rtol=1e-4, atol= 
> my.atol) # Initial values differ from help page
> matplot(out[,1],out[,2:5], type="l")
> out[dim(out)[1],] # The intent of my could was to cause population 1 
> to fall to zero as soon as it reached < 0.4
>
> Any thoughts would be appreciated. Thanks!
> Hank Stevens
>
>
> Dr. Martin Henry H. Stevens, Assistant Professor
> 338 Pearson Hall
> Botany Department
> Miami University
> Oxford, OH 45056
>
> Office: (513) 529-4206
> Lab: (513) 529-4262
> FAX: (513) 529-4243
> http://www.cas.muohio.edu/botany/bot/henry.html
> http://www.muohio.edu/ecology/
> http://www.muohio.edu/botany/
> "E Pluribus Unum"
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
>
>
Dr. Martin Henry H. Stevens, Assistant Professor
338 Pearson Hall
Botany Department
Miami University
Oxford, OH 45056

Office: (513) 529-4206
Lab: (513) 529-4262
FAX: (513) 529-4243
http://www.cas.muohio.edu/botany/bot/henry.html
http://www.muohio.edu/ecology/
http://www.muohio.edu/botany/
"E Pluribus Unum"