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 competition model by adding inter-specific terms to two coupled logistic equations and explain why every coefficient appears where it does.
Sketch the two nullclines $N_1 + \alpha_{12} N_2 = K_1$ and $N_2 + \alpha_{21} N_1 = K_2$ directly from their intercepts.
List the four equilibria $(0,0),\ (K_1,0),\ (0,K_2)$ and the coexistence point, and find each by inspection.
Compute the Jacobian at any equilibrium and classify it from the signs of the eigenvalues.
Predict which of four regimes will play out (1 wins, 2 wins, coexistence, bistability) from the two comparisons $K_1\ \text{vs.}\ \alpha_{12} K_2$ and $K_2\ \text{vs.}\ \alpha_{21} K_1$ .
Implement the model with an RK4 integrator in plain Python and again as a batched nn.Module in PyTorch.

The Question Behind the Math

"Two species share one bowl of food. Can they share — or does one drive the other out?"

Predator–prey models in the previous section described one species eating another. The dynamic was inherently asymmetric. Here we tackle the more common situation in ecology: two species who do not eat each other but who both want the same things — nutrients, sunlight, shelter, nesting sites. Their interaction is symmetric in flavour but its long-time outcome can be wildly different, and the difference is entirely captured by two numbers we will call competition coefficients.

The question that organises everything below is: given the growth rates, the carrying capacities, and the strengths of competition, who survives? The answer turns out to be a clean four-way classification you can read off a single picture — and the picture is two intersecting straight lines in the $(N_1, N_2)$ plane.

From Logistic Growth to Competition

For a single species the logistic equation says

\displaystyle \frac{dN}{dt} = r\,N\left(1 - \frac{N}{K}\right).

Read this in two steps. The factor $r N$ is the unrestrained exponential growth rate. The factor in parentheses is a brake: when the population $N$ equals the carrying capacity $K$ , the bracket is zero and growth stops. Halfway to $K$ the brake is at 50%, and at $N = 2K$ the brake is $-1$ — the population is overshooting and actually shrinking.

Now put a second species in the same bowl. The cleanest possible modification is: let species 2 count too, but in species-1 units. One individual of species 2 occupies as much of the niche as $\alpha_{12}$ individuals of species 1. So replace $N$ in the brake by $N_1 + \alpha_{12} N_2$ . Do the symmetric thing for species 2 — replace its $N$ by $N_2 + \alpha_{21} N_1$ . That is the entire derivation.

Pedagogical note. The coefficient

\alpha_{12}

means "what species 2 does to species 1". The first subscript is the SUFFERER, the second is the OFFENDER. This is the opposite of matrix indexing, and almost every student gets it wrong the first time. Read the index twice before plugging numbers in.

The Competition Model

The full system is

\displaystyle \frac{dN_1}{dt} = r_1\,N_1\left(1 - \frac{N_1 + \alpha_{12} N_2}{K_1}\right),

\displaystyle \frac{dN_2}{dt} = r_2\,N_2\left(1 - \frac{N_2 + \alpha_{21} N_1}{K_2}\right).

Six parameters, two state variables. Every term has a direct ecological meaning:

Symbol	Meaning	Units
r₁, r₂	Intrinsic per-capita growth rate of each species when alone and rare	1/time
K₁, K₂	Carrying capacity if the other species is absent (the solo logistic asymptote)	individuals
α₁₂	Per-capita effect of species 2 on species 1 (in units of species 1)	dimensionless
α₂₁	Per-capita effect of species 1 on species 2 (in units of species 2)	dimensionless
N₁, N₂	Population sizes at time t	individuals

Notice that putting $\alpha_{12} = \alpha_{21} = 0$ gives back two decoupled logistic equations — the species ignore each other and grow toward their own carrying capacities. The full model only differs from this trivial limit through the two cross-terms $\alpha_{12} N_2 / K_1$ and $\alpha_{21} N_1 / K_2$ . All of the rich behaviour below lives inside those two terms.

Nullclines: Where Each Species Stops Changing

The geometric backbone of any 2D ODE system is its nullclines: the curves on which one of the rates vanishes. Set $dN_1/dt = 0$ :

r_1 N_1 \left(1 - \frac{N_1 + \alpha_{12} N_2}{K_1}\right) = 0.

A product is zero iff a factor is zero, so either $N_1 = 0$ (the vertical axis) or the bracket is zero. The bracket gives a straight line:

