Python Complex Simulation Code Examples

Simulation lets you model uncertainty, test theories, and explore outcomes that would be impossible or dangerous to observe in the real world. Python's ecosystem of scientific libraries makes it one of the strongest languages for building simulations that range from simple coin flips to full-scale physics engines. This article walks through four distinct categories of simulation with complete, runnable code examples you can adapt to your own projects.

Simulation is one of those areas where Python truly shines. The language's readable syntax, combined with powerful libraries like NumPy, SimPy, and Pymunk, means you can go from a rough idea to a working model in a single afternoon. Whether you are estimating financial risk, modeling a hospital waiting room, simulating how a disease spreads through a population, or testing how objects collide under gravity, Python has the tools to make it happen.

Each example in this article is self-contained and runs with standard Python packages. Copy the code, experiment with the parameters, and watch the behavior change. That hands-on process is the fastest way to build intuition about complex systems.

Monte Carlo Simulations

Monte Carlo simulation uses repeated random sampling to estimate outcomes in systems with uncertainty. The name comes from the famous casino in Monaco, and the core idea is straightforward: if you run enough random trials, the average of your results will converge on the true answer. This technique is widely used in finance, engineering, physics, and risk analysis.

Estimating Pi with Random Points

The classic introductory Monte Carlo example involves throwing darts at a square and counting how many land inside an inscribed circle. The ratio of hits to total throws approximates pi divided by four. Here is the full implementation:

import numpy as np

def estimate_pi(num_samples: int) -> float:
    """Estimate pi using Monte Carlo random sampling."""
    # Generate random (x, y) points in the unit square [0, 1)
    x = np.random.uniform(0, 1, num_samples)
    y = np.random.uniform(0, 1, num_samples)

    # Check which points fall inside the quarter circle
    inside_circle = (x**2 + y**2) <= 1.0

    # Ratio of points inside circle to total points approximates pi/4
    pi_estimate = 4.0 * np.sum(inside_circle) / num_samples
    return pi_estimate

# Run the simulation at increasing sample sizes
for n in [1_000, 10_000, 100_000, 1_000_000]:
    result = estimate_pi(n)
    print(f"Samples: {n:>10,} | Pi estimate: {result:.6f}")

When you run this, you will notice that accuracy improves as sample size increases. With one million samples, the estimate typically lands within a few thousandths of the true value of pi. This demonstrates a fundamental property of Monte Carlo methods: convergence improves proportionally to the square root of the number of samples.

Stock Portfolio Risk Assessment

A more practical Monte Carlo application is modeling potential portfolio outcomes. This simulation generates thousands of possible future scenarios based on historical return distributions, giving you a range of expected results rather than a single number:

import numpy as np

def portfolio_monte_carlo(
    initial_investment: float,
    annual_return_mean: float,
    annual_return_std: float,
    years: int,
    num_simulations: int
) -> dict:
    """Simulate portfolio growth across multiple scenarios."""
    np.random.seed(42)
    all_final_values = np.zeros(num_simulations)

    for i in range(num_simulations):
        portfolio_value = initial_investment
        for year in range(years):
            # Each year's return is drawn from a normal distribution
            annual_return = np.random.normal(annual_return_mean, annual_return_std)
            portfolio_value *= (1 + annual_return)
        all_final_values[i] = portfolio_value

    return {
        "mean": np.mean(all_final_values),
        "median": np.median(all_final_values),
        "std": np.std(all_final_values),
        "percentile_5": np.percentile(all_final_values, 5),
        "percentile_95": np.percentile(all_final_values, 95),
        "probability_of_loss": np.mean(all_final_values < initial_investment)
    }

# Simulate a $100,000 portfolio over 10 years
results = portfolio_monte_carlo(
    initial_investment=100_000,
    annual_return_mean=0.08,
    annual_return_std=0.15,
    years=10,
    num_simulations=50_000
)

print("Portfolio Simulation Results (10 Years, 50,000 Scenarios)")
print(f"  Mean final value:       ${results['mean']:>12,.2f}")
print(f"  Median final value:     ${results['median']:>12,.2f}")
print(f"  5th percentile (worst): ${results['percentile_5']:>12,.2f}")
print(f"  95th percentile (best): ${results['percentile_95']:>12,.2f}")
print(f"  Probability of loss:    {results['probability_of_loss']:.2%}")

Discrete-Event Simulations with SimPy

Discrete-event simulation (DES) models systems where state changes happen at specific points in time rather than continuously. Think of a hospital emergency room, a network router processing packets, or customers moving through a checkout line. SimPy is a process-based DES framework built on Python generators, and it remains one of the go-to libraries for this type of modeling.

