Introduction to R deSolve package

package_deSolve.Rmd

---
title: "R Notebook"
output: html_notebook
---

# Solving ODEs with the package `DeSolve`

## A first example

Load the library `deSolve`.

```{r}
library(deSolve)
```

Define the ODE system to solve: equations, parameter values, initial condition, and time span.  Equations are defined as a `function`, the function must be defined as

```
func <- function(t,y,pars,...)
```

where `t` is the time, `y` is the vector of the state variables and `pars` is a vector or list of parameters. 
The return value of the function should be a list, whose first element is a vector of the derivative of y with 
respect to time: dy/dt. The next elements are used to pass global values. For example, the equations for the Lorenz system for atmoshperic chaos are defined in the function `Lorenz`:

```{r}
Lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dX <-  a * X + Y * Z
    dY <-  b * (Y - Z)
    dZ <- -X * Y + c * Y - Z
    list(c(dX, dY, dZ))
  })
}
```

Here, the argument `state` expects an array with names `X`, `Y` and `Z`. The initial conditions both define the names and the initial values of the variables.

```{r}
vars <- c(X = 1, Y = 1, Z = 1)
```

The parameter values are set in a named list:

```{r}
params <- c(a = -8/3, b = -10, c = 28)
```

The time points at which the solution musst be computed

```{r}
tim <- seq(0, 100, by = 0.01)
```

The main interface for solving ODEs is the function `ode`. It takes following arguments:

```
ode(y, times, func, parms)
```

For the Lorenz system we have just defined, the call to  `ode` is

```{r}
sol <- ode(y = vars, times = tim, func = Lorenz, parms = params)
```

The output is a deSolve matrix. The first colmun is time, the other columns are the solution for each of the variables of the equations. The is a plot method (`plot.deSolve`) for plotting the solution, with can be called just with `plot`.

```{r}
plot(sol)
```


The solution can also  be plotted with the default plot command

```{r}
plot(sol[,'Y'], sol[,'Z'], ty = 'l')
```

Or the argument `which` can be used to select which variables to plot.

```{r}
plot(sol, which = c('Y','Z'))
```

The optional argument `method` selects the numerical integrator among several. The default integrator is `lsoda`.
LSODA is an integrator for solving stiff and non-stiff systems of ordinary differential equations. It was written in FORTRAN by Linda Petzold and Alan Hindmarsh. It can solve systems with dense or banded Jacobian when the problem is stiff. LSODA autamatically selects between a stiff solver (backward differencing formulae, BDF) and a non-stiff solver (Adams). This and other integrators can be used directly by their names. LSODA should be the first integrator to try when solving a new system. 

```{r}
sol.lsoda <- lsoda(y = vars, times = tim, func = Lorenz, parms = params)
plot(sol.lsoda, which = 'X')
```

# Stiff equations

The stiffness of an equation has many definitions. Mathematically, a system 

$$ \frac{dx}{dt} = F(x)$$

is stiff if some of the eigenvalues $\lambda_1, \lambda_2$, of the Jacobian matrix of $F$ have much different orders of magnitudes 

$$|\lambda_1| \ll |\lambda_2|.$$

Numerically, stiff systems are best solved with *implicit* methods. Solving stiff equations with explicit methods while maintaining a good accuracy is very slow. Let's use a prototype stiff system, the van der Pol oscillator.
The van der Pol oscillator is a second-order equation 

$$\frac{d^2x}{dt} = \mu \Bigl( (1-x^2) \frac{dx}{dt} - x \Bigr) + F(t).$$ 
The term $F(t)$  is a force term that we set to zero. The van der Pol oscilator is a nonlinear oscillator. Superficially, it looks like a damped harmonic oscillator, except that the damping force is *negative* when $x$ is small, so that instead of losing energy, the system gains in energy. The system is stiff when the paramter value is large; we will set $\mu = 1e6$.

Solve the van der Pol equation with an explicit integrator first, for times between 0 and 6.3, with $y(0)=0.0$ and
$x(0) = 2.0$. First, express the equation as a system of two first-order equations

$$\frac{dy}{dt} = \mu \Bigl( (1-x^2) y - x \Bigr),$$
$$\frac{dx}{dt} = y.$$

We use the Runge-Kutta method `rk45ck`. This is an explicit, variable step method of the fourth and fifth orders. 