N_1 + \alpha_{12} N_2 = K_1.

It crosses the $N_1$ -axis at $N_1 = K_1$ (set $N_2 = 0$ ) and the $N_2$ -axis at $N_2 = K_1/\alpha_{12}$ (set $N_1 = 0$ ). Identically for species 2:

N_2 + \alpha_{21} N_1 = K_2,

with intercepts $(K_2/\alpha_{21},\,0)$ and $(0,\,K_2)$ .

The geometric punchline

On its own nullcline, a species is momentarily frozen: its rate is exactly zero. Above it the bracket is negative and the species shrinks; below it the bracket is positive and the species grows. Combined, the two nullclines partition the positive quadrant into 3 or 4 regions, and the SIGN of $(dN_1/dt,\,dN_2/dt)$ in each region is enough to predict the long-time outcome.

The Four Equilibria

An equilibrium is a point where both rates are zero simultaneously — geometrically, an intersection of an $N_1$ -nullcline with an $N_2$ -nullcline.

$E_0 = (0, 0)$ — both extinct. The axes $N_1 = 0$ and $N_2 = 0$ intersect at the origin. Both rates vanish there.
$E_1 = (K_1, 0)$ — only species 1. The species-1 sloped nullcline meets the $N_2 = 0$ axis. Species 1 is at its solo carrying capacity.
$E_2 = (0, K_2)$ — only species 2. Symmetric.
$E^* = (N_1^*, N_2^*)$ — coexistence. The two sloped nullclines meet (if they meet inside the positive quadrant) at $\displaystyle N_1^* = \frac{K_1 - \alpha_{12} K_2}{1 - \alpha_{12}\alpha_{21}},\quad N_2^* = \frac{K_2 - \alpha_{21} K_1}{1 - \alpha_{12}\alpha_{21}}.$

The coexistence formulas come from solving the two linear nullcline equations $N_1 + \alpha_{12} N_2 = K_1$ and $\alpha_{21} N_1 + N_2 = K_2$ — straight Cramer's rule. The denominator $1 - \alpha_{12}\alpha_{21}$ appearing in both expressions is exactly the determinant of the 2×2 coefficient matrix. When that determinant is zero the lines are parallel and coexistence either doesn't exist or fills a whole line — a degenerate case.

Stability via the Jacobian

Near any equilibrium $(\bar{N}_1, \bar{N}_2)$ we linearise. Write $f_1(N_1, N_2) = r_1 N_1 (1 - (N_1 + \alpha_{12} N_2)/K_1)$ and similarly $f_2$ . The Jacobian is

\displaystyle J = \begin{pmatrix} \dfrac{\partial f_1}{\partial N_1} & \dfrac{\partial f_1}{\partial N_2} \\[6pt] \dfrac{\partial f_2}{\partial N_1} & \dfrac{\partial f_2}{\partial N_2} \end{pmatrix} = \begin{pmatrix} r_1\left(1 - \dfrac{2N_1 + \alpha_{12} N_2}{K_1}\right) & -\dfrac{r_1\,\alpha_{12} N_1}{K_1} \\[10pt] -\dfrac{r_2\,\alpha_{21} N_2}{K_2} & r_2\left(1 - \dfrac{2N_2 + \alpha_{21} N_1}{K_2}\right) \end{pmatrix}.

At each equilibrium plug in $(\bar{N}_1, \bar{N}_2)$ and read off the eigenvalues from the trace and determinant. Their signs determine the type:

Eigenvalues	Type	Long-time behaviour near the point
Both negative	Stable node	All nearby trajectories converge to the point
Both positive	Unstable node	All nearby trajectories run away
Opposite signs	Saddle	Most trajectories run away; one special direction comes in
Complex with Re<0	Stable spiral	Trajectories spiral inward (rare in pure competition)

Quick check at the boundary equilibria. At

E_1 = (K_1, 0)

the Jacobian becomes upper triangular and you can read the eigenvalues straight off the diagonal:

\lambda_1 = -r_1

and

\lambda_2 = r_2 (1 - \alpha_{21} K_1 / K_2)

. The first is always negative, so

E_1

is stable iff

\alpha_{21} K_1 > K_2

— species 1 wins against any small invasion by species 2.

The Four Regimes — Geometry Predicts the Future

Compare $K_1$ with $\alpha_{12} K_2$ , and $K_2$ with $\alpha_{21} K_1$ . These two comparisons encode which nullcline lies above which on each axis, and they yield exactly four cases.

