| lsoda {odesolve} | R Documentation |
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.
lsoda(y, times, func, parms, rtol, atol, tcrit=NULL,jac=NULL)
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). |
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:
mathbf{rtol} times mathrm{abs}(mathbf{y}) + mathbf{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.
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.
R. Woodrow Setzer setzer.woodrow@epa.gov, fixups by Martin Maechler maechler@stat.math.ethz.ch
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: http://www.netlib.org
### 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))
}