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).\]

Transformation for numerical stability

The codes below integrate the GLV equations numerically. To increase numerical stability, they do 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.

Computing the flow

The following function integrates the Lotka-Volterra vectorfield. Specifically, if x contains a given set of initial states (at time 0), r and A contain the vector \(r\) and matrix \(A\) as above, and \(t\) is a vector of times at which the flow is to be evaluated, then lv_flow(t,x,r,A) returns the evaluated flow.

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), # cf equation for dy/dt above
        NULL
      )
    },
    method="lsoda"
  ) |>
    as.matrix() -> out
  ## Transform back onto the natural scale.
  out[,-1] <- exp(out[,-1])
  cleanup(out,x)
}

Integrator details

Internally, the numerical integration is done using the ode function from the deSolve package. The call sequence for ode is

ode(y, times, func, parms, method, ...)

Here, y is the initial condition (i.e., the state at time 0), times are the times at which solutions are desired, func is a function that returns the vectorfield at a point, parms is a vector that contains any constant parameters upon which the vectorfield depends, and method is a string that specifies which of several integration algorithms is to be used. You can view the ode manual page by executing ?ode at the R prompt. By specifying method="lsoda", we cause ode to use the Livermore Solver of Ordinary Differential equations, with Automatic method switching.

Cleanup function

The cleanup function reshapes the output and adds labels if needed.

cleanup <- function (out, x) {
  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

lv_flow(
  t=seq(0,100,by=0.01),
  x=c(1,0.1),
  r=c(1,-1),
  A=matrix(c(0,-1,1,0),2,2,byrow=TRUE)
) -> dat

Plot 1

dat |>
  ggplot(aes(x=x_1,y=x_2))+
  geom_path()+
  labs(x=expression(x[1]),y=expression(x[2]))

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

Note that here we integrate three trajectories simultaneously. Each column of matrix x is a distinct initial state.

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

lv_flow(
  t=seq(1000,10000,by=0.02),
  x=c(0.1,0.1,0.3),
  r=c(1,1,-1),
  A=matrix(
    c(
      -1.0, -1.0, -10,
      -1.5, -1.0,  -1,
       5.0,  0.5,  -0.01
    ),
    3,3,byrow=TRUE
  )
) -> dat

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

lv_flow(
  t=seq(0,10000,by=0.05),
  x=matrix(
    c(
      0.1, 0.1, 0.3,
      0.1, 0.1, 0.3001
    ),
    3,2
  ),
  r=c(1,1,-1),
  A=matrix(
    c(
      -1.0, -1.0, -10,
      -1.5, -1.0,  -1,
       5.0,  0.5,  -0.01
    ),
    3,3,byrow=TRUE
  )
) -> dat

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:

circulant(1,2,3,4)
\[\begin{bmatrix} 1 &4 &3 &2 \\2 &1 &4 &3 \\3 &2 &1 &4 \\4 &3 &2 &1 \\ \end{bmatrix}\]

Triangle plot

The following codes project the positive cone \(\mathbb{R}^3_+\) onto the 3-simplex \(S_3=\left\{x\in\mathbb{R}^3_+\middle\vert x_1+x_2+x_3=1\right\}\).

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_coords

Stable interior rest point

With \(r=\begin{bmatrix} 1 &1 &1 \\ \end{bmatrix}\) 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\).

Code

lv_flow(
  t=seq(0,200,by=0.01),
  x=c(1,1,0.1),
  r=1,
  A=-circulant(1,beta=0.4,alpha=0.1)
) -> dat

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(data=.vertices)+
  geom_path()+
  geom_text(
    data=.label_coords,
    mapping=aes(label=c("x[1]","x[2]","x[3]")),
    parse=TRUE,
    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

lv_flow(
  t=seq(0,400,by=0.01),
  x=c(0.1,0.1,0.1),
  r=1,
  A=-circulant(1,beta=1.9,alpha=0.2)
) -> dat

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("x[1]","x[2]","x[3]")),
    parse=TRUE,
    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

lv_flow(
  t=seq(0,1200,by=0.03),
  x=c(0.1,0.1,0.099),
  r=1,
  A=-circulant(1,beta=1.9,alpha=0.2)
) -> dat

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("x[1]","x[2]","x[3]")),
    parse=TRUE,
    size=16,
    size.unit="pt"
  )+
  theme_minimal()+
  theme(
    axis.text=element_blank(),
    axis.title=element_blank(),
    panel.grid=element_blank()
  )

