stode | R Documentation |
Estimates the steady-state condition for a system of ordinary differential equations (ODE) written in the form:
dy/dt = f(t,y)
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, ...)
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 |
atol |
absolute error tolerance, either a scalar or a vector, one
value for each |
ctol |
if between two iterations, the maximal change in |
jacfunc |
if not If the Jacobian is a full matrix, If the Jacobian is banded, |
jactype |
the structure of the Jacobian, one of "fullint", "fullusr", "bandusr", or "bandint" - either full or banded and estimated internally or by the user. |
verbose |
if |
bandup |
number of non-zero bands above the diagonal, in case the Jacobian is banded. |
banddown |
number of non-zero bands below the diagonal, in case the jacobian is banded. |
positive |
either a logical or a vector with indices of the state
variables that have to be non-negative; if |
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 and the model is
defined in compiled code: the number of output variables calculated in
the compiled function |
outnames |
only used if ‘dllname’ is specified and
|
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. 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 )
\leq
1, where ewt is a vector of positive error
weights. The values of rtol
and atol
should all be
non-negative.
The form of ewt is:
\mathbf{rtol} \times \mathrm{abs}(\mathbf{y}) + \mathbf{atol}
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))
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).
The implementation of stode
and substantial parts of the help file
is similar to the implementation of the integration routines (e.g.
lsode
) from package deSolve.
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.
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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.