rootSolve (version 1.8.2.1)

steady.band: Steady-state solver for ordinary differential equations; assumes a banded jacobian

Description

Estimates the steady-state condition for a system of ordinary differential equations.

Assumes a banded jacobian matrix.

Usage

steady.band(y, time = 0, func, parms = NULL, 
            nspec = NULL, bandup = nspec, banddown = nspec, 
            times = time, ...)

Arguments

y

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, times

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').

func

either an R-function that computes the values of the derivatives in the ode system (the model defininition) at time time, or a character string giving the name of a compiled function in a dynamically loaded shared library. If func is an R-function, it must be defined as: yprime = func(t, y, parms,...). t is the current time point in the integration, 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 or list of parameters; ... (optional) are any other arguments passed to the function.

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 are global values whose steady-state value is also required.

The derivatives should be specified in the same order as the state variables y.

parms

parameters passed to func.

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 stode.

Value

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 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 the precision attained during each iteration.

Details

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

See Also

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

Examples

# NOT RUN {
## =======================================================================
## 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"))

# }