stodes | R Documentation |
Estimates the steady-state 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 newton-raphson 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.
stodes(y, time = 0, func, parms = NULL, rtol = 1e-6, atol = 1e-8,
ctol = 1e-8, 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, times = time, ...)
y |
the initial guess of (state) values for the ode system, a vector.
If |
time, times |
time for which steady-state is wanted; the default is
|
func |
either a user-supplied 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, steady-state is assumed to be reached. |
sparsetype |
the sparsity structure of the Jacobian, one of
"sparseint", "sparseusr", "sparsejan", "sparsereturn", ...,
The sparsity can be estimated internally by stodes (first and last option)
or given by the user (other two). See details.
Note: setting |
verbose |
if TRUE: full output to the screen, e.g. will output the steady-state 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 non-negative; if TRUE, the state variables are forced to be non-negative 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 dll-functions whose names are specified
by |
ipar |
only when ‘dllname’ is specified: a vector with integer
values passed to the dll-functions 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 two-columned 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 |
The work is done by a Fortran 77 routine that implements the Newton-Raphson 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 non-zero 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 attribute
dims
).
In this case, inz
need not be specified.
sparsetype
= "sparseusr"
. The sparsity is determined by
the user.
In this case, inz
should be a matrix
, containing indices
(row, column) to the nonzero elements in the Jacobian matrix.
The number of nonzeros nnz
will be set equal to the number of rows
in inz
.
sparsetype
= "sparsejan"
. The sparsity is also determined by
the user.
In this case, inz
should be a vector
, containting the ian
and
jan
elements of the sparse storage format, as used in the sparse solver.
Elements of ian
should be the first n+1
elements of this vector, and
contain the starting locations in jan
of columns 1.. n.
jan
contains the row indices of the nonzero locations of
the jacobian, reading in columnwise order.
The number of nonzeros nnz
will be set equal to the length of inz
- (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
functions ode.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 fill-in 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.
A list containing
y |
a vector with the state variable values from the last iteration
during estimation of steady-state condition of the system of equations.
If |
... |
the number of "global" values returned. |
The output will have the attribute steady
, which returns TRUE
,
if steady-state 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).
In case the argument sparsetype
is set to "sparsereturn",
then two extra attributes will be returned, i.e. ian
and jan
.
These can then be used to speed up subsequent calculations - see last example.
Karline Soetaert <karline.soetaert@nioz.nl>
For a description of the Newton-Raphson 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, 1145-1151.
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.
steady
, for a general interface to most of the steady-state
solvers
steady.band
, to find the steady-state of ODE models with a
banded Jacobian
steady.1D
, steady.2D
,
steady.3D
, steady-state solvers for 1-D, 2-D and 3-D
partial differential equations.
stode
, iterative steady-state solver for ODEs with full
or banded Jacobian.
runsteady
, steady-state solver by dynamically running to
steady-state
## =======================================================================
## 1000 simultaneous equations
## =======================================================================
model <- function (time, OC, parms, decay, ing)
{
# model describing C in a sediment,
# Upper boundary = imposed flux, lower boundary = zero-gradient
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 first-order 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
print(system.time(
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)
## =======================================================================
## deriving sparsity structure and speeding up calculations
## =======================================================================
sparse <- stodes(runif(N), func = model, parms = NULL,
sparsetype = "sparsereturn", decay = 5, ing = 20)
ian <- attributes(sparse)$ian
jan <- attributes(sparse)$jan
nnz <- length(jan)
inz <- c(ian, jan)
print(system.time(
s2 <- stodes(runif(N), func = model, parms = NULL, positive = TRUE,
sparsetype = "sparsejan", inz = inz, decay = 5, ing = 20)
))
# Can also be used with steady.1D, by setting jactype = "sparsejan".
# The advantage is this allows easy plotting...
print(system.time(
s2b <- steady.1D(runif(N), func = model, parms = NULL, method = "stodes",
nspec = 1, jactype = "sparsejan", inz = inz,
decay = 5, ing = 20, verbose = FALSE)
))
plot(s2b, grid = dist, xyswap = TRUE, type = "l", lwd = 2,
xlab = "Nonlocal exchange", ylab = "sediment depth",
main = "steady 1-D, sparse jacobian imposed")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.