# steady.band: Steady-state solver for ordinary differential equations;... In rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations

## Description

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

Assumes a banded jacobian matrix.

## Usage

 ```1 2 3``` ```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`.

## 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`

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

## Author(s)

Karline Soetaert <[email protected]>

`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

 ``` 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 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")) ```

rootSolve documentation built on May 29, 2017, 3:28 p.m.