bvptwpC | R Documentation |
Interface to bvptwp()
bvptwpC( yini = NULL, x, func, yend = NULL, parms, xguess = NULL, yguess = NULL, ... )
yini |
named vector of type numeric. Initial values to be overwritten. |
x |
vector of type numeric. Integration times |
func |
return value from funC() with a boundary argument. |
yend |
named vector of type numeric. End values to be overwritten. |
parms |
named vector of type numeric. The dynamic parameters. |
xguess |
vector of type numeric, the x values |
yguess |
matrix with as many rows as variables and columns as x values |
... |
further arguments going to |
See bvpSolve-package for a full description of possible arguments
matrix with times and states
## Not run: ###################################################################### ## Boundary value problem: Ozon formation with fixed ozon/oxygen ratio ## at final time point ###################################################################### library(bvpSolve) # O2 + O <-> O3 # diff = O2 - O3 # build_O3 = const. f <- c( O3 = " build_O3 * O2 * O - decay_O3 * O3", O2 = "-build_O3 * O2 * O + decay_O3 * O3", O = "-build_O3 * O2 * O + decay_O3 * O3", diff = "-2 * build_O3 * O2 * O + 2 * decay_O3 * O3", build_O3 = "0" ) bound <- data.frame( name = names(f), yini = c(0, 3, 2, 3, NA), yend = c(NA, NA, NA, 0, NA) ) # Generate ODE function func <- funC(f, jacobian="full", boundary = bound, modelname = "example4") # Initialize times, states, parameters and forcings times <- seq(0, 15, by = .1) pars <- c(decay_O3 = .1) xguess <- times yguess <- matrix(c(1, 1, 1, 1, 1), ncol=length(times), nrow = length(f)) # Solve BVP out <- bvptwpC(x = times, func = func, parms = pars, xguess = xguess, yguess = yguess) # Solve BVP for different ini values, end values and parameters yini <- c(O3 = 2) yend <- c(diff = 0.2) pars <- c(decay_O3 = .01) out <- bvptwpC(yini = yini, yend = yend, x = times, func = func, parms = pars, xguess = xguess, yguess = yguess) # Plot solution par(mfcol=c(1,2)) t <- out[,1] M1 <- out[,2:5] M2 <- cbind(out[,6], pars) matplot(t, M1, type="l", lty=1, col=1:4, xlab="time", ylab="value", main="states") legend("topright", legend = c("O3", "O2", "O", "O2 - O3"), lty=1, col=1:4) matplot(t, M2, type="l", lty=1, col=1:2, xlab="time", ylab="value", main="parameters") legend("right", legend = c("build_O3", "decay_O3"), lty=1, col=1:2) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.