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.0000000 1.000000 1.000000
## [2,] 0.01 0.9848912 1.012567 1.259918
## [3,] 0.02 0.9731148 1.048823 1.523999
## [4,] 0.03 0.9651593 1.107207 1.798314
## [5,] 0.04 0.9617377 1.186866 2.088545
## [6,] 0.05 0.9638068 1.287555 2.400161
summary(out)
##                    X             Y             Z
## Min.    9.617377e-01   -18.4060950   -24.9372365
## 1st Qu. 1.762149e+01    -6.8491758    -6.4923103
## Median  2.348657e+01    -0.5046839    -0.6147006
## Mean    2.388163e+01    -0.4399923    -0.4512985
## 3rd Qu. 3.014746e+01     5.2467173     4.7473544
## Max.    4.783396e+01    19.5550956    27.1835768
## N       1.000100e+04 10001.0000000 10001.0000000
## sd      8.300157e+00     7.9778965     8.9761641
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)

Plotting in 3D:

library(plot3D)

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

# 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

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.000000  0.00000000
## [2,] 0.01 1.999901 -0.01970323
## [3,] 0.02 1.999608 -0.03882242
## [4,] 0.03 1.999127 -0.05737271
## [5,] 0.04 1.998462 -0.07537148
## [6,] 0.05 1.997621 -0.09283422
plot(out, xlab = "time", ylab = "-")

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

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.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [2,]    1 1.667194e-55 9.555028e-52 2.555091e-48 4.943131e-45
## [3,]    2 3.630860e-41 4.865105e-39 5.394287e-37 5.053775e-35
## [4,]    3 2.051210e-34 9.207997e-33 3.722714e-31 1.390691e-29
## [5,]    4 1.307456e-30 3.718598e-29 9.635350e-28 2.360716e-26
## [6,]    5 6.839152e-28 1.465288e-26 2.860056e-25 5.334391e-24
summary(out)
##                Aphid
## Min.    0.000000e+00
## 1st Qu. 1.007562e-02
## Median  1.119879e-01
## Mean    2.075457e-01
## 3rd Qu. 3.087324e-01
## Max.    1.193554e+00
## N       1.806000e+04
## sd      2.500619e-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")

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)