Hospital Emergency Room

This simulation models patients arriving at an emergency room with limited doctors. Patients are triaged into priority levels, and higher-priority patients are seen first. The simulation tracks wait times across all priority levels:

import simpy
import random

class EmergencyRoom:
    """Simulates an ER with priority-based triage."""

    def __init__(self, env, num_doctors: int):
        self.env = env
        # PriorityResource lets higher-priority patients jump the queue
        self.doctors = simpy.PriorityResource(env, capacity=num_doctors)
        self.wait_times = []

    def treat_patient(self, patient_id: int, severity: int):
        """Treatment takes longer for more severe cases."""
        treatment_time = random.uniform(10, 30) + (severity * 5)
        yield self.env.timeout(treatment_time)

def patient_arrival(env, er: EmergencyRoom, patient_id: int):
    """A single patient arrives and waits for treatment."""
    arrival_time = env.now

    # Severity: 1 = critical, 2 = serious, 3 = minor
    severity = random.choices([1, 2, 3], weights=[15, 35, 50])[0]

    # Request a doctor with priority (lower number = higher priority)
    with er.doctors.request(priority=severity) as req:
        yield req
        wait_time = env.now - arrival_time
        er.wait_times.append({"patient": patient_id, "severity": severity,
                              "wait": wait_time})

        yield env.process(er.treat_patient(patient_id, severity))

def run_er_simulation(
    num_doctors: int = 4,
    sim_duration: int = 480,
    arrival_rate: float = 5.0
):
    """Run the full ER simulation."""
    random.seed(42)
    env = simpy.Environment()
    er = EmergencyRoom(env, num_doctors)

    def generate_patients():
        patient_id = 0
        while True:
            yield env.timeout(random.expovariate(1.0 / arrival_rate))
            patient_id += 1
            env.process(patient_arrival(env, er, patient_id))

    env.process(generate_patients())
    env.run(until=sim_duration)

    # Analyze results by severity level
    print(f"ER Simulation: {num_doctors} doctors, {sim_duration} min shift")
    print(f"Total patients seen: {len(er.wait_times)}")

    for sev in [1, 2, 3]:
        label = {1: "Critical", 2: "Serious", 3: "Minor"}[sev]
        waits = [w["wait"] for w in er.wait_times if w["severity"] == sev]
        if waits:
            avg_wait = sum(waits) / len(waits)
            max_wait = max(waits)
            print(f"  {label:>10}: avg wait {avg_wait:.1f} min, "
                  f"max wait {max_wait:.1f} min (n={len(waits)})")

run_er_simulation(num_doctors=4, sim_duration=480, arrival_rate=5.0)

Try changing the number of doctors or the arrival rate. You will quickly see how sensitive wait times are to small changes in capacity. This kind of sensitivity analysis is one of the primary reasons simulation is so valuable for operational planning.

Network Packet Processing Queue

Here is a second SimPy example that models a network switch handling incoming packets with a limited buffer. Packets that arrive when the buffer is full are dropped:

import simpy
import random

def packet_generator(env, switch_buffer, stats, packet_rate):
    """Generate packets at a given rate."""
    packet_id = 0
    while True:
        yield env.timeout(random.expovariate(packet_rate))
        packet_id += 1
        size = random.choice([64, 128, 256, 512, 1024, 1500])

        if switch_buffer.level + size <= switch_buffer.capacity:
            yield switch_buffer.put(size)
            stats["received"] += 1
        else:
            stats["dropped"] += 1

def packet_processor(env, switch_buffer, stats, processing_rate):
    """Process packets from the buffer."""
    while True:
        # Wait until there is something to process
        yield env.timeout(1.0 / processing_rate)
        if switch_buffer.level > 0:
            # Process one packet worth of data
            amount = min(switch_buffer.level, 1500)
            yield switch_buffer.get(amount)
            stats["processed"] += 1

def run_network_sim(buffer_size=65536, packet_rate=500, processing_rate=450):
    """Run the network buffer simulation."""
    random.seed(7)
    env = simpy.Environment()
    switch_buffer = simpy.Container(env, capacity=buffer_size, init=0)
    stats = {"received": 0, "dropped": 0, "processed": 0}

    env.process(packet_generator(env, switch_buffer, stats, packet_rate))
    env.process(packet_processor(env, switch_buffer, stats, processing_rate))

    env.run(until=60)  # Simulate 60 seconds

    total = stats["received"] + stats["dropped"]
    drop_rate = stats["dropped"] / total * 100 if total > 0 else 0
    print(f"Network Simulation (60 seconds)")
    print(f"  Buffer size:    {buffer_size:,} bytes")
    print(f"  Packets in:     {stats['received']:,}")
    print(f"  Packets dropped:{stats['dropped']:,}")
    print(f"  Packets out:    {stats['processed']:,}")
    print(f"  Drop rate:      {drop_rate:.2f}%")

