lsodar {deSolve}  R Documentation 
Solver for Ordinary Differential Equations (ODE), Switching Automatically Between Stiff and Nonstiff Methods and With Root Finding
Description
Solving initial value problems for stiff or nonstiff systems of firstorder ordinary differential equations (ODEs) and including rootfinding.
The R function lsodar
provides an interface to the FORTRAN ODE
solver of the same name, written by Alan C. Hindmarsh and Linda
R. Petzold.
The system of ODE's is written as an R function or be defined in
compiled code that has been dynamically loaded.  see description of
lsoda
for details.
lsodar
differs from lsode
in two respects.
It switches automatically between stiff and nonstiff methods (similar as lsoda).
It finds the root of at least one of a set of constraint functions g(i) of the independent and dependent variables.
Two uses of lsodar
are:
To stop the simulation when a certain condition is met
To trigger events, i.e. sudden changes in one of the state variables when a certain condition is met.
when a particular condition is met.
Usage
lsodar(y, times, func, parms, rtol = 1e6, atol = 1e6,
jacfunc = NULL, jactype = "fullint", rootfunc = NULL,
verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0,
hmax = NULL, hini = 0, ynames = TRUE, maxordn = 12,
maxords = 5, bandup = NULL, banddown = NULL, maxsteps = 5000,
dllname = NULL, initfunc = dllname, initpar = parms,
rpar = NULL, ipar = NULL, nout = 0, outnames = NULL, forcings=NULL,
initforc = NULL, fcontrol=NULL, events=NULL, lags = NULL, ...)
Arguments
y 
the initial (state) values for the ODE system. If 
times 
times at which explicit estimates for 
func 
either an Rfunction that computes the values of the derivatives in the ODE system (the model definition) at time t, or a character string giving the name of a compiled function in a dynamically loaded shared library. If The return value of If 
parms 
vector or list of parameters used in 
rtol 
relative error tolerance, either a scalar or an array as
long as 
atol 
absolute error tolerance, either a scalar or an array as
long as 
jacfunc 
if not In some circumstances, supplying
If the Jacobian is a full matrix, If the Jacobian is banded, 
jactype 
the structure of the Jacobian, one of

rootfunc 
if not 
verbose 
a logical value that, when 
nroot 
only used if ‘dllname’ is specified: the number of
constraint functions whose roots are desired during the integration;
if 
tcrit 
if not 
hmin 
an optional minimum value of the integration stepsize. In
special situations this parameter may speed up computations with the
cost of precision. Don't use 
hmax 
an optional maximum value of the integration stepsize. If
not specified, 
hini 
initial step size to be attempted; if 0, the initial step size is determined by the solver. 
ynames 
logical, if 
maxordn 
the maximum order to be allowed in case the method is
nonstiff. Should be <= 12. Reduce 
maxords 
the maximum order to be allowed in case the method is stiff. Should be <= 5. Reduce maxord to save storage space. 
bandup 
number of nonzero bands above the diagonal, in case the Jacobian is banded. 
banddown 
number of nonzero bands below the diagonal, in case the Jacobian is banded. 
maxsteps 
maximal number of steps per output interval taken by the solver. 
dllname 
a string giving the name of the shared library
(without extension) that contains all the compiled function or
subroutine definitions refered to in 
initfunc 
if not 
initpar 
only when ‘dllname’ is specified and an
initialisation function 
rpar 
only when ‘dllname’ is specified: a vector with
double precision values passed to the dllfunctions whose names are
specified by 
ipar 
only when ‘dllname’ is specified: a vector with
integer values passed to the dllfunctions whose names are specified
by 
nout 
only used if 
outnames 
only used if ‘dllname’ is specified and

