plot.bvpSolve | R Documentation |
Plot the output of boundary value solver routines.
## S3 method for class 'bvpSolve'
plot(x, ..., which = NULL, ask = NULL,
obs = NULL, obspar= list())
## S3 method for class 'bvpSolve'
print(x, ...)
x |
the output of It is allowed to pass several objects of class
|
which |
the name(s) or the index to the variables that should be plotted. Default = all variables, except the first column. |
ask |
logical; if |
obs |
a By default the first column of an observed data set should contain
the If the first column of If |
obspar |
additional graphics arguments passed to |
... |
additional arguments. The graphical arguments are passed to The dots may also contain other objects of class
The arguments after ... must be matched exactly. |
print.bvpSolve
prints the matrix without the attributes.
plot.bvpSolve
plots multiple figures on a page.
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.
diagnostics.bvpSolve
, for a description of diagnostic messages.
## =============================================================================
## The example MUSN from Ascher et al., 1995.
## Numerical Solution of Boundary Value Problems for Ordinary Differential
## Equations, SIAM, Philadelphia, PA, 1995.
##
## The problem is
## u' = 0.5*u*(w - u)/v
## v' = -0.5*(w - u)
## w' = (0.9 - 1000*(w - y) - 0.5*w*(w - u))/z
## z' = 0.5*(w - u)
## y' = -100*(y - w)
## on the interval [0 1] and with boundary conditions:
## u(0) = v(0) = w(0) = 1, z(0) = -10, w(1) = y(1)
## =============================================================================
musn <- function(t, Y, pars) {
with (as.list(Y),
{
du <- 0.5 * u * (w-u)/v
dv <- -0.5 * (w-u)
dw <- (0.9 - 1000 * (w-y) - 0.5 * w * (w-u))/z
dz <- 0.5 * (w-u)
dy <- -100 * (y-w)
return(list(c(du, dv, dw, dz, dy)))
})
}
#--------------------------------
# Boundaries
#--------------------------------
bound <- function(i,y,pars) {
with (as.list(y), {
if (i ==1) return (u-1)
if (i ==2) return (v-1)
if (i ==3) return (w-1)
if (i ==4) return (z+10)
if (i ==5) return (w-y)
})
}
xguess <- seq(0, 1, len = 5)
yguess <- matrix(ncol = 5, (rep(c(1, 1, 1, -10, 0.91), times = 5)) )
rownames(yguess) <- c("u", "v", "w", "z", "y")
sol <- bvpcol (bound = bound, x = seq(0, 1, by = 0.05),
leftbc = 4, func = musn, xguess = xguess, yguess = yguess)
mf <- par("mfrow")
plot(sol)
par(mfrow = mf)
## =============================================================================
## Example 2. Example Problem 31 from Jeff Cash's website
## =============================================================================
Prob31 <- function(t, Y, pars) {
with (as.list(Y), {
dy <- sin(Tet)
dTet <- M
dM <- -Q/xi
T <- 1/cos (Tet) +xi*Q*tan(Tet)
dQ <- 1/xi*((y-1)*cos(Tet)-M*T)
list(c( dy, dTet, dM, dQ))
})
}
ini <- c(y = 0, Tet = NA, M = 0, Q = NA)
end <- c(y = 0, Tet = NA, M = 0, Q = NA)
# run 1
xi <-0.1
twp <- bvptwp(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
func = Prob31, atol = 1e-10)
# run 2
xi <- 0.05
twp2 <- bvptwp(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
func = Prob31, atol = 1e-10)
# run 3
xi <- 0.01
twp3 <- bvptwp(yini = ini, yend = end, x = seq(0, 1, by = 0.01),
func = Prob31, atol = 1e-10)
# print all outputs at once
plot(twp, twp2, twp3, xlab = "x", ylab = names(ini))
# change type, colors, ...
plot(twp, twp2, twp3, type = c("l", "b", "p"),
main = paste ("State Variable", names(ini)),
col = c("red", "blue", "orange"), cex = 2)
## =============================================================================
## Assume we have two 'data sets':
## =============================================================================
# data set in 'wide' format
obs1 <- cbind(time = c(0, 0.5, 1), Tet = c(0.4, 0.0, -0.4))
# data set in 'long' format
obs2 <- data.frame(name = "Tet", time = c(0, 0.5, 1), value = c(0.35, 0.0, -0.35))
plot(twp, twp2, obs = obs1, obspar = list(pch = 16, cex = 1.5))
plot(twp, twp2, obs = list(obs1, obs2),
obspar = list(pch = 16, cex = 1.5))
plot(twp, twp2, obs = list(obs1, obs2), which = c("Tet", "Q"),
obspar = list(pch = 16:17, cex = 1.5, col = c("red", "black"))
)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.