Math, Nature

The Lotka-Volterra equations

November 6, 2025November 15, 2025 by Jeff Wollen

The Lotka-Volterra equations are a pair of differential equations used to model population dynamics. They assume an idealized environment with two species—one predator and one prey—and describe how the two populations change over time.

The equations were developed simultaneously by two mathematicians in the 1920s, Alfred Lotka and Vito Volterra. Lotka theorized about oscillating chemical reaction while Volterra tried to make sense of changing fish populations off the coast of Italy during World War I. The equations’ versatility is something I’ll circle back to later.

Predator-prey models can never perfectly describe reality but they still provide a solid foundation for understanding real-world systems. Let’s take a closer look.

1. The equations.

The equations describe the rate of change of predator and prey populations with respect to time.

x — prey population
y — predator population

The key idea is that the rate of change (dx/dt) depends on the current value (x). In the real world, it means that when a population spikes, it sets itself up for a big contraction. And vice versa. The environment settles into a self-correcting cycle of booms and busts.

Let’s look at the coefficients and think about what they mean.

α (alpha) — Birth rate of prey. This obviously puts upward pressure on prey population (x) so the term is positive.
β (beta) — Death rate of prey. The term puts downward pressure on prey population. βxy describes interaction between predator and prey. Either population growing will lead to more encounters. A larger β implies that encounters are more likely to end in a kill.
δ (delta) — δxy is similar in that it describes predator-prey interaction. This time it puts upward pressure on predator population. A higher δ means prey are more efficiently converted into births.
γ (gamma) — Natural death rate of predators. This term is independent of prey.

2. The solution.

The equations are considered a system of Ordinary Differential Equations (ODEs). Sometimes ODEs can be solved analytically. You can manipulate the terms and find that y(t)=3t²+4 or something like that.

In this case, we’ll need to use numerical methods to approximate a solution. Euler’s Method is the most intuitive to me. To briefly summarize:

We start with a given point.
We have an equation for the derivative, i.e. slope, at any given point.
So calculate the slope and draw a tiny line segment.
Now we have a new point. Repeat the process from there.
Eventually we’ll have a series of line segments that approximate the solution.

It’s important to use a small step size (the difference between successive values of t) because error builds up over time. Since we’re using a modern computer, that’s not a problem. This technique is much less labor-intensive than it was in Euler’s day.

And since we’re approximating curves, we will of course plot them.

3. The code.

First let’s define coefficients. For our purposes, the values are mostly arbitrary. I think the ones below provide a nice visualization so that’s what we’ll do.

alpha = 1.0
beta = 0.1
delta = 0.02
gamma = 0.5

Next, prepare the time variable, t.

dt is the time step discussed above. It’s how much we increment t before calculating a new estimated value of the dependent variable. Visually, it represents the horizontal width of the tiny line segments.

You won’t be able to see these tiny line segments on the plot. The output will look like a normal smooth curve. That’s the idea!

t_max is simply the largest value of t that we’ll consider. With a step size of 0.001, there will be 40,000 (or really 40,001) elements in t_list.

from numpy import arange

dt = 0.001
t_max = 40.0
num_steps = int(t_max / dt)
t_list = arange(0, t_max + dt, dt)

To use Euler’s Method, we need an initial point where we can begin drawing tiny line segments. Create x_list and y_list and give them initial values.

Remember, x (prey population) and y (predator population) are both functions of time. They could be written as x(t) and y(t). Usually y is a function of x but differential equations throw that convention out the window.

x_list = [10]
y_list = [5]

Now we’ll step through a for loop to generate line segments.

i will correspond to the index of x_list and y_list. We can use it to access the most recent estimated value.

for i in range(1, num_steps + 1):
    prev_x = x_list[i - 1]
    prev_y = y_list[i - 1]

    [...]

From that point, we calculate instantaneous slope using the equations above. The Lotka-Volterra equations isolate dx/dt and dy/dt on the left-hand side, so the right-hand side represents slope.

The next value will be slope (rise/run) multiplied by step size (run). That leaves us with rise, i.e. the vertical change. Add it to the previous value, append to the list, and return to the top of the loop.

Now the previous value is the one we just calculated. Draw a line segment beginning there and repeat the process 40,000 times.

for i in range(1, num_steps + 1):
    prev_x = x_list[i - 1]
    prev_y = y_list[i - 1]

    dx_dt = alpha * prev_x - beta * prev_x * prev_y
    next_x = prev_x + dx_dt * dt
    x_list.append(next_x)

    dy_dt = delta * prev_x * prev_y - gamma * prev_y
    next_y = prev_y + dy_dt * dt
    y_list.append(next_y)

4. Plot the curves.

I’ll use a custom Matplotlib style that will be linked at the bottom of this post.

Create an Axes instance (ax) for plotting.

import matplotlib.pyplot as plt

plt.style.use("wollen_blue.mplstyle")

fig, ax = plt.subplots()

Now plot the two curves, x(t) and y(t).

ax.plot(t_list, x_list, label="Prey (x)")

