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.

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

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

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

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

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

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

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