Steadystate solver for ordinary differential equations (ODE) with a sparse jacobian.
Description
Estimates the steadystate condition for a system of ordinary differential equations (ODE) in the form:
dy/dt = f(t,y)
and where the jacobian matrix df/dy has an arbitrary sparse structure.
Uses a newtonraphson method, implemented in Fortran.
The system of ODE's is written as an R function or defined in compiled code that has been dynamically loaded.
Usage
1 2 3 4 5 6 7 8  stodes(y, time = 0, func, parms = NULL, rtol = 1e6, atol = 1e8,
ctol = 1e8, sparsetype = "sparseint", verbose = FALSE,
nnz = NULL, inz = NULL, lrw = NULL, ngp = NULL,
positive = FALSE, maxiter = 100, ynames = TRUE,
dllname = NULL, initfunc = dllname, initpar = parms,
rpar = NULL, ipar = NULL, nout = 0, outnames = NULL,
forcings = NULL, initforc = NULL, fcontrol = NULL,
spmethod = "yale", control = NULL,...)

Arguments
y 
the initial guess of (state) values for the ode system, a vector.
If 
time 
time for which steadystate is wanted; the default is

func 
either a usersupplied function that computes the values of the
derivatives in the ode system (the model definition) at time
If The return value of The derivatives
should be specified in the same order as the state variables If 
parms 
other parameters passed to 
rtol 
relative error tolerance, either a scalar or a vector, one value for each y. 
atol 
absolute error tolerance, either a scalar or a vector, one value for each y. 
ctol 
if between two iterations, the maximal change in y is less than this amount, steadystate is assumed to be reached. 
sparsetype 
the sparsity structure of the Jacobian, one of "sparseint" or "sparseusr", "sparsejan", ..., The sparsity can be estimated internally by stodes (first option) or given by the user (last two). See details. 
verbose 
if TRUE: full output to the screen, e.g. will output the steadystate settings. 
nnz 
the number of nonzero elements in the sparse Jacobian (if
this is unknown, use an estimate); If NULL, a guess will be made, and
if not sufficient, If a solution is found, the minimal value of 
inz 
if 
lrw 
the length of the work array of the solver; due to the sparsicity,
this cannot be readily predicted. If If a solution is found, the minimal value of In case of an error induced by a too small value of 
ngp 
number of groups of independent state variables. Due to the
sparsicity, this cannot be readily predicted. If NULL, a guess will be
made, and if not sufficient, If a solution is found, the minimal value of 
positive 
either a logical or a vector with indices of the state variables that have to be nonnegative; if TRUE, the state variables are forced to be nonnegative numbers. 
maxiter 
maximal number of iterations during one call to the solver. 
ynames 
if FALSE: names of state variables are not passed to
function 
dllname 
a string giving the name of the shared library (without
extension) that contains all the compiled function or subroutine
definitions referred to in 
initfunc 
if not NULL, the name of the initialisation function (which initialises values of parameters), as provided in ‘dllname’. See details. 
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 ‘dllname’ is specified: the number of output
variables calculated in the compiled function 
outnames 
only used if ‘dllname’ is specified and

