plot.steady: Plot and Summary Method for steady1D, steady2D and steady3D...

In rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations

Description

Plot the output of steady-state solver routines.

Usage

## S3 method for class 'steady1D'
plot(x, ..., which = NULL, grid = NULL, 
              xyswap = FALSE, ask = NULL, 
              obs = NULL, obspar = list(), vertical = FALSE)
## S3 method for class 'steady2D'
image(x, which = NULL, add.contour = FALSE, 
              grid = NULL, ask = NULL, 
              method = "image", legend = FALSE, ...)
## S3 method for class 'steady2D'
subset(x, which = NULL, ...)
## S3 method for class 'steady3D'
image(x, which = NULL, dimselect = NULL, 
              add.contour = FALSE, grid = NULL, ask = NULL, 
              method = "image", legend = FALSE, ...)
## S3 method for class 'rootSolve'
summary(object, ...)

Arguments

x	an object of class steady1D, or steady2D as returned by the solvers steady.1D and steady.2D, and to be plotted. For steady1D objects, it is allowed to pass several objects after x (unnamed) - see second example.
which	the name(s) or the index to the variables that should be plotted. Default = all variables.
grid	For 1-D plots of output generated with steady.1D, a vector of values against which the 1-D steady-state solution has to be plotted. If NULL, then steady-state solutions are plotted against the index. for image plots of output generated with steady.2D or steady.3D: the x- and y-grid, as a list.
ask	logical; if TRUE, the user is asked before each plot, if NULL the user is only asked if more than one page of plots is necessary and the current graphics device is set interactive, see par(ask=.) and dev.interactive.
xyswap	if TRUE, then x-and y-values are swapped and the y-axis is from top to bottom. Useful for drawing vertical profiles.
vertical	if TRUE, then 1. x-and y-values are swapped, the y-axis is from top to bottom, the x-axis is on top, margin 3 and the main title gets the value of the x-axis. Useful for drawing vertical profiles; see example 2.
obs	a data.frame or matrix with "observed data" that will be added as points to the plots. obs can also be a list with multiple data.frames and/or matrices containing observed data. The first column of obs should contain the grid-variables as specified in argument grid. The other columns contain the observed values and they should have names that are known in x. If the first column of obs consists of factors, or characters (strings), then it is assumed that the data are presented in long (database) format, where the first three columns contain (name, grid, value). If obs is not NULL and which is NULL, then the variables, common to both obs and x will be plotted.
obspar	additional graphics arguments passed to points, for plotting the observed data. If obs is a list containing multiple observed data sets, then the graphics arguments can be a vector or a list (e.g. for xlim, ylim), specifying each data set separately.
dimselect	a list or NULL. The dimension over which the 3-D image loops. If NULL, will loop over the 3rd (z) dimension. This is similar as setting dimselect = list(z = 1:Nz) where Nz is the number of grid cells in the 3rd dimension; setting dimselect = list(z = seq(1, Nz, by =2)) will loop over the 3rd dimension, but every 2nd cell; dimselect = list(x = ...) or dimselect = list(y = ...) will loop over the x respectively y-dimension. See steady.3D for some examples.
add.contour	if TRUE, will add contours to the image plot.
method	the name of the plotting function to use, one of "image", "filled.contour", "contour" or "persp".
legend	if TRUE, a color legend will be drawn next to the "image", or "persp" plot.
object	object of class rootSolve whose summary has to be calculated.
...	additional arguments passed to the methods. The graphical arguments are passed to plot.default (for 1D) or image (for 2D, 3D) For plot.steady1D, the dots may contain other objects of class steady1D, as returned by steady.1D, and to be plotted on the same graphs as x - see second example. x and and these other objects should be compatible, i.e. the column names should be the same. For plot.steady1D, the arguments after ... must be matched exactly.

Details

The number of panels per page is automatically determined up to 3 x 3 (par(mfrow=c(3, 3))). This default can be overwritten by specifying user-defined settings for mfrow or mfcol. Set mfrow equal to NULL to avoid the plotting function to change user-defined mfrow or mfcol settings

Other graphical parameters can be passed as well. Parameters are vectorized, either according to the number of plots (xlab, ylab, main, sub, xlim, ylim, log, asp, ann, axes, frame.plot,panel.first,panel.last, cex.lab,cex.axis,cex.main) or according to the number of lines within one plot (other parameters e.g. col, lty, lwd etc.) so it is possible to assign specific axis labels to individual plots, resp. different plotting style. Plotting parameter ylim, or xlim can also be a list to assign different axis limits to individual plots.

Similarly, the graphical parameters for observed data, as passed by obspar can be vectorized, according to the number of observed data sets.

For steady3D objects, 2-D images are generated by looping over one of the axies; by default the third axis. See steady.3D.

See Also

steady.1D, steady.2D, steady.3D

Examples

## =======================================================================
##  EXAMPLE 1: 1D model, BOD + O2                                
## =======================================================================
## Biochemical Oxygen Demand (BOD) and oxygen (O2) dynamics
## in a river

