Lotka-Volterra equations

Uses Derivatives of real functions of one variable

The Lotka-Volterra equations, also known as the predator-prey equations, are a model for the dynamics of two interacting species in an ecosystem. This model takes the form of two coupled non-linear differential equations.

The number of prey
 x : ℝ → ℝ
and the number of predators
 y : ℝ → ℝ
are functions of time. For mathematical convience (the model goes back to the early 20th century, before digital computers), they are represented by real numbers rather than integers, as this allows the use of differential equations instead of the clumsier difference equations.

The time evolution of prey and predator populations is given by
 lvx: 𝒟(x) = (α × x) - (β × x × y),
 lvy: 𝒟(y) = (δ × x × y) - (γ × y),
where α : ℝ.p, β : ℝ.p, γ : ℝ.p, and δ : ℝ.p are constants describing the interactions of the two species.

The term α × x describes the exponential growth of the prey population x in the absence of predators. The model assumes no other contraints (such as limited food supply or limited space) on the growth of the prey population.

The term β × x × y describes the decrease of the prey population x due to predation, which is assumed to be proportional to the number of encounters between individuals of the two species, which in turn is proportional to both x and y.

The term γ × y describes an exponential decay of the predator population y in the absence of prey, the assumption being that predators will either starve or leave the region of interest.

The term δ × x × y describes the increase of the predator population, which is assumed to be proportional to the amount of food, given by the prey population x.

Numerical simulation: Euler method

A simple numerical simulation can be performed by applying the Euler method. It computes the values of approximations x≅e : (ℕ → ℝ) and y≅e : (ℕ → ℝ) to the populations
n:ℕ, x[n × Δt] and n:ℕ, y[n × Δt] on a time grid with spacing Δt : ℝ.p.

Given initial values x≅e[0] and y≅e[0], values for later times are computed as
 lvx/euler: n:ℕ.nz, x≅e[n] ⇒ x≅e[n - 1] + (Δt × ((α × x≅e) - (β × x≅e × y≅e))[n - 1])
 lvy/euler: n:ℕ.nz, y≅e[n] ⇒ y≅e[n - 1] + (Δt × ((δ × x≅e × y≅e) - (γ × y≅e))[n - 1])

Example

Note: this example is taken from the Wikipedia entry cited above.

Suppose there are two species of animals, a baboon (prey) and a cheetah (predator). If the initial conditions are _example x≅e[0] ⇒ 10 baboons and _example y≅e[0] ⇒ 10 cheetahs, one can plot the progression of the two species over time; given the parameters that the growth and death rates of baboon are _example α ⇒ 11/10 and _example β ⇒ 2/5 while that of cheetahs are _example δ ⇒ 1/10 and _example γ ⇒ 2/5 respectively.

The Euler method yield two recurrence equations, lvx/euler and lvy/euler. They can be solved using the solver implemented in RESequence Object subclass: #RESequence instanceVariableNames: 'context rules variables sequence stepSize indexVariable parameters iterationCode' classVariableNames: '' package: 'RecurrenceEquations' (see Solving recurrence equations):

sequence1 := RESequence
	context: (thisSnippet page lzSubcontext: #example)
	rules: #('lvx/euler' 'lvy/euler')
	parameters: { 'Δt ⇒ 1/10' }
	initialValues: (Dictionary
						with: 'x≅e' -> #(10.0)
						with: 'y≅e' -> #(10.0))
	indexVariable: 't'
	stepSize: 'Δt'.
sequence1 addSteps: 500.
  

Before trying to interpret those results, let's be cautious. These curves have sharp peaks, which are always suspicious in the context of numerical solutions. Is our time step small enough? Is our discretization (the Euler method) precise enough? If the answer to one of these questions is "no", then we could be looking at numerical artifacts. Let's try again with a smaller step size:

sequence2 := RESequence
	context: (thisSnippet page lzSubcontext: #example)
	rules: #('lvx/euler' 'lvy/euler')
	parameters: { 'Δt ⇒ 1/100' }
	initialValues: (Dictionary
						with: 'x≅e' -> #(10.0)
						with: 'y≅e' -> #(10.0))
	indexVariable: 't'
	stepSize: 'Δt'.
sequence2 addSteps: 5000.
  

This looks rather different! The increase of the population peaks from one cycle to the next looks like a numerical artifact. This is in fact a well-known issue with the Lotka-Volterra equations, as explained e.g. in this blog post.

Numerical simulation: the symplectic Euler method

A common strategy for stabilizing numerical solutions is to switch to an implicit integration scheme, in which the quantities on the right-hand side are evaluated at the end of the step rather than at the beginning . Doing this would result in the implicit Euler method . For the specific case of the Lotka-Volterra equations, which are a Hamiltonian system, the best approach is the semi-implicit Euler method, also known as the symplectic Euler method, which consists of treating one variable (we choose x) implicly and the other one explicitly.

We define the solutions obtained with the symplectic Euler method as x≅s : (ℕ → ℝ) and y≅s : (ℕ → ℝ) with the rules:
lvx/symplectic: n:ℕ.nz, x≅s[n] ⇒ x≅s[n - 1] ÷ (1 - (Δt × (α - (β × y≅s[n - 1]))))
 lvy/symplectic: n:ℕ.nz, y≅s[n] ⇒ y≅s[n - 1] + (Δt × ((δ × x≅s[n] × y≅s[n - 1]) - (γ × y≅s[n - 1])))

Back to the example

These equations yield the expected periodic cycles:

sequence3 := RESequence
	context: (thisSnippet page lzSubcontext: #example)
	rules: #('lvx/symplectic' 'lvy/symplectic')
	parameters: { 'Δt ⇒ 1/10' }
	initialValues: (Dictionary
						with: 'x≅s' -> #(10.0)
						with: 'y≅s' -> #(10.0))
	indexVariable: 't'
	stepSize: 'Δt'.
sequence3 addSteps: 500.
  

It is also interesting to plot x≅s against y≅s:

sequence3 parametricPlot: 'x≅s' with: 'y≅s'