```{r}
vanDerPol <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx <-  y
    dy <-  mu*((1-x*x)*y - x)
    list(c(dx, dy))
  })
}
vars <- c(x = 2.0, y = 0.0)
params <- c(mu = 1e6)
tim <- seq(0, 6.3, by = 0.01)

start_time <- Sys.time()
sol.rk <- rk(y = vars, times = tim, func = vanDerPol, parms = params, method = 'rk45ck', rtol = 1e-6, atol = 1e-6)
end_tiem <- Sys.time()
end_time - start_time
```

So it doesn't work at all. Let's try with a mininal time step, to speed things up, at the expense of accuracy,  of course. The maximal number of steps also need to be increased. 

```{r}
tim <- seq(0, 6.3, by = 0.01)
start_time <- Sys.time()
sol.rk45 <- rk(y = vars, times = tim, func = vanDerPol, parms = params, method = 'rk45ck', hmin = 1e-7, rtol = 1e-6, atol = 1e-6, maxsteps = 1e8)
end_time <- Sys.time()
end_time - start_time
```

```{r}
plot(sol.rk45, which = 'y')
```

Now, we use a backward-differentiation formula `BDF` to solve the system

```{r}
tim <- seq(0, 6.3, by = 0.01)
start_time <- Sys.time()
sol.bdf <- ode(y = vars, times = tim, func = vanDerPol, parms = params, method = 'bdf', rtol = 1e-6, atol = 1e-6)
end_time <- Sys.time()
end_time - start_time
```

```{r}
plot(sol.bdf, which = 'y')
```

Plot both solutions on the same axes

```{r}
plot(sol.rk45, which = 'y', ylim = c(-4,4))
lines(sol.bdf[,'time'],sol.bdf[,'y'], col = 'blue')
```

The implicit BDF method is *much* faster than the explicit RK method. 

Non-stiff systems are much easier to solve, and in most cases, explicit integrators are faster. If low accuracy is acceptable (for a non-stiff system), lower-order integrators such as `ode23`/`rk23bs`. Using explicit methods for stiff system leads either very small step sizes or to blow up of the solutions. In the worst case, the step size becomes so small it is smaller than espilon machine, so that t+dt = t, and the integrator becomes stuck. Blow up of solutions occurs when the time step is too large and solutions go out of their range of definition. For example, solution may become negative if the derivative is negative. Nonlinear rhs functions may not be appropriately defined for negative solutions, take $1/(1+x)$, which is bounded for $x>0$ but unbounded near $x=-1$.

# Non-stiff equations

Subtrate Producer Consumer model (3D) with a nonlinear function passed as argument to the model

```{r}
# define the ODE equations
SPCmodel <- function(t, x, parms, input)  {
  with(as.list(c(parms, x)), {
    s0 <- input(t)                      # substrate production is time-dependent
    dS <- s0 - b*S*P + g*C              # substrate
    dP <- c*S*P  - d*C*P                # producer
    dC <- e*P*C  - f*C                  # consumer
    res <- c(dS, dP, dC)
    list(res)
  })
}

# set the parameter names and values
pars <- c(b = 0.001, c = 0.1, d = 0.1, e = 0.1, f = 0.1, g = 0.0)

# set time points
times <- seq(0, 200, by = 1.0)

# define the substrate production function
subs <- data.frame(times = times,
                     prod = rep(0, length(times)))

# = 0.2 between t = 10 and t = 11
subs$prod[subs$times >= 10 & subs$times <= 11] <- 0.2

# function interpolating the vector subs$prod
subsprod <- approxfun(subs$times, subs$prod, rule = 2)

# set initial conditions
init_conds <- c(S = 1, P = 1, C = 1)

# solve the ODE system
sol <- ode(y = init_conds, times = times, func = SPCmodel, parms = pars, input = subsprod)

# plot solutions
plot(sol)
```

# Error control

The integrators `lsoda` and other routines of the `deSolve` package have *error control*. 
The numerical methods compute approximate solution at the next time step at two different order of accuracy. The difference between the two solutions is a good approximation of the
size of the true error. This error is local, i.e. it is calculated assuming that the 
solution at time $t$ is exact. Error is controled by setting *tolerance* options; if
the error is larger than the tolerance, the integrator will reject the new solution and pick a smaller time step. If the error is less than the tolerance, the integrator will accept
the new solution and increase the next time step to speed up the computation.

