Solving Differential Equations in R

1 CONTRIBUTEDRESEARCHARTICLES5 Solving Differential Equations in Rby Karline Soetaert, Thomas Petzoldt and R. WoodrowSetzer1 AbstractAlthough R is still predominantly ap-plied for statistical analysis and graphical repre-sentation, it is rapidly becoming more suitablefor mathematical computing. One of the fieldswhere considerable progress has been made re-cently is the solution of Differential we give a brief overview of differentialequations that can now be solved by Equations describe exchanges of matter,energy, information or any other quantities, often asthey vary in time and/or space. Their thorough ana-lytical treatment forms the basis of fundamental the-ories in mathematics and physics, and they are in-creasingly applied in chemistry, life sciences and Equations are solved by integration,but unfortunately, for many practical applicationsin science and engineering, systems of differentialequations cannot be integrated to give an analyticalsolution, but rather need to be solved advanced numerical algorithms that solvedifferential Equations are available as (open-source)computer codes, written in programming languageslike FORTRAN or C and that are available fromrepositories like GAMS ( ) orNETLIB ( ).

2 Depending on the problem, mathematical for-malisations may consist of ordinary differentialequations (ODE), partial Differential Equations (PDE), Differential algebraic Equations (DAE), or de-lay Differential Equations (DDE). In addition, a dis-tinction is made between initial value problems (IVP)and boundary value problems (BVP).With the introduction of R-packageodesolve(Setzer, 2001), it became possible to use R (R Devel-opment Core Team, 2009) for Solving very simple ini-tial value problems of systems of ordinary differen-tial Equations , using thelsodaalgorithm of Hind-marsh (1983) and Petzold (1983). However, manyreal-life applications, including physical transportmodeling, equilibrium chemistry or the modeling ofelectrical circuits, could not be solved with this , much effort has been made toimprove R s capabilities to handle Differential equa-tions, mostly by incorporating published and welltested numerical codes, such that now a much morecomplete repertoire of Differential Equations can benumerically specifically, the following types of differen-tial Equations can now be handled with add-on pack-ages in R: Initial value problems (IVP) of ordinary differ-ential Equations (ODE), using packagedeSolve(Soetaert et al.)

3 , 2010b). Initial value Differential algebraic Equations (DAE), packagedeSolve. Initial value partial Differential Equations (PDE),packagesdeSolveandReacTra n(Soetaert and Meysman, 2010). boundary value problems (BVP) of ordinarydifferential Equations , using packagebvpSolve(Soetaert et al., 2010a), orReacTranandroot-Solve(Soetaert, 2009). Initialvaluedelaydifferentialequations(D DE), using packagesdeSolveorPBSddes-olve(Couture-Be il et al., 2010). Stochastic Differential Equations (SDE), usingpackagessde(Iacus, 2008) andpomp(Kinget al., 2008).In this short overview, we demonstrate how tosolve the first four types of Differential equationsin R. It is beyond the scope to give an exhaustiveoverview about the vast number of methods to solvethese Differential Equations and their theory, so thereader is encouraged to consult one of the numer-ous textbooks ( , Ascher and Petzold, 1998; Presset al., 2007; Hairer et al., 2009; Hairer and Wanner,2010; LeVeque, 2007, and many others).

4 In addition, a large number of analytical and nu-merical methods exists for the analysis of bifurca-tions and stability properties of deterministic sys-tems, the efficient simulation of stochastic differen-tial Equations or the estimation of parameters. Wedo not deal with these methods of Differential equationsOrdinary Differential equationsOrdinary Differential Equations describe the changeof astate variable yas a functionfof oneindependentvariable t( , time or space), ofyitself, and, option-ally, a set of other variablesp, often calledparameters:y =dydt=f(t,y,p)1 The views expressed in this paper are those of the authors and do not necessarily reflect the views or policies of the EnvironmentalProtection AgencyThe R Journal Vol. 2/2, December 2010 ISSN 2073-48596 CONTRIBUTEDRESEARCHARTICLESIn many cases, Solving Differential Equations re-quires the introduction of extra conditions. In the fol-lowing, we concentrate on the numerical treatmentof two classes of problems, namely initial value prob-lems and boundary value value problemsIf the extra conditions are specified at the initial valueof the independent variable, the Differential equa-tions are calledinitial value problems(IVP).

5 There exist two main classes of algorithms to nu-merically solve such problems, so-calledRunge-Kuttaformulas andlinear multistepformulas (Hairer et al.,2009; Hairer and Wanner, 2010). The latter containstwo important families, the Adams family and thebackward differentiation formulae (BDF).Another important distinction is betweenexplicitandimplicitmethods, where the latter methods cansolve a particular class of Equations (so-called stiff Equations ) where explicit methods have problemswith stability and efficiency. Stiffness occurs for in-stance if a problem has components with differentrates of variation according to the independent vari-able. Very often there will be a tradeoff between us-ing explicit methods that require little work per inte-gration step and implicit methods which are able totake larger integration steps, but need (much) morework for one R, initial value problems can be solved withfunctions from packagedeSolve(Soetaert et al.)