ax.plot(t_list, y_list, label="Predator (y)")

Since we’re already hard-coding several values in this script, I won’t feel guilty about hard-coding ticks and window limits as well. We can combine them into a single set method along with labels.

ax.set(xticks=range(0, 45, 5),
       yticks=range(0, 80, 10),
       xlim=(-0.5, 40.5),
       ylim=(0, 70.5),
       xlabel="Time (t)",
       ylabel="Population",
       title="Lotka-Volterra Predator-Prey Model")

I want to write the Lotka-Volterra equations on the plot using Matplotlib’s built-in LaTeX rendering. Activate it by putting $ symbols at the beginning and end of a string. Let’s locate text halfway between the two rightmost grid lines.

ax.text(37.5, 57.5, r'$\frac{dx}{dt} = \alpha x - \beta x y$', ha="center", va="center")
ax.text(37.5, 53.0, r'$\frac{dy}{dt} = \delta x y - \gamma y$', ha="center", va="center")
ax.text(37.5, 48.7, r'$\alpha = ' + f'{alpha:.02}' + '$', size=9, ha="center", va="center")
ax.text(37.5, 46.2, r'$\beta = ' + f'{beta:.02}' + '$', size=9, ha="center", va="center")
ax.text(37.5, 43.7, r'$\delta = ' + f'{delta:.02}' + '$', size=9, ha="center", va="center")
ax.text(37.5, 41.2, r'$\gamma = ' + f'{gamma:.02}' + '$', size=9, ha="center", va="center")

Place a legend directly above the equations. bbox_to_anchor allows for more precise control of the legend’s location. (0, 0) corresponds to the plot’s bottom-left and (1, 1) is top-right.

ax.legend(bbox_to_anchor=(0.93, 0.9), loc="center")

Finally, save the figure with an increased dpi.

plt.savefig("predator_prey.png", dpi=150)

5. The output.

I think that’s a cool result! The model shows that when there are few predators around, prey population booms. But that means plenty of food for predators. They reproduce faster and eventually catch up, leading to a bust in prey population. That means less food for predators, and so on.

Our approximation is pretty good with a time step of dt=0.001. We see a change of about 0.2% from peak to peak. If you scaled back step size to 0.01, you would easily see values drift with your naked eye.

Maybe theoretical rabbits and foxes aren’t that exciting to you, but the Lotka-Volterra equations can be applied in other fields as well. For example, they’re often used in economics. You can imagine large firms as predators and smaller firms as prey.

Or sellers and buyers in financial markets. Day trading stocks is hard because some of the smartest people in the world are in those markets. Many of them probably build models based on these concepts (and a lot more).

If nothing else, it’s nice to know that oftentimes you can brute force your way to a differential equation solution with a little code.

A cute bunny picture. Photo credit: instagram.com/mister.carrots.

Download the Matplotlib style.

Full code:

from numpy import arange
import matplotlib.pyplot as plt


alpha = 1.0
beta = 0.1
delta = 0.02
gamma = 0.5

dt = 0.001
t_max = 40.0
num_steps = int(t_max / dt)
t_list = arange(0, t_max + dt, dt)

x_list = [10]
y_list = [5]

for i in range(1, num_steps + 1):
    prev_x = x_list[i - 1]
    prev_y = y_list[i - 1]

    dx_dt = alpha * prev_x - beta * prev_x * prev_y
    next_x = prev_x + dx_dt * dt
    x_list.append(next_x)

    dy_dt = delta * prev_x * prev_y - gamma * prev_y
    next_y = prev_y + dy_dt * dt
    y_list.append(next_y)

plt.style.use("wollen_blue.mplstyle")

fig, ax = plt.subplots()

ax.plot(t_list, x_list, label="Prey (x)")

ax.plot(t_list, y_list, label="Predator (y)")

ax.set(xticks=range(0, 45, 5),
       yticks=range(0, 80, 10),
       xlim=(-0.5, 40.5),
       ylim=(0, 70.5),
       xlabel="Time (t)",
       ylabel="Population",
       title="Lotka-Volterra Predator-Prey Model")

ax.text(37.5, 57.5, r'$\frac{dx}{dt} = \alpha x - \beta x y$', ha="center", va="center")
ax.text(37.5, 53.0, r'$\frac{dy}{dt} = \delta x y - \gamma y$', ha="center", va="center")
ax.text(37.5, 48.7, r'$\alpha = ' + f'{alpha:.02}' + '$', size=9, ha="center", va="center")
ax.text(37.5, 46.2, r'$\beta = ' + f'{beta:.02}' + '$', size=9, ha="center", va="center")
ax.text(37.5, 43.7, r'$\delta = ' + f'{delta:.02}' + '$', size=9, ha="center", va="center")
ax.text(37.5, 41.2, r'$\gamma = ' + f'{gamma:.02}' + '$', size=9, ha="center", va="center")

ax.legend(bbox_to_anchor=(0.93, 0.9), loc="center")

plt.savefig("predator_prey.png", dpi=150)

wollen.org