steady.band | R Documentation |
Estimates the steady-state condition for a system of ordinary differential equations.
Assumes a banded jacobian matrix.
steady.band(y, time = 0, func, parms = NULL,
nspec = NULL, bandup = nspec, banddown = nspec,
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 an R-function 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 |
This is the method of choice for single-species 1-D models.
For multi-species 1-D 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 1-D *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
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 the precision attained during each iteration.
Karline Soetaert <karline.soetaert@nioz.nl>
steady
, for a general interface to most of the steady-state
solvers
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.
stodes
, iterative steady-state solver for ODEs with arbitrary
sparse Jacobian.
runsteady
, steady-state solver by dynamically running to
steady-state
## =======================================================================
## 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 first-order 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 + (1-1/(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+1-2yi+yi-1)/dx^2
d2y <- (c(y[-1],y6) -2*y + c(y1,y[-N])) /dx/dx
# dy/dx = (yi+1-yi-1)/(2dx)
dy <- (c(y[-1],y6) - c(y1,y[-N])) /2/dx
res <- d2y+dy/x+(1-1/(4*x*x))*y-sqrt(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"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.