lsoda {deSolve}  R Documentation 
Solving initial value problems for stiff or nonstiff systems of firstorder 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) or be defined in
compiled code that has been dynamically loaded. 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 nonlinear models in R such as
optim
, nls
, nlm
or
nlme
lsoda
differs from the other integrators (except lsodar
)
in that it switches automatically between stiff and nonstiff methods.
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.
lsoda(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,...)
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, or a list of symbols returned by
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 
if 
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 
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/opkdmain.f’). The implementation
is based on the 12 November 2003 version of lsoda, from Netlib.
lsoda
switches automatically between stiff and nonstiff
methods. 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.
The form of the Jacobian can be specified by jactype
which can
take the following values:
a full Jacobian, calculated internally by lsoda, the default,
a full Jacobian, specified by user function jacfunc
,
a banded Jacobian, specified by user function jacfunc
the size of the bands specified by bandup
and banddown
,
banded Jacobian, calculated by lsoda; the size of the bands
specified by bandup
and banddown
.
If jactype
= "fullusr" or "bandusr" then the user must supply a
subroutine jacfunc
.
The following description of error control is adapted from the
documentation of the lsoda source code
(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 maxnorm of ( e/ewt ) \leq
1, where ewt is a vector of positive error weights. The
values of rtol
and atol
should all be nonnegative. The
form of ewt is:
\mathbf{rtol} \times \mathrm{abs}(\mathbf{y}) +
\mathbf{atol}
where multiplication of two vectors is elementbyelement.
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.
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.
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 ‘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 ‘demo’ directory contains some examples of using
gnls
to estimate parameters in a
dynamic model.
R. Woodrow Setzer <setzer.woodrow@epa.gov>
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, NorthHolland, 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. doi:10.1137/0904010
Netlib: https://netlib.org
lsode
, which can also find a root
lsodes
, lsodar
, 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.
## =======================================================================
## Example 1:
## A simple resource limited LotkaVolterraModel
##
## Note:
## 1. parameter and state variable names made
## accessible via "with" function
## 2. function sigimp accessible through lexical scoping
## (see also ode and rk examples)
## =======================================================================
SPCmod < function(t, x, parms) {
with(as.list(c(parms, x)), {
import < sigimp(t)
dS < import  b*S*P + g*C #substrate
dP < c*S*P  d*C*P #producer
dC < e*P*C  f*C #consumer
res < c(dS, dP, dC)
list(res)
})
}
## Parameters
parms < c(b = 0.0, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)
## vector of timesteps
times < seq(0, 100, length = 101)
## external signal with rectangle impulse
signal < as.data.frame(list(times = times,
import = rep(0,length(times))))
signal$import[signal$times >= 10 & signal$times <= 11] < 0.2
sigimp < approxfun(signal$times, signal$import, rule = 2)
## Start values for steady state
y < xstart < c(S = 1, P = 1, C = 1)
## Solving
out < lsoda(xstart, times, SPCmod, parms)
## Plotting
mf < par("mfrow")
plot(out, main = c("substrate", "producer", "consumer"))
plot(out[,"P"], out[,"C"], type = "l", xlab = "producer", ylab = "consumer")
par(mfrow = mf)
## =======================================================================
## Example 2:
## 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(1e6, 1e10, 1e6)
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, yd1yd3, yd3), c(massbalance = sum(y)))
}
exampjac < function(t, y, p) {
matrix(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
), 3, 3)
}
## measure speed (here and below)
system.time(
out < lsoda(c(1, 0, 0), times, lsexamp, parms, rtol = 1e4,
atol = my.atol, hmax = Inf)
)
out
## This is what the authors of lsoda got for the example:
## the output of this program (on a cdc7600 in single precision)
## is as follows..
##
## at t = 4.0000e01 y = 9.851712e01 3.386380e05 1.479493e02
## at t = 4.0000e+00 y = 9.055333e01 2.240655e05 9.444430e02
## at t = 4.0000e+01 y = 7.158403e01 9.186334e06 2.841505e01
## at t = 4.0000e+02 y = 4.505250e01 3.222964e06 5.494717e01
## at t = 4.0000e+03 y = 1.831975e01 8.941774e07 8.168016e01
## at t = 4.0000e+04 y = 3.898730e02 1.621940e07 9.610125e01
## at t = 4.0000e+05 y = 4.936363e03 1.984221e08 9.950636e01
## at t = 4.0000e+06 y = 5.161831e04 2.065786e09 9.994838e01
## at t = 4.0000e+07 y = 5.179817e05 2.072032e10 9.999482e01
## at t = 4.0000e+08 y = 5.283401e06 2.113371e11 9.999947e01
## at t = 4.0000e+09 y = 4.659031e07 1.863613e12 9.999995e01
## at t = 4.0000e+10 y = 1.404280e08 5.617126e14 1.000000e+00
## Using the analytic Jacobian speeds up execution a little :
system.time(
outJ < lsoda(c(1, 0, 0), times, lsexamp, parms, rtol = 1e4,
atol = my.atol, jacfunc = exampjac, jactype = "fullusr", hmax = Inf)
)
all.equal(as.data.frame(out), as.data.frame(outJ)) # TRUE
diagnostics(out)
diagnostics(outJ) # shows what lsoda did internally