\documentclass[article,nojss]{jss} \DeclareGraphicsExtensions{.pdf,.eps} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Add-on packages and fonts \usepackage{graphicx} \usepackage{amsmath} \newcommand{\noun}[1]{\textsc{#1}} %% Bold symbol macro for standard LaTeX users \providecommand{\boldsymbol}[1]{\mbox{\boldmath $#1$}} %% Because html converters don't know tabularnewline \providecommand{\tabularnewline}{\\} \usepackage{array} % tabel commands \setlength{\extrarowheight}{0.1cm} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands. \newcommand{\ls}{\textbf{\textsf{limSolve }}} \newcommand{\rs}{\textbf{\textsf{rootSolve }}} \newcommand{\rt}{\textbf{\textsf{ReacTran }}} \newcommand{\bs}{\textbf{\textsf{bvpSolve }}} \newcommand{\R}{\proglang{R }} \newcommand{\ds}{\textbf{\textsf{deSolve }}} \newcommand{\rb}[1]{\raisebox{1.5ex}{#1}} \title{Package \rs: roots, gradients and steady-states in \proglang{R}} \Plaintitle{Package rootSolve: roots, gradients and steady-states in R} \Keywords{roots of nonlinear equations, gradient, Jacobian, Hessian, steady-state, boundary value ODE, method of lines, \proglang{R}} \Plainkeywords{roots of nonlinear equations, gradient, Jacobian, Hessian, steady-state, boundary value ODE, method of lines, R} \author{Karline Soetaert\\ Royal Netherlands Institute of Sea Research (NIOZ)\\ Yerseke\\ The Netherlands} \Plainauthor{Karline Soetaert} \Abstract{ \R package \rs \citep{rootSolve} includes \begin{itemize} \item root-finding algorithms to solve for the roots of n nonlinear equations, using a Newton-Raphson method. \item An extension of \R function \code{uniroot} \item Functions that find the steady-state condition of a set of ordinary differential equations (ODE). These functions are compatible with the solvers in package \ds \citep{deSolve}, \citep{deSolve_jss}, which solve initial value differential equatons. Separate solvers for full, banded and generally sparse problems are included. These allow to estimate the steady-state of 1-D, 2-D and 3-D partial differential equations (PDE) that have been rewritten as a set of ODEs by numerical differencing, e.g. using the functions available in package \rt \citep{ReacTran}. \item Functions that calculate the Hessian and Jaobian matrix or - more general - the gradient of functions with respect to independent variables. \end{itemize} } %% The address of (at least) one author should be given %% in the following format: \Address{ Karline Soetaert\\ Royal Netherlands Institute of Sea Research (NIOZ)\\ 4401 NT Yerseke, Netherlands E-mail: \email{karline.soetaert@nioz.nl}\\ URL: \url{http://www.nioz.nl}\\ } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R/Sweave specific LaTeX commands. %% need no \usepackage{Sweave} %\VignetteIndexEntry{roots, gradients and steady-states in R} %\VignetteKeywords{nonlinear equations, differential equations, roots, Jacobian, Hessian} %\VignettePackage{rootSolve} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin of the document \begin{document} \SweaveOpts{engine=R,eps=FALSE} \SweaveOpts{keep.source=TRUE} <>= library("rootSolve") options(prompt = "> ") options(width=70) @ \maketitle The root of a function f(x) is the value of x for which f(x) = 0. Package \rs deals with finding the roots of \textbf{n} nonlinear (or linear) equations. \\This is, it finds the values \[ x_i ^* \quad (i=1,n) \] for which \[f_j ({\mathbf{x}}^* ) = \mathbf{0} \quad(j=1,n) \] Package \rs serves several purposes: \begin{itemize} \item it extends the root finding capabilities of \R for non-linear functions. \item it includes functions for finding the steady-state of systems of ordinary differential equations (ODE) and partial differential equations (PDE). \item it includes functions to numerically estimate gradient matrices. \end{itemize} The package was created to solve the steady-state and stability analysis examples in the book of \cite{Soetaert08}. Please cite this work if you use the package. The various functions in \rs are given in table (1). \section{Finding roots of nonlinear equations in R and rootSolve} The root-finding functions in \R are: \begin{itemize} \item \emph{uniroot}. Finds one root of one equation. \item \emph{polyroot}. Finds the complex roots of a polynomial. \end{itemize} \subsection{One equation} To find the root of function: \[ f(x)=\cos^3 (2x) \] in the interval [0,8] and plot the curve, we write: <>= fun <- function (x) cos(2*x)^3 curve(fun(x), 0, 8) abline(h = 0, lty = 3) uni <- uniroot(fun, c(0, 8))$root points(uni, 0, pch = 16, cex = 2) @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Root found with \emph{uniroot}} \label{fig:one} \end{figure} Although the graph (figure \ref{fig:one}) clearly demonstrates the existence of many roots in the interval [0,8] \R function \emph{uniroot} extracts only one. \rs function \emph{uniroot.all} is a simple extension of \emph{uniroot} which extracts many (presumably *all*) roots in the interval. <>= curve(fun(x), 0, 8) abline(h = 0, lty = 3) All <- uniroot.all(fun, c(0, 8)) points(All, y = rep(0, length(All)), pch = 16, cex = 2) @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Roots found with uniroot.all} \label{fig:two} \end{figure} \emph{uniroot.all} does that by first subdividing the interval into small sections and, for all sections where the function value changes sign, invoking \emph{uniroot} to locate the root. Note that this is not a full-proof method: in case subdivision is not fine enough some roots will be missed. Also, in case the curve does not cross the X-axis, but just "touches" it, the root will not be retrieved; (but neither will it be located by \emph{uniroot}). \subsection{n equations in n unknowns} Except for polynomial root finding, to date \R has no functions that retrieve the roots of multiple nonlinear equations. Function \emph{multiroot} in \rs implements the Newton-Raphson method (e.g. \cite{Press92}) to solve this type of problem. As the Newton-Raphson method locates the root iteratively, the result will depend on the initial guess of the root. Also, it is not guaranteed that the root will actually be found (i.e. the method may fail). The example below finds two different roots of a three-valued function: \begin{eqnarray*} f_1=x_1+x_2+x_3^2-12\\ f_2=x_1^2-x_2+x_3-2\\ f_3=2 \cdot x_1-x_2^2+x_3-1 \end{eqnarray*} %%\begin{verbatim} <<>>= model <- function(x) { F1 <- x[1] + x[2] + x[3]^2 -12 F2 <- x[1]^2 - x[2] + x[3] -2 F3 <- 2*x[1] - x[2]^2 + x[3] -1 c(F1 = F1, F2 = F2, F3 = F3) } # first solution (ss <- multiroot(f = model, start = c(1, 1, 1))) # second solution; use different start values (ss2 <- multiroot(model, c(0, 0, 0)))$root model(ss2$root) # the function value at the root @ As another example, we seek the 5x5 matrix \textbf{X} for which \[ \mathbf{X \cdot X \cdot X} = \left[ \begin{array}{*{20}l} 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 \\ \end{array} \right] \] <<>>= f2<-function(x) { X <- matrix(nr = 5, x) X %*% X %*% X - matrix(nrow = 5, data = 1:25, byrow = TRUE) } print(system.time( x<-multiroot(f2, start = 1:25)$root )) (X<-matrix(nrow = 5, x)) X%*%X%*%X @ \subsection{n equations in n unknowns with known Jacobian} If the Jacobian is known, OR it has a known sparsity structure, then it is much more efficient to take that into account; As an example, a set of linear equations, comprising 500 unknowns are solved. Of course, one would not do that using a nonlinear equation solver, but rather by using \code{solve}. However, the example is included here to show how to make good use of the flexibility at which the Jacobian can be specified, and to show that the methods are quite fast! We start by defining the linear system of equations: a (500x500) matrix \code{A} and a vector \code{B} of length 500. It is solved with \code{solve}, \R's default solver for systems of linear equations. The time it takes to do this is printed (in seconds): <<>>= A <- matrix(nrow = 500, ncol = 500, runif(500*500)) B <- runif(500) print(system.time(X1 <- solve(A, B))) @ To use \code{multiroot} the nonlinear equation solver, it is noted that the Jacobian is equal to matrix \code{A}. A function is created that returns the Jacobian (or matrix A); the calling sequence of the Jacobian function is \code{function(x)} <<>>= jfun <- function (x) A @ The function whose root is to be solved receives the current estimate of \code{x}, and returns the difference \code{A x-B}: <<>>= fun <- function(x) A %*%x - B @ Next \code{multiroot} is called with jactype = "fullusr" (it is a full Jacobian, and specified by the user): <<>>= print(system.time( X <- multiroot(start = 1:500, f = fun, jactype = "fullusr", jacfunc = jfun) )) @ Finally, both solutions are the same <<>>= sum( (X1 - X$y)^2) @ \section{Steady-state analysis} Ordinary differential equations are a special case of nonlinear equations, where \textbf{x} are called the state variables and the functions \textbf{f(x)} specify the derivatives of x with respect to some independent variable. This is: \[ f(\textbf{x},t) = \frac{{d\textbf{x}}}{{dt}} \] If "t", the independent variable is "time", then the root of the ODE system \[\frac{{d\textbf{x}}}{{dt}}=0\] is often referred to as the "steady-state" condition. Within \R, package \ds \citep{deSolve_jss} is designed to solve so-called initial value problems (IVP) of ODEs and PDEs - partial differential equations by integration. \ds includes integrators that deal efficiently with sparse and banded Jacobians or that are especially designed to solve initial value problems resulting from 1-Dimensional and 2-Dimensional partial differential equations. The latter are first written as ODEs using the method-of-lines approach. To ensure compatibility, \rs offers the same functionalities as \ds, and requires the ODE's to be similarly specified. The function specifying the ordinary differential equations should thus be defined as: \begin{verbatim} deriv = function(x,t,parms,...) \end{verbatim} where \emph{parms} are the ODE parameters, \emph{x} the state variables, \emph{t }the independent variable and \emph{...} are any other arguments passed to the function (optional). The return value of the function should be a list, whose first element is a vector containing the derivatives of x with respect to time. Two different approaches are used to solve for the \emph{steady-state} condition of ODE's: \begin{itemize} \item by dynamically running to steady-state. \item by solving for the root of the ODE using the Newton-Raphson method. \end{itemize} \subsection{Running dynamically to steady-state} Function \emph{runsteady} finds the steady-state condition by dynamically running (integrating) the ODE until the derivatives stop changing. This solves a particular case of an IVP, where the time instance for which the value of the state variable is sought equals infinity. The implementation is based on \ds solver function \emph{lsode} (\cite{Hindmarsh83}). Consider the following simple sediment biogeochemical model: \begin{eqnarray*} \frac{dOM}{dt}&=&Flux-r \cdot OM \cdot \frac{O_2}{O_2+ks} -r \cdot OM \cdot (1-\frac{O_2}{O_2+ks}) \cdot \frac{SO_4}{SO_4+ks2}\\ \frac{dO_2}{dt}&=& -r \cdot OM \cdot \frac{O_2}{O_2+ks} -2 rox \cdot HS \cdot \frac{O_2}{O_2+ks} +D\cdot (BO2-O_2)\\ \frac{dSO_4}{dt}&=& -0.5\cdot r \cdot OM \cdot (1-\frac{O_2}{O_2+ks}) \cdot \frac{SO_4}{SO_4+ks2}+rox \cdot HS \cdot \frac{O_2}{O_2+ks} +D\cdot (BSO4-SO_4)\\ \frac{dHS}{dt}&=& 0.5\cdot r \cdot OM \cdot (1-\frac{O_2}{O_2+ks}) \cdot \frac{SO_4}{SO_4+ks2}-rox \cdot HS \cdot \frac{O_2}{O_2+ks} +D\cdot (BHS-HS) \end{eqnarray*} In \R this model is defined as: <<>>= model <- function(t, y, pars) { with (as.list(c(y, pars)),{ oxicmin = r*OM*(O2/(O2+ks)) anoxicmin = r*OM*(1-O2/(O2+ks))* SO4/(SO4+ks2) dOM = Flux - oxicmin - anoxicmin dO2 = -oxicmin -2*rox*HS*(O2/(O2+ks)) + D*(BO2-O2) dSO4 = -0.5*anoxicmin +rox*HS*(O2/(O2+ks)) + D*(BSO4-SO4) dHS = 0.5*anoxicmin -rox*HS*(O2/(O2+ks)) + D*(BHS-HS) list(c(dOM, dO2, dSO4, dHS), SumS = SO4+HS) }) } @ After defining the value of the parameters (\code{pars}) and the initial values (\code{y}), the model can be run to steady-state (\code{runsteady}). We specify the maximal length of time the simulation can take ($1e^5$) <<>>= pars <- c(D = 1, Flux = 100, r = 0.1, rox = 1, ks = 1, ks2 = 1, BO2 = 100, BSO4 = 10000, BHS = 0) y <- c(OM = 1, O2 = 1, SO4 = 1, HS = 1) print(system.time( RS <- runsteady(y = y, fun = model, parms = pars, times = c(0, 1e5)) )) RS @ The output shows the steady-state concentrations (\code{ST$y}), and the output variable (\code{ST$SumS}), the time and the number of timesteps it took to reach steady-state (attribute "time", "steps"). \subsection{Using the Newton-Raphson method} Functions \code{stode}, \code{stodes}, \code{steady}, \code{steady.1D}, \code{steady.2D}, \code{steady.3D}, and \code{steady.band} find the steady-state by the iterative Newton-Raphson method (e.g. \cite{Press92}. The same model as above can also be solved using \code{stode}. This is faster than running dynamically to steady-state, but not all models can be solved this way <<>>= stode(y = y, fun = model, parms = pars, pos = TRUE) @ Note that we set \code{pos=TRUE} to ensure that only positive values are found. Thus the outcome will be biologically realistic (negative concentrations do not exist). \subsection{Steady-state of 1-D models} Two special-purpose functions solve for the steady-state of 1-D models. \begin{itemize} \item Function \code{steady.band} efficiently estimates the steady-state condition for 1-D models that comprise one species only. \item Function \code{steady-1D} finds the steady-state for multi-species 1-D problems. \end{itemize} \subsubsection{1-D models of 1 species} Consider the following 2nd order differential equation whose steady-state should be estimated: \[ \frac{ \partial {y}}{\partial {t}}= 0 = \frac{\partial ^2{y}}{\partial {dx^2}} + \frac{1}{x}\frac{\partial {y}}{\partial {x}}+ (1-\frac{1}{4\cdot x ^2})\cdot y - \sqrt(x) \cdot\ cos(x) \] over the interval [1,6] and with boundary conditions: $y(1) = 1$ and $y(6) = -0.5$ The spatial derivatives are approximated using centred differences \footnote{in a later section, an alternative approximation is used}: \[ \frac{\partial ^2{y}}{\partial {x^2}}\approx \frac{y_{i+1}-2 \cdot y_i + y_{i-1}}{\Delta x ^2} \] and \[ \frac{\partial {y}}{\partial {x}}\approx \frac{y_{i+1}-y_{i-1}}{2 \cdot \Delta x} \] First the model function is defined: <<>>= derivs <- function(t, y, parms, x, dx, N, y1, y6) { d2y <- (c(y[-1],y6) -2*y + c(y1,y[-N])) /dx/dx dy <- (c(y[-1],y6) - c(y1,y[-N])) /2/dx res <- d2y + dy/x + (1-1/(4*x*x))*y-sqrt(x)*cos(x) return(list(res)) } @ Then the interval [1,6] is subdivided in 5001 boxes (\code{x}) and the steady-state condition estimated, using \code{steady.band}; we specify that there is only one species (\code{nspec=1}). <<>>= dx <- 0.001 x <- seq(1, 6, by = dx) N <- length(x) print(system.time( y <- steady.band(y = rep(1, N), time = 0, func = derivs, x = x, dx = dx, N = N, y1 = 1, y6 = -0.5, nspec = 1)$y )) @ The steady-state of this system of 5001 nonlinear equations is retrieved in about 0.03 seconds \footnote{on my computer that dates from 2008}. The analytical solution of this equation is known; after plotting the numerical approximation, it is added to the figure (figure \ref{fig:two}): <>= plot(x, y, type = "l", main = "5001 nonlinear equations - banded Jacobian") curve(0.0588713*cos(x)/sqrt(x)+1/4*sqrt(x)*cos(x)+ 0.740071*sin(x)/sqrt(x)+1/4*x^(3/2)*sin(x),add=TRUE,type="p") legend("topright", pch = c(NA, 1), lty = c(1, NA), c("numeric", "analytic")) @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the 2nd order differential equation - see text for explanation} \label{fig:two} \end{figure} \subsubsection{1-D models of many species} In the following model, dynamics of BOD (biochemical oxygen demand) and oxygen is modeled in a river. Both are transported downstream (velocity \code{v}) \begin{eqnarray*} \frac{\partial BOD}{\partial t}&=&0=-\cdot \frac{\partial v \cdot BOD }{\partial x} - r \cdot BOD \cdot \frac{O_2}{O_2+ks}\\ \frac{\partial O_2}{\partial t}&=&0=-\cdot \frac{\partial v \cdot O_2 }{\partial x} - r \cdot BOD \cdot \frac{O_2}{O_2+ks} +p \cdot (O_2sat-O_2) \end{eqnarray*} subject to the boundary conditions $BOD(x=0)=BOD_0$ and $O_2(x=0)=O2_0$ First the advective fluxes (transport with velocity v) are calculated, taking into account the upstream concentrations (\code{FluxBOD}, \code{FluxO2}); then the rate of change is written as the sum of -1*Flux gradient and the consumption and production rate: <<>>= O2BOD <- function(t, state, pars) { BOD <- state[1:N] O2 <- state[(N+1):(2*N)] 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), BODrate = BODrate)) } @ After assigning values to the parameters, and setting up the computational grid (\code{x}), steady-state is estimated with function \code{steady.1D}; there are 2 species (BOD, O2) (\code{nspec=2}); we force the result to be positive (\code{pos=TRUE}). <<>>= dx <- 10 # grid size, meters v <- 1e2 # velocity, m/day 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 x <- seq(dx/2, 10000, by = dx) # m, distance from river N <- length(x) state <- c(rep(200, N), rep(200, N)) # initial guess of state variables: print(system.time( out <- steady.1D (y = state, func = O2BOD, parms = NULL, nspec = 2, pos = TRUE) )) @ Although this model consists of 2000 nonlinear equations, it takes only 0.09 seconds to solve it \footnote{on my computer that dates from 2008}. Finally the results are plotted (figure \ref{fig:bod}), using \code{rootSolve}s plot functions: <>= mf <- par(mfrow = c(2, 2)) plot(out, grid = x, xlab = "Distance from river", mfrow = NULL, ylab = "mmol/m3", main = c("Oxygen", "BOD"), type = "l") plot(out, which = "BODrate", grid = x, mfrow = NULL, xlab = "Distance from river", ylab = "mmol/m3/d", main = "BOD decay rate", type = "l") par(mfrow = mf) @ \setkeys{Gin}{width=0.8\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Steady-state solution of the BOD-$O_2$ model. See text for explanation} \label{fig:bod} \end{figure} Note: 1-D problems can also be run dynamically to steady-state. For some models this is the only way. See the help file of \code{steady.1D} for an example. \subsection{Steady-state solution of 2-D PDEs} Function \code{steady.2D} efficiently finds the steady-state of 2-dimensional problems. In the following model \[ \frac{\partial C}{\partial t} = D_x \cdot \frac{\partial ^2 C}{\partial x^2} + D_y \cdot \frac{\partial ^2 C}{\partial y^2} -r \cdot C^2 +p_{xy} \] a substance C is consumed at a quadratic rate ($r \cdot C^2$), while dispersing in X- and Y-direction. At certain positions (x,y) the substance is produced (rate \code{p}). The model is solved on a square (100*100) grid. There are zero-flux boundary conditions at the 4 boundaries. The term $D_x \cdot \frac{\partial ^2 C}{\partial x^2}$ is in fact shorthand for: \[ - \frac{\partial Flux}{\partial x} \] where \[ Flux = -D_x \cdot \frac{\partial C}{\partial x} \] i.e. it is the negative of the flux gradient, where the flux is due to diffusion. In the numerical approximation fo the flux, the concentration gradient is approximated as the subtraction of two matrices, with the columns or rows shifted (e.g. \code{Conc[2:n,]-Conc[1:(n-1),]}). The flux gradient is then also approximated by subtracting entire matrices (e.g. \code{Flux[2:(n+1),]-Flux[1:(n),]}). This is very fast. The zero-flux at the boundaries is imposed by binding a column or row with 0-s. <<>>= diffusion2D <- function(t, conc, par) { Conc <- matrix(nrow = n, ncol = n, data = conc) # vector to 2-D matrix dConc <- -r*Conc*Conc # consumption BND <- rep(1, n) # boundary concentration # constant production in certain cells dConc[ii]<- dConc[ii]+ p #diffusion in X-direction; boundaries=imposed concentration Flux <- -Dx * rbind(rep(0, n), (Conc[2:n,]-Conc[1:(n-1),]), rep(0, n) )/dx dConc <- dConc - (Flux[2:(n+1),] - Flux[1:n,])/dx #diffusion in Y-direction Flux <- -Dy * cbind(rep(0, n), (Conc[,2:n]-Conc[,1:(n-1)]), rep(0, n))/dy dConc <- dConc - (Flux[,2:(n+1)]-Flux[,1:n])/dy return(list(as.vector(dConc))) } @ After specifying the values of the parameters, 10 cells on the 2-D grid where there will be substance produced are randomly selected (\code{ii}). <<>>= # parameters dy <- dx <- 1 # grid size Dy <- Dx <- 1.5 # diffusion coeff, X- and Y-direction r <- 0.01 # 2-nd-order consumption rate (/time) p <- 20 # 0-th order production rate (CONC/t) n <- 100 # 10 random cells where substance is produced at rate p ii <- trunc(cbind(runif(10)*n+1, runif(10)*n+1)) @ The steady-state is found using function \code{steady.2D}. It takes as arguments a.o. the dimensionality of the problem (\code{dimens}) and \code{lrw=1000000}, the length of the work array needed by the solver. If this value is set too small, the solver will return with the size needed. It takes about 0.5 second to solve this 10000 state variable model. <<>>= Conc0 <- matrix(nrow = n, ncol = n, 10.) print(system.time( ST3 <- steady.2D(Conc0, func = diffusion2D, parms = NULL, pos = TRUE, dimens = c(n, n), lrw = 1000000, atol = 1e-10, rtol = 1e-10, ctol = 1e-10) )) @ The S3 image method is used to generate the steady-state plot. <>= image(ST3, main = "2-D diffusion+production", xlab = "x", ylab = "y", legend = TRUE) @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Steady-state solution of the nonlinear 2-Dimensional model} \label{fig:twoD} \end{figure} \subsection{Steady-state solution of 3-D PDEs} Function \code{steady.3D} estimates the steady-state of 3-dimensional problems. We repeat the example from its help file, which models diffusion in 3-D, and with imposed boundary values. <<>>= diffusion3D <- function(t, Y, par) { yy <- array(dim=c(n,n,n),data=Y) # vector to 3-D array dY <- -r*yy # consumption BND <- rep(1,n) # boundary concentration for (i in 1:n) { y <- yy[i,,] #diffusion in X-direction; boundaries=imposed concentration Flux <- -Dy * rbind(y[1,]-BND,(y[2:n,]-y[1:(n-1),]),BND-y[n,])/dy dY[i,,] <- dY[i,,] - (Flux[2:(n+1),]-Flux[1:n,])/dy #diffusion in Y-direction Flux <- -Dz * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dz dY[i,,] <- dY[i,,] - (Flux[,2:(n+1)]-Flux[,1:n])/dz } for (j in 1:n) { y <- yy[,j,] #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[,j,] <- dY[,j,] - (Flux[2:(n+1),]-Flux[1:n,])/dx #diffusion in Y-direction Flux <- -Dz * cbind(y[,1]-BND,(y[,2:n]-y[,1:(n-1)]),BND-y[,n])/dz dY[,j,] <- dY[,j,] - (Flux[,2:(n+1)]-Flux[,1:n])/dz } for (k in 1:n) { y <- yy[,,k] #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[,,k] <- dY[,,k] - (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[,,k] <- dY[,,k] - (Flux[,2:(n+1)]-Flux[,1:n])/dy } return(list(as.vector(dY))) } # parameters dy <- dx <- dz <-1 # grid size Dy <- Dx <- Dz <-1 # diffusion coeff, X- and Y-direction r <- 0.025 # consumption rate n <- 10 y <- array(dim=c(n, n, n), data = 10.) print(system.time( ST3 <- steady.3D(y, func =diffusion3D, parms = NULL, pos = TRUE, dimens = c(n, n, n), lrw = 100000, atol = 1e-10, rtol = 1e-10, ctol = 1e-10, verbose = TRUE) )) @ The S3 image method is used to generate the steady-state plot. We can loop over either one of the dimensions. Here we loop over the first dimenstion, selecting a subset of the 10 cells (\code{dimselect}). We add contours and a legend. <>= image(ST3, mfrow = c(2, 2), add.contour = TRUE, legend = TRUE, dimselect = list(x = c(1, 4, 8, 10))) @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Steady-state solution of the 3-Dimensional model} \label{fig:threeD} \end{figure} \section{Boundary value problems} The previous example from functions \code{steady.1D, steady.2D and steady.3D} solved bundary value problems, in a way that the function call is compatible with initial value problem solvers from package \ds. Function \code{multiroot.1D} can also be used to solve boundary value problems of ordinary differential equations, but has a simpler function interface. It also uses the method-of-lines approach. Another package, \bs provides two totally different methods to solve boundary values problems \citep{bvpSolve}. The following differential equation: \[ 0=f(x,y,\frac{dy}{dx},\frac{d^2y}{dx^2}) \] with boundary conditions $y_{x=a} = ya$, at the start and $y_{x=b}=yb$ at the end of the integration interval [a,b] is solved as follows: \begin{enumerate} \item First the integration interval x is discretized, \begin{verbatim} dx <- 0.01 x <- seq(a,b,by=dx) \end{verbatim} where \code{dx} should be small enough, such as to keep the numerical discretisation error reasonable. \item Then the first- and second-order derivatives are differenced on this numerical grid. R's \code{diff} function is very efficient in taking numerical differences, so it is used to approximate the first-, and second-order derivates as follows. A \emph{first-order derivative y'} can be approximated either as: \begin{itemize} \item y'=\code{diff(c(ya,y))/dx} if only the initial condition ya is prescribed, \item y'=\code{diff(c(y,yb))/dx} if only the final condition, yb is prescribed, \item y'=\code{0.5*(diff(c(ya,y))/dx+diff(c(y,yb))/dx)} if initial, ya, and final condition, yb are prescribed. \end{itemize} The latter (centered differences) is to be preferred. A \emph{second-order derivative y''} can be approximated by differencing twice. y''=\code{diff(diff(c(ya,y,yb))/dx)/dx} \item Finally, function \code{multiroot.1D} is used to locate the root. \end{enumerate} \subsection{test problem 22} As an example, the following boundary value problem will be solved: \begin{eqnarray*} \xi y'' + y' + y^2 = 0 \\ y_{x=0} = 0 \\ y_{x=1} = 1/2 \end{eqnarray*} This is problem number 22 from a set of test boundary value problems which can be found at: \url{http://www.ma.ic.ac.uk/~jcash/BVP_software/PROBLEMS.PDF}. First the function whose root has to solved is implemented: <<>>= bvp22 <- function (y, xi) { dy2 <- diff(diff(c(ya, y, yb))/dx)/dx dy <- 0.5*(diff(c(ya, y))/dx + diff(c(y, yb))/dx) return(xi*dy2+dy+y^2) } @ Then the grid [0,1] is discretised (\code{x}) and the boundary values (\code{ya,yb}) defined <<>>= dx <- 0.001 x <- seq(0, 1, by = dx) N <- length(x) ya <- 0 yb <- 0.5 @ The model is solved for different values of $\xi$ and the output plotted. With the settings of dx, the root of 1001 equations needs to be found; the time it takes (in \emph{mili}seconds) is printed for the first application. <<>>= print(system.time( Y1<- multiroot.1D(f = bvp22, start = runif(N), nspec = 1, xi = 0.1) )*1000) Y2<- multiroot.1D(f = bvp22, start = runif(N), nspec = 1, xi = 0.05) print(system.time( Y3<- multiroot.1D(f = bvp22, start = runif(N), nspec = 1, xi = 0.01) )*1000) @ <>= plot(x, Y3$root, type = "l", col = "green", lwd = 2, main = "bvp test problem 22" , ylab = "y") lines(x, Y2$root, col = "red", lwd = 2) lines(x, Y1$root, col = "blue", lwd = 2) @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the boundary value problem, for three values of $\xi$} \label{fig:bnd} \end{figure} \section{writing functions in compiled code} Similarly as for the models that are solved with integration routines from package \ds, the models solved by the steady-state routines (\code{stode, stodes, steady, steady.1D, steady.2D, steady.3D}) can be written in compiled code (C or Fortran). A vignette ("compiledCode") from package \ds can be consulted for how to do that \citep{compiledCode}. Here the simple sediment biogeochemical model from chapter 2.1 is implemented in \code{C} and \code{Fortran}. \subsection{main function in C-code} For \emph{code written in C}, the calling sequence for \code{func} must be as follows: \begin{verbatim} void anoxmod(int *neq, double *t, double *y, double *ydot, double *yout, int *ip) double OM, O2, SO4, HS; double Min, oxicmin, anoxicmin; if (ip[0] <1) error("nout should be at least 1"); OM = y[0]; O2 = y[1]; SO4 = y[2]; HS = y[3]; Min = r*OM; oxicmin = Min*(O2/(O2+ks)); anoxicmin = Min*(1-O2/(O2+ks))* SO4/(SO4+ks2); ydot[0] = Flux - oxicmin - anoxicmin; ydot[1] = -oxicmin -2*rox*HS*(O2/(O2+ks)) + D*(BO2-O2); ydot[2] = -0.5*anoxicmin +rox*HS*(O2/(O2+ks)) + D*(BSO4-SO4); ydot[3] = 0.5*anoxicmin -rox*HS*(O2/(O2+ks)) + D*(BHS-HS); yout[0] = SO4+HS; \end{verbatim} where \code{*neq} is the number of equations, \code{*t} is the value of the independent variable, \code{y} points to a double precision array of length \code{*neq} that contains the current value of the state variables, and \code{ydot} points to an array that will contain the calculated derivatives. \code{yout} points to a double precision vector whose first \code{nout} values are other output variables (different from the state variables \code{y}), and the next values are double precision values as passed by parameter \code{rpar} when calling the steady-state solver. The key to the elements of \code{yout} is set in \code{*ip}. \code{*ip} points to an integer vector whose length is at least 3; the first element contains the number of output values (which should be equal to \code{nout}), its second element contains the length of \code{*yout}, and the third element contains the length of \code{*ip}; next are integer values, as passed by parameter \code{ipar} when calling the steady-state solver. \subsection{main function in FORTRAN-code} For \emph{code written in Fortran}, the calling sequence for \code{func} must be as in the following example: \begin{verbatim} subroutine model (neq, t, y, ydot, yout, ip) double precision t, y(4), ydot(4), yout(*) double precision OM,O2,SO4,HS double precision min, oxicmin, anoxicmin integer neq, ip(*) double precision D, Flux, r, rox, ks, ks2, BO2, BSO4, BHS common /myparms/D, Flux, r, rox, ks, ks2, BO2, BSO4, BHS IF (ip(1) < 1) call rexit("nout should be at least 1") OM = y(1) O2 = y(2) SO4 = y(3) HS = y(4) Min = r*OM oxicmin = Min*(O2/(O2+ks)) anoxicmin = Min*(1-O2/(O2+ks))* SO4/(SO4+ks2) ydot(1) = Flux - oxicmin - anoxicmin ydot(2) = -oxicmin -2*rox*HS*(O2/(O2+ks)) + D*(BO2-O2) ydot(3) = -0.5*anoxicmin +rox*HS*(O2/(O2+ks)) + D*(BSO4-SO4 ydot(4) = 0.5*anoxicmin -rox*HS*(O2/(O2+ks)) + D*(BHS-HS) yout(1) = SO4+HS return end \end{verbatim} Note that we start by checking whether enough room is allocated for the output variables, else an error is passed to R (\code{rexit}) and the integration is stopped. In this example, parameters are kept in a common block (called \code{myparms}) in the Fortran code \subsection{initialisation subroutine} In order to put parameters in the common block from the calling \R code, an \emph{initialisation subroutine} as specified in \code{initfunc} should be defined. This function has as its sole argument a function \code{steadyparms} that fills a double array with double precision values. In the example here, the initialisation subroutine is called \code{myinit}: \begin{verbatim} subroutine myinit(steadyparms) external steadyparms double precision parms(9) common /myparms/parms call steadyparms(9, parms) return end \end{verbatim} Here \code{myinit} just calls \code{steadyparms} with the dimension of the parameter vector, and the array \code{parms} that will contain the parameter values. The corresponding C-code is: \begin{verbatim} void initanox (void (* steadyparms)(int *, double *)) int N = 9; steadyparms(&N, parms); \end{verbatim} \subsection{jacobian subroutine} If it is desired to supply a Jacobian to the solver, then the Jacobian must be defined in compiled code if the ode system is. The C function call for such a function must be as follows: \begin{verbatim} void myjac(int *neq, double *t, double *y, int *ml, int *mu, double *pd, int *nrowpd, double *yout, int *ip) \end{verbatim} The corresponding subroutine call in Fortran is: \begin{verbatim} subroutine myjac (neq, t, y, ml, mu, pd, nrowpd, yout, ip) integer neq, ml, mu, nrowpd, ip(*) double precision y(*), pd(nrowpd,*), yout(*) \end{verbatim} \subsection{Estimating steady-state for models written in compiled code} To run the model using e.g. the Fortran code, the code in anoxmod.f must first be compiled. This can be done in R itself: \code{system("R CMD SHLIB anoxmod.f")} which will create file anoxmod.dll After loading the DLL, the model can be solved: \begin{verbatim} dyn.load("anoxmod.dll") ST2 <- stode(y = y, func = "model", parms = pars, dllname = "anoxmod", initfunc = "myinit", pos = TRUE, nout = 1) \end{verbatim} Examples in both C and Fortran are in the \code{dynload} subdirectory of the \code{rootSolve} package directory. \section{Gradients, Jacobians and Hessians} \subsection{Gradient and Hessian matrices} Function \code{gradient} returns a forward difference approximation for the derivative of the function \code{f(y,...)} evaluated at the point specified by \code{x}. Function \code{hessian} returns a forward difference approximation of the hessian matrix. In the example below, the root of the "banana function" is first estimated (using R-function \code{nlm}), after which the gradient and the hessian at this point are taken. All this can also be achieved using function \code{nlm}. Note that, as \code{hessian} returns a (forward or centered) difference approximation of the gradient, which itself is also estimated by differencing, it is not very precise. <<>>= # the banana function fun <- function(x) 100*(x[2] - x[1]^2)^2 + (1 - x[1])^2 # the minimum mm <-nlm(fun, p=c(0,0))$estimate # the Hessian (Hes <- hessian(fun,mm)) # the gradient (grad <- gradient(fun,mm,centered=TRUE)) # Hessian and gradient can also be estimated by nlm: nlm(fun, p=c(0,0), hessian=TRUE) @ The inverse of the Hessian gives an estimate of parameter uncertainty <<>>= solve(Hes) @ \subsection{Jacobian matrices} Function \code{jacobian.full} and \code{jacobian.band} returns a forward difference approximation of the jacobian (the gradient matrix, where the function f is the derivative) for full and banded problems. <<>>= mod <- function (t=0,y, parms=NULL,...) { dy1 <- y[1] + 2*y[2] dy2 <- 3*y[1] + 4*y[2] + 5*y[3] dy3 <- 6*y[2] + 7*y[3] + 8*y[4] dy4 <- 9*y[3] +10*y[4] return(as.list(c(dy1, dy2, dy3, dy4))) } jacobian.full(y = c(1, 2, 3, 4), func = mod) jacobian.band(y = c(1, 2, 3, 4), func = mod) @ \clearpage \section{Finally} This vignette was created using Sweave \citep{Leisch02}. It takes 2.2 seconds to "Sweave" the file; this is about the time needed to solve all the examples. The computer on which this is run is an Intel Core (TM)2 Duo CPU T9300 2.5 GHz pentium PC with 3 GB of RAM \begin{table*}[t] \caption{Summary of the functions in package \rs; values in \textbf{bold} are vectors}\label{tb:rs} \centering \begin{tabular}{p{.15\textwidth}p{.15\textwidth}p{.75\textwidth}}\hline \rule[-3mm]{0mm}{8mm} & Function &Description\\ \hline \hline $f(x) = 0$, &&\\ $a