6 ,2010b), which implements many solvers from ODE-PACK (Hindmarsh, 1983), the codevode(Brownet al., 1989), the Differential algebraic equation solverdaspk(Brenan et al., 1996), all belonging to the linearmultistep methods , and comprising Adams meth-ods as well as backward differentiation former methods are explicit, the latter addition, this package contains a de-novo imple-mentation of a rather general Runge-Kutta solverbased on Dormand and Prince (1980); Prince andDormand (1981); Bogacki and Shampine (1989); Cashand Karp (1990) and using ideas from Butcher (1987)and Press et al. (2007). Finally, the implicit Runge-Kutta methodradau(Hairer et al., 2009) has beenadded value problemsIf the extra conditions are specified at differentvalues of the independent variable, the differen-tial Equations are calledboundary value problems(BVP). A standard textbook on this subject is Ascheret al. (1995).PackagebvpSolve(Soetaert et al.

7 , 2010a) imple-ments three methods to solve boundary value simplest solution method is thesingleshooting method, which combines initial value prob-lem integration with a nonlinear root finding algo-rithm (Press et al., 2007).Two more stable solu-tion methods implement a mono implicit Runge-Kutta (MIRK) code, based on the FORTRAN codetwpbvpC(Cash and Mazzia, 2005), and the collocationmethod, based on the FORTRAN codecolnew(Baderand Ascher, 1987). Some boundary value problemscan also be solved with functions from packagesRe-acTranandrootSolve(see below).Partial Differential equationsIn contrast to ODEs where there is only one indepen-dent variable, partial Differential Equations (PDE)contain partial derivatives with respect to more thanone independent variable, for instancet(time) andx(a spatial dimension).To distinguish this typeof Equations from ODEs, the derivatives are repre-sented with the symbol, y t=f(t,x,y, y x,p)Partial Differential Equations can be solved by sub-dividing one or more of the continuous independentvariables in a number of grid cells, and replacing thederivatives by discrete, algebraic approximate equa-tions, so-called finite differences (cf.

8 LeVeque, 2007;Hundsdorfer and Verwer, 2003).For time-varying cases, it is customary to discre-tise the spatial coordinate(s) only, while time is left incontinuous form. This is called the method-of-lines,and in this way, one PDE is translated into a largenumber of coupled ordinary Differential Equations ,that can be solved with the usual initial value prob-lem solvers (cf. Hamdi et al., 2007). This applies toparabolic PDEs such as the heat equation, and to hy-perbolic PDEs such as the wave time-invariant problems, usually all indepen-dent variables are discretised, and the derivatives ap-proximated by algebraic Equations , which are solvedby root-finding techniques. This technique applies toelliptic functions to gener-ate finite differences on a structured grid. After that,the resulting time-varying cases can be solved withspecially-designed functions from packagedeSolve,while time-invariant cases can be solved with root- Solving methods from algebraic equationsDifferential-algebraic Equations (DAE) contain amixture of Differential (f) and algebraic Equations (g), the latter for maintaining mass-balance con-ditions:y =f(t,y,p)0=g(t,y,p)Important for the solution of a DAE is its index of a DAE is the number of differentiationsThe R Journal Vol.

9 2/2, December 2010 ISSN 2073-4859 CONTRIBUTEDRESEARCHARTICLES7needed until a system consisting only of ODEs is (Brenan et al., 1996) from pack-agedeSolvesolves (relatively simple) DAEs of indexat most 1, while functionradau(Hairer et al., 2009)solves DAEs of index up to detailsThe implemented solver functions are explained bymeans of theode-function, used for the solution ofinitial value problems. The interfaces to the othersolvers have an analogous definition:ode(y, times, func, parms, method = c("lsoda","lsode", "lsodes", "lsodar","vode", "daspk", "euler", "rk4","ode23", "ode45", "radau", "bdf","bdf_d", "adams", "impAdams","impAdams_d"), ..)To use this, the system of Differential equationscan be defined as an R-function (func) that computesderivatives in the ODE system (the model definition)according to the independent variable ( timet).funccan also be a function in a dynamically loadedshared library (Soetaert et al., 2010c) and, in addition,some solvers support also the supply of an analyti-cally derived function of partial derivatives (Jacobianmatrix).

10 Iffuncis an R-function, it must be defined as:func <- function(t, y, parms, ..)wheretis the actual value of the independent vari-able ( the current time point in the integration),yis the current estimate of the variables in the ODEsystem,parmsis the parameter vector beused to pass additional arguments to the return value offuncshould be a list, whosefirst element is a vector containing the derivativesofywith respect tot, and whose next elements areoptional global values that can be recorded at eachpoint intimes. The derivatives must be specified inthe same order as the state on the algorithm specified in argu-mentmethod, numerical simulation proceeds eitherexactly at the time steps specified intimes, or us-ing time steps that are independent fromtimesandwhere the output is generated by interpolation. Withthe exception of methodeulerand several fixed-stepRunge-Kutta methods all algorithms have automatictime stepping, which can be controlled by setting ac-curacy requirements (see below) or by using optionalarguments likehini(initial time step),hmin(minimaltime step) andhmax(maximum time step).