Steadystate solver for ordinary differential equations; assumes a banded jacobian
Description
Estimates the steadystate condition for a system of ordinary differential equations.
Assumes a banded jacobian matrix.
Usage
1 2  steady.band(y, time = 0, func, parms = NULL,
nspec = NULL, bandup = nspec, banddown = nspec, ...)

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 an Rfunction that computes the values of the
derivatives in the ode system (the model defininition) at time The return value of The derivatives
should be specified in the same order as the state variables 
parms 
parameters passed to 
nspec 
the number of *species* (components) in the model. 
bandup 
the number of nonzero bands above the Jacobian diagonal. 
banddown 
the number of nonzero bands below the Jacobian diagonal. 
... 
additional arguments passed to function 
Details
This is the method of choice for singlespecies 1D models.
For multispecies 1D models, this method can only be used if the state variables are arranged per box, per species (e.g. A[1],B[1],A[2],B[2],A[3],B[3],.... for species A, B).
Usually a 1D *model* function will have the species arranged as
A[1],A[2],A[3],....B[1],B[2],B[3],....
in this case, use steady.1D
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 the precision attained during each iteration.
Author(s)
Karline Soetaert <karline.soetaert@nioz.nl>
See Also
steady
, for a general interface to most of the steadystate
solvers
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.
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  ## =======================================================================
## 1000 simultaneous equations, solved 6 times for different
## values of parameter "decay"
## =======================================================================
model < function (time, OC, parms, decay) {
# model of particles (OC) that sink out of the water and decay
# exponentially declining sinking rate, maximal 100 m/day
sink < 100 * exp(0.01*dist)
# boundary flux; upper boundary=imposed concentration (100)
Flux < sink * c(100 ,OC)
# Rate of change= Flux gradient and firstorder consumption
dOC < diff(Flux)/dx  decay*OC
list(dOC, maxC = max(OC))
}
dx < 1 # thickness of boxes
dist < seq(0, 1000, by = dx) # water depth at each modeled box interface
ss < NULL
for (decay in seq(from = 0.1, to = 1.1, by = 0.2))
ss < cbind(ss, steady.band(runif(1000), func = model,
parms = NULL, nspec = 1, decay = decay)$y)
matplot(ss, 1:1000, type = "l", lwd = 2, main = "steady.band",
ylim=c(1000, 0), ylab = "water depth, m",
xlab = "concentration of sinking particles")
legend("bottomright", legend = seq(from = 0.1, to = 1.1, by = 0.2),
lty = 1:10, title = "decay rate", col = 1:10, lwd = 2)
## =======================================================================
## 5001 simultaneous equations: solve
## dy/dt = 0 = d2y/dx2 + 1/x*dy/dx + (11/(4x^2)y  sqrx(x)*cos(x),
## over the interval [1,6], with boundary conditions: y(1)=1, y(6)=0.5
## =======================================================================
derivs < function(t, y, parms, x, dx, N, y1, y6) {
# Numerical approximation of derivates:
# d2y/dx2 = (yi+12yi+yi1)/dx^2
d2y < (c(y[1],y6) 2*y + c(y1,y[N])) /dx/dx
# dy/dx = (yi+1yi1)/(2dx)
dy < (c(y[1],y6)  c(y1,y[N])) /2/dx
res < d2y+dy/x+(11/(4*x*x))*ysqrt(x)*cos(x)
return(list(res))
}
dx < 0.001
x < seq(1,6,by=dx)
N < length(x)
y < steady.band(y = rep(1, N), time = 0, func = derivs, x = x, dx = dx,
N = N, y1 = 1, y6 = 0.5, nspec = 1)$y
plot(x, y, type = "l", main = "5001 nonlinear equations  banded Jacobian")
# add the analytic solution for comparison:
xx < seq(1, 6, by = 0.1)
ana < 0.0588713*cos(xx)/sqrt(xx)+1/4*sqrt(xx)*cos(xx)+
0.740071*sin(xx)/sqrt(xx)+1/4*xx^(3/2)*sin(xx)
points(xx, ana)
legend("topright", pch = c(NA, 1), lty = c(1, NA),
c("numeric", "analytic"))
