Solving Diﬀerential Equations in R : Package deSolve

In this paper we present the R package deSolve to solve initial value problems (IVP) written as ordinary diﬀerential equations (ODE), diﬀerential algebraic equations (DAE) of index 0 or 1 and partial diﬀerential equations (PDE), the latter solved using the method of lines approach. The diﬀerential equations can be represented in R code or as compiled code. In the latter case, R is used as a tool to trigger the integration and post-process the results, which facilitates model development and application, whilst the compiled code sig-niﬁcantly increases simulation speed. The methods implemented are eﬃcient, robust, and well documented public-domain Fortran routines. They include four integrators from the ODEPACK package (LSODE, LSODES, LSODA, LSODAR), DVODE and DASPK2.0. In addition, a suite of Runge-Kutta integrators and special-purpose solvers to eﬃciently integrate 1, 2-and 3-dimensional partial diﬀerential equations are available. The routines solve both stiﬀ and non-stiﬀ systems, and include many options, e.g., to deal in an eﬃcient way with the sparsity of the Jacobian matrix, or ﬁnding the root of equations. In this article, our objectives are threefold: (1) to demonstrate the potential of using R for dynamic modeling, (2) to highlight typical uses of the diﬀerent methods implemented and (3) to compare the performance of models speciﬁed in R code and in compiled code for a number of test cases. These comparisons demonstrate that, if the use of loops is avoided, R code can eﬃciently integrate problems comprising several thousands of state variables. Nevertheless, the same problem may be solved from 2 to more than 50 times faster by using compiled code compared to an implementation using only R code. Still, amongst the beneﬁts of R are a more ﬂexible and interactive implementation, better readability of the code, and access to R ’s high-level procedures. deSolve is the successor of package odesolve which will be deprecated in the future; it is free software and distributed under the GNU General Public License, as part of the R software project.


Introduction
Many phenomena in science and engineering can be mathematically represented as initial value problems (IVP) of ordinary differential equations (ODE, Asher and Petzold 1998).ODEs describe how a certain quantity changes as a function of time or space, or some other variable (called the independent variable).They can be mathematically represented as: where y are the differential variables, y are the derivatives, v are other variables, and t is the independent variable.For the remainder, we will assume that the independent variable is "time".For this equation to have a solution, an extra condition is required.Here we deal only with models where some initial condition (at t = t 0 ) is specified: These are called initial value problems (IVP).The formalism above provides an explicit expression for y as a function of y, x and t.A more general mathematical form is the implicit expression: If, in addition to the ordinary differential equations, the differential variables obey some algebraic constraints at each time point: then we obtain a set of differential algebraic equations (DAE).The two previous functions G (eq. 1) and g (eq.2) can be combined to a new function F : 0 = F (y , y, v, t) which is the formalism that we will use in this paper.Solving a DAE is more complex than solving an ODE.For instance, the initial conditions for a DAE must be chosen to be consistent.This is, the initial values of t, y and y , must obey: 0 = F (y (t 0 ), y(t 0 ), v, t 0 ) DAEs are commonly encountered in a number of scientific and engineering disciplines, e.g., in the modelling of electrical circuits or mechanical systems, in constrained variational problems, or in equilibrium chemistry (e.g., Brenan, Campbell, and Petzold 1996).
Most of the ODEs and DAEs are complicated enough to preclude finding an analytical solution, and therefore they are solved by numerical techniques, which calculate the solution only at a limited number of values of the independent variable (t).
A common theme in many of the numerical solvers, are their capabilities to solve "stiff " ODE or DAE problems.Formally, if the eigenvalue spectrum of the ODE system (i.e., of its Jacobian, see below) is large, the ODE system is said to be stiff (Hairer and Wanner 1980).As a less formal definition, an ODE system is called stiff if the problem changes on a wide variety of time scales, i.e., it contains both very rapidly and very slowly changing terms.Unless these stiff problems are solved with especially-designed methods, they require an excessive amount of computing time, as they need to use very small time steps to satisfy stability requirements (Press, Teukolsky, Vetterling, and Flannery 2007, p. 931).
Very often, stiff systems are most efficiently solved with implicit methods, which require the creation of a Jacobian matrix ( ∂f ∂y ) and the solution of a system of equations involving this Jacobian.As we will see below, there is much to be gained by taking advantage of the sparsity of the Jacobian matrix.Except for the Runge-Kutta methods, all solvers implemented in deSolve are variable order, variable step methods, that use the backward differentiation formulas and Adams methods, two important families of multistep methods (Asher and Petzold 1998).
The remainder of the paper is organized as follows.In Section 2, the different solvers are briefly discussed and some implementation issues noted.Section 3 gives some example implementations of ODE, PDE and DAE systems in R (R Development Core Team 2009).In Section 4, we demonstrate how to implement the models in a compiled language.Numerical benchmarks of computational performance are conducted in Section 5. Finally, concluding remarks are given in Section 6.
The package is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=deSolve.