spmethod 
the sparse method to be used, one of 
control 
only used if 
forcings 
only used if ‘dllname’ is specified: a vector with the
forcing function values, or a list with the forcing function data sets,
each present as a twocolumned matrix, with (time,value); interpolation
outside the interval [min( This feature is here for compatibility with models defined in compiled code
from package deSolve; see deSolve's package vignette 
initforc 
if not 
fcontrol 
A list of control parameters for the forcing functions.
See deSolve's package vignette 
... 
additional arguments passed to 
Details
The work is done by a Fortran 77 routine that implements the NewtonRaphson method.
stodes
is to be used for problems, where the Jacobian has a sparse
structure.
There are several choices for the sparsity specification, selected by
argument sparsetype
.

sparsetype
="sparseint"
. The sparsity is estimated by the solver, based on numerical differences. In this case, it is advisable to provide an estimate of the number of nonzero elements in the Jacobian (nnz
). This value can be approximate; upon return the number of nonzero elements actually required will be known (1st element of attributedims
). In this case,inz
need not be specified. 
sparsetype
="sparseusr"
. The sparsity is determined by the user. In this case,inz
should be amatrix
, containing indices (row, column) to the nonzero elements in the Jacobian matrix. The number of nonzerosnnz
will be set equal to the number of rows ininz
. 
sparsetype
="sparsejan"
. The sparsity is also determined by the user. In this case,inz
should be avector
, containting theian
andjan
elements of the sparse storage format, as used in the sparse solver. Elements ofian
should be the firstn+1
elements of this vector, and contain the starting locations injan
of columns 1.. n.jan
contains the row indices of the nonzero locations of the jacobian, reading in columnwise order. The number of nonzerosnnz
will be set equal to the length ofinz
 (n+1). 
sparsetype
="1D"
,"2D"
,"3D"
. The sparsity is estimated by the solver, based on numerical differences. Assumes finite differences in a 1D, 2D or 3D regular grid  used by functionsode.1D
,ode.2D
,ode.3D
. Similar are"2Dmap"
, and"3Dmap"
, which also include a mapping variable (passed in nnz).
The Jacobian itself is always generated by the solver (i.e. there is no provision to provide an analytic Jacobian).
This is done by perturbing simulataneously a combination of state variables that do not affect each other.
This significantly reduces computing time. The number of groups with
independent state variables can be given by ngp
.
The input parameters rtol
, atol
and ctol
determine
the error control performed by the solver. See help for stode
for details.
Models may be defined in compiled C or Fortran code, as well as in R. See package vignette for details on how to write models in compiled code.
When the spmethod
equals ilut
or ilutp
, a number of parameters
can be specified in argument control
. They are:
fillin, the fillin parameter. Each row of L and each row of U will have a maximum of lfil elements (excluding the diagonal element). lfil must be >= 0.
droptol, sets the threshold for dropping small terms in the factorization.
When ilutp
is chosen the following arguments can also be specified:
permtol = tolerance ratio used to determne whether or not to permute two columns. At step i columns i and j are permuted when abs(a(i,j))*permtol .gt. abs(a(i,i)) [0 –> never permute; good values 0.1 to 0.01]
lenplufac = sets the working array  increase its value if a warning.
Value
A list containing
y 
a vector with the state variable values from the last iteration
during estimation of steadystate condition of the system of equations.
If 
... 
the number of "global" values returned. 
The output will have the attribute steady
, which returns TRUE
,
if steadystate has been reached and the attribute precis
with an
estimate of the precision attained during each iteration, the mean absolute
rate of change (sum(abs(dy))/n).
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
References
For a description of the NewtonRaphson method, e.g.
Press, WH, Teukolsky, SA, Vetterling, WT, Flannery, BP, 1996. Numerical Recipes in FORTRAN. The Art of Scientific computing. 2nd edition. Cambridge University Press.
When spmethod = "yale" then the algorithm uses linear algebra routines from the Yale sparse matrix package:
Eisenstat, S.C., Gursky, M.C., Schultz, M.H., Sherman, A.H., 1982. Yale Sparse Matrix Package. i. The symmetric codes. Int. J. Num. meth. Eng. 18, 11451151.
else the functions ilut and ilutp from sparsekit package are used:
Yousef Saad, 1994. SPARSKIT: a basic tool kit for sparse matrix computations. VERSION 2
Yousef Saad, 2003. Iterative methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics.
See Also
steady
, for a general interface to most of the steadystate
solvers
steady.band
, to find the steadystate of ODE models with a
banded Jacobian
steady.1D
, steady.2D
,
steady.3D
, steadystate solvers for 1D, 2D and 3D
partial differential equations.
stode
, iterative steadystate solver for ODEs with full
or banded Jacobian.
runsteady
, steadystate solver by dynamically running to
steadystate
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  ## =======================================================================
## 1000 simultaneous equations
## =======================================================================
model < function (time, OC, parms, decay, ing)
{
# model describing C in a sediment,
# Upper boundary = imposed flux, lower boundary = zerogradient
Flux < v * c(OC[1] ,OC) + # advection
Kz*diff(c(OC[1],OC,OC[N]))/dx # diffusion;
Flux[1]< flux # imposed flux
# Rate of change= Flux gradient and firstorder consumption
dOC < diff(Flux)/dx  decay*OC
# Fraction of OC in first 5 layers is translocated to mean depth
# (layer N/2)
dOC[1:5] < dOC[1:5]  ing*OC[1:5]
dOC[N/2] < dOC[N/2] + ing*sum(OC[1:5])
list(dOC)
}
v < 0.1 # cm/yr
flux < 10
dx < 0.01
N < 1000
dist < seq(dx/2, by = dx, len = N)
Kz < 1 #bioturbation (diffusion), cm2/yr
ss < stodes(runif(N), func = model, parms = NULL,
positive = TRUE, decay = 5, ing = 20, verbose = TRUE)
plot(ss$y[1:N], dist, ylim = rev(range(dist)), type = "l", lwd = 2,
xlab = "Nonlocal exchange", ylab = "sediment depth",
main = "stodes, sparse jacobian")
# the size of lrw is in the attributes()$dims vector.
attributes(ss)