forcings 
only used if ‘dllname’ is specified: a list with
the forcing function data sets, each present as a twocolumned matrix,
with (time,value); interpolation outside the interval
[min( See forcings or package vignette 
initforc 
if not 
fcontrol 
A list of control parameters for the forcing functions.
See forcings or vignette 
events 
A matrix or data frame that specifies events, i.e. when the value of a state variable is suddenly changed. See events for more information. 
lags 
A list that specifies timelags, i.e. the number of steps that has to be kept. To be used for delay differential equations. See timelags, dede for more information. 
... 
additional arguments passed to 
Details
The work is done by the FORTRAN subroutine lsodar
, whose
documentation should be consulted for details (it is included as
comments in the source file ‘src/opkdmain.f’). The
implementation is based on the November, 2003 version of lsodar, from
Netlib.
lsodar
switches automatically between stiff and nonstiff
methods (similar as lsoda
). This means that the user does not
have to determine whether the problem is stiff or not, and the solver
will automatically choose the appropriate method. It always starts
with the nonstiff method.
lsodar
can find the root of at least one of a set of constraint functions
rootfunc
of the independent and dependent variables. It then returns the
solution at the root if that occurs sooner than the specified stop
condition, and otherwise returns the solution according the specified
stop condition.
Caution: Because of numerical errors in the function
rootfun
due to roundoff and integration error, lsodar
may
return false roots, or return the same root at two or more
nearly equal values of time
.
The form of the Jacobian can be specified by jactype
which can take the following values:
 jactype = "fullint":
a full Jacobian, calculated internally by lsodar, the default,
 jactype = "fullusr":
a full Jacobian, specified by user function
jacfunc
, jactype = "bandusr":
a banded Jacobian, specified by user function
jacfunc
; the size of the bands specified bybandup
andbanddown
, jactype = "bandint":
banded Jacobian, calculated by lsodar; the size of the bands specified by
bandup
andbanddown
.
If jactype
= "fullusr" or "bandusr" then the user must supply a
subroutine jacfunc
.
The input parameters rtol
, and atol
determine the
error control performed by the solver. See lsoda
for details.
The output will have the attribute iroot, if a root was found iroot is a vector, its length equal to the number of constraint functions it will have a value of 1 for the constraint function whose root that has been found and 0 otherwise.
The diagnostics of the integration can be printed to screen
by calling diagnostics
. If verbose
= TRUE
,
the diagnostics will written to the screen at the end of the integration.
See vignette("deSolve") for an explanation of each element in the vectors containing the diagnostic properties and how to directly access them.
Models may be defined in compiled C or FORTRAN code, as well as
in an Rfunction. See package vignette "compiledCode"
for details.
More information about models defined in compiled code is in the package vignette ("compiledCode"); information about linking forcing functions to compiled code is in forcings.
Examples in both C and FORTRAN are in the ‘dynload’ subdirectory
of the deSolve
package directory.
Value
A matrix of class deSolve
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 next elements 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 ‘lsodar’
returns with an unrecoverable error. If y
has a names
attribute, it will be used to label the columns of the output value.
If a root has been found, the output will have the attribute
iroot
, an integer indicating which root has been found.
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
References
Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.), NorthHolland, Amsterdam, 1983, pp. 5564.
Linda R. Petzold, Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations, Siam J. Sci. Stat. Comput. 4 (1983), pp. 136148. doi:10.1137/0904010
Kathie L. Hiebert and Lawrence F. Shampine, Implicitly Defined Output Points for Solutions of ODEs, Sandia Report SAND800180, February 1980.
Netlib: https://netlib.org
See Also

roots
for more examples on roots and events 
lsoda
,lsode
,lsodes
,vode
,daspk
for other solvers of the Livermore family, 
ode
for a general interface to most of the ODE solvers, 
ode.band
for solving models with a banded Jacobian, 
ode.1D
for integrating 1D models, 
ode.2D
for integrating 2D models, 
ode.3D
for integrating 3D models,
diagnostics
to print diagnostic messages.
Examples
## =======================================================================
## Example 1:
## from lsodar source code
## =======================================================================
Fun < function (t, y, parms) {
ydot < vector(len = 3)
ydot[1] < .04*y[1] + 1.e4*y[2]*y[3]
ydot[3] < 3.e7*y[2]*y[2]
ydot[2] < ydot[1]  ydot[3]
return(list(ydot, ytot = sum(y)))
}
rootFun < function (t, y, parms) {
yroot < vector(len = 2)
yroot[1] < y[1]  1.e4
yroot[2] < y[3]  1.e2
return(yroot)
}
y < c(1, 0, 0)
times < c(0, 0.4*10^(0:8))
out < lsodar(y = y, times = times, fun = Fun, rootfun = rootFun,
rtol = 1e4, atol = c(1e6, 1e10, 1e6), parms = NULL)
print(paste("root is found for eqn", which(attributes(out)$iroot == 1)))
print(out[nrow(out),])
diagnostics(out)
## =======================================================================
## Example 2:
## using lsodar to estimate steadystate conditions
## =======================================================================
## Bacteria (Bac) are growing on a substrate (Sub)
model < function(t, state, pars) {
with (as.list(c(state, pars)), {
## substrate uptake death respiration
dBact < gmax*eff*Sub/(Sub+ks)*Bact  dB*Bact  rB*Bact
dSub < gmax *Sub/(Sub+ks)*Bact + dB*Bact + input
return(list(c(dBact,dSub)))
})
}
## root is the condition where sum of rates of change
## is very small
rootfun < function (t, state, pars) {
dstate < unlist(model(t, state, pars)) # rate of change vector
return(sum(abs(dstate))  1e10)
}
pars < list(Bini = 0.1, Sini = 100, gmax = 0.5, eff = 0.5,
ks = 0.5, rB = 0.01, dB = 0.01, input = 0.1)
tout < c(0, 1e10)
state < c(Bact = pars$Bini, Sub = pars$Sini)
out < lsodar(state, tout, model, pars, rootfun = rootfun)
print(out)
## =======================================================================
## Example 3:
## using lsodar to trigger an event
## =======================================================================
## a state variable is decaying at a firstorder rate.
## when it reaches the value 0.1, a random amount is added.
derivfun < function (t,y,parms)
list (0.05 * y)
rootfun < function (t,y,parms)
return(y  0.1)
eventfun < function(t,y,parms)
return(y + runif(1))
yini < 0.8
times < 0:200
out < lsodar(func=derivfun, y = yini, times=times,
rootfunc = rootfun, events = list(func=eventfun, root = TRUE))
plot(out, type = "l", lwd = 2, main = "lsodar with event")