Case	Conditions	Stable equilibrium	Ecology
1. Species 1 wins	K₁ > α₁₂ K₂ and K₂ < α₂₁ K₁	(K₁, 0)	Competitive exclusion — species 2 driven extinct
2. Species 2 wins	K₁ < α₁₂ K₂ and K₂ > α₂₁ K₁	(0, K₂)	Symmetric — species 1 driven extinct
3. Stable coexistence	K₁ > α₁₂ K₂ and K₂ > α₂₁ K₁	(N₁, N₂)	Weak inter-specific competition. Both survive forever
4. Founder control	K₁ < α₁₂ K₂ and K₂ < α₂₁ K₁	(K₁, 0) or (0, K₂)	Strong competition both ways — initial conditions decide who wins

The intuition behind the conditions: the inequality $K_1 > \alpha_{12} K_2$ says that, even when species 2 is at its solo carrying capacity, the brake on species 1 is still less than 100%, so species 1 can invade. If both invasion criteria hold, both species can come back from rare, and the only place they can co-rest is $E^*$ .

Case 4 is the most counter-intuitive: both single-species equilibria are stable. Whichever species starts higher pushes the other below its invasion threshold and wins. The coexistence point $E^*$ still exists but is a saddle; its stable manifold is the boundary between the two basins of attraction.

Interactive Phase-Plane Explorer

Below is the phase plane with vector field, nullclines, and equilibria for any parameter setting you choose. Use the four presets to flip between the four regimes; drag the sliders to deform the picture; and click anywhere in the field to drop a new trajectory. Notice how the same geometric ingredients (two lines + four dots) tell the whole story.

Loading competition explorer…

Take three minutes to play. The single most useful experiment is to load Founder control, drop a trajectory at $(60, 20)$ , then drop another at $(20, 60)$ . The two trajectories live on opposite sides of an invisible line — the stable manifold of the saddle — and end up at different attractors.

Worked Example (Step-by-Step)

Take the canonical "stable coexistence" setting and grind through every step by hand.

Click to expand the full hand calculation

Step 1. Parameters

$r_1 = r_2 = 1,\quad K_1 = 100,\quad K_2 = 80,\quad \alpha_{12} = 0.5,\quad \alpha_{21} = 0.3$ . Time is measured in whatever units make $r = 1$ ; populations in head count.

Step 2. Which regime?

Check the two comparisons:

$K_1$ vs $\alpha_{12} K_2 = 0.5 \times 80 = 40$ : $100 > 40$ , so species 1 CAN invade when species 2 sits at $K_2$ .
$K_2$ vs $\alpha_{21} K_1 = 0.3 \times 100 = 30$ : $80 > 30$ , so species 2 can also invade when species 1 sits at $K_1$ .

Both invasion criteria hold ⇒ regime 3, stable coexistence.

Step 3. Locate the coexistence equilibrium

Apply the formula directly. The shared denominator is

1 - \alpha_{12}\alpha_{21} = 1 - 0.5 \times 0.3 = 1 - 0.15 = 0.85.

Then

N_1^* = \frac{100 - 0.5 \times 80}{0.85} = \frac{60}{0.85} \approx 70.588,

N_2^* = \frac{80 - 0.3 \times 100}{0.85} = \frac{50}{0.85} \approx 58.824.

Step 4. Jacobian at $E^*$

Plug $(N_1^*, N_2^*) \approx (70.588, 58.824)$ into the Jacobian. With $r_1 = r_2 = 1$ :

$J_{11} = 1 - (2 \times 70.588 + 0.5 \times 58.824)/100 = 1 - 1.41176 - 0.29412 = -0.706$ .
$J_{12} = -0.5 \times 70.588 / 100 = -0.353$ .
$J_{21} = -0.3 \times 58.824 / 80 = -0.221$ .
$J_{22} = 1 - (2 \times 58.824 + 0.3 \times 70.588)/80 = 1 - 1.47059 - 0.26471 = -0.735$ .

J^* \approx \begin{pmatrix} -0.706 & -0.353 \\ -0.221 & -0.735 \end{pmatrix}.

Step 5. Eigenvalues

Trace and determinant first:

$\mathrm{tr}(J^*) = -0.706 + (-0.735) = -1.441$ .
$\det(J^*) = (-0.706)(-0.735) - (-0.353)(-0.221) = 0.519 - 0.078 = 0.441$ .