run_network_sim()

This type of simulation is a staple in network engineering. By varying buffer sizes and processing rates, you can identify the point at which packet loss becomes unacceptable and right-size your hardware accordingly.

Agent-Based Modeling from Scratch

Agent-based modeling (ABM) simulates individual actors (agents) that follow simple rules and interact with each other and their environment. Complex, emergent behaviors arise from these simple interactions. ABM is used extensively in epidemiology, economics, ecology, and social science research.

Disease Spread in a Population

This example simulates a classic SIR (Susceptible-Infected-Recovered) epidemic model. Each person in the population is an individual agent who can move, encounter other agents, and potentially transmit or recover from an infection:

import random
from dataclasses import dataclass, field
from enum import Enum

class State(Enum):
    SUSCEPTIBLE = "S"
    INFECTED = "I"
    RECOVERED = "R"

@dataclass
class Person:
    id: int
    x: float
    y: float
    state: State = State.SUSCEPTIBLE
    days_infected: int = 0

class EpidemicSimulation:
    """SIR agent-based epidemic model."""

    def __init__(
        self,
        population_size: int = 500,
        grid_size: float = 100.0,
        infection_radius: float = 3.0,
        infection_prob: float = 0.3,
        recovery_days: int = 14,
        initial_infected: int = 5,
    ):
        self.grid_size = grid_size
        self.infection_radius = infection_radius
        self.infection_prob = infection_prob
        self.recovery_days = recovery_days

        # Create population with random positions
        self.people = []
        for i in range(population_size):
            p = Person(
                id=i,
                x=random.uniform(0, grid_size),
                y=random.uniform(0, grid_size),
            )
            self.people.append(p)

        # Infect a few people to start
        for p in random.sample(self.people, initial_infected):
            p.state = State.INFECTED
            p.days_infected = 1

        self.history = []

    def _move_agents(self):
        """Each person takes a random step."""
        for p in self.people:
            p.x = max(0, min(self.grid_size,
                             p.x + random.gauss(0, 2.0)))
            p.y = max(0, min(self.grid_size,
                             p.y + random.gauss(0, 2.0)))

    def _spread_infection(self):
        """Infected agents can transmit to nearby susceptible agents."""
        infected = [p for p in self.people if p.state == State.INFECTED]
        susceptible = [p for p in self.people if p.state == State.SUSCEPTIBLE]

        for s in susceptible:
            for i in infected:
                dist = ((s.x - i.x)**2 + (s.y - i.y)**2) ** 0.5
                if dist <= self.infection_radius:
                    if random.random() < self.infection_prob:
                        s.state = State.INFECTED
                        s.days_infected = 1
                        break  # Only need one successful transmission

    def _update_recovery(self):
        """Infected agents recover after enough days."""
        for p in self.people:
            if p.state == State.INFECTED:
                p.days_infected += 1
                if p.days_infected > self.recovery_days:
                    p.state = State.RECOVERED

    def _record_state(self, day: int):
        """Snapshot current counts."""
        counts = {"day": day, "S": 0, "I": 0, "R": 0}
        for p in self.people:
            counts[p.state.value] += 1
        self.history.append(counts)

    def run(self, days: int = 100):
        """Run the simulation for a given number of days."""
        for day in range(days):
            self._record_state(day)
            self._move_agents()
            self._spread_infection()
            self._update_recovery()

            # Stop early if no one is infected
            if all(p.state != State.INFECTED for p in self.people):
                self._record_state(day + 1)
                break

        return self.history

# Run the epidemic simulation
random.seed(123)
sim = EpidemicSimulation(
    population_size=500,
    infection_radius=3.0,
    infection_prob=0.3,
    recovery_days=14,
    initial_infected=5
)
history = sim.run(days=120)

# Print summary at key intervals
print("Day |   S   |   I   |   R")
print("-" * 33)
for record in history:
    if record["day"] % 10 == 0 or record == history[-1]:
        print(f"{record['day']:>3} | {record['S']:>5} | "
              f"{record['I']:>5} | {record['R']:>5}")

What makes this different from a simple mathematical SIR model is that each agent has a unique position and moves independently. This spatial component can produce realistic effects like localized outbreaks and geographic clustering that differential equation models cannot capture.

Physics-Based Simulations with Pymunk

