Data engineering and BI courses are free!
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 Rfunction (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 non-linear models in Rsuch 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=1e-6, atol=1e-6,
jacfunc=NULL, jactype="fullint", verbose=FALSE, 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, ...)y has a name attribute, the names will be used to label the output matrix.y are desired. The first value in times must be the initial time.func or jacfunc.y. See details.y. See details.NULL, an Rfunction that computes
the jacobian of the system of differential equations
dydot(i)/dy(j), or a string giving the name of a function or
subroutine in NULL, then lsoda cannot integrate past tcrit. The Fortran routine lsoda overshoots its targets
(times points in the vector times), and interpolates values
for the desired hmin if you don't know why!hmax is set to the largest difference in times, to avoid that the simulation possibly ignores short-term events. If 0, no maximal size is specifiedfunc ; this may speed up the simulation especially for large modelsfunc and jacfunc. See package vignette.initfunc is in the dll: the parameters passed to the initialiser, to initialise the common blocks (fortran) or global variables (C, C++)func and jacfuncfunc and jacfuncdllname is specified and the model is defined in compiled code: the number of output variables calculated in the compiled function func, present in the shared library. Note:
it is not automatically checked whetnout > 0: the names of output variables calculated in the compiled function func, present in the shared libraryfunc and jacfunc allowing this to be a generic functiontimes 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 output will have the attributes istate, and rstate, two vectors with several useful elements.
See details.
The first element of istate returns the conditions under which the last call to lsoda returned. Normal is istate[1] = 2.
If verbose = TRUE, the settings of istate and rstate will be written to the screenlsoda,
whose documentation should be consulted for details (it is included as
comments in the source file 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:
jacfunc
jacfunc; the size of the bands specified by bandup and banddown
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 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 )
$\leq$ 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:
lsoda will attempt a reasonable
reduction of precision in order to get an answer. It will write a
warning if it does so.
Models may be defined in compiled C or Fortran code, as well as in R. See package vignette for details.
Examples in both C and Fortran are in the deSolve package directory.
The output will have the attributes *istate*, and *rstate*, two vectors with several useful elements.
if verbose = TRUE, the settings of istate and rstate will be written to the screen.
the following elements of istate are meaningful:
ode, lsode, lsodes,
lsodar, vode, daspk, rk.gnls to estimate parameters in a
dynamic model.#########################################
## Example 1: Lotka-volterra model
#########################################
## A simple resource limited Lotka-Volterra-Model
## Note:
## 1. parameter and state variable names made
## accessible via "with" statement
## 2. function sigimp accessible through lexical scoping
## (see also ode and rk examples)
lvmodel <-function(t, x, parms) {
with(as.list(c(parms,x)), {
import <- sigimp(t)
dS <- import - b*S*P + g*K #substrate
dP <- c*S*P - d*K*P #producer
dK <- e*P*K - f*K #consumer
res<-c(dS, dP, dK)
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, K=1)
## Solving
out2 <- as.data.frame(lsoda(xstart, times, lvmodel, parms))
mf <- par(mfrow=c(2,2))
plot (out2$time, out2$S, type="l", ylab="substrate")
plot (out2$time, out2$P, type="l", ylab="producer")
plot (out2$time, out2$K, type="l", ylab="consumer")
plot (out2$P, out2$K, 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(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
)
}
## 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,
jacfunc = exampjac)
)
all.equal(out, outJ) # TRUERun the code above in your browser using DataLab