#==================#
# Model equations  #
#==================#
O2BOD <- function(t, state, pars) {
  BOD <- state[1:N]
  O2  <- state[(N+1):(2*N)]

# BOD dynamics
  FluxBOD <-  v * c(BOD_0, BOD)  # fluxes due to water transport
  FluxO2  <-  v * c(O2_0, O2)
  
  BODrate <- r*BOD*O2/(O2+10)  # 1-st order consumption, Monod in oxygen

#rate of change = flux gradient - consumption  + reaeration (O2)
  dBOD         <- -diff(FluxBOD)/dx  - BODrate
  dO2          <- -diff(FluxO2)/dx   - BODrate + p*(O2sat-O2)

  return(list(c(dBOD = dBOD, dO2 = dO2)))
}    # END O2BOD
 
 
#==================#
# Model application#
#==================#
# parameters
dx      <- 100       # grid size, meters
v       <- 1e2       # velocity, m/day
x       <- seq(dx/2,10000,by=dx)  # m, distance from river
N       <- length(x)
r       <- 0.1       # /day, first-order decay of BOD
p       <- 0.1       # /day, air-sea exchange rate
O2sat   <- 300       # mmol/m3 saturated oxygen conc
O2_0    <- 50        # mmol/m3 riverine oxygen conc
BOD_0   <- 1500      # mmol/m3 riverine BOD concentration

# initial guess:
state <- c(rep(200,N), rep(200,N))

# running the model
out   <- steady.1D (y = state, func = O2BOD, parms = NULL,
                    nspec = 2, pos = TRUE, 
                    names = c("BOD", "O2"))

summary(out)

# output
plot(out, grid = x, type = "l", lwd = 2, 
     ylab = c("mmol/m3", "mmol O2/m3"))

# observations
obs <- matrix (ncol = 2, data = c(seq(0, 10000, 2000),
                                c(1400, 900,400,100,10,10)))

colnames(obs) <- c("Distance", "BOD")

# plot with observations
plot(out, grid = x, type = "l", lwd = 2, ylab = "mmol/m3", obs = obs, 
     pch = 16, cex = 1.5)

# similar but data in "long" format
OUT <- data.frame(name = "BOD", obs)
## Not run: 
plot(out, grid = x, type = "l", lwd = 2, ylab = "mmol/m3", obs = OBS, 
     pch = 16, cex = 1.5)

## End(Not run)

## =======================================================================
##  EXAMPLE 2: 1D model, BOD + O2 - second run                               
## =======================================================================
# new runs with different v
v       <- 50       # velocity, m/day

# running the model a second time
out2   <- steady.1D (y = state, func = O2BOD, parms = NULL,
                     nspec = 2, pos = TRUE, names = c("BOD", "O2"))


v       <- 200       # velocity, m/day

# running the model a second time
out3   <- steady.1D (y = state, func = O2BOD, parms = NULL,
                     nspec = 2, pos = TRUE, names = c("BOD", "O2"))

# output of all three scenarios at once
plot(out, out2, out3, type = "l", lwd = 2, 
     ylab = c("mmol/m3", "mmol O2/m3"), grid = x,
     obs = obs, which = c("BOD", "O2"))
  
# output of all three scenarios at once, and using vertical style
plot(out, out2, out3, type = "l", lwd = 2, vertical = TRUE,
     ylab = "Distance [m]",
     main = c("BOD [mmol/m3]", "O2 [mmol O2/m3]"), grid = x,
     obs = obs, which = c("BOD", "O2"))

# change plot pars
plot(out, out2, out3, type = "l", lwd = 2, 
     ylab = c("mmol/m3", "mmol O2/m3"), 
     grid = x, col = c("blue", "green"), log = "y",  
     obs = obs, obspar = list(pch = 16, col = "red", cex = 2))

## =======================================================================
## EXAMPLE 3: Diffusion in 2-D; zero-gradient boundary conditions
## =======================================================================

diffusion2D <- function(t,Y,par)  {
   y    <- matrix(nr=n,nc=n,data=Y)  # vector to 2-D matrix
   dY   <- -r*y        # consumption
   BND   <- rep(1,n)   # boundary concentration 

   #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   <- dY - (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    <- dY - (Flux[ ,2:(n+1)]-Flux[ ,1:n])/dy
                                              
   return(list(as.vector(dY)))
}

  # parameters
dy    <- dx <- 1   # grid size
Dy    <- Dx <- 1   # diffusion coeff, X- and Y-direction
r     <- 0.025     # consumption rate

n  <- 100
Y  <- matrix(nrow = n, ncol = n, 10.)

ST <- steady.2D(Y, func = diffusion2D, parms = NULL, pos = TRUE,
                dimens = c(n, n), lrw = 1000000,
                atol = 1e-10, rtol = 1e-10, ctol = 1e-10)
grid <- list(x = seq(dx/2, by = dx, length.out = n), 
             y = seq(dy/2, by = dy, length.out = n))
image(ST, grid = grid)
summary(ST)

rootSolve documentation built on Sept. 23, 2021, 3 a.m.