Iterative steadystate solver for ordinary differential equations (ODE) and a full or banded Jacobian.
Description
Estimates the steadystate 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 newtonraphson 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.
Usage
1 2 3 4 5 6 7 8  stode(y, time = 0, func, parms = NULL,
rtol = 1e6, atol = 1e8, ctol = 1e8,
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, ...)

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 
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 nonzero bands above the diagonal, in case the Jacobian is banded. 
banddown 
number of nonzero 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 nonnegative; 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 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 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 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. 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 bybandup
andbanddown
.jactype = "bandint" : a banded jacobian, calculated by the solver; the size of the bands specified by
bandup
andbanddown
.
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 maxnorm of ( e/ewt )
<= 1, where ewt is a vector of positive error
weights. The values of rtol
and atol
should all be
nonnegative.
The form of ewt is:
\bold{rtol} * abs(\bold{y}) + \bold{atol}
where multiplication of two vectors is elementbyelement.
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 usersupplied Rfunction, it must be called as:
yprime = func(t, y, parms,...).
t is the time
at which the steadystate 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
steadystate 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*(1O2/(O2+ks))* SO4/(SO4+ks2
dOM = Flux  oxicmin  anoxicmin
dO2 = oxicmin 2*rox*HS*(O2/(O2+ks)) + D*(BO2O2)
dSO4 = 0.5*anoxicmin +rox*HS*(O2/(O2+ks)) + D*(BSO4SO4)
dHS = 0.5*anoxicmin rox*HS*(O2/(O2+ks)) + D*(BHSHS)
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))
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).
Note
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.
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.
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.
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.
stodes
, iterative steadystate solver for ODEs with arbitrary
sparse 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 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98  ## =======================================================================
## 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*(1O2/(O2+ks))* SO4/(SO4+ks2)
dOM = Flux  oxicmin  anoxicmin
dO2 = oxicmin 2*rox*HS*(O2/(O2+ks)) + D*(BO2O2)
dSO4 = 0.5*anoxicmin +rox*HS*(O2/(O2+ks)) + D*(BSO4SO4)
dHS = 0.5*anoxicmin rox*HS*(O2/(O2+ks)) + D*(BHSHS)
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 = 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
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%*%xb, as a list:
fun < function (t, x, p) # t and p are dummies here...
list(A%*%xB)
# 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((X1X$y)^2)