Implementation issues
The R package deSolve (Soetaert, Petzoldt, and Setzer 2009) is the successor of package odesolve (Setzer 2001), which will be deprecated in the future.Compared to odesolve, it includes a more complete set of integrators, a more extensive set of options to tune the integration routines, and provides more complete output.For instance, there was no provision to specify the structure of the Jacobian in odesolve's routine lsoda, whereas this is now fully supported in deSolve.Moreover, as the applicability domain of the new package now includes DAEs and PDEs, the name odesolve was considered too narrow, warranting a new one (deSolve).Several integration methods in deSolve implement efficient, robust, and frequentlyused open source routines.These routines are similar to one another and well documented, both in the source code, or by separate works (e.g., Brenan et al. 1996, for DASSL).The lsoda (Petzold 1983) implementation in deSolve is fully compatible with that in odesolve for systems coded fully in R.However, the calling sequence for systems using native language calls has changed between odesolve and deSolve.In the current version, the solvers take care of forcing function interpolation in compiled code.They also support events and time lags.
The collaboration between the three authors was greatly facilitated by the use of R-Forge (Theußl and Zeileis 2009, http://R-Forge.R-project.org/), the framework for R project developers, based on GForge (Copeland, Mas, McCullagh, Perdue, Smet, and Spisser 2006).

Integration options
When calling the integration routines, many options can be specified.We have tried, as far as possible, to retain the flexibility of the original codes; for most applications the defaults will do.
The sparsity structure of the Jacobian is crucial to efficiently solve (moderately) large systems of stiff equations.Sparse Jacobians can not only be generated much faster (fewer function calls), storing only the nonzero elements greatly reduces memory requirements.In addition, the resulting equations can be much more efficiently solved if the sparsity is taken into account.
Therefore, users have the option to specify whether the Jacobian has certain particular properties.By default it is considered to be a full matrix which is calculated by the solver, by numerical differencing (i.e., where the function gradient is estimated by successive perturbation of the state variables).To take advantage of a sparse Jacobian, the solver can be informed of any sparsity patterns.Thus, it is possible to specify that the Jacobian has a banded structure (vode, daspk, lsode, lsoda), or to use a more general sparse Jacobian (lsodes).In the latter case, lsodes can by itself determine the sparsity, but the user can also provide a matrix with row and column positions of its nonzero elements.In addition, ode.1D, ode.2D and ode.3D are specially designed to deal efficiently with the sparsity that arises in (PDE) models described in 1-, 2-and 3 (spatial) dimensions.With the exception of the Runge-Kutta solvers, all integration methods also provide the specification of an analytical Jacobian as an option, which may improve performance.Note, though, that a clear advantage of the finite difference approximation is that it is simpler.
Other important options are rtol and atol, relative and absolute tolerances that define error control.They are important because they not only affect the integration step taken, but also the numerical differencing used in the creation of the Jacobian.

A short description of the integrators
All lsode-type codes, vode and daspk use the variable-step, variable-order backward differentiation formula (BDF), suitable for solving stiff problems (order 1-5).The lsode family and vode also contain variable-step, variable-order Adams methods (order 1-12), which are well suited for nonstiff problems (Asher and Petzold 1998).
In detail: ode, ode.1D, ode.2D and ode.3D are wrappers around the integration routines described below.The latter three are especially designed to solve partial differential equations, where, in addition to the time derivative, the components also change in one, two or three (spatial) dimensions.lsoda automatically selects a stiff or nonstiff method.It may switch between the two methods during the simulation, in case the stiffness of the system changes.This is the default method used in ode and especially well suited for simple problems.
lsodar is similar to lsoda but includes a method to find the root of a function.lsode, and vode also solve stiff and nonstiff problems.However, the user must decide whether a stiff or nonstiff method is most suited for a particular problem and select an appropriate solution method.zvode is a variant of vode that solves equations involving variables that are complex numbers.lsode is the default solver used in ode.1D lsodes exploits the sparsity in the Jacobian matrix by using linear algebra routines from the Yale sparse matrix package (Eisenstat, Gursky, Schultz, and Sherman 1982).It can determine the sparsity on its own or take it as input.Especially for large stiff problems with few interactions between state variables (leading to sparse Jacobians), dramatic savings in computing time can be achieved when using lsodes.It is the solver used in ode.2D and ode.3D daspk is the only integrator in the package that also solves differential algebraic equations of index zero and one.It can also solve ODEs.
Finally, the package also includes solvers for several methods of the Runge-Kutta family (rk), with variable or fixed time steps.This includes the classical 4th order Runge-Kutta and the Euler method (rk4, euler).

Output
All solvers return an array that contains, in its columns, the time values (1 st column) and the values of all state variables (subsequent columns) followed by the ordinary output variables (if any).This format is particularly suited for graphical routines of R (e.g., matplot).In addition, a plot method is included which, for models with not too many state variables, plots all output in one figure.
All Fortran codes have in common that they monitor essential properties of the integration, such as the number of Jacobian evaluations, the number of time steps performed, the number of integration error test failures, the stepsize to be attempted on the next step and so on.These performance indicators can be called upon by a method called diagnostics.

Examples: Model implementations in R
In this section we first implement a simple biological model, the Lotka-Volterra consumerprey model, which is solved with the integration routine ode (which uses method lsoda).This model is then extended with a root function that stops the simulation at steady-state and which uses routine lsodar.An implementation in a 1-D and 2-D setting, demonstrates the capabilities of ode.1D and ode.2D.Finally, we end with a simple DAE model.Many more examples can be found in the deSolve package example files.Each model was run on a 2.5 GHz portable pc with an Intel Core 2 Duo T9300 processor and 3 GB of RAM.The CPU-times reported were estimated as the mean of 10 runs as follows: R> print(system.time(for(i in 1:10) + out <-ode(func = LVmod, y = yini, parms = pars, times = times) + )/10)

A simple consumer-prey model
We start with a simple Lotka-Volterra type of model (Lotka 1925;Volterra 1926) describing consumer-prey interactions: Where P and C are prey and consumer concentrations, r G is the growth rate of prey, K the carrying capacity, r I the ingestion rate of the consumer, k AE its assimilation efficiency and r M the consumer's mortality rate.The implementation in R consists of three parts: First the model function is defined.Here this function is called LVmod0D.It takes as input the current simulation time, the values of the state variables, and the model parameters.These three arguments have to be always present, and in this order; other arguments can be added after them.The solver will call this function at each time step during the integration process; at which Time contains the current simulation time; State the current value of the state variables and Pars the values of the parameters as passed to the solver.
Both State and Pars are a vector, with named elements; the with statement, allows using their names within the function.
The model returns a list, where the first element contains the derivatives, concatenated.Note that Time is not used here, but in many models, it is used, e.g., when there are external variables that depend on time.

R>
Then the parameters are given a name and a value (pars), the state variables initialized (yini) and the time points at which we want output specified (times).Based on these inputs, the model is simulated.Here we use the default integration function ode, which is based on the lsoda method; its returned model output is written in a matrix called out.We print the first part of this matrix (head(out)).Matrix out has in its first column the time sequence, and in its next columns the prey and consumer concentrations.

The consumer-prey model with stopping criterion
Supposing that we are interested in the initial phase only, we use the root finding functionality of function lsodar to halt the simulation after the state variables change less than some predefined amount (here 10 −4 ).Below, we define the root function (rootfun), which first estimates the rate of change and then calculates the difference between the sum of absolute values and a given tolerance (10 −4 ).If we run lsodar, the integration will stop if the sum of absolute values equals 10 −4 ; this is after 93.5 days.The results are depicted in Figure 1B; it takes 0.05 seconds to complete.This is slightly longer than lsoda takes to simulate the system for 200 days (0.04 seconds), since LVmod0D is called twice as often: once to evaluate the derivative and once to evaluate the root.

Consumer and prey dispersing on a 1-D grid
Now the same consumer-prey dynamics is implemented on a 1-dimensional grid, and including dispersion with dispersion coefficient Da (Crank 1975).The extended formulations are: These partial differential equations (PDE) are solved by discretizing in space first and integrating the resulting initial value ODE; this is a technique called the "method of lines" (Schiesser 1991).It is beyond the scope of this paper to derive how the spatial derivative is numerically approximated and implemented in R; interested readers may refer to Soetaert and Herman (2009) where this is discussed.Also, package ReacTran The matrix out has in its first column the time sequence, followed by 1000 columns with prey concentrations, one for each box, followed by 1000 columns with consumer concentrations.We plot only the prey concentrations (Figure 2).Function ode.1D was run using either lsode, vode, lsoda and lsodes as the integrator; it took 0.8 (lsode), 0.85 (vode), 1.2 (lsoda) and 0.65 (lsodes) seconds to finish the run.

Consumer and prey dispersing on a 2-D grid
Finally, we also implement the same consumer-prey dynamics on a 2-dimensional grid.The extended formulations now include dispersion in the x-and y-direction (dispersion coefficient Da): The function below implements this 2-D model.Note that, as for the 1-D case, the use of explicit looping is avoided: to estimate the gradient, we just subtract two matrices shifted with one row (x-direction) or one column (y-direction).The zero-fluxes at the boundaries are implemented by binding a row or column of 0-values (zero).Results are in Figure 3 for the initial condition, and after 20, 30 and 40 days.

R>
Note: with the initial conditions used (nonzero concentration in the centre), this is not a particularly clever way of modeling dispersion on a 2-D surface.For this special case it is much more efficient to represent these dynamics in a 1-dimensional model, and using cylindrical coordinates; this model is included as an example model in the help file of ode.1D.The 2-D implementation here was added just for illustrative purposes.

A chemical example: DAE
In chemistry, stiffness frequently arises from the fact that some reactions occur much faster than others.One way to deal with this stiffness is to reformulate the ODE as a DAE (e.g., Hofmann, Meysman, Soetaert, and Middelburg 2008).Consider the following model: Three chemical species A, B, D are kept in a vessel; the following reversible reaction occurs: In addition, D is produced at a constant rate, k prod , while B is consumed at a 1 st order rate, r.Implementing this model as an ODE system: The ODEs are now reformulated as a DAE.If the reversible reactions (involving k 1 and k 2 ) are much faster compared to the other rates (k prod , r) then the three quantities D, A and B can be assumed to be in local equilibrium.Thus, at all times, the following relationship exists between the concentrations of A, B and D (here concentration of x is denoted with [x]): where K = k 2 /k 1 is the so-called equilibrium constant.The equilibrium description is complete by taking linear combinations of rates of changes, such that the fast reversible reactions vanish.
In this DAE, the fast equilibrium reactions (involving k 1 and k 2 ) have been removed.
DAEs are specified in implicit form (see Section 1): In R we define a function that takes as input time (t), the values of the state variables (y) and their derivatives (yprime) and the parameter vector (pars), and that returns the results as a list; the first element of this list contains the implicit form of the differential equations (res1, res2) and of the algebraic equation (eq), concatenated.Additionally, other quantities can be returned as well (here CONC).Note that y, yprime and pars are vectors with named elements; their names are made available through the with statement.

Model implementation in a compiled language
This final model is an example of writing the 0-dimensional consumer-prey model equations in a compiled language such as Fortran, C, or C++, and embedding it in R code.Implementations in compiled languages (cf.Tables 1 and 2) consist of: an initializer function, which sets the values of the parameters (initparms() in the example), if forcing functions are to be used, an initializer for the data for the forcing function (here absent), the model function which calculates the rate of change and output variables (derivs() in the example), and -optionally-a function that calculates an analytic Jacobian (here absent).
Each function has a standard calling sequence (see the package vignette compiledCode in the package deSolve for the details).The initializer subroutines serve just to link data from the R side of things with memory accessible to the native code, and will rarely be more complicated than the example shown here.
The bulk of the computation is carried out in the subroutine that defines the system of differential equations.In the initializer routine, parameters are passed to the native programs as one vector (containing 5 values).In Fortran, parameters are stored in a common block, in which the values are given a name (rI,..) in the model function to make it easier to understand the code, while it is a vector in the initializer routine.In the C code names can be assigned to these parameters as well as state variables and their derivatives via #define statements that make the code more readable.The variable yout(1) holds the total concentration of consumer and prey at any given time.The model function must check whether enough memory is allocated for the output variable(s), to prevent a fatal error if the memory allocation is inadequate.The subroutines rexit() and error() are provided by R to gracefully exit with a message back to the R command prompt.
To run this model, the code must first be compiled.Given that the appropriate toolset is installed this can be done in R itself, using the system command: system("R CMD SHLIB LVmod0D.f") or system("R CMD SHLIB LVmod0D.c").The compiler will create a file that can be linked dynamically into an R session, for example the dynamic link library (DLL) LVmod0D.dll on Windows respectively a shared library LVmod0D.so on other operating systems.
The dynamic library is loaded into the current R process using a call to dyn.load().
The gain is much more pronounced when the model is implemented in Fortran rather than in R: here the Fortran implementation executes 2 to 20 (1-D, 2-D) to 66 times (0-D) faster (than R code 1).Finally, the difference between a Fortran model triggered by R, or a model completely implemented in Fortran is very small; part of the difference is due to the checking for illegal input values in the R integration routines.

Concluding remarks
The software R is rapidly gaining in popularity among scientists.With the launch of package odesolve (Setzer 2001), it became possible to use R as a tool to solve initial value problems of ordinary differential equations.The integration routines in this package opened up an entirely new field of application, although it took a while before this was acknowledged.More recent packages (rootSolve, bvpSolve) (Soetaert 2009;Soetaert, Cash, and Mazzia 2010) offer to solve boundary value problems of differential equations.
The paper in R News by Petzoldt (2003), demonstrated the suitability of R for running dynamic (ecological) simulations.More recently, a specially designed framework for ecological modelling in R, simecol, emerged (Petzoldt and Rinke 2007); packages for inverse modelling, (FME) and reactive transport modelling (ReacTran) (Soetaert and Petzoldt 2010;Soetaert and Meysman 2009) were created, while a framework for more general continuous dynamic modeling, Rdynamic (Setzer in prep.) is under construction.An increasing number of textbooks deal with the subject (Ellner and Guckenheimer 2006;Bolker 2008;Soetaert and Herman 2009;Stevens 2009).
In order to efficiently solve a variety of differential equation models, a flexible set of integration routines is required.It is with this goal in mind that the integration routines in deSolve were selected.Whereas the original integration routine in odesolve only efficiently solved relatively simple ODE systems, the suite of routines now also includes a solver for differential algebraic equations and methods to solve partial differential equations.
In this paper we have shown that thanks to these new functions, R can now more efficiently run 0-dimensional, 1-2-, and even 3-dimensional models of small to moderate size.Apart from models implemented in pure R, it is possible to specify model functions in compiled code written in any higher-level language that can produce shared libraries (resp.DLLs).The integration routines then communicate directly with this code, without passing arguments to and from R, so R is used just to trigger the integration and post-process the results.As the entire simulation occurs in compiled code, there is no loss in execution speed compared to a model that is fully implemented in the higher level language.But even in this case, all the power of R as a pre-and post processing environment as well as its graphical and statistical facilities are immediately available -no need to import the model output from an external source.
In the examples, linking compiled code to the integrator, indeed made the model run faster with a factor 2 (for the 10000 state variable 1-D model) up to more than 50 times for the smallest (2 state variable) model application, than when implemented as an R function.There was only a small difference when running the model entirely in compiled code.
There are several reasons why compiled code is faster.First of all, R is an interpreted language, and therefore processes the program at runtime.Every line is interpreted multiple times at each time step.This makes interpreted code significantly slower than compiled code, which transform programs directly into machine code, before running.Note though that R is a vectorized language and, compared to some other interpreted languages, less performance is lost if R's high level functions, based on optimized machine code, are efficiently exploited (Ligges and Fox 2008).In our 1-D example, we used R function diff to take numerical differences, whilst in the 2-D model, entire matrices were subtracted.Because of this use of high-level functions, the simulation speed of these models, entirely specified in R, was quite impressive, approaching the implementation in Fortran.Performance of R code especially deteriorates when using loops.For instance, if the 2-D model is implemented by looping over all rows, then the simulation time increases tenfold; when looping over rows and columns, computation speed drops with 2 orders of magnitude!There are also trade-offs in using complex variable types of R, especially if R performs extensive copying or internal data conversion.For instance, the use of named variables and parameters introduced a computational overhead of around 70% in our simplest model example.However, the effect was relatively less significant, in the order of 10-20%, in more demanding models.
The use of code in a dynamically linked library also has its drawbacks.First of all, it is less flexible.Whereas it is simple to interact with models specified in R code, this is not at all the case for compiled code: before the model code can be executed, it has to be formally compiled, and the DLL loaded.Secondly, errors may be particularly hard to trace and may even cause R to terminate.The lack of easy access to R's high-level procedures is another drawback of using compiled code, where much more has to be hand-coded.Note though that, as from deSolve version 1.5, the interpolation of external signals (also called forcing functions) to the current timepoints is taken care of by the integration routines; the compiled-code equivalent of R function approxfun.
Putting these pros and cons together, the optimal approach is probably to use pure R for the initial model development (rapid prototyping).In case the model executes too slowly, or when a large number of simulations are performed, implementing the model in C, C++ or Fortran may be considered.
Finally, the creation and solution of a mathematical model is never a goal in itself.Models are used, amongst other things to challenge our understanding of a natural system, to make budgets or to quantify immeasurable processes or rates.When used in this way, the interaction with data is crucial, as is statistical treatment and graphical representation of the model outcome and the data.We hope that R's excellence in these fields, and the fact that it is entirely free, will give impetus to also using R as a modelling platform.

deSolve:Figure 2 :
Figure 2: Results of the Lotka-Volterra model on a one-dimensional grid.

Figure 3 :
Figure 3: Results of the Lotka-Volterra model in a two-dimensional grid.

Figure 4 :
Figure 4: Results of the DAE chemical model.CONC = summed concentration of A, B and D.