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
\itemjac = "fullint" : a full jacobian, calculated internally by vode, corresponds to mf=22
\itemjac = "fullusr" : a full jacobian, specified by user function jacfunc, corresponds to mf=21
\itemjac = "bandusr" : a banded jacobian, specified by user function jacfunc; the size of the bands specified by bandup and banddown, corresponds to mf=24
\itemjac = "bandint" : a banded jacobian, calculated by vode; the size of the bands specified by 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
\itemJSV = SIGN(mf) indicates the Jacobian-saving strategy: JSV = 1 means a copy of the Jacobian is saved for reuse
in the corrector iteration algorithm. JSV = -1 means a copy of the Jacobian is not saved.
\itemMETH indicates the basic linear multistep method: METH = 1 means the implicit Adams method. METH = 2 means the method based on backward
differentiation formulas (BDF-s).
\itemMITER indicates the corrector iteration method: MITER = 0 means functional iteration (no Jacobian matrix is involved).
MITER = 1 means chord iteration with a user-supplied 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 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:
\itemel 1 : returns the conditions under which the last call to vode returned.
2 if DVODE was successful, -1 if excess work done, -2 means excess accuracy requested. (Tolerances too small),
-3 means illegal input detected. (See printed message.), -4 means repeated error test failures. (Check all input),
-5 means repeated convergence failures. (Perhaps bad Jacobian supplied or wrong choice of MF or tolerances.),
-6 means error weight became zero during problem. (Solution component i vanished, and atol or atol(i) = 0.)
\itemel 12 : The number of steps taken for the problem so far.
\itemel 13 : The number of function evaluations for the problem so far.",
\itemel 14 : The number of Jacobian evaluations so far.,
\itemel 15 : The method order last used (successfully).,
\itemel 16 : The order to be attempted on the next step.,
\itemel 17 : if el 1 =-4,-5: the largest component in the error vector,
\itemel 20 : The number of matrix LU decompositions so far.,
\itemel 21 : The number of nonlinear (Newton) iterations so far.,
\itemel 22 : The number of convergence failures of the solver so far ,
\itemel 23 : The number of error test failures of the integrator so far.)
rstate contains the following:
\item1: The step size in t last used (successfully).
\item2: The step size to be attempted on the next step.
\item3: The current value of the independent variable which the solver has actually reached, i.e. the
current internal mesh point in t.
\item4: A tolerance scale factor, greater than 1.0, computed when a request for too much accuracy was detected.
For more information, see the comments in the original code dvode.fode, 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