Vandermeer case

On the \(\alpha+\beta=2\) surface, we see neutrally stable orbits that fill the simplex.

Code

lv_flow(
  t=seq(0,1200,by=0.03),
  x=matrix(
    c(
      0.33,0.33,0.34,
      0.25,0.35,0.40,
      0.20,0.30,0.50,
      0.15,0.15,0.70
    ),
    nrow=3
  ),
  r=1,
  A=-circulant(1,beta=0,alpha=2)
) -> dat

Neutral case

dat |>
  pivot_longer(-time) |>
  separate(name,into=c("x","var","ic")) |>
  unite(col="name",x,var) |>
  pivot_wider() -> dat

plot_grid(
  dat |>
    ggplot(aes(x=x_1,y=x_2,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[1]),y=expression(x[2])),
  dat |>
    ggplot(aes(x=x_1,y=x_3,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[1]),y=expression(x[3])),
  dat |>
    ggplot(aes(x=x_2,y=x_3,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[2]),y=expression(x[3])),
  dat |>
    filter(time > 1100) |>
    pivot_longer(-c(time,ic)) |>
    separate(name,into=c(NA,"index")) |>
    ggplot(aes(x=time,y=value,group=interaction(index,ic),color=ic))+
    geom_line()+
    guides(color="none")+
    facet_grid(index~.,scales="free_y",
      labeller=label_bquote(rows=x[.(index)])
    )+
    labs(y="")
)

Just off the neutral case (1)

Here, we use \(\alpha=1.85\) and \(\beta=0.1\).

lv_flow(
  t=seq(0,1200,by=0.03),
  x=matrix(
    c(
      0.33,0.33,0.34,
      0.34,0.33,0.33
    ),
    nrow=3
  ),
  r=1,
  A=-circulant(1,alpha=1.85,beta=0.1)
) |>
  pivot_longer(-time) |>
  separate(name,into=c("x","var","ic")) |>
  unite(col="name",x,var) |>
  pivot_wider() -> dat

plot_grid(
  dat |>
    ggplot(aes(x=x_1,y=x_2,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[1]),y=expression(x[2])),
  dat |>
    ggplot(aes(x=x_1,y=x_3,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[1]),y=expression(x[3])),
  dat |>
    ggplot(aes(x=x_2,y=x_3,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[2]),y=expression(x[3])),
  dat |>
    filter(time > 1000) |>
    pivot_longer(-c(time,ic)) |>
    separate(name,into=c(NA,"index")) |>
    ggplot(aes(x=time,y=value,group=interaction(index,ic),color=ic))+
    geom_line()+
    guides(color="none")+
    facet_grid(index~.,scales="free_y",
      labeller=label_bquote(rows=x[.(index)])
    )+
    labs(y="")
)

Just off the neutral case (2)

Here, we use \(\alpha=1.85\) and \(\beta=0.2\).

lv_flow(
  t=seq(0,1200,by=0.03),
  x=matrix(
    c(
      0.33,0.33,0.34,
      0.34,0.33,0.33
    ),
    nrow=3
  ),
  r=1,
  A=-circulant(1,alpha=1.85,beta=0.2)
) |>
  pivot_longer(-time) |>
  separate(name,into=c("x","var","ic")) |>
  unite(col="name",x,var) |>
  pivot_wider() -> dat

plot_grid(
  dat |>
    ggplot(aes(x=x_1,y=x_2,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[1]),y=expression(x[2])),
  dat |>
    ggplot(aes(x=x_1,y=x_3,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[1]),y=expression(x[3])),
  dat |>
    ggplot(aes(x=x_2,y=x_3,color=ic))+
    geom_point(size=0.1)+
    guides(color="none")+
    labs(x=expression(x[2]),y=expression(x[3])),
  dat |>
    filter(time > 1000) |>
    pivot_longer(-c(time,ic)) |>
    separate(name,into=c(NA,"index")) |>
    ggplot(aes(x=time,y=value,group=interaction(index,ic),color=ic))+
    geom_line()+
    guides(color="none")+
    facet_grid(index~.,scales="free_y",
      labeller=label_bquote(rows=x[.(index)])
    )+
    labs(y="")
)

Download codes