Two errors are computed, the *relative error* and the *absolute error*. In `lsoda` the solution is accepted if  

$$err = \max |y_{high} - y_{low}| < \mathrm{rtol} \times |y| + \mathrm{atol}$$
The option `rtol` and `atol` are the *relative tolerance* and the *absolute tolerance* respectively. The smaller the tolerance the more stringent the condition for accpeting the solution. The relative tolerance control the precision of the solution when far away from zero. A relative tolerance of 1e-6 specifies 6 digits of accuracy. The absolute tolerance 
specifies the accuracy when solutions are close to zero. An absolute tolerance of 1e-6 states that any solution smaller than 1e-6 should be regarded as effectively zero. 
for system of equations, the error is set to the max of the errors on each coefficients.

*Remark* The error control expression involves the sum of the tolerances. This is equivalent to an OR: if the error is smaller than either tolerances, the solution will be accepted. Therefore, it is absolutely possible to set rtol or atol to zero (but not both).

In summary, rtol controls the solution away from zero and atol controls the solution close to zero.

# Non-negative solutions

It can often inconvenient to have negative a solution to an initial value problem to ODEs that theoretically should have only positive solution. 
Problems in biology and chemistry involving concentrations or population densities should admit only positive solutions, the RHS of the ODEs are often not
defined outside the range of positive states. The initial value problem (IVP) for a system of ODE is

$$\frac{dy}{dt} = f(y,t), \; y(0) = y_0.$$
For example, the IVP

$$\frac{dy}{dt} = -|y|, \; y(0) = 1.$$
The solution is $y(t) = \exp(-t)$ is positive, and decays to zero. Suppose that the numerical solution reaches a value $y^* < 0$ at a time $t^*$. 
Then the solution, starting at $t^*$, is $u(t) = y^*\exp(t-t^*)$. This solution is decreasing (towards $-\infty$). Thus the IVP is *unstable*, in the sense that small 
pertubation of the solution lead to large deviation from the true solution. This IVP will be the first test problem, solved on $[0,40]$

```{r}
# define the ODE equations
nonStiffTest1 <- function(t, x, parms)  {
  with(as.list(c(parms, x)), {
    dy <- - abs(y)
    list(dy)
  })
}

# set time points
times <- seq(0, 40)

# set initial conditions
y0 <- c(y = 1)

# solve the ODE system
sol <- ode(y = y0, times = times, parms = NULL, func = nonStiffTest1)

# plot solutions
plot(sol)
```

The second test is the Robertson problem. It is a semi-stable chemical reaction system

$$
\frac{dy_1}{dt}  = -0.04 y_1 + 10^4 y_2 y_3, \\
\frac{dy_2}{dt}  =  0.04  y_1 - 10^4 y_2 y_3 - 3 \times 10^7 y_2^2, \\
\frac{dy_3}{dt}  =  3 \times 10^7 y_2^2,
$$

