Explaining the options of Grind’s function

run()

A model can be solved numerically from the initial state by calling run(), and the output will be a timeplot, trajectory or table. Next to the graphics output, run() returns the final state attained by the simulation (or all data when table=TRUE). The former can be helpful if one wants to continue from the previous state (e.g., f<-run(); f<-run(state=f)). The full definition of run() is:

run <- function(tmax=100, tstep=1, state=s, parms=p, odes=model, ymin=0, ymax=NULL, log="", xlab="Time", ylab="Density", tmin=0, draw=lines, times=NULL, show=NULL, arrest=NULL, events=NULL, after=NULL, tweak=NULL, timeplot=TRUE, traject=FALSE, table=FALSE, add=FALSE, legend=TRUE, solution=FALSE, delay=FALSE, lwd=2, col="black", pch=20, ...) {}

run() calls the ode() function from the deSolve package. Additional arguments, like main and sub, are passed on to ode() and plot() via the ....

plane()

The phase plane function plane() sets up a space with the first variable on the horizontal axis, and the second on the vertical axis. The full definition of plane() is:

plane <- function(xmin=-0.001, xmax=1.05, ymin=-0.001, ymax=1.05, xlab="", ylab="", log="", npixels=500, state=s, parms=p, odes=model, x=1, y=2, time=0, grid=5, show=NULL, addone=FALSE, portrait=FALSE, vector=FALSE, add=FALSE, legend=TRUE, zero=TRUE, lwd=2, col="black", pch=20, ...) {}

Additional arguments ... are passed on to run() (for the phase portrait) and to plot() (e.g., for main and sub). Note that plane() calls the “vectorized” R-function outer(), which implies that if one calls functions in the ODEs they should also be vectorized, e.g., one should use pmax() instead of max(). There is an extension, cube.R, for 3-dimensional nullclines.

newton()

The function newton() finds a steady state from a nearby initial state, and can report the Jacobi matrix with its eigenvalues and eigenvectors. The full definition of newton() is:

newton <- function(state=s, parms=p, odes=model, time=0, positive=FALSE, jacobian=FALSE, vector=FALSE, plot=FALSE, silent=FALSE, addone=FALSE, ...) {  }

newton() calls the function steady() from the rootSolve package (which calls stode()). Additional arguments, ..., are passed on to both of them. newton() needs an initial state close to an equilibrium point.

continue()

The function continue() continues a steady state by changing a “bifurcation" parameter defined by the horizontal axis of the bifurcation diagram. The full definition of continue() is:

continue <- function(state=s, parms=p, odes=model, step=0.01, x=1, y=2, time=0, xmin=0, xmax=1,ymin=0, ymax=1.05, xlab="", ylab="", log="", col=c("red","black","blue"), lwd= c(2,1,1), addone=FALSE, positive=FALSE, nvar=FALSE, add=FALSE, ...) { }

continue() calls the function steady()from the rootSolve package (additional arguments (...) are passed on), and needs an initial state close to an equilibrium point. Note that there is much more proper software for bifurcation analysis like XPPAUT or MatCont, which reports the type of bifurcations encountered, and automatically continues all branches of branch points. Additionally, deBif is a Grind-like R-package allowing for interactive bifurcation analysis using a shiny-based interface.

fit()

The function fit() fits a model to data by non-linear parameter estimation. The output is an object containing the estimated parameters, the summed squared residuals, confidence ranges, and correlations (see the modFit() manual). The data and the model behavior for its best fit parameters are depicted in a timeplot. Its full definition is:

fit <- function(datas=data, state=s, parms=p, odes=model, free=NULL, who=NULL, differ=NULL, fixed=NULL, tmin=0, tmax=NULL, ymin=NULL, ymax=NULL, log="", xlab="Time", ylab="Density", bootstrap=0, show=NULL, fun=NULL, costfun=cost, logpar=FALSE, lower=-Inf, upper=Inf, initial=FALSE, add=FALSE, timeplot=TRUE, legend=TRUE, main=NULL, sub=NULL, pchMap=NULL, ...) { }

fit() calls the function modFit() from the FME package (which calls modCost()). Additional arguments, ..., are passed on to both of them, and to run() and ode().

timePlot()

Finally the internal function timePlot() can be used to plot data and is defined as:

timePlot <- function(data, tmin=0, tmax=NULL, ymin=0, ymax=NULL, log="", xlab="Time", ylab="Density", show=NULL, legend=TRUE, draw=lines, lwd=2, add=FALSE, main=NULL, sub=NULL, colMap=NULL, pchMap=NULL, ...) { }

These functions have many options (arguments). Fortunately most of them have good default values, and can typically be omitted. Since many options are the same in the five Grind functions, we list them alphabetically, giving their default values, and a short explanation.

Options and their defaults

More information

Check Grind’s homepage for more info.