Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section you will be able to:

Derive the Lotka–Volterra equations $x' = \alpha x - \beta x y$ and $y' = -\gamma y + \delta x y$ from one ecological idea — coupling through encounter terms $xy$ .
Find all equilibria of the system and identify the coexistence fixed point $(x^*, y^*) = (\gamma/\delta,\, \alpha/\beta)$ .
Read a phase portrait — closed orbits, the direction of motion, the location of the center — without doing any algebra.
Prove that $H(x,y) = \delta x - \gamma \ln x + \beta y - \alpha \ln y$ is conserved, and use it to explain WHY the orbits are closed loops.
Linearize about the coexistence point, compute the Jacobian, get purely imaginary eigenvalues $\pm i\sqrt{\alpha \gamma}$ , and conclude the period of small oscillations $T \approx 2\pi / \sqrt{\alpha \gamma}$ .
Implement the system in plain Python with RK4 and verify the conserved quantity with PyTorch's autograd to machine precision.
Explain Volterra's paradox — why predators always lag prey by a quarter-period — and connect it to data from Hudson's Bay fur returns and the Adriatic fishery.

The Story: Fish, Foxes and Vito Volterra

"Why, during the war years when fishing stopped in the Adriatic, did the proportion of predatory fish — sharks, skates, rays — actually INCREASE?" — the question Umberto D'Ancona put to his father-in-law Vito Volterra in 1925.

Volterra was a mathematician, not an ecologist. But he noticed that the question had nothing to do with the biology of any particular species — it was purely about the coupling between two populations. He wrote down two equations, solved them on the back of an envelope, and produced a prediction so counter-intuitive it became one of the founding results of mathematical biology: if you stop fishing, the predators benefit MORE than the prey, because the prey were suppressing them by attracting other predators.

At almost the same moment, Alfred Lotka — a Polish-American chemist working at Johns Hopkins — wrote down the identical equations while thinking about autocatalytic chemical reactions. Two people, two continents, two fields, one system. The equations have been called Lotka–Volterra ever since.

Why this section is the heart of the chapter

Every other system in this chapter is linear (we just have to find eigenvalues). Lotka–Volterra is the first nonlinear system you can't solve with eigenvalues alone — and yet it has a complete qualitative theory, a conserved quantity, and an analytic period formula. It is the proof-of-concept that nonlinear systems can be understood without closed-form solutions.

Building the Model from One Idea

Imagine just prey, no predators. Rabbits in a meadow with infinite grass and no foxes. Birth rate proportional to the existing population, so

\frac{dx}{dt} = \alpha\, x \quad \Longrightarrow \quad x(t) = x_0 e^{\alpha t}.

Pure exponential growth. Rabbits everywhere.

Now imagine just predators, no prey. Foxes starving in an empty meadow. Death rate proportional to existing predators:

\frac{dy}{dt} = -\gamma\, y \quad \Longrightarrow \quad y(t) = y_0 e^{-\gamma t}.

Pure exponential decay. Foxes vanish.

The interesting physics is in the coupling. Suppose prey and predators wander randomly. Then the rate of encounters between them is proportional to the product of densities: $\text{encounters per unit time} \propto x\,y$ . This is the same mass-action law that governs chemical reactions (rate of $A + B \to C$ scales with $[A][B]$ ) — which is exactly why Lotka came across the same equations.

Each encounter removes some prey (predation) and feeds the predators. Add those two coupling terms to the previous equations and the closed system pops out:

\frac{dx}{dt} = \alpha\,x \;-\; \beta\,x\,y \qquad \text{(prey)}

\frac{dy}{dt} = -\gamma\,y \;+\; \delta\,x\,y \qquad \text{(predator)}

Read it in plain English

Prey: "I grow at rate α per individual, but predators eat me at rate β per predator-prey pair."
Predator: "I starve at rate γ per individual, but every prey I eat converts (with efficiency δ/β) into more of me."

The minus sign on the $\beta x y$ term is where the entire ecological tension lives. Without it, prey explode. Add it back, and the prey grow only when predators are scarce — which automatically makes predators grow — which suppresses prey again. The system "breathes."

The Four Parameters and What They Really Mean

Every Lotka–Volterra investigation hinges on four positive constants. Their dimensions matter — units make the formulas honest.

Symbol	Name	Units	Mental picture
α	Prey birth rate	1 / time	Rabbits with infinite grass double every ln 2 / α time units.
β	Predation rate	1 / (predator · time)	One predator eats β · x prey per unit time when prey density is x.
γ	Predator death rate	1 / time	Predators starve to half-population in ln 2 / γ time units without prey.
δ	Conversion efficiency	1 / (prey · time)	Each prey eaten produces δ / β fraction of a new predator. Always δ < β in reality.

