\documentclass[article,nojss]{jss} \DeclareGraphicsExtensions{.pdf, .eps, .png, .jpeg} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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} % table commands \setlength{\extrarowheight}{0.1cm} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands. \newcommand{\R}{\proglang{R }} \newcommand{\ds}{\textbf{\textsf{deSolve }}} \newcommand{\bs}{\textbf{\textsf{bvpSolve }}} \newcommand{\rt}{\textbf{\textsf{ReacTran }}} \newcommand{\rb}[1]{\raisebox{1.5ex}{#1}} \title{Package \pkg{deSolve}: Solving Initial Value Differential Equations in \proglang{R}} \Plaintitle{Package deSolve: Solving Initial Value Differential Equations in R} \Keywords{differential equations, ordinary differential equations, differential algebraic equations, partial differential equations, initial value problems, \proglang{R}} \Plainkeywords{differential equations, ordinary differential equations, differential algebraic equations, partial differential equations, initial value problems, R} \author{ Karline Soetaert\\ Royal Netherlands Institute\\ of Sea Research (NIOZ)\\ Yerseke, The Netherlands \And Thomas Petzoldt\\ Technische Universit\"at \\ Dresden\\ Germany \And R. Woodrow Setzer\\ National Center for\\ Computational Toxicology\\ US Environmental Protection Agency } \Plainauthor{Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer} \Abstract{ \R package \ds \citep{deSolve_jss,deSolve} the successor of \proglang{R} package \pkg{odesolve} is a package to solve initial value problems (IVP) of: \begin{itemize} \item ordinary differential equations (ODE), \item differential algebraic equations (DAE), \item partial differential equations (PDE) and \item delay differential equations (DeDE). \end{itemize} The implementation includes stiff and nonstiff integration routines based on the \pkg{ODEPACK} \proglang{FORTRAN} codes \citep{Hindmarsh83}. It also includes fixed and adaptive time-step explicit Runge-Kutta solvers and the Euler method \citep{Press92}, and the implicit Runge-Kutta method RADAU \citep{Hairer2}. In this vignette we outline how to implement differential equations as \R-functions. Another vignette (``compiledCode'') \citep{compiledCode}, deals with differential equations implemented in lower-level languages such as \proglang{FORTRAN}, \proglang{C}, or \proglang{C++}, which are compiled into a dynamically linked library (DLL) and loaded into \proglang{R} \citep{Rcore}. Note that another package, \bs provides methods to solve boundary value problems \citep{bvpSolve}. } %% The address of (at least) one author should be given %% in the following format: \Address{ Karline Soetaert\\ Centre for Estuarine and Marine Ecology (CEME)\\ Royal Netherlands Institute of Sea Research (NIOZ)\\ 4401 NT Yerseke, Netherlands \\ E-mail: \email{karline.soetaert@nioz.nl}\\ URL: \url{https://www.nioz.nl}\\ \\ Thomas Petzoldt\\ Institut f\"ur Hydrobiologie\\ Technische Universit\"at Dresden\\ 01062 Dresden, Germany\\ E-mail: \email{thomas.petzoldt@tu-dresden.de}\\ URL: \url{https://tu-dresden.de/Members/thomas.petzoldt/}\\ \\ R. Woodrow Setzer\\ National Center for Computational Toxicology\\ US Environmental Protection Agency } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R/Sweave specific LaTeX commands. %% need no \usepackage{Sweave} %\VignetteIndexEntry{R Package deSolve: Solving Initial Value Differential Equations in R} %\VignetteKeywords{differential equations, ordinary differential equations, differential algebraic equations, partial differential equations, initial value problems, R} %\VignettePackage{deSolve} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin of the document \begin{document} \SweaveOpts{engine=R,eps=FALSE} \SweaveOpts{keep.source=TRUE} <>= library("deSolve") options(prompt = "> ") options(width=70) @ \maketitle \section{A simple ODE: chaos in the atmosphere} The Lorenz equations (Lorenz, 1963) were the first chaotic dynamic system to be described. They consist of three differential equations that were assumed to represent idealized behavior of the earth's atmosphere. We use this model to demonstrate how to implement and solve differential equations in \proglang{R}. The Lorenz model describes the dynamics of three state variables, $X$, $Y$ and $Z$. The model equations are: \begin{align*} \frac{dX}{dt} &= a \cdot X + Y \cdot Z \\ \frac{dY}{dt} &= b \cdot (Y - Z) \\ \frac{dZ}{dt} &= - X \cdot Y + c \cdot Y - Z \end{align*} with the initial conditions: \[ X(0) = Y(0) = Z(0) = 1 \] Where $a$, $b$ and $c$ are three parameters, with values of -8/3, -10 and 28 respectively. Implementation of an IVP ODE in \R can be separated in two parts: the model specification and the model application. Model specification consists of: \begin{itemize} \item Defining model parameters and their values, \item Defining model state variables and their initial conditions, \item Implementing the model equations that calculate the rate of change (e.g. $dX/dt$) of the state variables. \end{itemize} The model application consists of: \begin{itemize} \item Specification of the time at which model output is wanted, \item Integration of the model equations (uses R-functions from \pkg{deSolve}), \item Plotting of model results. \end{itemize} Below, we discuss the \proglang{R}-code for the Lorenz model. \subsection{Model specification} \subsubsection{Model parameters} There are three model parameters: $a$, $b$, and $c$ that are defined first. Parameters are stored as a vector with assigned names and values: <<>>= parameters <- c(a = -8/3, b = -10, c = 28) @ \subsubsection{State variables} The three state variables are also created as a vector, and their initial values given: <<>>= state <- c(X = 1, Y = 1, Z = 1) @ \subsubsection{Model equations} The model equations are specified in a function (\code{Lorenz}) that calculates the rate of change of the state variables. Input to the function is the model time (\code{t}, not used here, but required by the calling routine), and the values of the state variables (\code{state}) and the parameters, in that order. This function will be called by the \R routine that solves the differential equations (here we use \code{ode}, see below). The code is most readable if we can address the parameters and state variables by their names. As both parameters and state variables are `vectors', they are converted into a list. The statement \code{with(as.list(c(state, parameters)), {...})} then makes available the names of this list. The main part of the model calculates the rate of change of the state variables. At the end of the function, these rates of change are returned, packed as a list. Note that it is necessary \textbf{to return the rate of change in the same ordering as the specification of the state variables. This is very important.} In this case, as state variables are specified $X$ first, then $Y$ and $Z$, the rates of changes are returned as $dX, dY, dZ$. <<>>= Lorenz<-function(t, state, parameters) { with(as.list(c(state, parameters)),{ # rate of change dX <- a*X + Y*Z dY <- b * (Y-Z) dZ <- -X*Y + c*Y - Z # return the rate of change list(c(dX, dY, dZ)) }) # end with(as.list ... } @ \subsection{Model application} \subsubsection{Time specification} We run the model for 100 days, and give output at 0.01 daily intervals. R's function \code{seq()} creates the time sequence: <<>>= times <- seq(0, 100, by = 0.01) @ \subsubsection{Model integration} The model is solved using \ds function \code{ode}, which is the default integration routine. Function \code{ode} takes as input, a.o. the state variable vector (\code{y}), the times at which output is required (\code{times}), the model function that returns the rate of change (\code{func}) and the parameter vector (\code{parms}). Function \code{ode} returns an object of class \code{deSolve} with a matrix that contains the values of the state variables (columns) at the requested output times. <<>>= library(deSolve) out <- ode(y = state, times = times, func = Lorenz, parms = parameters) head(out) @ \subsubsection{Plotting results} Finally, the model output is plotted. We use the plot method designed for objects of class \code{deSolve}, which will neatly arrange the figures in two rows and two columns; before plotting, the size of the outer upper margin (the third margin) is increased (\code{oma}), such as to allow writing a figure heading (\code{mtext}). First all model variables are plotted versus \code{time}, and finally \code{Z} versus \code{X}: <>= par(oma = c(0, 0, 3, 0)) plot(out, xlab = "time", ylab = "-") plot(out[, "X"], out[, "Z"], pch = ".") mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5) @ \setkeys{Gin}{width=0.8\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the ordinary differential equation - see text for R-code} \label{fig:dae} \end{figure} \clearpage \section{Solvers for initial value problems of ordinary differential equations} Package \ds contains several IVP ordinary differential equation solvers, that belong to the most important classes of solvers. Most functions are based on original (\proglang{FORTRAN}) implementations, e.g. the Backward Differentiation Formulae and Adams methods from \pkg{ODEPACK} \citep{Hindmarsh83}, or from \citep{Brown89,Petzold1983}, the implicit Runge-Kutta method RADAU \citep{Hairer2}. The package contains also a de novo implementation of several Runge-Kutta methods \citep{Butcher1987, Press92, Hairer1}. All integration methods\footnote{except \code{zvode}, the solver used for systems containing complex numbers.} can be triggered from function \code{ode}, by setting \code{ode}'s argument \code{method}), or can be run as stand-alone functions. Moreover, for each integration routine, several options are available to optimise performance. For instance, the next statements will use integration method \code{radau} to solve the model, and set the tolerances to a higher value than the default. Both statements are the same: <<>>= outb <- radau(state, times, Lorenz, parameters, atol = 1e-4, rtol = 1e-4) outc <- ode(state, times, Lorenz, parameters, method = "radau", atol = 1e-4, rtol = 1e-4) @ The default integration method, based on the \proglang{FORTRAN} code LSODA is one that switches automatically between stiff and non-stiff systems \citep{Petzold1983}. This is a very robust method, but not necessarily the most efficient solver for one particular problem. See \citep{deSolve_jss} for more information about when to use which solver in \pkg{deSolve}. For most cases, the default solver, \code{ode} and using the default settings will do. Table \ref{tb:rs} also gives a short overview of the available methods. To show how to trigger the various methods, we solve the model with several integration routines, each time printing the time it took (in seconds) to find the solution: <<>>= print(system.time(out1 <- rk4 (state, times, Lorenz, parameters))) print(system.time(out2 <- lsode (state, times, Lorenz, parameters))) print(system.time(out <- lsoda (state, times, Lorenz, parameters))) print(system.time(out <- lsodes(state, times, Lorenz, parameters))) print(system.time(out <- daspk (state, times, Lorenz, parameters))) print(system.time(out <- vode (state, times, Lorenz, parameters))) @ \subsection{Runge-Kutta methods and Euler} The explicit Runge-Kutta methods are de novo implementations in \proglang{C}, based on the Butcher tables \citep{Butcher1987}. They comprise simple Runge-Kutta formulae (Euler's method \code{euler}, Heun's method \code{rk2}, the classical 4th order Runge-Kutta, \code{rk4}) and several Runge-Kutta pairs of order 3(2) to order 8(7). The embedded, explicit methods are according to \citet{Fehlberg1967} (\code{rk..f}, \code{ode45}), \citet{Dormand1980,Dormand1981} (\code{rk..dp.}), \citet{Bogacki1989} (\code{rk23bs}, \code{ode23}) and \citet{Cash1990} (\code{rk45ck}), where \code{ode23} and \code{ode45} are aliases for the popular methods \code{rk23bs} resp. \code{rk45dp7}. With the following statement all implemented methods are shown: <<>>= rkMethod() @ This list also contains implicit Runge-Kutta's (\code{irk..}), but they are not yet optimally coded. The only well-implemented implicit Runge-Kutta is the \code{radau} method \citep{Hairer2} that will be discussed in the section dealing with differential algebraic equations. The properties of a Runge-Kutta method can be displayed as follows: <<>>= rkMethod("rk23") @ Here \code{varstep} informs whether the method uses a variable time-step; \code{FSAL} whether the first same as last strategy is used, while \code{stage} and \code{Qerr} give the number of function evaluations needed for one step, and the order of the local truncation error. \code{A, b1, b2, c} are the coefficients of the Butcher table. Two formulae (\code{rk45dp7, rk45ck}) support dense output. It is also possible to modify the parameters of a method (be very careful with this) or define and use a new Runge-Kutta method: <<>>= func <- function(t, x, parms) { with(as.list(c(parms, x)),{ dP <- a * P - b * C * P dC <- b * P * C - c * C res <- c(dP, dC) list(res) }) } rKnew <- rkMethod(ID = "midpoint", varstep = FALSE, A = c(0, 1/2), b1 = c(0, 1), c = c(0, 1/2), stage = 2, Qerr = 1 ) out <- ode(y = c(P = 2, C = 1), times = 0:100, func, parms = c(a = 0.1, b = 0.1, c = 0.1), method = rKnew) head(out) @ \subsubsection{Fixed time-step methods} There are two explicit methods that do not adapt the time step: the \code{euler} method and the \code{rk4} method. They are implemented in two ways: \begin{itemize} \item as a \code{rkMethod} of the \textbf{general} \code{rk} solver. In this case the time step used can be specified independently from the \code{times} argument, by setting argument \code{hini}. Function \code{ode} uses this general code. \item as \textbf{special} solver codes \code{euler} and \code{rk4}. These implementations are simplified and with less options to avoid overhead. The timestep used is determined by the time increment in the \code{times} argument. \end{itemize} For example, the next two statements both trigger the Euler method, the first using the ``special'' code with a time step = 1, as imposed by the \code{times} argument, the second using the generalized method with a time step set by \code{hini}. Unsurprisingly, the first solution method completely fails (the time step $= 1$ is much too large for this problem). \begin{verbatim} out <- euler(y = state, times = 0:40, func = Lorenz, parms = parameters) outb <- ode(y = state, times = 0:40, func = Lorenz, parms = parameters, method = "euler", hini = 0.01) \end{verbatim} \subsection{Model diagnostics and summaries} Function \code{diagnostics} prints several diagnostics of the simulation to the screen. For the Runge-Kutta and \code{lsode} routine called above they are: <<>>= diagnostics(out1) diagnostics(out2) @ There is also a \code{summary} method for \code{deSolve} objects. This is especially handy for multi-dimensional problems (see below) <<>>= summary(out1) @ \clearpage \section{Partial differential equations} As package \ds includes integrators that deal efficiently with arbitrarily sparse and banded Jacobians, it is especially well suited to solve initial value problems resulting from 1, 2 or 3-dimensional partial differential equations (PDE), using the method-of-lines approach. The PDEs are first written as ODEs, using finite differences. This can be efficiently done with functions from R-package \rt \citep{ReacTran}. However, here we will create the finite differences in R-code. Several special-purpose solvers are included in \pkg{deSolve}: \begin{itemize} \item \code{ode.band} integrates 1-dimensional problems comprizing one species, \item \code{ode.1D} integrates 1-dimensional problems comprizing one or many species, \item \code{ode.2D} integrates 2-dimensional problems, \item \code{ode.3D} integrates 3-dimensional problems. \end{itemize} As an example, consider the Aphid model described in \citet{Soetaert08}. It is a model where aphids (a pest insect) slowly diffuse and grow on a row of plants. The model equations are: \[ \frac{{\partial N}}{{\partial t}} = - \frac{{\partial Flux}}{{\partial {\kern 1pt} x}} + g \cdot N \] and where the diffusive flux is given by: \[ Flux = - D\frac{{\partial N}}{{\partial {\kern 1pt} x}} \] with boundary conditions \[ N_{x=0}=N_{x=60}=0 \] and initial condition \begin{center} $N_x=0$ for $x \neq 30$ $N_x=1$ for $x = 30$ \end{center} In the method of lines approach, the spatial domain is subdivided in a number of boxes and the equation is discretized as: \[ \frac{{dN_i }}{{dt}} = - \frac{{Flux_{i,i + 1} - Flux_{i - 1,i} }}{{\Delta x_i }} + g \cdot N_i \] with the flux on the interface equal to: \[ Flux_{i - 1,i} = - D_{i - 1,i} \cdot \frac{{N_i - N_{i - 1} }}{{\Delta x_{i - 1,i} }} \] Note that the values of state variables (here densities) are defined in the centre of boxes (i), whereas the fluxes are defined on the box interfaces. We refer to \citet{Soetaert08} for more information about this model and its numerical approximation. Here is its implementation in \proglang{R}. First the model equations are defined: <<>>= Aphid <- function(t, APHIDS, parameters) { deltax <- c (0.5, rep(1, numboxes - 1), 0.5) Flux <- -D * diff(c(0, APHIDS, 0)) / deltax dAPHIDS <- -diff(Flux) / delx + APHIDS * r # the return value list(dAPHIDS ) } # end @ Then the model parameters and spatial grid are defined <<>>= D <- 0.3 # m2/day diffusion rate r <- 0.01 # /day net growth rate delx <- 1 # m thickness of boxes numboxes <- 60 # distance of boxes on plant, m, 1 m intervals Distance <- seq(from = 0.5, by = delx, length.out = numboxes) @ Aphids are initially only present in two central boxes: <<>>= # Initial conditions: # ind/m2 APHIDS <- rep(0, times = numboxes) APHIDS[30:31] <- 1 state <- c(APHIDS = APHIDS) # initialise state variables @ The model is run for 200 days, producing output every day; the time elapsed in seconds to solve this 60 state-variable model is estimated (\code{system.time}): <<>>= times <-seq(0, 200, by = 1) print(system.time( out <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid") )) @ Matrix \code{out} consist of times (1st column) followed by the densities (next columns). <<>>= head(out[,1:5]) @ The \code{summary} method gives the mean, min, max, ... of the entire 1-D variable: <<>>= summary(out) @ Finally, the output is plotted. It is simplest to do this with \pkg{deSolve}'s \proglang{S3}-method \code{image} %% Do this offline %%<>= \begin{verbatim} image(out, method = "filled.contour", grid = Distance, xlab = "time, days", ylab = "Distance on plant, m", main = "Aphid density on a row of plants") \end{verbatim} %%@ \setkeys{Gin}{width=0.8\textwidth} \begin{figure} \begin{center} %%<>= %%<> %%@ \includegraphics{aphid.png} \end{center} \caption{Solution of the 1-dimensional aphid model - see text for \R-code} \label{fig:aphid} \end{figure} As this is a 1-D model, it is best solved with \ds function \code{ode.1D}. A multi-species IVP example can be found in \citet{Soetaert08}. For 2-D and 3-D problems, we refer to the help-files of functions \code{ode.2D} and \code{ode.3D}. The output of one-dimensional models can also be plotted using S3-method \code{plot.1D} and \code{matplot.1D}. In both cases, we can simply take a \code{subset} of the output, and add observations. <<>>= data <- cbind(dist = c(0,10, 20, 30, 40, 50, 60), Aphid = c(0,0.1,0.25,0.5,0.25,0.1,0)) @ <>= par (mfrow = c(1,2)) matplot.1D(out, grid = Distance, type = "l", mfrow = NULL, subset = time %in% seq(0, 200, by = 10), obs = data, obspar = list(pch = 18, cex = 2, col="red")) plot.1D(out, grid = Distance, type = "l", mfrow = NULL, subset = time == 100, obs = data, obspar = list(pch = 18, cex = 2, col="red")) @ \setkeys{Gin}{width=1.0\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the Aphid model - plotted with matplot.1D, plot.1D - see text for R-code} \label{fig:matplot1d} \end{figure} \clearpage \section{Differential algebraic equations} Package \ds contains two functions that solve initial value problems of differential algebraic equations. They are: \begin{itemize} \item \code{radau} which implements the implicit Runge-Kutta RADAU5 \citep{Hairer2}, \item \code{daspk}, based on the backward differentiation code DASPK \citep{Brenan96}. \end{itemize} Function \code{radau} needs the input in the form $M y' = f(t,y,y')$ where $M$ is the mass matrix. Function \code{daspk} also supports this input, but can also solve problems written in the form $F(t, y, y') = 0$. \code{radau} solves problems up to index 3; \code{daspk} solves problems of index $\leq$ 1. \subsection{DAEs of index maximal 1} Function \code{daspk} from package \ds solves (relatively simple) DAEs of index\footnote{note that many -- apparently simple -- DAEs are higher-index DAEs} maximal 1. The DAE has to be specified by the \emph{residual function} instead of the rates of change (as in ODE). Consider the following simple DAE: \begin{eqnarray*} \frac{dy_1}{dt}&=&-y_1+y_2\\ y_1 \cdot y_2 &=& t \end{eqnarray*} where the first equation is a differential, the second an algebraic equation. To solve it, it is first rewritten as residual functions: \begin{eqnarray*} 0&=&\frac{dy_1}{dt}+y_1-y_2\\ 0&=&y_1 \cdot y_2 - t \end{eqnarray*} In \R we write: <<>>= daefun <- function(t, y, dy, parameters) { res1 <- dy[1] + y[1] - y[2] res2 <- y[2] * y[1] - t list(c(res1, res2)) } library(deSolve) yini <- c(1, 0) dyini <- c(1, 0) times <- seq(0, 10, 0.1) ## solver system.time(out <- daspk(y = yini, dy = dyini, times = times, res = daefun, parms = 0)) @ <>= matplot(out[,1], out[,2:3], type = "l", lwd = 2, main = "dae", xlab = "time", ylab = "y") @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the differential algebraic equation model - see text for R-code} \label{fig:dae2} \end{figure} \subsection{DAEs of index up to three} Function \code{radau} from package \ds can solve DAEs of index up to three provided that they can be written in the form $M dy/dt = f(t,y)$. Consider the well-known pendulum equation: \begin{eqnarray*} x' &=& u\\ y' &=& v\\ u' &=& -\lambda x\\ v' &=& -\lambda y - 9.8\\ 0 &=& x^2 + y^2 - 1 \end{eqnarray*} where the dependent variables are $x, y, u, v$ and $\lambda$. Implemented in \R to be used with function \code{radau} this becomes: <<>>= pendulum <- function (t, Y, parms) { with (as.list(Y), list(c(u, v, -lam * x, -lam * y - 9.8, x^2 + y^2 -1 )) ) } @ A consistent set of initial conditions are: <<>>= yini <- c(x = 1, y = 0, u = 0, v = 1, lam = 1) @ and the mass matrix $M$: <<>>= M <- diag(nrow = 5) M[5, 5] <- 0 M @ Function \code{radau} requires that the index of each equation is specified; there are 2 equations of index 1, two of index 2, one of index 3: <<>>= index <- c(2, 2, 1) times <- seq(from = 0, to = 10, by = 0.01) out <- radau (y = yini, func = pendulum, parms = NULL, times = times, mass = M, nind = index) @ <>= plot(out, type = "l", lwd = 2) plot(out[, c("x", "y")], type = "l", lwd = 2) @ \setkeys{Gin}{width=0.8\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Solution of the pendulum problem, an index 3 differential algebraic equation using \code{radau} - see text for \proglang{R}-code} \label{fig:pendulum} \end{figure} \clearpage \section{Integrating systems containing complex numbers, function zvode} Function \code{zvode} solves ODEs that are composed of complex variables. We use \code{zvode} to solve the following system of 2 ODEs: \begin{align*} \frac{dz}{dt} &= i \cdot z\\ \frac{dw}{dt} &= -i \cdot w \cdot w \cdot z\\ \intertext{where} w(0) &= 1/2.1 \\ z(0) &= 1 \end{align*} on the interval $t = [0, 2 \pi]$ <<>>= ZODE2 <- function(Time, State, Pars) { with(as.list(State), { df <- 1i * f dg <- -1i * g * g * f return(list(c(df, dg))) }) } yini <- c(f = 1+0i, g = 1/2.1+0i) times <- seq(0, 2 * pi, length = 100) out <- zvode(func = ZODE2, y = yini, parms = NULL, times = times, atol = 1e-10, rtol = 1e-10) @ The analytical solution is: \begin{align*} f(t) &= \exp (1i \cdot t) \intertext{and} g(t) &= 1/(f(t) + 1.1) \end{align*} The numerical solution, as produced by \code{zvode} matches the analytical solution: <<>>= analytical <- cbind(f = exp(1i*times), g = 1/(exp(1i*times)+1.1)) tail(cbind(out[,2], analytical[,1])) @ \clearpage \section{Making good use of the integration options} The solvers from \pkg{ODEPACK} can be fine-tuned if it is known whether the problem is stiff or non-stiff, or if the structure of the Jacobian is sparse. We repeat the example from \code{lsode} to show how we can make good use of these options. The model describes the time evolution of 5 state variables: <<>>= f1 <- function (t, y, parms) { ydot <- vector(len = 5) ydot[1] <- 0.1*y[1] -0.2*y[2] ydot[2] <- -0.3*y[1] +0.1*y[2] -0.2*y[3] ydot[3] <- -0.3*y[2] +0.1*y[3] -0.2*y[4] ydot[4] <- -0.3*y[3] +0.1*y[4] -0.2*y[5] ydot[5] <- -0.3*y[4] +0.1*y[5] return(list(ydot)) } @ and the initial conditions and output times are: <<>>= yini <- 1:5 times <- 1:20 @ The default solution, using \code{lsode} assumes that the model is stiff, and the integrator generates the Jacobian, which is assummed to be \emph{full}: <<>>= out <- lsode(yini, times, f1, parms = 0, jactype = "fullint") @ It is possible for the user to provide the Jacobian. Especially for large problems this can result in substantial time savings. In a first case, the Jacobian is written as a full matrix: <<>>= fulljac <- function (t, y, parms) { jac <- matrix(nrow = 5, ncol = 5, byrow = TRUE, data = c(0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1, -0.2, 0 , 0 , 0 , -0.3, 0.1)) return(jac) } @ and the model solved as: <<>>= out2 <- lsode(yini, times, f1, parms = 0, jactype = "fullusr", jacfunc = fulljac) @ The Jacobian matrix is banded, with one nonzero band above (up) and one below(down) the diagonal. First we let \code{lsode} estimate the banded Jacobian internally (\code{jactype = "bandint"}): <<>>= out3 <- lsode(yini, times, f1, parms = 0, jactype = "bandint", bandup = 1, banddown = 1) @ It is also possible to provide the nonzero bands of the Jacobian in a function: <<>>= bandjac <- function (t, y, parms) { jac <- matrix(nrow = 3, ncol = 5, byrow = TRUE, data = c( 0 , -0.2, -0.2, -0.2, -0.2, 0.1, 0.1, 0.1, 0.1, 0.1, -0.3, -0.3, -0.3, -0.3, 0)) return(jac) } @ in which case the model is solved as: <<>>= out4 <- lsode(yini, times, f1, parms = 0, jactype = "bandusr", jacfunc = bandjac, bandup = 1, banddown = 1) @ Finally, if the model is specified as ``non-stiff'' (by setting \code{mf=10}), there is no need to specify the Jacobian: <<>>= out5 <- lsode(yini, times, f1, parms = 0, mf = 10) @ \clearpage \section{Events and roots} As from version 1.6, \code{events} are supported. Events occur when the values of state variables are instantaneously changed. They can be specified as a \code{data.frame}, or in a function. Events can also be triggered by a root function. Several integrators (\code{lsoda}, \code{lsodar}, \code{lsode}, \code{lsodes} and \code{radau}) can estimate the root of one or more functions. For the first 4 integration methods, the root finding algorithm is based on the algorithm in solver LSODAR, and implemented in FORTRAN. For \code{radau}, the root solving algorithm is written in C-code, and it works slightly different. Thus, some problems involving roots may be more efficient to solve with either \code{lsoda, lsode}, or \code{lsodes}, while other problems are more efficiently solved with \code{radau}. If a root is found, then the integration will be terminated, unless an event function is defined. A help file with information on roots and events can be opened by typing \code{?events} or \code{?roots}. \subsection{Event specified in a data.frame} In this example, two state variables with constant decay are modeled: <<>>= eventmod <- function(t, var, parms) { list(dvar = -0.1*var) } yini <- c(v1 = 1, v2 = 2) times <- seq(0, 10, by = 0.1) @ At time 1 and 9 a value is added to variable \code{v1}, at time 1 state variable \code{v2} is multiplied with 2, while at time 5 the value of \code{v2} is replaced with 3. These events are specified in a \code{data.frame}, eventdat: <<>>= eventdat <- data.frame(var = c("v1", "v2", "v2", "v1"), time = c(1, 1, 5, 9), value = c(1, 2, 3, 4), method = c("add", "mult", "rep", "add")) eventdat @ The model is solved with \code{ode}: <<>>= out <- ode(func = eventmod, y = yini, times = times, parms = NULL, events = list(data = eventdat)) @ <>= plot(out, type = "l", lwd = 2) @ \setkeys{Gin}{width=0.8\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{A simple model that contains events} \label{fig:event1} \end{figure} \subsection{Event triggered by a root function} This model describes the position (\code{y1}) and velocity (\code{y2}) of a bouncing ball: <<>>= ballode<- function(t, y, parms) { dy1 <- y[2] dy2 <- -9.8 list(c(dy1, dy2)) } @ An event is triggered when the ball hits the ground (height = 0) Then velocity (\code{y2}) is reversed and reduced by 10 percent. The root function, \code{y[1] = 0}, triggers the event: <<>>= root <- function(t, y, parms) y[1] @ The event function imposes the bouncing of the ball <<>>= event <- function(t, y, parms) { y[1]<- 0 y[2]<- -0.9 * y[2] return(y) } @ After specifying the initial values and times, the model is solved, here using \code{lsode}. <<>>= yini <- c(height = 0, v = 20) times <- seq(from = 0, to = 20, by = 0.01) out <- lsode(times = times, y = yini, func = ballode, parms = NULL, events = list(func = event, root = TRUE), rootfun = root) @ <>= plot(out, which = "height", type = "l",lwd = 2, main = "bouncing ball", ylab = "height") @ \setkeys{Gin}{width=0.4\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{A model, with event triggered by a root function} \label{fig:event2} \end{figure} \subsection{Events and time steps} The use of events requires that all event times are contained in the output time steps, otherwise such events would be skipped. This sounds easy but sometimes problems can occur due to the limited accuracy of floating point arithmetics of the computer. To make things work as excpected, two requirements have to be fulfilled: \begin{enumerate} \item all event times have to be contained \textbf{exactly} in times, i.e. with the maximum possible accuracy of floating point arithmetics. \item two time steps should not be too close together, otherwise numerical problems would occur during the integration. \end{enumerate} Starting from version 1.10 of \pkg{deSolve} this is now checked (and if necessary also fixed) automatically by the solver functions. A warning is issued to inform the user about possible problems, especially that the output time steps were now adjusted and therefore different from the ones originally specified by the user. This means that all values of \code{eventtimes} are now contained but only the subset of times that have no exact or ``rather close'' neighbors in \code{eventtimes}. Instead of relying on this automatism, matching times and eventtimes can also be managed by the user, either by appropriate rounding or by using function \code{cleanEventTimes} shown below. Let's assume we have a vector of time steps \code{times} and another vector of event times \code{eventtimes}: <<>>= times <- seq(0, 1, 0.1) eventtimes <- c(0.7, 0.9) @ If we now check whether the \code{eventtimes} are in \code{times}: <<>>= eventtimes %in% times @ we get the surprising answer that this is only partly the case, because \code{seq} made small numerical errors. The easiest method to get rid of this is rounding: <<>>= times2 <- round(times, 1) times - times2 @ The last line shows us that the error was always smaller than, say $10^{-15}$, what is typical for ordinary double precision arithmetics. The accuracy of the machine can be determined with \code{.Machine\$double.eps}. To check if all \code{eventtimes} are now contained in the new times vector \code{times2}, we use: <<>>= eventtimes %in% times2 @ or <<>>= all(eventtimes %in% times2) @ and see that everything is o.k. now. In few cases, rounding may not work properly, for example if a pharmacokinetic model is simulated with a daily time step, but drug injection occurs at precisely fixed times within the day. Then one has to add all additional event times to the ordinary time stepping: <<>>= times <- 1:10 eventtimes <- c(1.3, 3.4, 4, 7.9, 8.5) newtimes <- sort(unique(c(times, eventtimes))) @ If, however, an event and a time step are almost (but not exactly) the same, then it is more safe to use: <<>>= times <- 1:10 eventtimes <- c(1.3, 3.4, 4, 7.9999999999999999, 8.5) newtimes <- sort(c(eventtimes, cleanEventTimes(times, eventtimes))) @ because \code{cleanEventTimes} removes not only the doubled 4 (like \code{unique}, but also the ``almost doubled'' 8, while keeping the exact event time. The tolerance of \code{cleanEventTimes} can be adjusted using an optional argument \code{eps}. As said, this is normally done automatically by the differential equation solvers and in most cases appropriate rounding will be sufficient to get rid of the warnings. \clearpage \section{Delay differential equations} As from \pkg{deSolve} version 1.7, time lags are supported, and a new general solver for delay differential equations, \code{dede} has been added. We implement the lemming model, example 6 from \citep{ST2000}. Function \code{lagvalue} calculates the value of the state variable at \code{t - 0.74}. As long a these lag values are not known, the value 19 is assigned to the state variable. Note that the simulation starts at \code{time = - 0.74}. <<>>= library(deSolve) #----------------------------- # the derivative function #----------------------------- derivs <- function(t, y, parms) { if (t < 0) lag <- 19 else lag <- lagvalue(t - 0.74) dy <- r * y * (1 - lag/m) list(dy, dy = dy) } #----------------------------- # parameters #----------------------------- r <- 3.5; m <- 19 #----------------------------- # initial values and times #----------------------------- yinit <- c(y = 19.001) times <- seq(-0.74, 40, by = 0.01) #----------------------------- # solve the model #----------------------------- yout <- dede(y = yinit, times = times, func = derivs, parms = NULL, atol = 1e-10) @ <>= plot(yout, which = 1, type = "l", lwd = 2, main = "Lemming model", mfrow = c(1,2)) plot(yout[,2], yout[,3], xlab = "y", ylab = "dy", type = "l", lwd = 2) @ \setkeys{Gin}{width=\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{A delay differential equation model} \label{fig:dde} \end{figure} \clearpage \section{Discrete time models, difference equations} There is one special-purpose solver, triggered with \code{method = "iteration"} which can be used in cases where the new values of the state variables are directly estimated by the user, and need not be found by numerical integration. This is for instance useful when the model consists of difference equations, or for 1-D models when transport is implemented by an implicit or a semi-implicit method. We give here an example of a discrete time model, represented by a difference equation: the Teasel model as from \citet[p287]{Soetaert08}. The dynamics of this plant is described by 6 stages and the transition from one stage to another is in a transition matrix: We define the stages and the transition matrix first: <<>>= Stages <- c("DS 1yr", "DS 2yr", "R small", "R medium", "R large", "F") NumStages <- length(Stages) # Population matrix A <- matrix(nrow = NumStages, ncol = NumStages, byrow = TRUE, data = c( 0, 0, 0, 0, 0, 322.38, 0.966, 0, 0, 0, 0, 0 , 0.013, 0.01, 0.125, 0, 0, 3.448 , 0.007, 0, 0.125, 0.238, 0, 30.170, 0.008, 0, 0.038, 0.245, 0.167, 0.862 , 0, 0, 0, 0.023, 0.75, 0 ) ) @ The difference function is defined as usual, but does not return the ``rate of change'' but rather the new relative stage densities are returned. Thus, each time step, the updated values are divided by the summed densities: <<>>= Teasel <- function (t, y, p) { yNew <- A %*% y list (yNew / sum(yNew)) } @ The model is solved using method ``iteration'': <<>>= out <- ode(func = Teasel, y = c(1, rep(0, 5) ), times = 0:50, parms = 0, method = "iteration") @ and plotted using R-function \code{matplot}: <>= matplot(out[,1], out[,-1], main = "Teasel stage distribution", type = "l") legend("topright", legend = Stages, lty = 1:6, col = 1:6) @ \setkeys{Gin}{width=0.6\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{A difference model solved with method = ``iteration''} \label{fig:difference} \end{figure} \section{Plotting deSolve Objects} There are \proglang{S3} \code{plot} and \code{image} methods for plotting 0-D (plot), and 1-D and 2-D model output (image) as generated with \code{ode}, \code{ode.1D}, \code{ode.2D}. How to use it and examples can be found by typing \code{?plot.deSolve}. \subsection{Plotting Multiple Scenario's} The \code{plot} method for \code{deSolve} objects can also be used to compare different scenarios, e.g from the same model but with different sets of parameters or initial values, with one single call to \code{plot}. % Scholarpedia not anymore well maintained; no https available As an example we implement the simple combustion model, which can be found on \texttt{http://www.scholarpedia.org/article/Stiff\_systems}: \[ y' = y^2 \cdot (1-y) \] The model is run with 4 different values of the initial conditions: $y = 0.01, 0.02, 0.03, 0.04$ and written to \code{deSolve} objects \code{out}, \code{out2}, \code{out3}, \code{out4}. <<>>= library(deSolve) combustion <- function (t, y, parms) list(y^2 * (1-y) ) @ <<>>= yini <- 0.01 times <- 0 : 200 @ <<>>= out <- ode(times = times, y = yini, parms = 0, func = combustion) out2 <- ode(times = times, y = yini*2, parms = 0, func = combustion) out3 <- ode(times = times, y = yini*3, parms = 0, func = combustion) out4 <- ode(times = times, y = yini*4, parms = 0, func = combustion) @ The different scenarios are plotted at once, and a suitable legend is written. <>= plot(out, out2, out3, out4, main = "combustion") legend("bottomright", lty = 1:4, col = 1:4, legend = 1:4, title = "yini*i") @ \setkeys{Gin}{width=\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Plotting 4 outputs in one figure} \label{fig:plotdeSolve} \end{figure} \subsection{Plotting Output with Observations} With the help of the optional argument \code{obs} it is possible to specify observed data that should be added to a \code{deSolve} plot. We exemplify this using the \code{ccl4model} in package \code{deSolve}. (see \code{?ccl4model} for what this is about). This model example has been implemented in compiled code. An observed data set is also available, called \code{ccl4data}. It contains toxicant concentrations in a chamber where rats were dosed with CCl4. <<>>= head(ccl4data) @ We select the data from animal ``A'': <<>>= obs <- subset (ccl4data, animal == "A", c(time, ChamberConc)) names(obs) <- c("time", "CP") head(obs) @ After assigning values to the parameters and providing initial conditions, the \code{ccl4model} can be run. We run the model three times, each time with a different value for the first parameter. Output is written to matrices \code{out} \code{out2}, and \code{out3}. <<>>= parms <- c(0.182, 4.0, 4.0, 0.08, 0.04, 0.74, 0.05, 0.15, 0.32, 16.17, 281.48, 13.3, 16.17, 5.487, 153.8, 0.04321671, 0.40272550, 951.46, 0.02, 1.0, 3.80000000) yini <- c(AI = 21, AAM = 0, AT = 0, AF = 0, AL = 0, CLT = 0, AM = 0) out <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = parms) par2 <- parms par2[1] <- 0.1 out2 <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = par2) par3 <- parms par3[1] <- 0.05 out3 <- ccl4model(times = seq(0, 6, by = 0.05), y = yini, parms = par3) @ We plot all these scenarios and the observed data at once: <>= plot(out, out2, out3, which = c("AI", "MASS", "CP"), col = c("black", "red", "green"), lwd = 2, obs = obs, obspar = list(pch = 18, col = "blue", cex = 1.2)) legend("topright", lty = c(1,2,3,NA), pch = c(NA, NA, NA, 18), col = c("black", "red", "green", "blue"), lwd = 2, legend = c("par1", "par2", "par3", "obs")) @ \setkeys{Gin}{width=\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Plotting output and observations in one figure} \label{fig:plotobs} \end{figure} If we do not select specific variables, then only the ones for which there are observed data are plotted. Assume we have measured the total mass at the end of day 6. We put this in a second data set: <<>>= obs2 <- data.frame(time = 6, MASS = 12) obs2 @ then we plot the data together with the three model runs as follows: <>= plot(out, out2, out3, lwd = 2, obs = list(obs, obs2), obspar = list(pch = c(16, 18), col = c("blue", "black"), cex = c(1.2 , 2)) ) @ \setkeys{Gin}{width=1.0\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Plotting variables in common with observations} \label{fig:plotobs2} \end{figure} \subsection{Plotting Summary Histograms} The \code{hist} function plots the histogram for each variable; all plot parameters can be set individually (here for \code{col}). To generate the next plot, we overrule the default \code{mfrow} setting which would plot the figures in 3 rows and 3 columns (and hence plot one figure in isolation) <>= hist(out, col = grey(seq(0, 1, by = 0.1)), mfrow = c(3, 4)) @ \setkeys{Gin}{width=1.0\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Plotting histograms of all output variables} \label{fig:plothist} \end{figure} \subsection{Plotting multi-dimensional output} The \code{image} function plots time versus x images for models solved with \code{ode.1D}, or generates x-y plots for models solved with \code{ode.2D}. \subsubsection{1-D model output} We exemplify its use by means of a Lotka-Volterra model, implemented in 1-D. The model describes a predator and its prey diffusing on a flat surface and in concentric circles. This is a 1-D model, solved in the cylindrical coordinate system. Note that it is simpler to implement this model in R-package \code{ReacTran} \citep{ReacTran}. <>= options(prompt = " ") options(continue = " ") @ We start by defining the derivative function <<>>= lvmod <- function (time, state, parms, N, rr, ri, dr, dri) { with (as.list(parms), { PREY <- state[1:N] PRED <- state[(N+1):(2*N)] ## Fluxes due to diffusion ## at internal and external boundaries: zero gradient FluxPrey <- -Da * diff(c(PREY[1], PREY, PREY[N]))/dri FluxPred <- -Da * diff(c(PRED[1], PRED, PRED[N]))/dri ## Biology: Lotka-Volterra model Ingestion <- rIng * PREY * PRED GrowthPrey <- rGrow * PREY * (1-PREY/cap) MortPredator <- rMort * PRED ## Rate of change = Flux gradient + Biology dPREY <- -diff(ri * FluxPrey)/rr/dr + GrowthPrey - Ingestion dPRED <- -diff(ri * FluxPred)/rr/dr + Ingestion * assEff - MortPredator return (list(c(dPREY, dPRED))) }) } @ <>= options(prompt = " ") options(continue = " ") @ Then we define the parameters, which we put in a list <<>>= R <- 20 # total radius of surface, m N <- 100 # 100 concentric circles dr <- R/N # thickness of each layer r <- seq(dr/2,by = dr,len = N) # distance of center to mid-layer ri <- seq(0,by = dr,len = N+1) # distance to layer interface dri <- dr # dispersion distances parms <- c(Da = 0.05, # m2/d, dispersion coefficient rIng = 0.2, # /day, rate of ingestion rGrow = 1.0, # /day, growth rate of prey rMort = 0.2 , # /day, mortality rate of pred assEff = 0.5, # -, assimilation efficiency cap = 10) # density, carrying capacity @ After defining initial conditions, the model is solved with routine \code{ode.1D} <<>>= state <- rep(0, 2 * N) state[1] <- state[N + 1] <- 10 times <- seq(0, 200, by = 1) # output wanted at these time intervals print(system.time( out <- ode.1D(y = state, times = times, func = lvmod, parms = parms, nspec = 2, names = c("PREY", "PRED"), N = N, rr = r, ri = ri, dr = dr, dri = dri) )) @ The \code{summary} method provides summaries for both 1-dimensional state variables: <<>>= summary(out) @ while the S3-method \code{subset} can be used to extract only specific values of the variables: <<>>= p10 <- subset(out, select = "PREY", subset = time == 10) head(p10, n = 5) @ We first plot both 1-dimensional state variables at once; we specify that the figures are arranged in two rows, and 2 columns; when we call \code{image}, we overrule the default mfrow setting (\code{mfrow = NULL}). Next we plot "PREY" again, once with the default xlim and ylim, and next zooming in. Note that xlim and ylim are a list here. When we call \code{image} for the second time, we overrule the default \code{mfrow} setting by specifying (\code{mfrow = NULL}). %% This is done offline. %%<>= \begin{verbatim} image(out, grid = r, mfrow = c(2, 2), method = "persp", border = NA, ticktype = "detailed", legend = TRUE) image(out, grid = r, which = c("PREY", "PREY"), mfrow = NULL, xlim = list(NULL, c(0, 10)), ylim = list(NULL, c(0, 5)), add.contour = c(FALSE, TRUE)) \end{verbatim} %%@ \setkeys{Gin}{width=1.0\textwidth} \begin{figure} \begin{center} %%<>= %%<> %%@ \includegraphics{image1D.png} \end{center} \caption{image plots} \label{fig:plotimg} \end{figure} \subsubsection{2-D model output} When using \code{image} with a 2-D model, then the 2-D values at all output times will be plotted. Sometimes we want only output at a specific time value. We then use \proglang{S3}-method \code{subset} to extract 2-D variables at suitable time-values and use \proglang{R}'s \code{image}, \code{filled.contour} or \code{contour} method to depict them. Consider the very simple 2-D model (100*100), containing just 1-st order consumption, at a rate \code{r_x2y2}, where \code{r_x2y2} depends on the position along the grid. First the derivative function is defined: <<>>= Simple2D <- function(t, Y, par) { y <- matrix(nrow = nx, ncol = ny, data = Y) # vector to 2-D matrix dY <- - r_x2y2 * y # consumption return(list(dY)) } @ Then the grid is created, and the consumption rate made a function of grid position (\code{outer}). <<>>= dy <- dx <- 1 # grid size nx <- ny <- 100 x <- seq (dx/2, by = dx, len = nx) y <- seq (dy/2, by = dy, len = ny) # in each grid cell: consumption depending on position r_x2y2 <- outer(x, y, FUN=function(x,y) ((x-50)^2 + (y-50)^2)*1e-4) @ After defining the initial values, the model is solved using solver \code{ode.2D}. We use Runge-Kutta method \code{ode45}. <<>>= C <- matrix(nrow = nx, ncol = ny, 1) ODE3 <- ode.2D(y = C, times = 1:100, func = Simple2D, parms = NULL, dimens = c(nx, ny), names = "C", method = "ode45") @ We print a summary, and extract the 2-D variable at \code{time = 50} <<>>= summary(ODE3) t50 <- matrix(nrow = nx, ncol = ny, data = subset(ODE3, select = "C", subset = (time == 50))) @ We use function \code{contour} to plot both the consumption rate and the values of the state variables at \code{time = 50}. <>= par(mfrow = c(1, 2)) contour(x, y, r_x2y2, main = "consumption") contour(x, y, t50, main = "Y(t = 50)") @ \setkeys{Gin}{width=1.0\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{Contour plot of 2-D variables} \label{fig:twoD} \end{figure} \clearpage \section{Troubleshooting} \subsection{Avoiding numerical errors} The solvers from \pkg{ODEPACK} should be first choice for any problem and the defaults of the control parameters are reasonable for many practical problems. However, there are cases where they may give dubious results. Consider the following Lotka-Volterra type of model: <<>>= PCmod <- function(t, x, parms) { with(as.list(c(parms, x)), { dP <- c*P - d*C*P # producer dC <- e*P*C - f*C # consumer res <- c(dP, dC) list(res) }) } @ and with the following (biologically not very realistic)% \footnote{they are not realistic because producers grow unlimited with a high rate and consumers with 100 \% efficiency} parameter values: <<>>= parms <- c(c = 10, d = 0.1, e = 0.1, f = 0.1) @ After specification of initial conditions and output times, the model is solved -- using \code{lsoda}: <<>>= xstart <- c(P = 0.5, C = 1) times <- seq(0, 200, 0.1) out <- ode(y = xstart, times = times, func = PCmod, parms = parms) tail(out) @ We see that the simulation was stopped before reaching the final simulation time and both producers and consumer values may have negative values. What has happened? Being an implicit method, \code{lsoda} generates very small negative values for producers, from day 40 on; these negative values, small at first grow in magnitude until they become infinite or even NaNs (not a number). This is because the model equations are not intended to be used with negative numbers, as negative concentrations are not realistic. A quick-and-dirty solution is to reduce the maximum time step to a considerably small value (e.g. \code{hmax = 0.02} which, of course, reduces computational efficiency. However, a much better solution is to think about the reason of the failure, i.e in our case the \textbf{absolute} accuracy because the states can reach very small absolute values. Therefore, it helps here to reduce \code{atol} to a very small number or even to zero: <<>>= out <- ode(y = xstart,times = times, func = PCmod, parms = parms, atol = 0) matplot(out[,1], out[,2:3], type = "l", xlab = "time", ylab = "Producer, Consumer") @ It is, of course, not possible to set both, \code{atol} and \code{rtol} simultaneously to zero. As we see from this example, it is always a good idea to test simulation results for plausibility. This can be done by theoretical considerations or by comparing the outcome of different ODE solvers and parametrizations. \subsection{Checking model specification} If a model outcome is obviously unrealistic or one of the \ds functions complains about numerical problems it is even more likely that the ``numerical problem'' is in fact a result of an unrealistic model or a programming error. In such cases, playing with solver parameters will not help. Here are some common mistakes we observed in our models and the codes of our students: \begin{itemize} \item The function with the model definition must return a list with the derivatives of all state variables in correct order (and optionally some global values). Check if the number and order of your states is identical in the initial states \code{y} passed to the solver, in the assignments within your model equations and in the returned values. Check also whether the return value is the last statement of your model definition. \item The order of function arguments in the model definition is \code{t, y, parms, ...}. This order is strictly fixed, so that the \ds solvers can pass their data, but naming is flexible and can be adapted to your needs, e.g. \code{time, init, params}. Note also that all three arguments must be given, even if \code{t} is not used in your model. \item Mixing of variable names: if you use the \code{with()}-construction explained above, you must ensure to avoid naming conflicts between parameters (\code{parms}) and state variables (\code{y}). \end{itemize} The solvers included in package \ds are thoroughly tested, however they come with \textbf{no warranty} and the user is solely responsible for their correct application. If you encounter unexpected behavior, first check your model and read the documentation. If this doesn't help, feel free to ask a question to an appropriate mailing list, e.g. \email{r-help@r-project.org} or, more specific, \email{r-sig-dynamic-models@r-project.org}. \subsection{Making sense of deSolve's error messages} As many of \pkg{deSolve}'s functions are wrappers around existing \proglang{FORTRAN} codes, the warning and error messages are derived from these codes. Whereas these codes are highly robust, well tested, and efficient, they are not always as user-friendly as we would like. Especially some of the warnings/error messages may appear to be difficult to understand. Consider the first example on the \code{ode} function: <<>>= LVmod <- function(Time, State, Pars) { with(as.list(c(State, Pars)), { Ingestion <- rIng * Prey * Predator GrowthPrey <- rGrow * Prey * (1 - Prey/K) MortPredator <- rMort * Predator dPrey <- GrowthPrey - Ingestion dPredator <- Ingestion * assEff - MortPredator return(list(c(dPrey, dPredator))) }) } pars <- c(rIng = 0.2, # /day, rate of ingestion rGrow = 1.0, # /day, growth rate of prey rMort = 0.2 , # /day, mortality rate of predator assEff = 0.5, # -, assimilation efficiency K = 10) # mmol/m3, carrying capacity yini <- c(Prey = 1, Predator = 2) times <- seq(0, 200, by = 1) out <- ode(func = LVmod, y = yini, parms = pars, times = times) @ This model is easily solved by the default integration method, \code{lsoda}. Now we change one of the parameters to an unrealistic value: \code{rIng} is set to $100$. This means that the predator ingests 100 times its own body-weight per day if there are plenty of prey. Needless to say that this is very unhealthy, if not lethal. Also, \code{lsoda} cannot solve the model anymore. Thus, if we try: <>= pars["rIng"] <- 100 out2 <- ode(func = LVmod, y = yini, parms = pars, times = times) @ A lot of seemingly incomprehensible messages will be written to the screen. We repeat the latter part of them: \begin{verbatim} DLSODA- Warning..Internal T (=R1) and H (=R2) are such that in the machine, T + H = T on the next step (H = step size). Solver will continue anyway. In above message, R1 = 53.4272, R2 = 2.44876e-15 DLSODA- Above warning has been issued I1 times. It will not be issued again for this problem. In above message, I1 = 10 DLSODA- At current T (=R1), MXSTEP (=I1) steps taken on this call before reaching TOUT In above message, I1 = 5000 In above message, R1 = 53.4272 Warning messages: 1: In lsoda(y, times, func, parms, ...) : an excessive amount of work (> maxsteps ) was done, but integration was not successful - increase maxsteps 2: In lsoda(y, times, func, parms, ...) : Returning early. Results are accurate, as far as they go \end{verbatim} The first sentence tells us that at T = 53.4272, the solver used a step size H = 2.44876e-15. This step size is so small that it cannot tell the difference between T and T + H. Nevertheless, the solver tried again. The second sentence tells that, as this warning has been occurring 10 times, it will not be outputted again. As expected, this error did not go away, so soon the maximal number of steps (5000) has been exceeded. This is indeed what the next message is about: The third sentence tells that at T = 53.4272, maxstep = 5000 steps have been done. The one before last message tells why the solver returned prematurely, and suggests a solution. Simply increasing maxsteps will not work and it makes more sense to first see if the output tells what happens: <>= plot(out2, type = "l", lwd = 2, main = "corrupt Lotka-Volterra model") @ You may, of course, consider to use another solver: <>= pars["rIng"] <- 100 out3 <- ode(func = LVmod, y = yini, parms = pars, times = times, method = "ode45", atol = 1e-14, rtol = 1e-14) @ but don't forget to think about this too and, for example, increase simulation time to 1000 and try different values of \code{atol} and \code{rtol}. We leave this open as an exercise to the reader. \setkeys{Gin}{width=\textwidth} \begin{figure} \begin{center} <>= <> @ \end{center} \caption{A model that cannot be solved correctly} \label{fig:err} \end{figure} \clearpage %\section{Function overview} \begin{table*}[b] \caption{Summary of the functions that solve differential equations}\label{tb:rs} \centering \begin{tabular}{p{.15\textwidth}p{.75\textwidth}}\hline \rule[-3mm]{0mm}{8mm} Function &Description\\ \hline \hline \code{ode} & integrates systems of ordinary differential equations, assumes a full, banded or arbitrary sparse Jacobian \\ \hline \code{ode.1D} & integrates systems of ODEs resulting from 1-dimensional reaction-transport problems \\ \hline \code{ode.2D} & integrates systems of ODEs resulting from 2-dimensional reaction-transport problems \\ \hline \code{ode.3D} & integrates systems of ODEs resulting from 3-dimensional reaction-transport problems \\ \hline \code{ode.band} & integrates systems of ODEs resulting from unicomponent 1-dimensional reaction-transport problems \\ \hline \code{dede} & integrates systems of delay differential equations \\ \hline \code{daspk} & solves systems of differential algebraic equations, assumes a full or banded Jacobian \\ \hline \code{radau} & solves systems of ordinary or differential algebraic equations, assumes a full or banded Jacobian; includes a root solving procedure \\ \hline \code{lsoda} & integrates ODEs, automatically chooses method for stiff or non-stiff problems, assumes a full or banded Jacobian \\ \hline \code{lsodar} & same as \code{lsoda}, but includes a root-solving procedure \\ \hline \code{lsode} or \code{vode} & integrates ODEs, user must specify if stiff or non-stiff assumes a full or banded Jacobian; Note that, as from version 1.7, \code{lsode} includes a root finding procedure, similar to \code{lsodar}. \\ \hline \code{lsodes} & integrates ODEs, using stiff method and assuming an arbitrary sparse Jacobian. Note that, as from version 1.7, \code{lsodes} includes a root finding procedure, similar to \code{lsodar} \\ \hline \code{rk} & integrates ODEs, using Runge-Kutta methods (includes Runge-Kutta 4 and Euler as special cases) \\ \hline \code{rk4} & integrates ODEs, using the classical Runge-Kutta 4th order method (special code with less options than \code{rk}) \\ \hline \code{euler} & integrates ODEs, using Euler's method (special code with less options than \code{rk}) \\ \hline \code{zvode} & integrates ODEs composed of complex numbers, full, banded, stiff or nonstiff \\ \hline \hline \end{tabular} \end{table*} \begin{table*}[b] \caption{Meaning of the integer return parameters in the different integration routines. If \code{out} is the output matrix, then this vector can be retrieved by function \code{attributes(out)\$istate}; its contents is displayed by function \code{diagnostics(out)}. Note that the number of function evaluations, is without the extra evaluations needed to generate the output for the ordinary variables. } \centering \begin{tabular}{p{.05\textwidth}p{.95\textwidth}}\hline \rule[-3mm]{0mm}{8mm} Nr &Description\\ \hline \hline 1 & the return flag; the conditions under which the last call to the solver returned. For \code{lsoda, lsodar, lsode, lsodes, vode, rk, rk4, euler} these are: 2: the solver was successful, -1: excess work done, -2: excess accuracy requested, -3: illegal input detected, -4: repeated error test failures, -5: repeated convergence failures, -6: error weight became zero \\ \hline 2 & the number of steps taken for the problem so far\\ \hline 3 & the number of function evaluations for the problem so far\\ \hline 4 & the number of Jacobian evaluations so far\\ \hline 5 & the method order last used (successfully)\\ \hline 6 & the order of the method to be attempted on the next step\\ \hline 7 & If return flag = -4,-5: the largest component in the error vector\\ \hline 8 & the length of the real work array actually required. (\proglang{FORTRAN} code)\\ \hline 9 & the length of the integer work array actually required. (\proglang{FORTRAN} code)\\ \hline 10 & the number of matrix LU decompositions so far\\ \hline 11 & the number of nonlinear (Newton) iterations so far\\ \hline 12 & the number of convergence failures of the solver so far\\ \hline 13 & the number of error test failures of the integrator so far\\ \hline 14 & the number of Jacobian evaluations and LU decompositions so far\\ \hline 15 & the method indicator for the last succesful step, 1 = adams (nonstiff), 2 = bdf (stiff)\\ \hline 17 & the number of nonzero elements in the sparse Jacobian\\ \hline 18 & the current method indicator to be attempted on the next step, 1 = adams (nonstiff), 2 = bdf (stiff)\\ \hline 19 & the number of convergence failures of the linear iteration so far\\ \hline \hline \end{tabular} \end{table*} \begin{table*}[b] \caption{Meaning of the double precision return parameters in the different integration routines. If \code{out} is the output matrix, then this vector can be retrieved by function \code{attributes(out)\$rstate}; its contents is displayed by function \code{diagnostics(out)}} \centering \begin{tabular}{p{.05\textwidth}p{.95\textwidth}}\hline \rule[-3mm]{0mm}{8mm} Nr &Description\\ \hline \hline 1 & the step size in t last used (successfully)\\ \hline 2 & the step size to be attempted on the next step\\ \hline 3 & the current value of the independent variable which the solver has actually reached\\ \hline 4 & a tolerance scale factor, greater than 1.0, computed when a request for too much accuracy was detected\\ \hline 5 & the value of t at the time of the last method switch, if any (only \code{lsoda, lsodar}) \\ \hline \hline \end{tabular} \end{table*} %% this adds References to the PDF-Index without adding an obsolete section \phantomsection \addcontentsline{toc}{section}{References} \bibliography{integration} \end{document}