lsoda                package:odesolve                R Documentation

_S_o_l_v_e _S_y_s_t_e_m _o_f _O_D_E (_o_r_d_i_n_a_r_y _d_i_f_f_e_r_e_n_t_i_a_l _e_q_u_a_t_i_o_n)_s.

_D_e_s_c_r_i_p_t_i_o_n:

     Solving initial value problems for stiff or non-stiff systems of
     first-order ordinary differential equations (ODEs), The R function
     `lsoda' provides an interface to the Fortran ODE solver of the
     same name, written by Linda R. Petzold and Alan C. Hindmarsh.  The
     system of ODE's is written as an R function (which may, of course,
     use `.C', `.Fortran', `.Call', etc., to call foreign code).  A
     vector of parameters is passed to the ODEs, so the solver may be
     used as part of a modeling package for ODEs, or for parameter
     estimation using any appropriate modeling tool for non-linear
     models in R such as `optim', `nls', `nlm' or `nlme'.

_U_s_a_g_e:

     lsoda(y, times, func, parms, rtol, atol, tcrit=NULL,jac=NULL)

_A_r_g_u_m_e_n_t_s:

       y: the initial values for the ode system.  If `y' has a name
          attribute, the names will be used to label the output matrix.

   times: times at which explicit estimates for `y' are desired.  The
          first value in `times' must be the initial time.

    func: a user-supplied function that computes the values of the
          derivatives in the ode system at time t.  It must be called
          as: `yprime = func(t, y, parms)'.  `t' is the current time
          point in the integration, `y' is the current estimate of the
          variables in the ode system, and `parms' is a vector of
          parameters (which may have a names attribute, desirable in a
          large system).

          The return value of func should be a list, whose first
          element is a vector containing the derivatives of `y' with
          respect to `time', and whose second element is a vector
          (possibly with a `names' attribute) of global values that are
          required at each point in `times'. 

   parms: any parameters used in `func' that should be modifiable
          without rewriting the function.

    rtol: relative error tolerance, either a scaler or an array as long
          as `y'.  See details.

    atol: absolute error tolerance, either a scaler or an array as long
          as `y'.  See details.

   tcrit: the Fortran routine `lsoda' overshoots its targets (times
          points in the vector `times'), and interpolates values for
          the desired time points.  If there is a time beyond which
          integration should not proceed (perhaps because of a
          singularity), that should be provided in `tcrit'.  Note that
          it does not make sense (though it is not an error) to include
          times in `times' past `tcrit', since the solver will stop and
          return at the last point in `times' that is earlier than
          `tcrit'.

     jac: an R function (or `NULL') that computes the jacobian of the
          system of differential equations: dydot(i)/dy(j).  In some
          circumstances, supplying `jac' can speed up the computations,
          if the system is stiff.  The calling sequence for `jac' is
          identical to that of `func'.  `jac' should return a vector
          whose ((i-1)*length(y) + j)th value is dydot(i)/dy(j). That
          is, return the matrix dydot/dy, where the ith row is the
          derivative of dy_i/dt with respect to y_j, by columns (the
          way R and Fortran store matrices).

_D_e_t_a_i_l_s:

     All the hard work is done by the Fortran subroutine `lsoda', whose
     documentation should be consulted for details (it is included as
     comments in the source file `src/lsoda.f').  This is based on the
     Feb 24, 1997 version of lsoda, from Netlib.  The following
     description of error control is adapted from that documentation
     (input arguments `rtol' and `atol', above):

     The input parameters `rtol', and `atol' determine the error
     control performed by the solver.  The solver will control the
     vector e of estimated local errors in y, according to an
     inequality of the form max-norm of ( e/ewt ) <= 1, where ewt is a
     vector of positive error weights.  The values of `rtol' and `atol'
     should all be non-negative. The form of ewt is:

                   rtol times mathrm{abs}(y) + atol

     {rtol * abs(y) + atol},  where multiplication of two vectors is
     element-by-element.

     If the request for precision exceeds the capabilities of the
     machine, the Fortran subroutine lsoda will return an error code;
     under some circumstances, the R function `lsoda' will attempt a
     reasonable reduction of precision in order to get an answer.  It
     will write a warning if it does so.

_V_a_l_u_e:

     A matrix with up to as many rows as elements in `times' and as
     many columns as elements in `y' plus the number of "global" values
     returned in the second element of the return from `func', plus and
     additional column for the time value.  There will be a row for
     each element in `times' unless the Fortran routine `lsoda' returns
     with an unrecoverable error.  If `y' has a names attribute, it
     will be used to label the columns of the output value. The output
     will have the attribute `istate' which returns the conditions
     under which the last call to lsoda returned.  See the source code
     for an explanation of those values: normal is `istate = 2'.

_A_u_t_h_o_r(_s):

     R. Woodrow Setzer setzer.woodrow@epa.gov, fixups by Martin
     Maechler maechler@stat.math.ethz.ch

_R_e_f_e_r_e_n_c_e_s:

     Hindmarsh, Alan C. (1983) ODEPACK, A Systematized Collection of
     ODE Solvers; in p.55-64 of  Stepleman, R.W. et al.[ed.] (1983)
     Scientific Computing, North-Holland, Amsterdam.

     Petzold, Linda R. (1983)  Automatic Selection of Methods for
     Solving Stiff and Nonstiff Systems of Ordinary Differential
     Equations. Siam J. Sci. Stat. Comput. 4, 136-148.

     Netlib: <URL: http://www.netlib.org>

_E_x_a_m_p_l_e_s:

     ### lsexamp -- example from lsoda source code

     ## names makes this easier to read, but may slow down execution.
     parms <- c(k1=0.04, k2=1e4, k3=3e7)
     my.atol <- c(1e-6,  1e-10,  1e-6)
     times <- c(0,4 * 10^(-1:10))
     lsexamp <- function(t, y, p)
       {
         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)))
       }
     exampjac <- function(t, y, p)
       {
         c(-p["k1"],  p["k1"],  0,

             p["k2"]*y[3],
           - p["k2"]*y[3] - 2*p["k3"]*y[2],
                            2*p["k3"]*y[2],

           p["k2"]*y[2],  -p["k2"]*y[2],  0
           )
       }

     require(odesolve)
     ## measure speed (here and below)
     system.time( 
     out <- lsoda(c(1,0,0),times,lsexamp, parms, rtol=1e-4, atol= my.atol)
     )  
     out

     ## This is what the authors of lsoda got for the example:

     ## the output of this program (on a cdc-7600 in single precision)
     ## is as follows..
     ##
     ##   at t =  4.0000e-01   y =  9.851712e-01  3.386380e-05  1.479493e-02
     ##   at t =  4.0000e+00   y =  9.055333e-01  2.240655e-05  9.444430e-02
     ##   at t =  4.0000e+01   y =  7.158403e-01  9.186334e-06  2.841505e-01
     ##   at t =  4.0000e+02   y =  4.505250e-01  3.222964e-06  5.494717e-01
     ##   at t =  4.0000e+03   y =  1.831975e-01  8.941774e-07  8.168016e-01
     ##   at t =  4.0000e+04   y =  3.898730e-02  1.621940e-07  9.610125e-01
     ##   at t =  4.0000e+05   y =  4.936363e-03  1.984221e-08  9.950636e-01
     ##   at t =  4.0000e+06   y =  5.161831e-04  2.065786e-09  9.994838e-01
     ##   at t =  4.0000e+07   y =  5.179817e-05  2.072032e-10  9.999482e-01
     ##   at t =  4.0000e+08   y =  5.283401e-06  2.113371e-11  9.999947e-01
     ##   at t =  4.0000e+09   y =  4.659031e-07  1.863613e-12  9.999995e-01
     ##   at t =  4.0000e+10   y =  1.404280e-08  5.617126e-14  1.000000e+00

     ## Using the analytic jacobian speeds up execution a little :

     system.time( 
     outJ <- lsoda(c(1,0,0),times,lsexamp, parms, rtol=1e-4, atol= my.atol,
                   jac = exampjac)
     )

     all.equal(out, outJ) # TRUE

     ### example of an application to nonlinear regression
     if (require(gnlm)) {
     t <- 1:10

     # known initial conditions
     # explicit regression formula with gamma variability
     fn1 <- ~exp(-a*t)
     y <- rgamma(10, 2, finterp(fn1)(0.1)/2)
     print(gnlr(y, dist="gamma", mu=fn1, pmu=0.2, pshape=0))

     # regression formula defined by differential equation
     dfn <- function(t, y, p) list(-p*y)
     fn2 <- ~lsoda(1, c(0,t), dfn, p)[2:11,2]
     print(finterp(fn1)(0.1))
     print(finterp(fn2)(0.1))
     print(gnlr(y, dist="gamma", mu=fn2, pmu=0.2, pshape=0))

     # estimated initial conditions
     # explicit regression formula with gamma variability
     fn3 <- ~b*exp(-a*t)
     y <- rgamma(10, 2, finterp(fn3)(c(2,0.1))/2)
     print(gnlr(y, dist="gamma", mu=fn3, pmu=c(1,0.2), pshape=0))

     # regression formula defined by differential equation
     fn4 <- ~lsoda(b, c(0,t), dfn, a)[2:11,2]
     print(finterp(fn3)(c(1,0.1)))
     print(finterp(fn4)(c(1,0.1)))
     print(gnlr(y, dist="gamma", mu=fn4, pmu=c(1,0.2), pshape=0))
     }

