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)

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.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(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.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")

## 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)