deSolve
is a package to solve initial value problems of several types of differential equations.
Refs:
library(deSolve)
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
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 = ".")
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")
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)