with initial conditions $(1,0,0)^t$ with time interval $[0,4\times 10^{11}$. This is a *stiff* problem that need to be solved over a very large time interval.

```{r}
Robertson <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dy1 <-  -a*y1 + b*y2*y3
    dy2 <-   a*y1 - b*y2*y3 - c*y2*y2
    dy3 <-   c*y2*y2
    list(c(dy1, dy2, dy3))
  })
}

# set the parameter names and values
pars <- c(a = 0.04, b = 1e4, c = 3e7)

# set time points
times <- seq(0, 4e11, length.out = 50)

# set initial conditions
y0 <- c(y1 = 1, y2 = 0, y3 = 0)

# solve the ODE system
sol <- ode(y = y0, times = times, parms = pars, func = Robertson, rtol = 1e-6, atol = 1e-6)
```

```{r}
# plot solutions
plot(sol, which = 'y1')
```

The Lotka-Volterra model for predator-prey dynamics admits oscillatory solutions 
that periodically come arbitrarily close to zero. The IVP is 

$$
\frac{dy_1}{dt}  =  0.5 y_1 \Bigl( 1 - \frac{y_1}{20} \Bigr) - 0.1 y_1 y_2, \\
\frac{dy_2}{dt}  =  0.01  y_1 y_2 - 0.001 y_2,
$$
with initial conditions $(25,5)^t$ over the time interval $[0,870]$. If at some time $t^*$ the numerical solution  $y_1^* < 0$, and the solution is then reset at $(0,y_2^*)^t$, the
solution has $y_1(t) = 0$ for all $t > t^*$.

```{r}
LotkaVolterra <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dy1 <-  a * y1 * ( 1 - y1/b ) - c * y1 * y2
    dy2 <-  d * y1*y2 - e * y2
    list(c(dy1, dy2))
  })
}

# set the parameter names and values
pars <- c(a = 0.5, b = 20, c = 0.1, d = 0.01, e = 0.001)

# set time points
times <- seq(0,1000,length.out = 2)

# set initial conditions
y0 <- c(y1 = 25, y2 = 5)

# solve the ODE system
sol1 <- ode(y = y0, times = c(0,1000), parms = pars, func = LotkaVolterra, atol = 1e-12)
sol2 <- ode(y = y0, times = seq(0,1000,length.out = 1000), parms = pars, func = LotkaVolterra, atol = 1e-12)
```

```{r}
# plot solutions
plot(sol1[,1],sol1[,'y1'])
lines(sol2[,1],sol2[,'y1'])
```

The IVP

$$\frac{dy}{dt} = \sqrt{1-y^2}, \; y(0) = 0.$$
The solution $y(t) = \sin(t)$ increases to 1 at $t = \pi /2$, and the rhs is not defined
for $y>1$ (it is not Lipschitz). The solution it not unique after $t^* = \pi/2$. 

```{r}
# define the ODE equations
Ball <- function(t, x, parms)  {
  with(as.list(c(parms, x)), {
    dy <- sqrt(1-y*y)
    list(dy)
  })
}

# set time points
times <- seq(0, 3, length.out = 100)

# set initial conditions
y0 <- c(y = 0)

# solve the ODE system
sol <- ode(y = y0, times = times, parms = NULL, func = Ball, method = 'bdf')

# plot solutions
plot(sol)
```


# Example of a Runge-Kutta pair of order 2 and 3: Bogacki-Shampine RK3(2)

The numerical scheme is a pair of Runge-Kutta formulas, one of order 3 and one of order 2. For the differential
equation

$$\frac{dy(t)}{dt} = f(t,y(t)), \; a \leq t \leq b,$$
$$y(a) = y_a.$$

The two formulas are used to advance the solution by a time step $h$, 
given a solution $y_n$ at time step $t_n$. The difference between the third order $\hat y_{n+1}$ and the second order step $y_{n+1}$ provides an estimate of the error $E$ to the true solution. 
Each formula requires several evaluations of the function $f$. Evaluation is done in $S$ stages:

$$y_{n+1} = \hat y_{n} + h \sum_{i=1}^S b_i k_i.$$

Stages $k_i$ are calculated iteratively: $k_i$ depends only on the previous $k_j$, $j=1,...,i-1$. The first stage is always

$$k_1 = f(t_n, y_n).$$

The other stages are

$$k_i = f \bigl( x_n + c_i h, y_n + h \sum_{j=1}^{i-1} a_{i,j} k_j \bigr),$$
and

$$c_i = \sum_{j=1}^{i-1} a_{i,j}$$

To reduce the total number of function evaluations, a strategy consists in sharing the stages $k_i$ between the two formulas, this is called embedding. Third order Runge-Kutta formulas need at least three stages, and second order formulas at least two stages. Each stage require one evalutation of the function $f$. The first step $k_1$ is always shared. If the second stage of the third order formula is also used in the second order formula, the total number of function evaluations is 3: we get the second order approximation almost for free. 

However, Bogacki and Shampine reason in a different way. They want to match as closely as possible the regions of stability of the two formulas. All two-stage, second order formulas have the same region of stability and have all three-stages, third order formulas. The only way to match the regions of stability is to increase the number of stages. Increasing the number of stages leaves quite a lot of room for choosing the coefficient, and Bogacki and Shampine took a four stages for the second order formula. This seems rather superfluous, but they also reuse the first three stages in the second-order formula. The fourth stage of the second order formula is 

$$k_4 = f(t_n + h, \hat y_n + h \sum_{j=1}^3 \hat b_j k_j).$$
Notice that $k_4$ is excatly the value of the function $f$ at $t_{n+1} = t_n + h$ and $\hat y_{n+1} = \hat y_n + h \sum_{j=1}^3 \hat b_j k_j$. This is the first stage $k_1$ of the next time step! Once we realize that most steps are successful, we see that the four-stage, second order formula come at little extra computational cost: $k_4$ is carried over to the next time step. This approach is call FSAL: First Same As Last.

The Butcher tableau for the pair of formula is

    0   | 
    1/2 | 1/2
    3/4 |   0   3/4
      1 | 2/9   1/3 4/9
    --- - ----  --- --- ---
        | 2/9   1/3 4/9   0
        | 7/24  1/4 1/3 1/8     


 # Event location
 
 Event location is an option to track the solution of a system of ODE while it is being
 solved for. This useful, but not limited to cases when the dynamic variable are
 suddenly changed, i.e. when the dynamical variables are discontinuous.
 Normal integration routines cannot deal very well with discontinuities. The `events` option specifies when the events or discontiuities should occur, so that the integrator can deal with it. Events are imposed either by a `data.frame` that specifies the times at which the variables have jumps, or by an `event function` that monitor the state of the variable. A `root function` will trigger an event when the function takes the value zero.
 
The event function has the same syntax as the ode funcion: `function(t,y,params, ...)`.
The `data.frame` should have the folowing columns, in that order: `var` the state variable
name or number affected by the event, `time` the time at which the event takes place,
`value` the magnitude of the event, `method` how to set the new values, one of "replace", "add", "multiply".

The row: 

    "v1" 10 2 "add"
    
would **add** 2 to the variable `v1` at time 10.

The routines `lsoda` (and some others) have root-finding abilities. If a root function is specified, the solver will stop at the first root. To monitor the location of many events, an event function must be set that will leave the state variables unaltered.

## Example 1 Event in a data.frame

```{r}
# ODE function
myOdefunc <- function(t, var, pars) 
{
  with(as.list(c(pars, var)), {
    dx <- 0
    dy <- - a * y
    list(c(dx,dy))
  })
}

y0 <- c(x = 0, y = 1) # initial condition (x,y must be the same x,y as inside myOdeFunc)
pars <- c(a = 1.2) # parameters (a must be the same a as inisde myOdeFunc)
times <- seq(0,10,by=0.1)

# define a data.frame called eventDataFrame
eventDataFrame <- data.frame(var = c("x", "y", "y", "x"),
                             time = c(1,1,5,9),
                             value = c(1,2,3,4),
                             method = c("add", "mult", "rep", "add"))

sol <- ode(y = y0, func = myOdefunc, times = times, parms = pars, 
           event = list(data = eventDataFrame))

plot(sol)
```
 

## Example 2 Event in a function

```{r}
myOdefunc <- function(t, var, pars) 
{
  with(as.list(c(pars, var)), {
    dx <- 0
    dy <- - a * y
    list(c(dx,dy))
  })
}



# events: add 1 to x, multiply y with a random number
eventFun <- function(t, var, parms){
  with (as.list(var),{
    x <- x + 1
    y <- 5 * runif(1)
    return(c(x, y))
  })
}

y0 <- c(x = 1, y = 2)
times <- seq(0, 10, by = 0.1)
pars <- c(a = 1.2)

sol <- ode(func = myOdefunc, y = y0, times = times, parms = pars,
           events = list(func = eventFun, time = c(1:9, 2.2, 2.4)) )
plot(sol, type = "l")
```

# Event triggered by a root function

```{r}
## ODE: simple first-order decay
firstOrderDecay <- function(t, var, pars) {
  with(as.list(c(pars, var)), {
    dy <- - a * y
    list(c(dy))
  })
}

## event triggered if state variable y = 0.5
rootFun <- function (t, var, pars) {
  with(as.list(c(pars, var)), {
    y - 0.5
  })
}

## sets state variable = 1
eventFun <- function(t, y, pars) {
    with (as.list(var),{
    y <- 1
    return(y)
  })
}

y0 <- c(y = 2)
times <- seq(0, 100, 0.1)
pars <- c(a = 0.1)

## uses ode to solve; root = TRUE specifies that the event is
## triggered by a root.
sol <- ode(times = times, y = y0, func = firstOrderDecay, parms = pars,
           events = list(func = eventFun, root = TRUE),
           rootfun = rootFun)

plot(sol, type = "l")

## time of the root:
times_root <- attributes(sol)$troot
points(times_root, rep(0.5, length(times_root)))
```