The model is invariant under $(x, y) \mapsto (cx, dy)$ with $(\beta, \delta) \mapsto (\beta/d, \delta/c)$ . That means you can choose units of prey and predators freely — the PHASE PORTRAIT (the shape of the orbits) only cares about the dimensionless products $\alpha\gamma$ and $\beta\delta$ .

Equilibrium Points: Where Nothing Changes

An equilibrium of a system $\mathbf{x}' = \mathbf{F}(\mathbf{x})$ is a state where the velocity vector is zero — pop. you in there and you stay forever. To find them, set both right-hand sides to zero simultaneously:

\alpha x - \beta x y \;=\; x(\alpha - \beta y) \;=\; 0

-\gamma y + \delta x y \;=\; y(-\gamma + \delta x) \;=\; 0

Each equation factors into two roots. We need a state that satisfies BOTH equations at once, so we cross-pair the factors:

Prey factor = 0	Predator factor = 0	Equilibrium
x = 0	y = 0	(0, 0) — extinction
x = 0	−γ + δ·x = 0 → x = γ/δ	contradiction (x can't be both)
α − β·y = 0 → y = α/β	y = 0	contradiction (y can't be both)
α − β·y = 0 → y = α/β	−γ + δ·x = 0 → x = γ/δ	(γ/δ, α/β) — coexistence

Two equilibria, both intuitive:

Trivial fixed point $(0, 0)$ — nothing in the meadow, nothing happens.
Coexistence fixed point $(x^*, y^*) = (\gamma/\delta,\; \alpha/\beta)$ — exactly enough prey to sustain the predators, exactly enough predators to balance prey growth.

Notice the swap: the equilibrium prey density depends on predator parameters (γ, δ) and the equilibrium predator density depends on prey parameters (α, β). That is what coupling looks like.

The Phase Plane: Following the State Vector

Time-series plots show $x(t)$ and $y(t)$ separately. But the system is 2-D: its complete state at any instant is a point $(x(t), y(t))$ in the positive quadrant of the plane. Watching THAT point move is the geometric picture.

At every state, the velocity vector $\mathbf{F}(x, y) = (\alpha x - \beta x y,\; -\gamma y + \delta x y)$ tells you which way the orbit is heading. Walking around the coexistence point you can read off four qualitative regions, each named for the direction of motion:

Region	Where	dx/dt	dy/dt	What's happening
Right of x, below y	x > γ/δ and y < α/β	+ (prey grow)	+ (predators grow)	Lots of prey, few predators — both populations rising. State moves up-right.
Right of x, above y	x > γ/δ and y > α/β	− (prey crash)	+ (predators still grow)	Predators take over. State moves up-left.
Left of x, above y	x < γ/δ and y > α/β	− (prey still falling)	− (predators starve)	Both crashing. State moves down-left.
Left of x, below y	x < γ/δ and y < α/β	+ (prey recovering)	− (predators still few)	Recovery. State moves down-right back to start.

Round and round. The orbit is a counter-clockwise closed loop around $(\gamma/\delta,\, \alpha/\beta)$ . The next section's interactive widget lets you watch this with sliders.

Interactive: Move the Knobs

The widget below integrates the Lotka–Volterra system with RK4 (3000 steps over $t \in [0, 30]$ ). The left panel is the time series of $x(t)$ (green) and $y(t)$ (pink). The right panel is the same trajectory drawn in the phase plane $(x, y)$ . The white dot is the current state — watch how it walks along the closed orbit while the time-series oscillations slosh up and down.

Loading predator–prey animator…

Try this: push β (predation) up to 1.0 and watch the orbit shrink towards the fixed point — heavy predation means smaller oscillations around a stable balance. Push α (prey birth) up to 2.0 and the orbit balloons outward — more food at the bottom makes the cycle more violent. The center of the orbit moves too: it sits at $(\gamma/\delta, \alpha/\beta)$ , so increasing α slides the center upward.

Lotka–Volterra orbits are marginally stable, not asymptotically. If your numerical integrator has any dissipative error (Forward-Euler does!), the orbit will spiral OUT over many cycles. RK4 is conservative enough at this step size that you won't see drift — try setting $h = 0.1$ in your own integrator and watch the cycle slowly grow.

Why the Orbits Close: A Conserved Quantity

Why CLOSED loops? In a general 2-D nonlinear system, orbits can spiral, wander chaotically, or blow up. Closed loops are the rare special case. For Lotka–Volterra they happen because there is a smooth function $H(x, y)$ that is constant along every orbit — the equation $H(x, y) = c$ picks out one loop for each value of $c$ .

Here is how to find it. Divide the two equations to eliminate $dt$ :

\frac{dy}{dx} = \frac{-\gamma y + \delta x y}{\alpha x - \beta x y} = \frac{y(-\gamma + \delta x)}{x(\alpha - \beta y)}.

This is now an ODE in $(x, y)$ with $t$ gone. Separate variables — get all $y$ s on one side and all $x$ s on the other:

\frac{\alpha - \beta y}{y}\, dy = \frac{-\gamma + \delta x}{x}\, dx.

Split the fractions and integrate term by term:

\int \left(\frac{\alpha}{y} - \beta\right) dy = \int \left(-\frac{\gamma}{x} + \delta\right) dx

\alpha \ln y - \beta y = -\gamma \ln x + \delta x + C.

Move everything to one side, define $H(x,y)$ :

\boxed{\; H(x, y) = \delta x - \gamma \ln x + \beta y - \alpha \ln y \;}

Then every orbit satisfies $H(x(t), y(t)) = \text{const}$ . Different orbits = different constants = different level curves.

Why H has a minimum exactly at the fixed point

Compute the partial derivatives: $\partial_x H = \delta - \gamma/x$ and $\partial_y H = \beta - \alpha/y$ . Setting both to zero gives $x = \gamma/\delta$ and $y = \alpha/\beta$ — exactly the coexistence fixed point. The Hessian there is positive-definite (both $\partial^2_x H = \gamma/x^2 > 0$ and $\partial^2_y H = \alpha/y^2 > 0$ ), so $H$ has a strict minimum there. That is why orbits curl around the fixed point — they are level curves of a bowl-shaped surface.

Linearization at the Coexistence Point

The coexistence equilibrium $(x^*, y^*)$ is interesting; the extinction point $(0, 0)$ is not (it's an unstable saddle — perturb either species upward and it grows). Near $(x^*, y^*)$ we can replace the nonlinear system with its linear approximation, the Jacobian matrix:

J(x, y) = \begin{pmatrix} \alpha - \beta y & -\beta x \\ \delta y & -\gamma + \delta x \end{pmatrix}.

Substitute $x = \gamma/\delta$ , $y = \alpha/\beta$ . The diagonal entries vanish identically:

J(x^*, y^*) = \begin{pmatrix} 0 & -\beta\gamma/\delta \\ \delta\alpha/\beta & 0 \end{pmatrix}.

Trace is 0, determinant is $\alpha\gamma > 0$ . The characteristic polynomial is therefore $\lambda^2 + \alpha\gamma = 0$ , giving eigenvalues

\lambda_{1,2} = \pm i\sqrt{\alpha\gamma}.

Purely imaginary eigenvalues. The linearization says the equilibrium is a center: small perturbations rotate around it with angular frequency $\omega = \sqrt{\alpha\gamma}$ and don't decay or grow. This is consistent with what we already proved using the conserved quantity $H$ : closed orbits.

Linearization caveat

For a center, the linearization is not always reliable as a guide to the nonlinear system — generic nonlinear perturbations turn centers into spirals (either inward or outward). Lotka–Volterra happens to have the conserved $H$ , which forces the orbits to remain closed loops even at large amplitudes. Without that conservation law, ANY small perturbation could destroy the oscillation.

Interactive: The Jacobian as a Linear System

The widget below shows phase portraits of a generic linear 2-D system $\mathbf{u}' = A \mathbf{u}$ where you can slide the four entries of $A$ . The defaults match the Lotka–Volterra Jacobian at the coexistence point for our parameters $(\alpha, \beta, \gamma, \delta) = (1, 0.4, 0.8, 0.25)$ : $a = 0,\; b = -\beta\gamma/\delta = -1.28,\; c = \delta\alpha/\beta = 0.625,\; d = 0$ . Verify that the eigenvalues you read off match $\pm i\sqrt{\alpha\gamma} = \pm i\sqrt{0.8} \approx \pm 0.894 i$ .

Loading phase-portrait explorer…

Experiment: set $a$ to a small NEGATIVE value (say −0.05). The center turns into a stable spiral — orbits decay to the origin. Set it slightly positive (+0.05) and the center becomes an unstable spiral — orbits explode outward. The boundary case a = 0 (where the trace vanishes) is the ONLY way to get a center. It is a knife-edge condition, and the Lotka–Volterra system sits exactly on it.

The Period of Small Oscillations

For small departures from the coexistence point, the linearization rotates with angular frequency $\omega = \sqrt{\alpha\gamma}$ . The period is therefore

T \approx \frac{2\pi}{\omega} = \frac{2\pi}{\sqrt{\alpha\gamma}}.

Plug in $\alpha = 1$ , $\gamma = 0.8$ : $T \approx 2\pi/\sqrt{0.8} \approx 7.02$ time units. Look at the interactive widget's time series — count the gap between successive prey peaks. You will measure exactly ≈ 7 even though the oscillation is far from infinitesimal (we run amplitudes of order 1). The period formula is exact in the linear limit and remarkably accurate even for large orbits — a useful empirical bonus of the Hamiltonian structure.

The period depends ONLY on $\alpha$ and $\gamma$ — the "diagonal" rates of pure birth and pure death. The coupling constants $\beta$ and $\delta$ only set the AMPLITUDE of the oscillation and the position of the equilibrium, not its tempo. This is the kind of insight you would NEVER spot from numerics alone.

Phase Lag: Predators Always Arrive Late

In the interactive time series, the pink predator peak comes AFTER the green prey peak by roughly a quarter of the period — about $T/4 \approx 1.75$ time units in the default case. This is not a numerical artifact. It is what center-style oscillation does: the two variables trace out a circle in the phase plane, and a circle has its x-coordinate maximum exactly $\pi/2$ radians before its y-coordinate maximum.

Ecologically: predators feast on the prey peak, but it takes time for their offspring to mature and contribute to $y$ . By the time the predator population catches up, the prey is already crashing — which then starves the predators a quarter-cycle later. Hudson's Bay Company fur returns from the 1840s show the snowshoe-hare / lynx cycle with this exact phase lag, repeating with a ≈ 10-year period for over a century.

Volterra's paradox: a constant proportional harvest (say, fishing both species at rate $h$ ) replaces $\alpha \to \alpha - h$ and $\gamma \to \gamma + h$ . The fixed point shifts to $(x^*, y^*) = ((\gamma + h)/\delta,\; (\alpha - h)/\beta)$ — MORE prey, FEWER predators. That is exactly D'Ancona's observation in reverse: stop fishing (h → 0) and predator equilibrium rises while prey equilibrium falls.

Worked Example (Step-by-Step, Collapsible)

Open the box to follow a complete pen-and-paper computation: find the equilibria, derive the period, compute one RK4 step by hand, and check that $H$ is conserved across that step to 4 decimal places.

Expand: full hand-computed walkthrough with

(\alpha, \beta, \gamma, \delta) = (1, 0.4, 0.8, 0.25)

and

(x_0, y_0) = (2.5, 1.0)

Step 1 — Find the equilibria

Set both RHS to zero: $x(\alpha - \beta y) = 0$ and $y(-\gamma + \delta x) = 0$ . Either $x = 0$ (then $y = 0$ ), giving the trivial point (0, 0); or

y = \alpha/\beta = 1/0.4 = 2.50, \quad x = \gamma/\delta = 0.8/0.25 = 3.20.

Coexistence at $(x^*, y^*) = (3.20,\, 2.50)$ .

Step 2 — Linearize there

$J = \begin{pmatrix} 0 & -\beta\gamma/\delta \\ \delta\alpha/\beta & 0 \end{pmatrix} = \begin{pmatrix} 0 & -1.28 \\ 0.625 & 0 \end{pmatrix}.$ Determinant $= (-1.28)(0.625)\cdot(-1) = 0.80$ ; trace = 0. Eigenvalues from $\lambda^2 + 0.80 = 0$ : $\lambda = \pm i\sqrt{0.80} \approx \pm 0.894\,i$ .

Step 3 — Period of small oscillations

T \;=\; \frac{2\pi}{\sqrt{\alpha\gamma}} \;=\; \frac{2\pi}{\sqrt{0.80}} \;\approx\; 7.025.

Step 4 — One RK4 step from (2.5, 1.0) with h = 0.1

Right-hand side helper: $f(x,y) = (x - 0.4xy,\; -0.8y + 0.25xy)$ .

k1 = f(2.5, 1.0): $(2.5 - 0.4\cdot 2.5 \cdot 1,\; -0.8\cdot 1 + 0.25\cdot 2.5 \cdot 1) = (1.5,\; -0.175).$

Midpoint with k1: $(x_m, y_m) = (2.5 + 0.05\cdot 1.5,\; 1.0 + 0.05\cdot(-0.175)) = (2.575,\; 0.99125)$ . k2 = f(2.575, 0.99125): $(2.575 - 0.4\cdot 2.575 \cdot 0.99125,\; -0.8\cdot 0.99125 + 0.25\cdot 2.575 \cdot 0.99125) \approx (1.5539,\; -0.1551).$

Midpoint with k2: $(2.5 + 0.05\cdot 1.5539,\; 1.0 + 0.05\cdot(-0.1551)) = (2.5777,\; 0.99224)$ . k3 = f(2.5777, 0.99224): $\approx (1.5547,\; -0.1542).$

Endpoint with k3: $(2.5 + 0.1\cdot 1.5547,\; 1.0 + 0.1\cdot(-0.1542)) = (2.6555,\; 0.98458)$ . k4 = f(2.6555, 0.98458): $\approx (1.6093,\; -0.1343).$

Weighted update: $x_1 = 2.5 + (0.1/6)(1.5 + 2\cdot 1.5539 + 2\cdot 1.5547 + 1.6093) \approx 2.6559.$

$y_1 = 1.0 + (0.1/6)(-0.175 - 2\cdot 0.1551 - 2\cdot 0.1542 - 0.1343) \approx 0.98446.$

Step 5 — Verify H is conserved across this step

Before: $H(2.5, 1.0) = 0.25 \cdot 2.5 - 0.8\ln 2.5 + 0.4\cdot 1 - \ln 1$ $= 0.625 - 0.7330 + 0.4 - 0 = 0.2920.$

After: $H(2.6559, 0.98446) = 0.25\cdot 2.6559 - 0.8\ln 2.6559 + 0.4\cdot 0.98446 - \ln 0.98446$ $\approx 0.6640 - 0.7811 + 0.3938 + 0.0157 = 0.2924.$

Drift: $\Delta H \approx 0.0004$ after one step of $h = 0.1$ . RK4 is 4th-order accurate so we expect $\Delta H \sim O(h^5) \approx 10^{-5}$ ; the observed $4 \times 10^{-4}$ is consistent with the 6-decimal rounding we did on the intermediate slopes. The PyTorch test below will show drift ≪ $10^{-6}$ at $h = 0.01$ .

Plain Python: A Hand-Built RK4 Integrator

Read this block top to bottom. Every line is annotated. The point is not the integrator itself — that's a textbook RK4 — but to see how each constant in the code corresponds to a physical rate and how the right-hand side is JUST a literal translation of the equations.

Lotka–Volterra with Pure-Python RK4

🐍lotka_volterra_rk4.py

Explanation(19)

Code(36)

1Why pure Python first

Before reaching for a solver, we hand-roll one. The Lotka–Volterra equations are coupled and nonlinear — there is no closed-form x(t) or y(t). Numerics is the only way through. Pure Python lets you see every multiplication.

4α = 1.0 (prey birth rate)

Units: 1/time. With α = 1 and no predators, the prey alone would satisfy x′ = x, so x(t) = x₀ eᵗ — population doubles roughly every 0.7 time units.

EXAMPLE

no predators → x(t) = x₀ · exp(α·t)

5β = 0.4 (predation rate per encounter)

Units: 1/(predator·time). Each unit of predator removes 0.4·x prey per unit time. A larger β means predators are deadlier — even few of them suppress prey hard.

6γ = 0.8 (predator natural death)

Units: 1/time. With no prey to eat, predators alone follow y′ = −0.8 y, so y(t) = y₀ e^(−0.8 t) — population halves about every 0.87 time units.

7δ = 0.25 (predator gain per prey eaten)

Units: 1/(prey·time). Of every prey killed (rate β·x·y), only a fraction δ/β = 0.625 turns into new predator biomass. This is the ecological inefficiency — eating a rabbit doesn't make a whole fox.

9f(x, y) returns the whole RHS as a tuple

A 2-D system has TWO derivatives at every point. The function maps the state vector (x, y) to the velocity vector (x′, y′). That velocity vector is exactly what we will integrate.

EXAMPLE

state (x,y) → velocity (x′,y′)

11Prey equation: x′ = α x − β x y

Plain-English read-out: 'prey grow at rate α (per prey) MINUS predator-removal at rate β (per predator-prey pair).' The minus sign couples it to y — without predators it would explode exponentially.

EXAMPLE

with x=2.5, y=1.0:  x′ = 1·2.5 − 0.4·2.5·1 = 1.5

12Predator equation: y′ = −γ y + δ x y

Plain-English: 'predators die at rate γ (per predator) PLUS they gain mass δ for every prey eaten (rate x·y).' Without prey, the +δ x y term vanishes and y collapses exponentially.

EXAMPLE

with x=2.5, y=1.0:  y′ = −0.8·1 + 0.25·2.5·1 = −0.175

15rk4_step — one Runge–Kutta–4 step

Euler's method (just x + h·f(x)) is too inaccurate for closed orbits — error accumulates and the orbit spirals out. RK4 evaluates the slope FOUR times per step and combines them so the local error is O(h⁵). For a 30-unit run with h = 0.01 the orbit closes to ~6 decimal places.

17k1 — slope at the start of the interval

First sample of the velocity field, taken exactly at the current (x, y). If we used only this slope we would get Euler's method.

EXAMPLE

k1 = f(x, y) = (1.5, −0.175)

18k2 — slope at the midpoint using k1 to estimate halfway

Take an Euler half-step using k1, then sample the field at THAT predicted midpoint. This is the trick that turns Euler into a midpoint method (already 2nd-order).

19k3 — slope at the midpoint AGAIN, but using k2

Re-estimate the midpoint slope using the better midpoint guess k2. Two looks at the middle of the interval, each more refined than the last.

20k4 — slope at the END of the interval, using k3 to predict it

Use k3 to take a FULL Euler step and then evaluate the field at the end point. So we have slopes at t, t+h/2, t+h/2, t+h — both ends and the middle (twice).

21Weighted average: (k1 + 2·k2 + 2·k3 + k4) / 6

Simpson's-rule weights: the middle samples get twice the weight of the endpoints. This combination cancels the lowest-order Taylor errors so the leftover error scales as h⁵.

EXAMPLE

RK4 local error ≈ C · h⁵

25x_next, y_next — the new state after one step h

Both x and y advance together, using their own weighted slopes. Coupling is automatic because every f(...) call sees the CURRENT (x, y) of the joint state.

28Initial condition (x₀, y₀) = (2.5, 1.0)

Start with 2.5 'units of prey' and 1.0 'unit of predator'. Units could be hundreds of fish or thousands of bacteria — Lotka–Volterra is scale-free in this respect.

29Step size h = 0.01

Small enough that the velocity field hardly changes across one step (prey period here is ≈ 7 time units, so 700 steps per period — way more than enough for RK4).

32The main loop — 3000 steps means t ranges over [0, 30]

Each iteration: advance the state by h, record the new (x, y, t). At the end we have three parallel lists that can be plotted as time-series OR as a phase-plane orbit by plotting xs vs ys.

36Sanity check — the fixed point is (γ/δ, α/β) = (3.20, 2.50)

We will derive this algebraically in the Equilibria section. Print it here so you can compare to the orbit's centroid — the orbit should circle exactly this point.

EXAMPLE

γ/δ = 0.8/0.25 = 3.20      α/β = 1.0/0.4 = 2.50

17 lines without explanation

1# Lotka-Volterra in pure Python — no NumPy, no SciPy.
2# Goal: build every step of an RK4 integrator yourself so the math has nowhere to hide.
3
4ALPHA = 1.0   # prey birth rate
5BETA  = 0.4   # predation rate (per encounter)
6GAMMA = 0.8   # predator death rate (in absence of prey)
7DELTA = 0.25  # predator gain (per prey eaten)
8
9def f(x, y):
10    """Right-hand side of the Lotka-Volterra ODE: returns (dx/dt, dy/dt)."""
11    dxdt =  ALPHA * x - BETA  * x * y
12    dydt = -GAMMA * y + DELTA * x * y
13    return dxdt, dydt
14
15def rk4_step(x, y, h):
16    """One classical fourth-order Runge-Kutta step of size h."""
17    k1x, k1y = f(x,                   y                  )
18    k2x, k2y = f(x + 0.5 * h * k1x,   y + 0.5 * h * k1y  )
19    k3x, k3y = f(x + 0.5 * h * k2x,   y + 0.5 * h * k2y  )
20    k4x, k4y = f(x +       h * k3x,   y +       h * k3y  )
21    x_next = x + (h / 6.0) * (k1x + 2 * k2x + 2 * k3x + k4x)
22    y_next = y + (h / 6.0) * (k1y + 2 * k2y + 2 * k3y + k4y)
23    return x_next, y_next
24
25# Initial condition: 2.5 prey, 1.0 predator.
26x, y = 2.5, 1.0
27t, h = 0.0, 0.01
28
29xs, ys, ts = [x], [y], [t]
30for _ in range(3000):              # integrate from t=0 to t=30
31    x, y = rk4_step(x, y, h)
32    t   += h
33    xs.append(x); ys.append(y); ts.append(t)
34
35print(f"final (t={t:.2f}):  x={x:.4f},  y={y:.4f}")
36print(f"fixed point (γ/δ, α/β) = ({GAMMA/DELTA:.2f}, {ALPHA/BETA:.2f})")

To check yourself: at $t = 30$ with the default parameters and initial condition, you should print $x \approx 0.49$ and $y \approx 3.30$ — the orbit has nearly closed (returning to a state very close to $(2.5, 1.0)$ would happen near $t \approx 4T \approx 28$ ).

PyTorch: Verifying the Conserved Quantity with Autograd

We claimed $H(x, y)$ is invariant along the flow. Numerics with finite step sizes can't prove it — the drift you saw above is partly the truth and partly numerical error. PyTorch autograd, on the other hand, computes the gradient of $H$ EXACTLY by the chain rule. If $\nabla H \cdot \mathbf{F}$ evaluates to zero at every state we pick — not just on one orbit — then $H$ really is conserved.

Autograd Sanity-Check: ∇H · v = 0 at every state

🐍conserved_quantity.py

Explanation(12)

Code(32)

1Why PyTorch instead of more NumPy

PyTorch's autograd gives us EXACT partial derivatives of H(x, y). No finite differences, no truncation error. We can numerically check the analytical claim ∇H · velocity = 0 to machine precision.

8Import torch

torch carries tensors that remember every operation done to them, so calling torch.autograd.grad(H, x) walks the recorded graph backwards using the chain rule.

12Define H(x, y) — the conserved quantity

Derived in the 'Why orbits close' section. Notice the algebraic mirror: the prey part (δ x − γ ln x) and the predator part (β y − α ln y) each look like a Lyapunov-style 'energy' term with a minimum at x* = γ/δ and y* = α/β.

EXAMPLE

H is a Hamiltonian for the conjugate variables (u, v) = (ln x, ln y).

14torch.log is element-wise natural log

Natural log appears because we are integrating dx / x and dy / y when separating variables to derive H. ln keeps the gradient autograd needs (d/dx ln x = 1/x).

16velocity(x, y) — the Lotka-Volterra RHS, again

Same formula as plain-Python, now built on tensors so autograd can track gradients through it if we ever need ∇·velocity (we won't here, but it's automatic).

21Sample 5 random positive states

We pick 5 unrelated points in the (x, y) > 0 quadrant. If H is really conserved, ∇H · velocity must vanish at EVERY one of them — not just on a single orbit.

EXAMPLE

xs ≈ tensor([2.49, 3.81, 0.55, 1.41, 1.85])

22requires_grad=True

Tells autograd to record all subsequent ops on these tensors. Without this flag, torch.autograd.grad would raise 'one of the differentiated Tensors does not require grad'.

26h_val = H(x, y)

A scalar tensor. Behind the scenes torch built a small graph: multiply, subtract, log, etc. — every op needed to differentiate H w.r.t. its inputs.

27grad_x, grad_y = torch.autograd.grad(h_val, (x, y))

Asks the engine: 'differentiate this scalar h_val with respect to BOTH x and y, give me the two partial derivatives.' The result is a tuple of two tensors.

EXAMPLE

∂H/∂x = δ − γ/x,   ∂H/∂y = β − α/y  (autograd derives both).

28vx, vy = velocity(x, y)

The right-hand side of the ODE evaluated at the same state. ∇H · v will tell us whether H is changing as we slide along the trajectory.

29dot = grad_x · vx + grad_y · vy

Compute (∂H/∂x)·(dx/dt) + (∂H/∂y)·(dy/dt). By the chain rule this is dH/dt along the flow. If H is conserved it must be 0 — and you'll see numbers like 1e-16, i.e. literally zero up to float-32 round-off.

30Print the verification line

For each of the 5 sampled points you should see a 'dot' value on the order of 1e-15 to 1e-7. That is autograd confirming, with no derivation tricks on our side, that H is exactly invariant under the Lotka–Volterra flow.

20 lines without explanation

1# PyTorch demo: USE AUTOGRAD to verify that
2#     H(x, y) = δ x − γ ln(x) + β y − α ln(y)
3# is conserved along the flow x' = αx − βxy,  y' = −γy + δxy.
4#
5# A function H is conserved iff its gradient is perpendicular to the
6# velocity field at every state — i.e. ∇H · (x', y') = 0 exactly.
7
8import torch
9
10ALPHA, BETA, GAMMA, DELTA = 1.0, 0.4, 0.8, 0.25
11
12def H(x, y):
13    """Lotka-Volterra Hamiltonian-like first integral."""
14    return DELTA * x - GAMMA * torch.log(x) + BETA * y - ALPHA * torch.log(y)
15
16def velocity(x, y):
17    return ALPHA * x - BETA * x * y, -GAMMA * y + DELTA * x * y
18
19# Sample 5 random states in the positive quadrant.
20torch.manual_seed(0)
21xs = torch.rand(5, requires_grad=True) * 4 + 0.5
22ys = torch.rand(5, requires_grad=True) * 4 + 0.5
23
24for x, y in zip(xs, ys):
25    h_val = H(x, y)
26    grad_x, grad_y = torch.autograd.grad(h_val, (x, y))
27    vx, vy = velocity(x, y)
28    dot = grad_x * vx + grad_y * vy
29    print(f"state=({x.item():.3f}, {y.item():.3f}) "
30          f"  ∇H={grad_x.item(): .4f},{grad_y.item(): .4f} "
31          f"  velocity={vx.item(): .4f},{vy.item(): .4f} "
32          f"  ∇H·v = {dot.item(): .2e}")

Expected printed output: five lines, each with $\nabla H \cdot \mathbf{v}$ on the order of $10^{-15}$ to $10^{-7}$ . That is autograd confirming the theorem at random points — a pen-and-paper proof would compute $(\delta - \gamma/x)(\alpha x - \beta x y) + (\beta - \alpha/y)(-\gamma y + \delta x y) = 0$ algebraically, which is exactly what this code verifies numerically.

Real-World Applications

Field	x = ?	y = ?	Real cycle observed
Boreal ecology	Snowshoe hare	Canada lynx	Hudson's Bay Company fur returns 1845–1935: ≈ 10-year cycle, predator lags prey by ≈ 2 yr.
Marine biology	Plankton	Predatory fish	D'Ancona's Adriatic data 1914–1923: predator fraction rose during WWI fishing pause.
Epidemiology	Susceptibles	Infectives	SIS / SIR variants are a special case — coupling term is contact rate β·S·I.
Chemical kinetics	Autocatalytic species	Inhibitor	Lotka's original 1910 motivation: oscillating Belousov–Zhabotinsky-style reactions.
Economics	Capital stock	Wage share	Goodwin's 1967 growth-cycle model is literally Lotka–Volterra in disguise.
Cybersecurity	Vulnerable hosts	Active worms	Code-Red and Slammer infection curves fit predator-prey before patches dampened them.

Every entry in that table came from someone noticing the same coupling structure — "a thing that grows when it eats another thing that regrows when not being eaten." Once you recognise the pattern you will see Lotka–Volterra in places that have nothing to do with biology.

Where the Model Breaks (and What to Do Next)

Honest accounting. The classical Lotka–Volterra model is the simplest possible coupling — it deliberately ignores almost every real biological complication. Each failure mode points to a more realistic extension you will meet in section 23.9 and beyond:

Unbounded prey in the absence of predators. The model says $x' = \alpha x$ grows forever. Reality has carrying capacity $K$ . Fix: Logistic prey term $\alpha x (1 - x/K)$ — gives the modified Rosenzweig–MacArthur model with limit cycles (asymptotically stable, NOT marginal).
Linear predator response to prey density. A real predator saturates — once it's eaten too many rabbits it slows down. Fix: Type-II functional response $\beta x y / (1 + h x)$ — captures handling time h.
No interspecific competition within prey or predator.Adds quadratic $-x^2$ / $-y^2$ terms (Lotka–Volterra competition — section 23.9).
Deterministic. Real populations are stochastic; for small populations the closed orbit can be lost to extinction in one unlucky cycle. Fix: stochastic differential equations driven by Brownian motion.
Spatially uniform. Real predators and prey occupy patches. Fix: reaction-diffusion PDEs — gives travelling waves and Turing patterns (Chapter 26).
Three or more species. Lotka–Volterra with three species can be CHAOTIC. The cycle structure that kept us safe with two species disappears.

A good rule of thumb: classical Lotka–Volterra is a cartoon. It exists not because it's realistic but because it's the ONLY nonlinear two-species model where a complete analytical theory exists. Every more accurate model uses numerical methods, the conserved quantity disappears, and you fall back on phase-plane geometry — exactly what we built here.

Summary

The Lotka–Volterra system $x' = \alpha x - \beta xy,\; y' = -\gamma y + \delta xy$ encodes mass-action coupling between a prey and a predator population. Every term has a direct ecological reading.
It has two equilibria: extinction $(0, 0)$ (saddle) and coexistence $(\gamma/\delta,\, \alpha/\beta)$ (center).
A conserved quantity $H(x,y) = \delta x - \gamma \ln x + \beta y - \alpha \ln y$ forces every orbit to be a closed loop around the coexistence point.
Linearizing at the coexistence point gives a Jacobian with purely imaginary eigenvalues $\pm i\sqrt{\alpha\gamma}$ , predicting period $T \approx 2\pi / \sqrt{\alpha\gamma}$ for small oscillations.
Predators always lag prey by $T/4$ — the geometric quarter-period of a center.
Volterra's paradox: harvesting both species at the same rate $h$ shifts equilibrium UP for prey and DOWN for predators. Stop harvesting, the reverse happens — exactly the 1920s Adriatic observation.
Pure-Python RK4 is enough to integrate this faithfully; PyTorch autograd confirms the conserved quantity to machine precision without a single hand-derivation.

The deepest message of Lotka–Volterra: nonlinear coupling between two variables, with no time delay, no noise, and no spatial structure, is already enough to produce sustained oscillation. Predator-prey cycles are not a quirk of biology — they are a consequence of one minus sign in the right place.