Solver for Ordinary Differential Equations (ODE)
Description
Solves the initial value problem for stiff or nonstiff systems of ordinary differential equations (ODE) in the form:
dy/dt = f(t,y)
.
The R function lsode
provides an interface to the FORTRAN ODE
solver of the same name, written by Alan C. Hindmarsh and Andrew
H. Sherman.
It combines parts of the code lsodar
and can thus find the root
of at least one of a set of constraint functions g(i) of the independent
and dependent variables. This can be used to stop the simulation or to
trigger events, i.e. a sudden change in one of the state variables.
The system of ODE's is written as an R function or be defined in compiled code that has been dynamically loaded.
In contrast to lsoda
, the user has to specify whether or
not the problem is stiff and choose the appropriate solution method.
lsode
is very similar to vode
, but uses a
fixedstepinterpolate method rather than the variablecoefficient
method in vode
. In addition, in vode
it is
possible to choose whether or not a copy of the Jacobian is saved for
reuse in the corrector iteration algorithm; In lsode
, a copy is
not kept.
Usage
1 2 3 4 5 6 7 8  lsode(y, times, func, parms, rtol = 1e6, atol = 1e6,
jacfunc = NULL, jactype = "fullint", mf = NULL, rootfunc = NULL,
verbose = FALSE, nroot = 0, tcrit = NULL, hmin = 0, hmax = NULL,
hini = 0, ynames = TRUE, maxord = NULL, 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 
time sequence for which output is wanted; the first
value of 
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,

jactype 
the structure of the Jacobian, one of

mf 
the "method flag" passed to function lsode  overrules

rootfunc 
if not 
verbose 
if TRUE: full output to the screen, e.g. will
print the 
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 
maxord 
the maximum order to be allowed. 
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 list 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 lsode
, 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 lsode, from Netlib.
Before using the integrator lsode
, the user has to decide
whether or not the problem is stiff.
If the problem is nonstiff, use method flag mf
= 10, which
selects a nonstiff (Adams) method, no Jacobian used.
If the problem
is stiff, there are four standard choices which can be specified with
jactype
or mf
.
The options for jactype are
 jactype = "fullint"
a full Jacobian, calculated internally by lsode, corresponds to
mf
= 22, jactype = "fullusr"
a full Jacobian, specified by user function
jacfunc
, corresponds tomf
= 21, jactype = "bandusr"
a banded Jacobian, specified by user function
jacfunc
; the size of the bands specified bybandup
andbanddown
, corresponds tomf
= 24, jactype = "bandint"
a banded Jacobian, calculated by lsode; the size of the bands specified by
bandup
andbanddown
, corresponds tomf
= 25.
More options are available when specifying mf directly.
The
legal values of mf
are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23,
24, 25.
mf
is a positive twodigit integer, mf
=
(10*METH + MITER), where
 METH
indicates the basic linear multistep method: METH = 1 means the implicit Adams method. METH = 2 means the method based on backward differentiation formulas (BDFs).
 MITER
indicates the corrector iteration method: MITER = 0 means functional iteration (no Jacobian matrix is involved). MITER = 1 means chord iteration with a usersupplied full (NEQ by NEQ) Jacobian. MITER = 2 means chord iteration with an internally generated (difference quotient) full Jacobian (using NEQ extra calls to
func
per df/dy value). MITER = 3 means chord iteration with an internally generated diagonal Jacobian approximation (using 1 extra call tofunc
per df/dy evaluation). MITER = 4 means chord iteration with a usersupplied banded Jacobian. MITER = 5 means chord iteration with an internally generated banded Jacobian (using ML+MU+1 extra calls tofunc
per df/dy evaluation).
If MITER = 1 or 4, the user must supply a subroutine jacfunc
.
Inspection of the example below shows how to specify both a banded and full Jacobian.
The input parameters rtol
, and atol
determine the
error control performed by the solver. See lsoda
for details.
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.
lsode
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, lsode
may
return false roots, or return the same root at two or more
nearly equal values of time
.
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 ‘lsode’
returns with an unrecoverable error. If y
has a names
attribute, it will be used to label the columns of the output value.
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.
See Also

rk
, 
rk4
andeuler
for RungeKutta integrators. 
lsoda
,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.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122  ## =======================================================================
## Example 1:
## Various ways to solve the same model.
## =======================================================================
## the model, 5 state variables
f1 < function (t, y, parms) {
ydot < vector(len = 5)
ydot[1] < 0.1*y[1] 0.2*y[2]
ydot[2] < 0.3*y[1] +0.1*y[2] 0.2*y[3]
ydot[3] < 0.3*y[2] +0.1*y[3] 0.2*y[4]
ydot[4] < 0.3*y[3] +0.1*y[4] 0.2*y[5]
ydot[5] < 0.3*y[4] +0.1*y[5]
return(list(ydot))
}
## the Jacobian, written as a full matrix
fulljac < function (t, y, parms) {
jac < matrix(nrow = 5, ncol = 5, byrow = TRUE,
data = c(0.1, 0.2, 0 , 0 , 0 ,
0.3, 0.1, 0.2, 0 , 0 ,
0 , 0.3, 0.1, 0.2, 0 ,
0 , 0 , 0.3, 0.1, 0.2,
0 , 0 , 0 , 0.3, 0.1))
return(jac)
}
## the Jacobian, written in banded form
bandjac < function (t, y, parms) {
jac < matrix(nrow = 3, ncol = 5, byrow = TRUE,
data = c( 0 , 0.2, 0.2, 0.2, 0.2,
0.1, 0.1, 0.1, 0.1, 0.1,
0.3, 0.3, 0.3, 0.3, 0))
return(jac)
}
## initial conditions and output times
yini < 1:5
times < 1:20
## default: stiff method, internally generated, full Jacobian
out < lsode(yini, times, f1, parms = 0, jactype = "fullint")
## stiff method, usergenerated full Jacobian
out2 < lsode(yini, times, f1, parms = 0, jactype = "fullusr",
jacfunc = fulljac)
## stiff method, internallygenerated banded Jacobian
## one nonzero band above (up) and below(down) the diagonal
out3 < lsode(yini, times, f1, parms = 0, jactype = "bandint",
bandup = 1, banddown = 1)
## stiff method, usergenerated banded Jacobian
out4 < lsode(yini, times, f1, parms = 0, jactype = "bandusr",
jacfunc = bandjac, bandup = 1, banddown = 1)
## nonstiff method
out5 < lsode(yini, times, f1, parms = 0, mf = 10)
## =======================================================================
## Example 2:
## diffusion on a 2D grid
## partially specified Jacobian
## =======================================================================
diffusion2D < function(t, Y, par) {
y < matrix(nrow = n, ncol = n, data = Y)
dY < r*y # production
## diffusion in Xdirection; boundaries = 0concentration
Flux < Dx * rbind(y[1,],(y[2:n,]y[1:(n1),]),y[n,])/dx
dY < dY  (Flux[2:(n+1),]Flux[1:n,])/dx
## diffusion in Ydirection
Flux < Dy * cbind(y[,1],(y[,2:n]y[,1:(n1)]),y[,n])/dy
dY < dY  (Flux[,2:(n+1)]Flux[,1:n])/dy
return(list(as.vector(dY)))
}
## parameters
dy < dx < 1 # grid size
Dy < Dx < 1 # diffusion coeff, X and Ydirection
r < 0.025 # production rate
times < c(0, 1)
n < 50
y < matrix(nrow = n, ncol = n, 0)
pa < par(ask = FALSE)
## initial condition
for (i in 1:n) {
for (j in 1:n) {
dst < (i  n/2)^2 + (j  n/2)^2
y[i, j] < max(0, 1  1/(n*n) * (dst  n)^2)
}
}
filled.contour(y, color.palette = terrain.colors)
## =======================================================================
## jacfunc need not be estimated exactly
## a crude approximation, with a smaller bandwidth will do.
## Here the halfbandwidth 1 is used, whereas the true
## halfbandwidths are equal to n.
## This corresponds to ignoring the ydirection coupling in the ODEs.
## =======================================================================
print(system.time(
for (i in 1:20) {
out < lsode(func = diffusion2D, y = as.vector(y), times = times,
parms = NULL, jactype = "bandint", bandup = 1, banddown = 1)
filled.contour(matrix(nrow = n, ncol = n, out[2,1]), zlim = c(0,1),
color.palette = terrain.colors, main = i)
y < out[2, 1]
}
))
par(ask = pa)
