Differential Equations in R(1)

• The model application consists of: • Specification of the time at which model output is wanted, • Integration of the model equations (uses Rfunctions from deSolve), • Plotting of model results.
• All integration methods can be triggered from function ode, by setting ode's argument method), or can be run as stand-alone functions. Moreover, for each integration routine, several options are available to optimise performance.
1.Solvers for initial value problems of ordinary dierential equations
• Implementation of an IVP ODE in R can be separated in two parts: the model specification and the model application. • Model specification consists of:
• • • • oute1 <- radau(state, times, expon, parameters1, atol = 1e-4, rtol = 1e-4) oute2<- ode(state, times, expon, parameters1, method = "radau", atol = 1e-4, rtol = 1e-4) outl1 <- radau(state, times, logistic, parameters2, atol = 1e-4, rtol = 1e-4) outl2<- ode(state, times, logistic, parameters2, method = "radau", atol = 1e-4, rtol = 1e-4)
specified in a function (Lorenz) that calculates the rate of change of the state variables
expon<-function(t, state, parameters1) { with(as.list(c(state, parameters1)),{ # rate of change As both parameters and state variables are `vectors', they are converted into a list. dR <- k*R # return the rate of change list(dR) }) # end with(as.list ... }
logistic<-function(t, state, parameters2) { with(as.list(c(state, parameters2)),{ # rate of change dR <- k*R*(1-R/b) # return the rate of change list(dR) }) # end with(as.list ... }
Differential Equations in R
outline
• • • • How to specify a model An overview of solver functions Plotting, scenario comparison Forcing function and events
Necessary packages
• deSolve: main integration package • rootSolve: steady-state solver • bvpSolve: boundary value problem solvers • ReacTran: partial differential equations • simecol: interactive environment for implementing models
Plotting results
• plot(out1,lwd=2,col=1,xlim=c(0,100),ylim=c(20 00,30000),xlab="time",ylab="number of rabbit",main="") • par(new=T) • plot(out2,lwd=2,col=2,xlim=c(0,100),ylim=c(20 00,30000),xlab="",ylab="",main="") • text(40,30000,"exponential model") • text(60,22000,"logistic model")
• • • • • • print(system.time(outl3 <- rk4 (state, times, logistic, parameters))) print(system.time(outl4 <- lsode (state, times, logistic, parameters))) print(system.time(outl5 <- lsoda (state, times, logistic, parameters))) print(system.time(out l6<- lsodes(state, times, logistic, parameters))) print(system.time(out l7<- daspk (state, times, logistic, parameters))) print(system.time(out l8<- voቤተ መጻሕፍቲ ባይዱe (state, times, logistic, parameters)))
Species model
• • • • Exponential model: dR/dt=KR,R(0)=2000,k=0.1 Logistic model dR/dt=KR(1-R/b),R(0)=2000,k=0.1,b=25000
Model parameters
library(deSolve) ##model specification ##parameters parameters1<-c(k=0.1) parameters2<-c(k=0.1,b=25000) ##state variable state<-c(R=2000) ##model equation
2.Solvers for initial value problems of ordinary differential equations
• Package deSolve contains several IVP ordinary differential equation solvers, that belong to the most important classes of solvers. • the Backward Differentiation Formulae and Adams methods from ODEPACK • The implicit Runge-Kutta method RADAU • The package contains also de novo implementation of several Runge-Kutta methods
Model integration
• The model is solved using deSolve function ode, which is the default integration routine. Function ode takes as input, a.o. the state variable vector (y), the times at which output is required (times), the model function that returns the rate of change (func) and the parameter vector (parms). • Function ode returns an object of class deSolve with a matrix that contains the values of the state variables (columns) at the requested output times. • out1 <- ode(y = state, times = times, func = expon, parms = parameters1) • out2 <- ode(y = state, times = times, func = logistic, parms = parameters2) • head(out1);head(out2)
• The default integration method, based on the FORTRAN code LSODA is one that switches automatically between stiff and non-stiff systems. • This is a very robust method, but not necessarily the most efficient solver for one particular problem.
Model application
• Time specification • We run the model for 100 months and give output at 0.01 intervals. R's function seq() creates the time sequence: • times <- seq(0, 100, by = 0.01)
Installing the R software and package
• Downloading R from the R-project website: http://www.r-project.org
• Packages can be installed from within the Rsoftware • Or via command line: • install.package(“deSolve”,dependencies=T)
• 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:
• Dening model parameters and their values, • Dening model state variables and their initial conditions, • Implementing the model equations that calculate the rate of change (e.g. dX=dt) of the state variables.