Integrating the General Lotka-Volterra (GLV) equations numerically in R
The R codes in this document are provided here as a script.
Numerical integration of the GLV vectorfield
The GLV equations for the dynamics of a community consisting of \(n\) populations, with densities \(x_i\), \(i=1,\dots,n\), are: \[\frac{dx_i}{dt} = x_i\,\left(r_i + \sum_j a_{ij}\,x_j\right).\]
The lv_flow function below integrates the GLV equations
numerically. To increase numerical stability, it does so in terms of
\(y=\log(x)\) instead of the original
variables \(x\). [Recall that we use
\(\log\) to refer to the
natural logarithm.]
In terms of \(y\), the GLV equations are: \[\frac{dy_i}{dt} = r_i + \sum_j a_{ij}\,\exp(y_j).\] Using matrix notation, this is \[\frac{dy}{dt} = r + A\cdot\exp(y),\] where the \(\cdot\) denotes matrix multiplication and the exponentiation is element-wise.
The numerical integration is done using the ode function
from the deSolve package. The lsoda option
tells ode to use the Livermore Solver of Ordinary
Differential equations, with Automatic method switching.
lv_flow <- function (t, x, r, A) {
x <- as.matrix(x)
ode(
y=log(x), # we transform from the natural to the log scale
times=c(0,t),
func=function (t, y, ...) {
dim(y) <- dim(x)
list(
r + A %*% exp(y),
NULL
)
},
method="lsoda"
) |>
as.matrix() -> out
## Transform back onto the natural scale.
out[,-1] <- exp(out[,-1])
if (is.null(rownames(x))) {
if (ncol(x)==1) {
colnames(out) <- c("time",paste0("x_",seq_len(nrow(x))))
} else {
colnames(out) <- c(
"time",
paste0(
"x_",rep(seq_len(nrow(x)),times=ncol(x)),
"_",rep(seq_len(ncol(x)),each=nrow(x))
)
)
}
}
as_tibble(out[-1,]) # a tibble is like a data frame, but better in some ways
}Classical predator prey model
With Malthusian parameters \(r=\begin{bmatrix} 1 &-1 \\ \end{bmatrix}\) and community matrix \[A=\begin{bmatrix} 0 &-1 \\1 &0 \\ \end{bmatrix},\] we have the classical Lotka-Volterra predator-prey model, with its constant of the motion and neutrally stable continuum of periodic orbits.
Code
Plot 2
These equations admit a first integral. Accordingly, the orbits lie on the level sets of this function.
lv_flow(
t=seq(0,100,by=0.01),
x=matrix(
c(1,0.1,
1,1.1,
1,2),
2,3),
r=c(1,-1),
A=matrix(c(0,-1,1,0),2,2,byrow=TRUE)
) |>
pivot_longer(-time) |>
separate(name,into=c(NA,"var","ic")) |>
pivot_wider(names_from="var") |>
rename(x1="1",x2="2") |>
ggplot(aes(x=x1,y=x2,group=ic))+
geom_path()+
labs(x=expression(x[1]),y=expression(x[2]))
Chaos
With three species and \[r=\begin{bmatrix} 1 \\1 \\-1 \\ \end{bmatrix}, \qquad A=\begin{bmatrix} -1 &-1 &-10 \\-1.5 &-1 &-1 \\5 &0.5 &-0.01 \\ \end{bmatrix},\] we obtain chaotic trajectories.
Code
Plot
plot_grid(
dat |>
ggplot(aes(x=x_1,y=x_2))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[2])),
dat |>
ggplot(aes(x=x_1,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[3])),
dat |>
ggplot(aes(x=x_2,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[2]),y=expression(x[3])),
dat |>
filter(time>5000,time<10000) |>
pivot_longer(-time) |>
separate(name,into=c(NA,"index")) |>
ggplot(aes(x=time,y=value,group=index))+
geom_line()+
facet_grid(
index~.,
scales="free_y",
labeller=label_bquote(rows=x[.(index)])
)+
labs(y="")
)
Sensitive dependence on initial conditions
At these parameters, the system is chaotic in the sense that it displays sensitive dependence on initial conditions. This means that trajectories starting from different initial conditions—no matter how close they are—diverge exponentially quickly. The result is that there is a fundamental, mathematical limit to the predictability of this system, even though it is deterministic.
Code
Plot
dat |>
pivot_longer(-time) |>
separate(name,into=c(NA,"index","ic")) -> dat
dat |>
ggplot(aes(x=time,y=value,color=ic))+
geom_line()+
scale_color_manual(values=c(`1`="red",`2`="blue"))+
facet_grid(index~.,scales="free_y",
labeller=label_bquote(rows=x[.(index)])
)+
labs(y="",color="initial condition")+
theme_bw()+
theme(legend.position="top") -> pl
pl %+% filter(dat,time<=1000) -> pl1
pl %+% filter(dat,time>7000,time<=8000) -> pl2
plot_grid(pl1,pl2)
Cyclic competition
With a circulant community matrix, we can obtain intransitive competitive loops.
Circulant matrices
The following function constructs a circulant matrix.
circulant <- function (...) {
x <- c(...)
n <- length(x)
a <- array(data=x,dim=c(2*n-1,n))
a[seq_len(n),]
}For example:
\[\begin{bmatrix} 1 &4 &3 &2 \\2 &1 &4 &3 \\3 &2 &1 &4 \\4 &3 &2 &1 \\ \end{bmatrix}\]Triangle plot
array(
c(
-sqrt(1/2),-sqrt(1/6),
sqrt(1/2) ,-sqrt(1/6),
0 , sqrt(2/3)
),
dim=c(2,3)
) -> .Qmat
triangularize <- function (data, x, y, z, scale = 1) {
x <- eval(substitute(x),data)
y <- eval(substitute(y),data)
z <- eval(substitute(z),data)
array(
c(x,y,z),
dim=c(length(x),3)
) -> coords
coords |>
sweep(1L,apply(coords,1L,sum)/scale,"/") |>
tcrossprod(y=.Qmat) -> coords
tibble(x=coords[,1L],y=coords[,2L])
}
tibble(x=c(1,0,0,1),y=c(0,1,0,0),z=c(0,0,1,0)) |>
triangularize(x,y,z) -> .vertices
tibble(x=c(1,0,0),y=c(0,1,0),z=c(0,0,1)) |>
triangularize(x,y,z,scale=1.03) -> .label_coordsStable interior rest point
With \(r=1\) and the circulant community matrix \[A=\begin{bmatrix} -1 &-0.1 &-0.4 \\-0.4 &-1 &-0.1 \\-0.1 &-0.4 &-1 \\ \end{bmatrix},\] we have a stable rest point in the interior of \(\mathbb{R}_+^3\).
Plot
plot_grid(
dat |>
ggplot(aes(x=x_1,y=x_2))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[2])),
dat |>
ggplot(aes(x=x_1,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[3])),
dat |>
ggplot(aes(x=x_2,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[2]),y=expression(x[3])),
dat |>
pivot_longer(-time) |>
separate(name,into=c(NA,"index")) |>
ggplot(aes(x=time,y=value,group=index))+
geom_line()+
facet_grid(index~.,scales="free_y",
labeller=label_bquote(rows=x[.(index)])
)+
labs(y="")
)
Triangle plot
dat |>
triangularize(x_1,x_2,x_3) |>
ggplot(aes(x=x,y=y))+
geom_path()+
geom_path(data=.vertices)+
geom_text(
data=.label_coords,
mapping=aes(label=c("R","S","P")),
size=13,
size.unit="pt"
)+
theme_minimal()+
theme(
axis.text=element_blank(),
axis.title=element_blank(),
panel.grid=element_blank()
)
Stable manifold
When \[A=\begin{bmatrix} -1 & -\alpha & -\beta \\ -\beta & -1 & -\alpha \\ -\alpha & -\beta & -1 \end{bmatrix},\] the unique interior rest point is a saddle whenever \(\alpha+\beta>2\). If we are careful, we can put the initial conditions onto its stable manifold.
Code
Plot
plot_grid(
dat |>
ggplot(aes(x=x_1,y=x_2))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[2])),
dat |>
ggplot(aes(x=x_1,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[3])),
dat |>
ggplot(aes(x=x_2,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[2]),y=expression(x[3])),
dat |>
pivot_longer(-time) |>
separate(name,into=c(NA,"index")) |>
ggplot(aes(x=time,y=value,group=index))+
geom_line()+
facet_grid(index~.,scales="free_y",
labeller=label_bquote(rows=x[.(index)])
)+
labs(y="")
)
Triangle plot
dat |>
triangularize(x_1,x_2,x_3) |>
ggplot(aes(x=x,y=y))+
geom_path()+
geom_path(data=.vertices)+
geom_text(
data=.label_coords,
mapping=aes(label=c("R","S","P")),
size=13,
size.unit="pt"
)+
theme_minimal()+
theme(
axis.text=element_blank(),
axis.title=element_blank(),
panel.grid=element_blank()
)
Stable heteroclinic cycle
Just off the stable manifold, we approach saddle point, lingering there for some time before leaving along its unstable manifold. As time goes on, we approach the stable heteroclinic cycle, with its characteristic dynamics.
Code
Plot
plot_grid(
dat |>
ggplot(aes(x=x_1,y=x_2))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[2])),
dat |>
ggplot(aes(x=x_1,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[1]),y=expression(x[3])),
dat |>
ggplot(aes(x=x_2,y=x_3))+
geom_point(size=0.1)+
labs(x=expression(x[2]),y=expression(x[3])),
dat |>
pivot_longer(-time) |>
separate(name,into=c(NA,"index")) |>
ggplot(aes(x=time,y=value,group=index))+
geom_line()+
facet_grid(index~.,scales="free_y",
labeller=label_bquote(rows=x[.(index)])
)+
labs(y="")
)
Triangle plot
dat |>
triangularize(x_1,x_2,x_3) |>
ggplot(aes(x=x,y=y))+
geom_path()+
geom_path(data=.vertices)+
geom_text(
data=.label_coords,
label=c(expression(x[1]),expression(x[2]),expression(x[3])),
size=16,
size.unit="pt"
)+
theme_minimal()+
theme(
axis.text=element_blank(),
axis.title=element_blank(),
panel.grid=element_blank()
)