Pymunk is a 2D rigid-body physics library built on top of the Chipmunk2D engine. It handles collision detection, gravity, friction, joints, and constraints out of the box. If you have ever wanted to simulate billiard balls bouncing, a ragdoll falling, or a chain of linked objects swinging under gravity, Pymunk is where you start.

Bouncing Balls in a Box

This example creates a bounded box and drops several balls into it. Each ball has mass, elasticity, and friction properties. The simulation tracks the kinetic energy of the system as it evolves over time:

import pymunk

def create_walls(space, width, height, thickness=5):
    """Create static walls forming a box."""
    walls = [
        # floor
        pymunk.Segment(space.static_body, (0, 0), (width, 0), thickness),
        # ceiling
        pymunk.Segment(space.static_body, (0, height), (width, height), thickness),
        # left wall
        pymunk.Segment(space.static_body, (0, 0), (0, height), thickness),
        # right wall
        pymunk.Segment(space.static_body, (width, 0), (width, height), thickness),
    ]
    for wall in walls:
        wall.elasticity = 0.8
        wall.friction = 0.5
    space.add(*walls)

def add_ball(space, position, radius=15, mass=1.0):
    """Add a dynamic ball to the simulation."""
    moment = pymunk.moment_for_circle(mass, 0, radius)
    body = pymunk.Body(mass, moment)
    body.position = position
    shape = pymunk.Circle(body, radius)
    shape.elasticity = 0.9
    shape.friction = 0.4
    space.add(body, shape)
    return body

def calculate_kinetic_energy(space):
    """Sum kinetic energy of all dynamic bodies."""
    total_ke = 0.0
    for body in space.bodies:
        if body.body_type == pymunk.Body.DYNAMIC:
            v = body.velocity
            speed_sq = v.x**2 + v.y**2
            total_ke += 0.5 * body.mass * speed_sq
    return total_ke

def run_physics_sim(num_balls=8, duration_seconds=10, dt=1/120):
    """Run a bouncing balls physics simulation."""
    space = pymunk.Space()
    space.gravity = (0, -981)  # Gravity pointing down

    width, height = 600, 400
    create_walls(space, width, height)

    # Add balls at staggered positions
    import random
    random.seed(99)
    balls = []
    for i in range(num_balls):
        x = random.uniform(50, width - 50)
        y = random.uniform(height * 0.5, height - 50)
        balls.append(add_ball(space, (x, y), radius=15, mass=1.0))

    steps = int(duration_seconds / dt)
    snapshots = []

    for step in range(steps):
        space.step(dt)
        if step % 120 == 0:  # Record every second
            ke = calculate_kinetic_energy(space)
            t = step * dt
            snapshots.append((t, ke))

    # Print energy over time
    print(f"Bouncing Balls Simulation ({num_balls} balls, {duration_seconds}s)")
    print(f"{'Time (s)':>10} | {'Kinetic Energy':>15}")
    print("-" * 30)
    for t, ke in snapshots:
        print(f"{t:>10.1f} | {ke:>15.1f}")

run_physics_sim(num_balls=8, duration_seconds=10)

You will notice that kinetic energy decreases over time as collisions dissipate energy (the elasticity is set below 1.0). Setting elasticity to exactly 1.0 creates perfectly elastic collisions where total kinetic energy is conserved, while values below 1.0 model real-world materials that absorb energy on impact.

Combining Techniques: A Network Intrusion Simulator

Real-world simulations often blend multiple techniques. This final example combines Monte Carlo randomness with agent-based logic to simulate network intrusion attempts against a defended system. Attackers probe for vulnerabilities while defenders detect and block threats. Each side has configurable parameters that affect the outcome:

import random
from dataclasses import dataclass

@dataclass
class AttackEvent:
    timestamp: float
    attack_type: str
    source_ip: str
    target_port: int
    detected: bool
    blocked: bool

