Plot the output of steady-state solver routines.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | ```
## 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, ...)
``` |

`x ` |
an object of class For |

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

`ask ` |
logical; if |

`xyswap ` |
if |

`vertical ` |
if |

`obs ` |
a The first column of If the first column of If |

`obspar ` |
additional graphics arguments passed to |

`dimselect ` |
a |

`add.contour ` |
if |

`method ` |
the name of the plotting function to use, one of "image", "filled.contour", "contour" or "persp". |

`legend ` |
if |

`object ` |
object of class |

`...` |
additional arguments passed to the methods. The graphical arguments are passed to
For For |

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`

.

`steady.1D`

, `steady.2D`

, `steady.3D`

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 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 | ```
## =======================================================================
## 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)
``` |

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.