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.

1 2 3 |

`y ` |
the initial guess of (state) values for the ODE system, a vector. |

`time, times ` |
time for which steady-state is wanted; the default is |

`func ` |
either an 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. |

`dimens ` |
a 3-valued vector with the dimensionality of the model, i.e. the number of *boxes* in x-, y- and z- direction. |

`method ` |
the solution method, one of "stodes", or "runsteady".
When |

`jactype ` |
the jacobian type - default is a regular 2-D structure where layers only interact with adjacent layers in both directions.
If the structure does not comply with this, the jacobian can be set equal to |

`cyclicBnd ` |
if not |

`names ` |
the names of the components; used to label the output, and for plotting. |

`... ` |
additional arguments passed to function stodes. See note. |

This is the method of choice to find the steady-state 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.

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.

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.

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.2D`

will return with an error message telling
that there was insufficient storage. In the second try then, increase
`lrw`

. you may need to experiment to find suitable value. The smaller the better.

The error message that says to increase `lrw`

may look like this:

1 2 3 |

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 "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.
Another attribute, called `dims`

returns a.o. the length of the work
array actually required.
This can be specified with input argument `lrw`

. See note and first example.

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.

It is likely that the estimated length of the work array (argument `lrw`

),
required for the solver stodes will be too small.
If that is the case, the solver will return with an error
saying to increase `lrw`

. The current value of the work array can be
found via the `attributes`

of the output. See first example.

Karline Soetaert <karline.soetaert@nioz.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

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 67 68 69 70 71 72 | ```
## =======================================================================
## 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.)
# stodes is used, so we should specify the size of the work array (lrw)
# We take a rather large value initially
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)
))
# the actual size of lrw is in the attributes()$dims vector.
# it is best to set lrw as small as possible
attributes(ST3)
# image plot
y <- array(dim=c(n, n, n), data = ST3$y)
filled.contour(y[ , ,n/2], color.palette = terrain.colors)
# rootSolve's image plot, looping over 3rd dimension
image(ST3, mfrow = c(4,3))
# loop over 1st dimension, contours, legends added
image(ST3, mfrow = c(2, 2), add.contour = TRUE, legend = TRUE,
dimselect = list(x = c(1, 4, 8, 10)))
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.