vode provides an interface to the
Fortran ODE solver of the same name, written by Peter N. Brown, Alan C. Hindmarsh and George D. Byrne.
The system of ODE's is written as an Rfunction 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.
vode is very similar to lsode, but uses a variable-coefficient method
rather than the fixed-step-interpolate methods in lsode.
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.vode(y, times, func, parms, rtol=1e-6, atol=1e-8,
jacfunc=NULL, jactype="fullint", mf=NULL, verbose=FALSE,
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, ...)y has a name attribute, the names will be used to label the output matrix.times must be the initial time; if only one step is to be taken; set times = NULLfunc 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 mf is not NULLjactype - provides more options than jactype - see detailsNULL, then vode cannot integrate past tcrit. The Fortran routine dvode overshoots its targets (times points in the vector times), and interpolates values
for the desired time potimes, 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 multi-D modelsfunc and jacfunc. See package vignetteinitfunc 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 functiony plus the number of "global" values returned
in the next elements of the return from func, plus an additional column (the first) for the time value.
There will be one row for each element in times unless the Fortran routine `vode' 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 screenvode, 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
mf=22
jacfunc, corresponds to mf=21
jacfunc; the size of the bands specified by bandup and banddown, corresponds to mf=24
bandup and banddown, corresponds to mf=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, -11, -12, -14, -15, -21, -22, -24, -25.
mf is a signed two-digit integer, mf = JSV*(10*METH + MITER),where
func per df/dy value).
MITER = 3 means chord iteration with an internally generated diagonal Jacobian approximation
(using 1 extra call to func per df/dy evaluation).
MITER = 4 means chord iteration with a user-supplied banded Jacobian.
MITER = 5 means chord iteration with an internally generated banded Jacobian (using ML+MU+1 extra calls to func per df/dy evaluation).
If MITER = 1 or 4, the user must supply a subroutine jacfunc.
The example for integrator lsode demonstrates how to specify both a banded and full jacobian.
The input parameters rtol, and atol determine the error
control performed by the solver.
If the request for precision exceeds the capabilities of the machine,
vode will return an error code. See lsoda for details.
Models may be defined in compiled C or Fortran code, as well as in an R-function. See package vignette for details.
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, lsoda, lsode, lsodes,
lsodar, daspk, rk.# The famous Lorenz equations: chaos in the earth's atmosphere
# Lorenz 1963. J. Atmos. Sci. 20, 130-141.
chaos<-function(t,state,parameters)
{
with(as.list(c(state)),{
dx <- -8/3*x+y*z
dy <- -10*(y-z)
dz <- -x*y+28*y-z
list(c(dx,dy,dz)) })
} # end of model
state <-c(x=1, y=1, z=1)
times <-seq(0,100,0.01)
out <-as.data.frame(vode(state,times,chaos,0))
plot(out$x,out$y,type="l",main="Lorenz butterfly")Run the code above in your browser using DataLab