steady.3D {rootSolve} | R Documentation |
Estimates the steady-state condition for a system of ordinary differential equations that result from 3-Dimensional partial differential equation models that have been converted to ODEs by numerical differencing.
It is assumed that exchange occurs only between adjacent layers.
steady.3D(y, time=0, func, parms=NULL, nspec=NULL, dimens, cyclicBnd = NULL, ...)
y |
the initial guess of (state) values for the ODE system, a vector. |
time |
time for which steady-state is wanted; the default is time=0. |
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. |
dimens |
a 3-valued vector with the dimensionality of the model, i.e. the number of *boxes* in x-, y- and z- direction. |
cyclicBnd |
if not NULL then a number or a 3-valued vector
with the dimensions where a cyclic boundary is used - 1 : x-dimension,
2 : y-dimension; 3 : z-dimension;see details.
|
... |
additional arguments passed to function stodes .
|
This is the method of choice for 3-dimensional models, that are only subjected to transport between adjacent layers.
Based on the dimension of the problem, the method first calculates the
sparsity pattern of the Jacobian, under the assumption
that transport is only occurring between adjacent layers.
Then stodes
is called to find the steady-state.
As stodes
is used, it will probably be necessary to specify the
length of the real work array, lrw
.
Although a reasonable guess of lrw
is made, it may occur that this
will be too low.
In this case, steady.3D
will return with an error message telling
the size of the work array actually needed. In the second try then, set
lrw
equal to this number.
In some cases, a cyclic boundary condition exists. This is when the first
boxes in x-, y-, or z-direction interact with the last boxes.
In this case, there will be extra non-zero fringes in the Jacobian which
need to be taken into account. The occurrence of cyclic boundaries can be
toggled on by specifying argument cyclicBnd
. For innstance,
cyclicBnd = 1
indicates that a cyclic boundary is required only for
the x-direction, whereas cyclicBnd = c(1,2)
imposes a cyclic boundary
for both x- and y-direction. The default is no cyclic boundaries.
See stodes
for the additional options.
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. |
... |
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.
It is advisable though not mandatory to specify BOTH nspec
and
dimens
. In this case, the solver can check whether the input makes
sense (as nspec*dimens[1]*dimens[2]*dimens[3] = length(y))
do NOT use this method for problems that are not 3D.
Karline Soetaert <k.soetaert@nioo.knaw.nl>
steady
, for a general interface to most of the steady-state
solvers
steady.band
, to find the steady-state of ODE models with a
banded Jacobian
steady.1D
, steady.2D
,
steady-state solvers for 1-D and 2-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
## ======================================================================= ## Diffusion in 3-D; imposed boundary conditions ## ======================================================================= diffusion3D <- function(t,Y,par) { yy <- array(dim=c(n,n,n),data=Y) # vector to 3-D array dY <- -r*yy # consumption BND <- rep(1,n) # boundary concentration for (i in 1:n) { y <- yy[i,,] #diffusion in X-direction; boundaries=imposed concentration Flux <- -Dy * rbind(y[1,]-BND,(y[2:n,]-y[1:(n-1),]),BND-y[n,])/dy dY[i,,] <- dY[i,,] - (Flux[2:(n+1),]-Flux[1:n,])/dy #diffusion in Y-direction Flux <- -Dz * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dz dY[i,,] <- dY[i,,] - (Flux[,2:(n+1)]-Flux[,1:n])/dz } for (j in 1:n) { y <- yy[,j,] #diffusion in X-direction; boundaries=imposed concentration Flux <- -Dx * rbind(y[1,]-BND,(y[2:n,]-y[1:(n-1),]),BND-y[n,])/dx dY[,j,] <- dY[,j,] - (Flux[2:(n+1),]-Flux[1:n,])/dx #diffusion in Y-direction Flux <- -Dz * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dz dY[,j,] <- dY[,j,] - (Flux[,2:(n+1)]-Flux[,1:n])/dz } for (k in 1:n) { y <- yy[,,k] #diffusion in X-direction; boundaries=imposed concentration Flux <- -Dx * rbind(y[1,]-BND,(y[2:n,]-y[1:(n-1),]),BND-y[n,])/dx dY[,,k] <- dY[,,k] - (Flux[2:(n+1),]-Flux[1:n,])/dx #diffusion in Y-direction Flux <- -Dy * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dy dY[,,k] <- dY[,,k] - (Flux[,2:(n+1)]-Flux[,1:n])/dy } return(list(as.vector(dY))) } # parameters dy <- dx <- dz <-1 # grid size Dy <- Dx <- Dz <-1 # diffusion coeff, X- and Y-direction r <- 0.025 # consumption rate n <- 10 y <- array(dim=c(n,n,n),data=10.) print(system.time( ST3 <- steady.3D(y, func=diffusion3D, parms=NULL, pos=TRUE, dimens=c(n,n,n), lrw=100000, atol=1e-10, rtol=1e-10, ctol=1e-10, verbose=TRUE) )) y <- array(dim=c(n,n,n),data=ST3$y) filled.contour(y[,,n/2],color.palette=terrain.colors)