Discriminant $\Delta = \mathrm{tr}^2 - 4\det = 2.077 - 1.764 = 0.313$ , so $\sqrt{\Delta} \approx 0.559$ , and

\lambda_{1,2} = \frac{-1.441 \pm 0.559}{2} \approx -0.441,\ -1.000.

Both eigenvalues are real and negative ⇒ stable node. Every trajectory that starts in the positive quadrant gets pulled into $(70.588, 58.824)$ .

Step 6. Sanity check at the other equilibria

At $E_1 = (100, 0)$ the Jacobian is upper-triangular: $\lambda_1 = -r_1 = -1$ and $\lambda_2 = 1 - \alpha_{21} K_1/K_2 = 1 - 30/80 = 0.625$ . One negative, one positive ⇒ saddle, as predicted by the regime table.

At $E_2 = (0, 80)$ the Jacobian is lower-triangular: $\lambda_1 = 1 - \alpha_{12} K_2/K_1 = 1 - 40/100 = 0.6$ and $\lambda_2 = -r_2 = -1$ . Another saddle.

At $E_0 = (0,0)$ both diagonal entries are $+1$ ⇒ unstable node.

Step 7. Picture

Three saddles/unstable points plus one stable node fully constrain the flow. Trajectories from anywhere in the positive quadrant curve into $(70.588, 58.824)$ . The geometry has just told us the future without us solving a single transcendental.

Time-Series Companion

The phase plane shows where the system goes. The time-series view shows how fast. Slide the initial populations and watch $N_1(t)$ (cyan) and $N_2(t)$ (rose) both bend toward the coexistence levels $(70.588, 58.824)$ . The dashed horizontal lines mark the solo carrying capacities — note that both species equilibrate below their solo $K$ , because each is making room for the other.

Loading time-series…

Plain Python Implementation

The cleanest way to internalise the model is to integrate it yourself. Below is a complete plain-Python program that reproduces the coexistence run with classical RK4. Every name in the code maps one-to-one onto the equations above.

Lotka–Volterra competition by hand

🐍python

Explanation(14)

Code(41)

1Import NumPy

We only need NumPy for array allocation (np.empty_like) and timestep grids (np.linspace). The ODE itself is written without any vectorised tricks so that the math is one-to-one with the equations on paper.

4Set the six parameters

These are exactly the values used in the worked example: weak inter-specific competition. r1=r2=1 fixes the time scale. K1=100, K2=80 are the solo carrying capacities. α12=0.5 says one extra individual of species 2 'eats' the same resources as half an individual of species 1; α21=0.3 is the symmetric statement.

EXECUTION STATE

r1, r2 = 1.0, 1.0

K1, K2 = 100.0, 80.0

a12, a21 = 0.5, 0.3

10The right-hand side rhs(N1, N2)

This is the heart of the model — one line per ODE. Read each as 'population × per-capita growth rate', and the per-capita rate is 1 minus how full the niche is. When N1 + α12·N2 hits K1, the bracket vanishes and N1 stops growing.

11Species-1 equation

f1 = r1·N1·(1 − (N1 + α12·N2)/K1). Notice α12·N2: species 2 is 'rented out' in species-1 units. If α12 < 1 the competitor counts for less than a conspecific (weak competition); if α12 > 1 it counts for more (strong competition).

12Species-2 equation

Symmetric. f2 = r2·N2·(1 − (N2 + α21·N1)/K2). Plug N1=0 here: you recover plain logistic growth toward K2. The coupling lives entirely inside the α21·N1 term.

16RK4 step

Classical 4-stage Runge-Kutta has local error O(h⁵) and global error O(h⁴) — vastly more accurate than Euler at the same h. Each (k_x, k_y) is a SLOPE estimate at some intermediate point in (t, t+h).

17k1: slope at the left endpoint

k1 is the derivative evaluated AT the current state (N1, N2). If we used Euler's method we'd stop here and step h·k1. RK4 will refine it three more times.

18k2: slope at the midpoint, using k1

Take a half-step with k1 to land at a tentative midpoint, then re-evaluate the slope there. If the field is curving, this midpoint slope is a much better representative than the endpoint slope k1.

19k3: another midpoint estimate, using k2

Same as k2 but using k2 to choose the midpoint. The double-midpoint redundancy is exactly what kills the third-order error term.

20k4: slope at the right endpoint, using k3

