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.

2 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 ( ).

3 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).

4 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.)

5 , 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).

6 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).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.

7 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.

8 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).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).

9 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.)

10 ,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.