Why avoid SciPy in this implementation?

The populations are stepped forward with Euler integration inside SimPy so it can run even when SciPy is unavailable.

How does the time_step affect accuracy?

A larger dt leads to bigger numerical errors since Euler integration is only first order.

Where is env.timeout used?

Each discrete step yields env.timeout(dt) to advance the SimPy clock.

Predator-Prey (Lotka–Volterra, SimPy)

Discrete-time SimPy approximation of the predator–prey model.

Level:Intermediate

population ecosystem nonlinear-dynamics

Stocks:prey, predator
Flows:prey_births, predations, predator_deaths
Feedback Loops:predation cycle
Probes:prey, predator

Run

FAQ

Why avoid SciPy in this implementation?: The populations are stepped forward with Euler integration inside SimPy so it can run even when SciPy is unavailable.
How does the time_step affect accuracy?: A larger dt leads to bigger numerical errors since Euler integration is only first order.
Where is env.timeout used?: Each discrete step yields env.timeout(dt) to advance the SimPy clock.

simulation.py

Predator and prey via SimPy

This example approximates the Lotka--Volterra equations in discrete time using SimPy's event loop. It intentionally avoids SciPy's ODE solver so it can run in environments where SciPy is unavailable.

Limitations: because we update the populations at a fixed time step, results depend heavily on the chosen time_step and may drift from the true continuous solution. A stiff system or a large step size can cause the simulation to diverge or behave unrealistically.


from tys import probe, progress


def simulate(cfg: dict):
    """Simulate predator and prey populations with SimPy."""

    import simpy
    env = simpy.Environment()

Parameters

    prey = float(cfg["initial_prey"])
    predator = float(cfg["initial_predator"])
    alpha = cfg["alpha"]
    beta = cfg["beta"]
    delta = cfg["delta"]
    gamma = cfg["gamma"]
    sim_time = float(cfg["sim_time"])
    dt = float(cfg.get("time_step", 0.01))

    assert all(p > 0 for p in (alpha, beta, delta, gamma)), "Rates must be positive"
    assert dt > 0 and sim_time > 0, "Time parameters must be positive"

    steps = int(sim_time / dt)
    done = env.event()

Discretely integrate the differential equations using Euler's method. This simple integrator is convenient but inaccurate for larger dt.

    def dynamics():
        nonlocal prey, predator
        for i in range(steps):
            prey_births = alpha * prey
            predations = beta * prey * predator
            predator_births = delta * prey * predator
            predator_deaths = gamma * predator

            prey += (prey_births - predations) * dt
            predator += (predator_births - predator_deaths) * dt

            prey = max(prey, 0.0)
            predator = max(predator, 0.0)

            probe("prey", env.now, prey)
            probe("predator", env.now, predator)
            progress(int(100 * (i + 1) / steps))
            yield env.timeout(dt)
        done.succeed({
            "final_prey": prey,
            "final_predator": predator,
        })

    env.process(dynamics())
    env.run(until=done)
    return done.value


def requirements():
    return {
        "builtin": ["micropip", "pyyaml"],
        "external": ["simpy==4.1.1"],
    }

Default.yaml

initial_prey: 10
initial_predator: 5
alpha: 1.5
beta: 1.0
delta: 1.0
gamma: 3.0
sim_time: 20
time_step: 0.01

Charts (Default)

prey

Samples	2000 @ 0.00–19.99
Values	min 0.05, mean 3.12, median 0.89, max 18.58, σ 4.37

predator

Samples	2000 @ 0.00–19.99
Values	min 0.00, mean 1.63, median 0.08, max 15.63, σ 3.22

Final Results (Default)

Metric	Value
final_prey	0.24
final_predator	10.79

Predator-Prey (Lotka–Volterra, SciPy)Machines Simulation

Open in Playground