Take a full step with k3 to reach the tentative right end and evaluate slope there. Combining k1, k2, k3, k4 with the magic weights (1, 2, 2, 1)/6 gives a fourth-order accurate step.

21Combine the four slopes

The weighted average (k1 + 2·k2 + 2·k3 + k4)/6 is the effective slope for the whole interval. Multiply by h and add to the current state — done. This is the standard RK4 recipe, identical to what SciPy's solve_ivp uses internally.

25simulate(N1_0, N2_0, …)

Pre-allocate the output arrays and run rk4_step in a tight loop. tmax=40 with h=0.02 gives 2000 steps — usually enough to reach the long-time attractor of any of the four regimes.

33Three starting points

We launch the same system from three corners of the (N1, N2) plane. With the chosen parameters every trajectory should converge to (N1*, N2*) ≈ (70.6, 58.8) — the stable-coexistence equilibrium. The printed end-states are the strongest possible confirmation that the geometry of nullclines correctly predicts the future.

35Printed result

Expected output (after the loop): start ( 10.0, 10.0) end ( 70.588, 58.824) start ( 20.0, 60.0) end ( 70.588, 58.824) start ( 80.0, 5.0) end ( 70.588, 58.824) All three start in very different places, all three settle at the SAME coexistence point.

27 lines without explanation

1import numpy as np
2
3# Parameters (the worked-example values)
4r1, r2  = 1.0, 1.0
5K1, K2  = 100.0, 80.0
6a12     = 0.5     # how much one unit of species 2 hurts species 1
7a21     = 0.3     # how much one unit of species 1 hurts species 2
8
9# 1.  The right-hand side of the ODE system.
10#     Returns dN1/dt and dN2/dt at the current (N1, N2).
11def rhs(N1, N2):
12    f1 = r1 * N1 * (1.0 - (N1 + a12 * N2) / K1)
13    f2 = r2 * N2 * (1.0 - (N2 + a21 * N1) / K2)
14    return f1, f2
15
16# 2.  One step of the classical 4-stage Runge-Kutta integrator.
17def rk4_step(N1, N2, h):
18    k1x, k1y = rhs(N1, N2)
19    k2x, k2y = rhs(N1 + 0.5 * h * k1x, N2 + 0.5 * h * k1y)
20    k3x, k3y = rhs(N1 + 0.5 * h * k2x, N2 + 0.5 * h * k2y)
21    k4x, k4y = rhs(N1 + h * k3x,        N2 + h * k3y)
22    N1_new = N1 + (h / 6.0) * (k1x + 2*k2x + 2*k3x + k4x)
23    N2_new = N2 + (h / 6.0) * (k1y + 2*k2y + 2*k3y + k4y)
24    return N1_new, N2_new
25
26# 3.  Integrate forward in time from an initial condition.
27def simulate(N1_0, N2_0, tmax=40.0, h=0.02):
28    steps = int(tmax / h)
29    ts = np.linspace(0.0, tmax, steps + 1)
30    N1 = np.empty_like(ts)
31    N2 = np.empty_like(ts)
32    N1[0], N2[0] = N1_0, N2_0
33    for i in range(steps):
34        N1[i+1], N2[i+1] = rk4_step(N1[i], N2[i], h)
35    return ts, N1, N2
36
37# 4.  Run it for three different starting points.
38for (N1_0, N2_0) in [(10.0, 10.0), (20.0, 60.0), (80.0, 5.0)]:
39    ts, N1, N2 = simulate(N1_0, N2_0)
40    print(f"start ({N1_0:5.1f}, {N2_0:5.1f})   "
41          f"end ({N1[-1]:6.3f}, {N2[-1]:6.3f})")

Run it and you should see all three trajectories — starting from very different corners — print the same final state $(70.588,\,58.824)$ . That is the computational confirmation of the geometric prediction.

Same Model in PyTorch

Once you have it in plain Python, it is one short hop to PyTorch: wrap the parameters as nn.Parameters and broadcast the right-hand side over a batch dimension. The payoff is parallel integration of many initial conditions on GPU and automatic differentiation through every coefficient — so you can fit the model to data with torch.optim.Adam instead of by hand.

Lotka–Volterra competition as an nn.Module

🐍python

Explanation(14)

Code(66)

1Imports

torch for tensors and autograd, torch.nn for the Module / Parameter machinery. We keep the model deliberately small — no neural net, just a tiny dynamical system written in the PyTorch idiom so it composes with everything else (CUDA, batches, autograd, optimisers).

