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"
Martin Henry H. Stevens
2004-Jun-09 17:43 UTC
[R] lsoda with arbitrary zero thresholds (with psuedo-solution)
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"