[R] lsoda

Fri May 11 17:14:18 CEST 2001

The return value is a list because, in the applications for which I
originally wrote the package, it is common to have state variables whose
values are determined by the ODEs, and others that are functions of the
state variables (the second component of the return value).   For example,
state variables for a physiologically-based pharmacokinetic model are
usually amounts of chemical in the various compartments, but the rate
equations usually depend upon concentrations, and usually it is
concentrations that we measure.  Also, we like to make sure that we have
really accounted for everything, so we calculate the mass balance at each
time step, which would just be the sum of all the amounts (assuming the
model tries to account for the fate of all the incoming chemical).
Since it is often the case that the odes themselves can be written to
depend upon those global variables, it seemed more efficient to output the
global variables with the solutions of the odes, rather than have to
recompute them as a function of the ode solutions.

If you do not have any such 'global' variables, you don't need to force it!
Just return a list with one component, the ode solution.   So, for your
problem, just return:
list(c(dN1.dt,dN2.dt))
However, make sure you are returning a list, and not a simple vector --
lsoda will bomb (I've done that too many times).

As to suggestions for atol and rtol, I've always used values on the order
of 1e-6, and compared the results of runs with different values to see if
it made a difference.  I'm not sure how low you can reasonably go.  No
doubt the literature on the performance of lsoda and similar solvers has
some better recommendations.

Finally, make sure you are using at least version 0.5-4, or you will crash
the first time you specify times as t <- 1:20.

R. Woodrow Setzer, Jr.                                            Phone:
(919) 541-0128
Experimental Toxicology Division                       Fax:  (919) 541-5394
Pharmacokinetics Branch
NHEERL MD-74; US EPA; RTP, NC 27711

                    "Martin Henry H.                                                                                           
                    Stevens"                    To:     r-help at stat.math.ethz.ch                                               
                    <hstevens at rci.rutger        cc:                                                                            
                    s.edu>                      Subject:     [R] lsoda                                                         
                    Sent by:                                                                                                   
                    owner-r-help at stat.ma                                                                                       
                    th.ethz.ch                                                                                                 

                    05/11/01 09:45 AM                                                                                          

I am running R 1.2.3 with ESS5.1.18 with Windows 98.

I am trying to use lsoda in the odesolve apckage and am having problems.
Question:
The return value of the function of the system of ode's has to be a list
that includes first, the ode's and second, "a vector
          (possibly with a `names' attribute) of global values that are
required at each point in `times'."
I didn't understand from the example why the massbalance = sum(y) was
needed
for that system of equations. Why is it needed?

For my function below ('lvcomp2'), I don't know what global values lsoda
needs at each point in time (see???), other than the previous time step's
values of y and the parms. Any help?

Also, any rule of thumb advice regarding rtol and atol would be
appreciated.

p <- c(r1=.01, r2=.01, a12=0.5, a21=0.6, k1=1,k2=1)
t <- 1:20
lvcomp2 <- function(y, t, p)
{ # Lotka - Volterra model of two competing populations
dN1.dt <- p["r1"] * y[1] * (1-(y[1] + p["a12"]*y[2])/p["k1"])
dN2.dt <- p["r2"] * y[2] * (1-(y[2] + p["a21"]*y[1])/p["k2"])
list(c(dN1.dt,dN2.dt),c(???))
}
outlv <- lsoda(c(.1,.1),t, lvcomp2, p,rtol=1e-4, atol=1e-6)
outlv

TIA,
Henry

Dr. M. Henry H. Stevens
Postdoctoral Associate
Department of Ecology, Evolution, & Natural Resources
14 College Farm Road
Cook College, Rutgers University
New Brunswick, NJ 08901-8551

email: hstevens at rci.rutgers.edu
phone: 732-932-9631
fax: 732-932-8746

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.
-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
_._._

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._