4Class CompetitionModel(nn.Module)

Wrapping the model as a Module gives us three free things: (a) every parameter is automatically registered for autograd, (b) .to(device) moves them to GPU with a single call, (c) torch.save / torch.load round-trips the full state.

11log-parameterise positive quantities

Growth rates and carrying capacities must stay positive. Storing log(r1), log(r2), log(K1), log(K2) and exp-ing on the way out is a standard reparameterisation trick: gradient descent sees a free real number, but the model only ever sees a positive value. No projection step needed.

16soft-positive parameterisation for the alphas

Same idea but with softplus(_a) = log(1 + exp(_a)). softplus(0)=ln 2≈0.693, so we initialise _a as log(exp(α)−1) — the inverse of softplus — to recover the requested α. This stays differentiable through 0 unlike a hard clamp.

27Property accessors

Exposing r1, r2, K1, K2, α12, α21 as @property means the rest of the code reads exactly like the math: model.r1, model.K1. The transformation back from the log/softplus representation is invisible at the call site.

33forward(N)

Vectorised in the LAST dimension only — N has shape (..., 2) so we can run a whole batch of initial conditions in parallel. f1 and f2 are computed with broadcasted scalars (model.r1, model.K1, …) and stacked back together as (..., 2).

34Unpack the state

N1 = N[..., 0], N2 = N[..., 1]. The Ellipsis '…' is the key — it preserves any leading batch dimensions. With N of shape (3, 2) we get N1, N2 of shape (3,), so f1, f2 of shape (3,), and torch.stack(..., dim=-1) gives back (3, 2).

35Species-1 ODE, broadcast version

Identical math to the plain Python version. Because every operand is a tensor, this line ALSO carries gradients: if you wanted to fit α12 to data, the chain rule from a loss through this line is computed automatically by autograd.

38torch.stack along dim=-1

Re-assemble the two channels into one (..., 2) tensor. dim=-1 means 'add the new axis at the very end', regardless of how many batch dimensions sat in front.

41rk4_step works on ANY tensor shape

Because model(N) is itself a (..., 2) tensor, the same RK4 weighted combination from the plain-Python code now updates 3 (or 3000) trajectories with one call. This is the parallelism that makes PyTorch attractive even for tiny dynamical systems.

50Instantiate the model

All six parameters are nn.Parameters — you could pass model.parameters() to any optimiser and use observed (N1(t), N2(t)) data to learn r1, r2, K1, K2, α12, α21 from a real ecosystem.

54Batch of three initial conditions

N has shape (3, 2). Row k is the (N1, N2) starting point of trajectory k. This is exactly the geometry that turns a serial Python loop into a single parallel call.

60Integration loop

2000 RK4 steps. On CPU this finishes in a fraction of a second; on GPU you could push the batch to a million initial conditions and still feel fast. The same code runs unchanged.

63Expected print

tensor([[70.5882, 58.8235], [70.5882, 58.8235], [70.5882, 58.8235]]). All three batched trajectories land at the same coexistence equilibrium — the basin of attraction for this regime is the whole positive quadrant.

52 lines without explanation

