
Estimates the steady-state condition for a system of
ordinary differential equations (ODE) written in the form:
i.e. finds the values of y
for which f(t,y) = 0.
Uses a newton-raphson method, implemented in Fortran 77.
The system of ODE's is written as an R function or defined in compiled code that has been dynamically loaded.
stode(y, time = 0, func, parms = NULL,
rtol = 1e-6, atol = 1e-8, ctol = 1e-8,
jacfunc = NULL, jactype = "fullint", verbose = FALSE,
bandup = 1, banddown = 1, 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,
times = time, ...)
A list containing
a vector with the state variable values from the last iteration
during estimation of steady-state condition of the system of equations.
If y
has a names attribute, it will be used to label the output
values.
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).
the initial guess of (state) values for the ode system, a vector.
If y
has a name attribute, the names will be used to label the
output matrix.
time for which steady-state is wanted; the default is
times
=0.
(note- since version 1.7, 'times' has been added as an alias to 'time').
either a user-supplied function that computes the values of the
derivatives in the ode system (the model definition) at time
time
, or a character string giving the name of a compiled function
in a dynamically loaded shared library.
If func
is a user-supplied function, it must be called as:
yprime = func(t, y, parms, ...)
. t
is the time point
at which the steady-state is wanted, y
is the current estimate of
the variables in the ode system. If the initial values y
has a names
attribute, the names will be available inside func
. parms
is
a vector of parameters (which may have a names attribute).
The return value of func
should be a list, whose first element is a
vector containing the derivatives of y
with respect to
time
, and whose next elements (possibly with a
names
attribute) are global values that are required as
output.
The derivatives
should be specified in the same order as the state variables y
.
If func
is a string, then dllname
must give the name
of the shared library (without extension) which must be loaded
before stode()
is called. see Details for more information.
other parameters passed to func
and jacfunc
.
relative error tolerance, either a scalar or a vector, one
value for each y
.
absolute error tolerance, either a scalar or a vector, one
value for each y
.
if between two iterations, the maximal change in y
is
less than this amount, steady-state is assumed to be reached.
if not NULL
, either a user-supplied R function that
estimates the Jacobian of the system of differential equations
dydot(i)/dy(j), or a character string giving the name of a compiled
function in a dynamically loaded shared library as provided in
dllname
. In some circumstances, supplying jacfunc
can speed up the computations. The R calling sequence for
jacfunc
is identical to that of func
.
If the Jacobian is a full matrix, jacfunc
should return a matrix
dydot/dy, where the ith row contains the derivative of
If the Jacobian is banded, jacfunc
should return a matrix containing
only the nonzero bands of the jacobian, (dydot/dy), rotated row-wise.
the structure of the Jacobian, one of "fullint", "fullusr", "bandusr", or "bandint" - either full or banded and estimated internally or by the user.
if TRUE
: full output to the screen, e.g. will output
the steady-state settings.
number of non-zero bands above the diagonal, in case the Jacobian is banded.
number of non-zero bands below the diagonal, in case the jacobian is banded.
either a logical or a vector with indices of the state
variables that have to be non-negative; if TRUE
, all state
variables y
are forced to be non-negative numbers.
maximal number of iterations during one call to the solver.\
if FALSE: names of state variables are not passed to function
func
; this may speed up the simulation especially for multi-D
models.
a string giving the name of the shared library (without
extension) that contains all the compiled function or subroutine
definitions referred to in func
and jacfunc
.
if not NULL, the name of the initialisation function
(which initialises values of parameters), as provided in dllname
.
See details.
only when dllname
is specified and an initialisation
function initfunc
is in the dll: the parameters passed to the
initialiser, to initialise the common blocks (FORTRAN) or global variables
(C, C++).
only when dllname
is specified: a vector with double
precision values passed to the dll-functions whose names are specified
by func
and jacfunc
.
only when dllname
is specified: a vector with integer
values passed to the dll-functions whose names are specified by func
and jacfunc
.
only used if dllname
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 whether this is indeed the number of
output variables calculated in the dll - you have to perform this check
in the code - see package vignette.
only used if dllname
is specified and
nout
> 0: the names of output variables calculated in the
compiled function func
, present in the shared library.
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(times
), max(times
)] is done by
taking the value at the closest data extreme.
This feature is here for compatibility with models defined in compiled code
from package deSolve; see deSolve's package vignette "compiledCode"
.
if not NULL
, the name of the forcing function
initialisation function, as provided in
dllname
. It MUST be present if forcings
has been given a
value.
See deSolve's package vignette "compiledCode"
.
A list of control parameters for the forcing functions.
See deSolve's package vignette "compiledCode"
.
additional arguments passed to func
and jacfunc
allowing this to be a generic function.
Karline Soetaert <karline.soetaert@nioz.nl>
The work is done by a Fortran 77 routine that implements the Newton-Raphson method. It uses code from LINPACK.
The form of the Jacobian can be specified by jactype
which can
take the following values:
jactype = "fullint" : a full jacobian, calculated internally by the solver, the default.
jactype = "fullusr" : a full jacobian, specified by user function
jacfunc
.
jactype = "bandusr" : a banded jacobian, specified by user function
jacfunc
; the size of the bands specified by bandup
and
banddown
.
jactype = "bandint" : a banded jacobian, calculated by the solver;
the size of the bands specified by bandup
and banddown
.
if jactype
= "fullusr" or "bandusr" then the user must supply a
subroutine jacfunc
.
The input parameters rtol
, atol
and ctol
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 )
rtol
and atol
should all be
non-negative.
The form of ewt is:
where multiplication of two vectors is element-by-element.
In addition, the solver will stop if between two iterations, the maximal
change in the values of y is less than ctol
.
Models may be defined in compiled C or Fortran code, as well as in R.
If func
or jacfunc
are a string, then they are assumed to be
compiled code.
In this case, dllname
must give the name of the shared library
(without extension) which must be loaded before stode()
is called.
See vignette("rooSolve") for how a model has to be specified in compiled code. Also, vignette("compiledCode") from package deSolve contains examples of how to do this.
If func
is a user-supplied R-function, it must be called as:
yprime = func(t, y, parms,...).
t is the time
at which the steady-state should be estimated,
y
is the current estimate of the variables in the ode system.
The return value of func
should be a list, whose first element is
a vector containing the derivatives of y
with respect to time
,
and whose next elements contains output variables whose values at
steady-state are also required.
An example is given below:
model<-function(t,y,pars)
{
with (as.list(c(y,pars)),{
Min = r*OM
oxicmin = Min*(O2/(O2+ks))
anoxicmin = Min*(1-O2/(O2+ks))* SO4/(SO4+ks2
dOM = Flux - oxicmin - anoxicmin
dO2 = -oxicmin -2*rox*HS*(O2/(O2+ks)) + D*(BO2-O2)
dSO4 = -0.5*anoxicmin +rox*HS*(O2/(O2+ks)) + D*(BSO4-SO4)
dHS = 0.5*anoxicmin -rox*HS*(O2/(O2+ks)) + D*(BHS-HS)
list(c(dOM,dO2,dSO4,dHS),SumS=SO4+HS)
})
}
This model can be solved as follows:
pars <- c(D=1,Flux=100,r=0.1,rox =1,
ks=1,ks2=1,BO2=100,BSO4=10000,BHS = 0)
y<-c(OM=1,O2=1,SO4=1,HS=1)
ST <- stode(y=y,func=model,parms=pars,pos=TRUE))
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.
The algorithm uses LINPACK code:
Dongarra, J.J., J.R. Bunch, C.B. Moler and G.W. Stewart, 1979. LINPACK user's guide, SIAM, Philadelphia.
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.
stodes
, iterative steady-state solver for ODEs with arbitrary
sparse Jacobian.
runsteady
, steady-state solver by dynamically running to
steady-state
## =======================================================================
## Example 1. A simple sediment biogeochemical model
## =======================================================================
model<-function(t, y, pars)
{
with (as.list(c(y, pars)),{
Min = r*OM
oxicmin = Min*(O2/(O2+ks))
anoxicmin = Min*(1-O2/(O2+ks))* SO4/(SO4+ks2)
dOM = Flux - oxicmin - anoxicmin
dO2 = -oxicmin -2*rox*HS*(O2/(O2+ks)) + D*(BO2-O2)
dSO4 = -0.5*anoxicmin +rox*HS*(O2/(O2+ks)) + D*(BSO4-SO4)
dHS = 0.5*anoxicmin -rox*HS*(O2/(O2+ks)) + D*(BHS-HS)
list(c(dOM, dO2, dSO4, dHS), SumS = SO4+HS)
})
}
# parameter values
pars <- c(D = 1, Flux = 100, r = 0.1, rox = 1,
ks = 1, ks2 = 1, BO2 = 100, BSO4 = 10000, BHS = 0)
# initial conditions
y<-c(OM = 1, O2 = 1, SO4 = 1, HS = 1)
# direct iteration - enforces positivitiy..
ST <- stode(y = y, func = model, parms = pars, pos = TRUE)
ST
## =======================================================================
## Example 2. 1000 simultaneous equations
## =======================================================================
model <- function (time, OC, parms, decay, ing) {
# model describing organic Carbon (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
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 <- stode(runif(N), func = model, parms = NULL, positive = TRUE,
decay = 5, ing = 20)))
plot(ss$y[1:N], dist, ylim = rev(range(dist)), type = "l", lwd = 2,
xlab = "Nonlocal exchange", ylab = "sediment depth",
main = "stode, full jacobian")
## =======================================================================
## Example 3. Solving a system of linear equations
## =======================================================================
# this example is included to demonstrate how to use the "jactype" option.
# (and that stode is quite efficient).
A <- matrix(nrow = 500, ncol = 500, runif(500*500))
B <- 1:500
# this is how one would solve this in R
print(system.time(X1 <- solve(A, B)))
# to use stode:
# 1. create a function that receives the current estimate of x
# and that returns the difference A%*%x-b, as a list:
fun <- function (t, x, p) # t and p are dummies here...
list(A%*%x-B)
# 2. jfun returns the Jacobian: here this equals "A"
jfun <- function (t, x, p) # all input parameters are dummies
A
# 3. solve with jactype="fullusr" (a full Jacobian, specified by user)
print (system.time(
X <- stode(y = 1:500, func = fun, jactype = "fullusr", jacfunc = jfun)
))
# the results are the same (within precision)
sum((X1-X$y)^2)
Run the code above in your browser using DataLab