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

Thu Jun 10 15:30:25 CEST 2004

Dear Hank,
Last question first:  really, only you can say for sure if 4e-281 and
5e-11 are small enough; it depends on the units you measure your state
variables in.  However, this strategy cannot get the state variables to
exactly 0.  Obviously, you could get closer to 0.0 faster by setting the
derivatives even larger in absolute value.  You may run into problems
with the solver when the derivatives are discontinuous functions of the
state variables.

There is a simple and elegant solution in theory, but not (yet)
available in odesolve.  soda has a variant called lsodar that returns
whenever a function of the state variables satisfies a given set of
conditions (in your case, you could tell lsodar to return whenever any
state variable drops below 0.4).  Once the call to lsodar returns, you'd
then reset all the state variables that were < 0.4 to 0, and restart the
integrator at that point.  I've been meaning to add lsodar to the
odesolve package for some time, but I never seem to have the week or so
of time I'd need to do it.  You can simulate the action of lsodar by
breaking your simulation in short sections, and doing a search to
identify time intervals where a state variable drops below its critical
value (that is, suppose you note that at t1 y[1] > 0.4, and at t3, y[1]
< 0.4]; then search the time interval between t1 and t3 for the value
where abs(y[1] - 0.4) < eps, say t2 ).  Stop the current integration at
t2, change the state variables, and restart.  For any given problem,
you'd probably have to experiment with time reporting intervals (the
intervals between points in the 'times' vector) so as not to miss any
important events.

Woody

On Jun 9, 2004, at 1:43 PM, Martin Henry H. Stevens wrote:

> 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. Woodrow Setzer, Jr.                        Phone: (919) 541-0128
Experimental Toxicology Division             Fax:  (919) 541-4284
Pharmacokinetics Branch
NHEERL B143-01; US EPA; RTP, NC 27711