1import torch
2from torch import nn
3
4class CompetitionModel(nn.Module):
5    """
6    Lotka-Volterra competition as a PyTorch Module.
7    The six parameters live as nn.Parameter so you can either
8      (a)  hold them fixed and only use forward() as the ODE RHS, or
9      (b)  optimise them to fit observed time-series data via autograd.
10    """
11    def __init__(self, r1=1.0, r2=1.0, K1=100.0, K2=80.0,
12                 a12=0.5, a21=0.3):
13        super().__init__()
14        self.log_r1 = nn.Parameter(torch.tensor(float(r1)).log())
15        self.log_r2 = nn.Parameter(torch.tensor(float(r2)).log())
16        self.log_K1 = nn.Parameter(torch.tensor(float(K1)).log())
17        self.log_K2 = nn.Parameter(torch.tensor(float(K2)).log())
18        # raw alphas live unconstrained; softplus keeps them non-negative
19        self._a12 = nn.Parameter(torch.tensor(float(a12)).log().expm1().log())
20        self._a21 = nn.Parameter(torch.tensor(float(a21)).log().expm1().log())
21
22    @property
23    def r1(self):  return self.log_r1.exp()
24    @property
25    def r2(self):  return self.log_r2.exp()
26    @property
27    def K1(self):  return self.log_K1.exp()
28    @property
29    def K2(self):  return self.log_K2.exp()
30    @property
31    def a12(self): return torch.nn.functional.softplus(self._a12)
32    @property
33    def a21(self): return torch.nn.functional.softplus(self._a21)
34
35    def forward(self, N):
36        """N has shape (..., 2). Returns dN/dt of the same shape."""
37        N1, N2 = N[..., 0], N[..., 1]
38        f1 = self.r1 * N1 * (1.0 - (N1 + self.a12 * N2) / self.K1)
39        f2 = self.r2 * N2 * (1.0 - (N2 + self.a21 * N1) / self.K2)
40        return torch.stack([f1, f2], dim=-1)
41
42
43def rk4_step(model, N, h):
44    k1 = model(N)
45    k2 = model(N + 0.5 * h * k1)
46    k3 = model(N + 0.5 * h * k2)
47    k4 = model(N + h * k3)
48    return N + (h / 6.0) * (k1 + 2*k2 + 2*k3 + k4)
49
50
51# 1.  Build the model with the worked-example parameters.
52model = CompetitionModel(r1=1.0, r2=1.0, K1=100.0, K2=80.0,
53                         a12=0.5, a21=0.3)
54
55# 2.  Batch of three different initial conditions, all integrated in parallel.
56N = torch.tensor([[10., 10.],
57                  [20., 60.],
58                  [80.,  5.]])                     # shape (3, 2)
59
60# 3.  Integrate for 40 units of time with step h = 0.02.
61h, tmax = 0.02, 40.0
62steps = int(tmax / h)
63for _ in range(steps):
64    N = rk4_step(model, N, h)
65
66print(N)        # shape (3, 2) — every row should be ~(70.59, 58.82)

The Competitive Exclusion Principle

Cases 1, 2 and 4 all end with at least one extinction. That is not an accident of the model — it is the mathematical face of a famous ecological idea due to G. F. Gause (1934): two species competing for exactly the same single resource cannot coexist indefinitely. In the limit $\alpha_{12} \to 1$ and $\alpha_{21} \to 1$ (the two species are functionally identical), the coexistence point either disappears or lies on a line of degenerate equilibria — exact coexistence is not robust to any perturbation.

Stable coexistence (regime 3) only exists when the two species specialise: each must compete more strongly with itself than with the other. Mathematically, that is $\alpha_{12} < K_1/K_2$ and $\alpha_{21} < K_2/K_1$ simultaneously. In ecology this is called niche partitioning, and it is why finches on different Galápagos islands evolved different beak shapes.

Common Pitfalls

Index direction. $\alpha_{12}$ is the effect OF species 2 ON species 1, not the other way round. The first index is the sufferer. Plug-and-chug errors here flip the whole prediction.
Confusing carrying capacity with coexistence. $(K_1, 0)$ is the carrying capacity of species 1 when alone; $N_1^*$ is its coexistence population. They are different and both can be below $K_1$ at the same time.
Picking the wrong rk step size. With $r = 1$ the characteristic time scale is ~1, so $h = 0.02$ is fine. If you scale $r$ up by 100 you must scale $h$ down by the same factor or RK4 will diverge.
Forgetting the positivity constraint. The model itself preserves $N_i \ge 0$ exactly (each rate is proportional to $N_i$ ), but numerical drift can push trajectories slightly negative. Either clip or use an adaptive solver.
Reading the saddle in case 4 as a coexistence equilibrium. The coordinates of $E^*$ are real and positive in case 4, but the equilibrium is unstable — almost every trajectory leaves it. Always classify before celebrating.

Summary

Two coupled logistic equations with cross-competition coefficients $\alpha_{12}, \alpha_{21}$ already produce a rich four-way ecology. The geometry is two straight nullclines whose relative positions are decided by two simple comparisons — and those two comparisons completely determine which of four futures the system will live: species 1 wins, species 2 wins, both coexist, or the winner is whoever happens to start ahead. The Jacobian gives a rigorous stability classification at every equilibrium, plain Python integrates the model in 20 lines, and PyTorch lets you scale it to batched and fitted versions without rewriting a single equation.

Next we will close out the chapter with the most famous nonlinear system of all — Lorenz's three-equation atmosphere model — where three coupled ODEs produce something that never settles to an equilibrium and yet is bounded forever: deterministic chaos.