class NetworkIntrusionSim:
    """Simulates attack and defense interactions on a network."""

    ATTACK_TYPES = {
        "port_scan":       {"frequency": 0.35, "stealth": 0.7, "severity": 2},
        "brute_force":     {"frequency": 0.25, "stealth": 0.3, "severity": 5},
        "sql_injection":   {"frequency": 0.15, "stealth": 0.5, "severity": 8},
        "phishing":        {"frequency": 0.15, "stealth": 0.8, "severity": 7},
        "zero_day":        {"frequency": 0.10, "stealth": 0.9, "severity": 10},
    }

    def __init__(
        self,
        detection_capability: float = 0.6,
        response_speed: float = 0.7,
        num_target_ports: int = 20,
    ):
        self.detection_capability = detection_capability
        self.response_speed = response_speed
        self.target_ports = list(range(1, num_target_ports + 1))
        self.events = []

    def _generate_source_ip(self) -> str:
        return ".".join(str(random.randint(1, 254)) for _ in range(4))

    def _choose_attack_type(self) -> str:
        types = list(self.ATTACK_TYPES.keys())
        weights = [self.ATTACK_TYPES[t]["frequency"] for t in types]
        return random.choices(types, weights=weights)[0]

    def _attempt_detection(self, attack_type: str) -> bool:
        stealth = self.ATTACK_TYPES[attack_type]["stealth"]
        detection_chance = self.detection_capability * (1 - stealth * 0.5)
        return random.random() < detection_chance

    def _attempt_block(self, detected: bool) -> bool:
        if not detected:
            return False
        return random.random() < self.response_speed

    def run(self, duration_hours: float = 24, attack_rate: float = 10):
        """Simulate attacks over a given duration."""
        current_time = 0.0
        duration_minutes = duration_hours * 60

        while current_time < duration_minutes:
            # Time until next attack follows exponential distribution
            interval = random.expovariate(attack_rate / 60)
            current_time += interval

            if current_time >= duration_minutes:
                break

            attack_type = self._choose_attack_type()
            source_ip = self._generate_source_ip()
            target_port = random.choice(self.target_ports)

            detected = self._attempt_detection(attack_type)
            blocked = self._attempt_block(detected)

            event = AttackEvent(
                timestamp=current_time,
                attack_type=attack_type,
                source_ip=source_ip,
                target_port=target_port,
                detected=detected,
                blocked=blocked,
            )
            self.events.append(event)

        return self.events

    def report(self):
        """Print a summary of the simulation results."""
        total = len(self.events)
        if total == 0:
            print("No attacks recorded.")
            return

        detected = sum(1 for e in self.events if e.detected)
        blocked = sum(1 for e in self.events if e.blocked)
        undetected = total - detected

        print(f"Network Intrusion Simulation Report")
        print(f"{'=' * 45}")
        print(f"  Total attacks:    {total}")
        print(f"  Detected:         {detected} ({detected/total:.1%})")
        print(f"  Blocked:          {blocked} ({blocked/total:.1%})")
        print(f"  Undetected:       {undetected} ({undetected/total:.1%})")
        print()

        print("  Per Attack Type:")
        for atype in self.ATTACK_TYPES:
            type_events = [e for e in self.events if e.attack_type == atype]
            if not type_events:
                continue
            count = len(type_events)
            det = sum(1 for e in type_events if e.detected)
            blk = sum(1 for e in type_events if e.blocked)
            sev = self.ATTACK_TYPES[atype]["severity"]
            print(f"    {atype:<18} count={count:<4} detected={det:<4} "
                  f"blocked={blk:<4} severity={sev}")

# Run the simulation
random.seed(2026)
sim = NetworkIntrusionSim(
    detection_capability=0.65,
    response_speed=0.75,
    num_target_ports=25
)
sim.run(duration_hours=24, attack_rate=12)
sim.report()

This simulator shows how you can layer randomness (Monte Carlo sampling for attack types, timing, and detection outcomes) on top of agent-like behavior (each attack follows a logical chain of detection and response). The configurable parameters let you model different security postures and see how changes to detection capability or response speed shift the overall risk profile.

Key Takeaways

1. Monte Carlo methods estimate uncertain outcomes through repeated random sampling. They are effective whenever you have a known or assumed probability distribution and want to understand the range of possible results. NumPy's vectorized operations make these simulations fast even at millions of iterations.
2. Discrete-event simulation models systems where things happen at specific moments. SimPy's generator-based processes and built-in resource types like PriorityResource and Container let you model queues, buffers, and shared resources without building queue logic from scratch.
3. Agent-based models produce emergent behavior from simple individual rules. By giving each agent a position, state, and set of interactions, you can observe macro-level phenomena like epidemic curves and clustering that arise naturally from the micro-level mechanics.
4. Physics engines handle the math of collisions, gravity, and constraints. Pymunk provides a tested, optimized foundation so you can focus on setting up the scenario rather than implementing collision detection algorithms from scratch.
5. Combining simulation techniques creates more realistic models. Real systems rarely fit neatly into one category. Blending Monte Carlo randomness with agent-based logic and event-driven scheduling produces simulations that capture the complexity of actual environments.

Simulation is a skill that compounds. Each model you build sharpens your ability to decompose complex systems into programmable components. Start with the examples here, adjust the parameters, break them, and rebuild them with your own variations. That iterative process is how you go from copying code to designing your own simulations from a blank file.