deSolve is a package to solve initial value problems of several types of differential equations.

Refs:

library(deSolve)

Ordinary Differential Equations (ODE)

An ordinary differential equation is an equation containing a function of one independent variable and its derivatives.

An eg solving the Lorentz system:

\[\frac{dX}{dt} = aX+Y+Z\] \[\frac{dY}{dt} = b(Y-Z)\] \[\frac{dZ}{dt} = -XY+cY-Z\]

with parameters \(a=-8/3, b=-10, c=28\) and initial position at \(X(0)=Y(0)=Z(0)=1\).

First we make the model specification, ie, the parameter values, initial position and a function that implement the model equations regarding their rate of change:

parameters    <- c(a = -8/3, b = -10, c =  28)
initial.state <- c(X = 1, Y = 1, Z = 1)

Lorenz<-function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    # rate of change
    dX <- a*X + Y*Z
    dY <- b * (Y-Z)
    dZ <- -X*Y + c*Y - Z
    
    list(c(dX, dY, dZ)) # return the rate of change
  })
}

Then we apply the model. For that we need to know what are the timestamps used:

times <- seq(0, 100, by = 0.01)

Finally we apply all into the ODE solver:

out <- ode(y = initial.state, times = times, func = Lorenz, parms = parameters)

Visualizing the results:

head(out)
##      time      X     Y     Z
## [1,] 0.00 1.0000 1.000 1.000
## [2,] 0.01 0.9849 1.013 1.260
## [3,] 0.02 0.9731 1.049 1.524
## [4,] 0.03 0.9652 1.107 1.798
## [5,] 0.04 0.9617 1.187 2.089
## [6,] 0.05 0.9638 1.288 2.400
summary(out)
##                 X         Y         Z
## Min.        0.962   -18.200   -24.600
## 1st Qu.    17.200    -6.660    -6.220
## Median     23.200    -0.524    -0.485
## Mean       23.800    -0.379    -0.388
## 3rd Qu.    30.200     5.150     4.660
## Max.       47.800    19.600    27.200
## N       10001.000 10001.000 10001.000
## sd          8.408     7.960     8.999
par(oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "X"], out[, "Z"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)

plot of chunk unnamed-chunk-5

Plotting in 3D:

library(plot3D)

points3D(out[, "X"], out[, "Y"], out[, "Z"], pch='.', colkey=F, colvar=out[, "Y"]) # color using Y values

plot of chunk unnamed-chunk-6

# zooming
out.zoom <- subset(out, select=c("X","Y","Z"), subset = Y < 8 & X > 10 & X < 30)
points3D(out.zoom[, "X"], out.zoom[, "Y"], out.zoom[, "Z"], 
         pch='.', colkey=F, colvar=out.zoom[, "Y"]) # color using Y values

plot of chunk unnamed-chunk-6

Second-Order ODE

This 2nd order ODE

\[y'' - \mu(1-y^2)y' + y=0\]

can be transformed into the following system of first order ODEs:

\[y_1' = y_2\] \[y_2'= \mu(1-y_1^2)y_2-y_1\]

vdpol <- function (t, y, mu) {
  list(c(y[2],
         mu * (1 - y[1]^2) * y[2] - y[1])
      )
}

y.init <- c(y1 = 2, y2 = 0)
out <- ode(y = y.init, func = vdpol, times = seq(0, 30, 0.01), parms = 1)
head(out)
##      time    y1       y2
## [1,] 0.00 2.000  0.00000
## [2,] 0.01 2.000 -0.01970
## [3,] 0.02 2.000 -0.03882
## [4,] 0.03 1.999 -0.05737
## [5,] 0.04 1.998 -0.07537
## [6,] 0.05 1.998 -0.09283
plot(out, xlab = "time", ylab = "-")

plot of chunk unnamed-chunk-7

plot(out[, "y1"], out[, "y2"], pch = ".")

plot of chunk unnamed-chunk-7

Partial Differential Equations (PDE)

A partial differential equation is a differential equation that contains unknown multivariable functions and their partial derivatives.

The functions to use are ode.1D, ode.2D, and ode.3D for problems in these respective dimensions.

A 1-D eg describing a model where aphids (a pest insect) slowly diffuse and grow on a row of plants:

\[\frac{\partial N}{\partial t} = - \frac{\partial Flux}{\partial x} + rN\]

\[Flux = -D \frac{\partial N}{\partial x}\]

with boundaries \(N_{x=0} = N_{x=60} = 0\) and initial state \(N_{30}=1\) and everything else at zero (ie, there’s one aphid at the 30th plant box).

The package uses the method of lines approach to split the spatial domain into a number of boxes and discretize \(dN/dt\).

parameters = list(D=0.3,       # diffusion rate; m^2/day
                  r=0.01,      # net growth rate; day^-1
                  numboxes=60, # number of boxes
                  delx=1)      # thickness of each box; m
  
Aphid <- function(t, N, parameters) {
  with(parameters,{

    deltax  <- c(0.5, rep(1, numboxes - 1), 0.5)
    Flux    <- -D * diff(c(0, N, 0)) / deltax
    dN      <- -diff(Flux) / delx + N * r
     
    list(dN)
  })
}

# initial condition
N <- rep(0, times =  parameters$numboxes)
N[30:31] <- 1
state <- c(N = N) # initialise state variables

# let's run for 300 days
times <- seq(0, 300, by = 1)
out <- ode.1D(state, times, Aphid, parms = parameters, nspec = 1, names = "Aphid")

head(out[,1:5])
##      time        N1        N2        N3        N4
## [1,]    0 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## [2,]    1 1.667e-55 9.555e-52 2.555e-48 4.943e-45
## [3,]    2 3.631e-41 4.865e-39 5.394e-37 5.054e-35
## [4,]    3 2.051e-34 9.208e-33 3.723e-31 1.391e-29
## [5,]    4 1.307e-30 3.719e-29 9.635e-28 2.361e-26
## [6,]    5 6.839e-28 1.465e-26 2.860e-25 5.334e-24
summary(out)
##             Aphid
## Min.    0.000e+00
## 1st Qu. 1.010e-02
## Median  1.120e-01
## Mean    2.080e-01
## 3rd Qu. 3.090e-01
## Max.    1.190e+00
## N       1.806e+04
## sd      2.501e-01
image(out, method = "filled.contour", 
      grid = seq(from = 0.5, by = parameters$delx, length.out = parameters$numboxes),
      xlab = "time, days", ylab = "Distance on plant, meters",
      main = "Aphid density on a row of plants")

plot of chunk unnamed-chunk-8

Differential algebraic equations (DAE)

A differential-algebraic equation is an equation involving an unknown function and its derivatives.

Eg:

\[\frac{dy_1}{dt} = y_2 - y_1\] \[y_1y_2 = t\]

To solve it, first rewrite the equation in their residual form:

\[\frac{dy_1}{dt} + y_1 - y_2 = 0\] \[y_1y_2 - t = 0\]

f <- function(t, y, dy, parameters) {
  res1 <- dy[1] + y[1] - y[2]
  res2 <- y[2] * y[1] - t

  list(c(res1, res2))
}

yini  <- c(2, 0) # initial conditions
dyini <- c(1, 0)

times <- seq(0, 20, 0.1)

out <- daspk(y = yini, dy = dyini, times = times, res = f, parms = 0)
matplot(out[,1], out[,2:3], type = "l", lwd = 2, col=c("red","blue"), lty=1, 
        main = "DAE", xlab = "time", ylab = "ys")
legend("bottomright",legend=c("y1","y2"), col=c("red","blue"), lty=1, lwd=2)

plot of chunk unnamed-chunk-9