Build Simulations in C

What is simulation?

A simulation reproduces the rules of a real-world system in a program and advances time step by step, letting you watch the state evolve. It lets you study behavior that is hard to solve analytically.

Stochastic

Random walk

A particle picks a direction at random each step. Models diffusion.

Physics

Projectile motion

Apply gravity and air drag at small time intervals (Euler integration).

Discrete

Game of Life

A 2D cellular automaton with simple birth/death rules.

Population

SIR epidemic

Track how susceptible, infected, and recovered people change over time.

Three common steps

Initialize: set up the state (positions, velocities, grid, populations, ...).
Update rule: a single function that advances the state by one step.
Loop: call the update a fixed number of times (or until some condition) and print intermediate values.

Prerequisites: for/while, arrays, rand(), math.h, functions.

Random walk

At each step the particle moves +1 or -1 with equal probability. It is the simplest stochastic simulation and the classic 1D model of diffusion.

1D version

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define STEPS 100

int main(void) {
    srand((unsigned)time(NULL));
    int x = 0;
    printf("step,x\n");
    for (int t = 0; t <= STEPS; t++) {
        printf("%d,%d\n", t, x);
        x += (rand() % 2 == 0) ? 1 : -1;
    }
    return 0;
}

Key point: rand() % 2 yields 0 or 1, which we map to +1 / -1. Without srand you will get the same trajectory every run.

Many trials and the √N rule

For a 1D random walk, the average distance from origin after N steps scales as √N. Let's verify with 10,000 trials.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define TRIALS 10000
#define STEPS  1000

int main(void) {
    srand((unsigned)time(NULL));
    double sum_sq = 0.0;
    for (int tr = 0; tr < TRIALS; tr++) {
        int x = 0;
        for (int t = 0; t < STEPS; t++) {
            x += (rand() % 2) ? 1 : -1;
        }
        sum_sq += (double)x * x;
    }
    double rms = sqrt(sum_sq / TRIALS);
    printf("RMS distance = %.3f (theory sqrt(%d) = %.3f)\n",
           rms, STEPS, sqrt((double)STEPS));
    return 0;
}

$ gcc walk_rms.c -o walk_rms -lm
$ ./walk_rms
RMS distance = 31.585 (theory sqrt(1000) = 31.623)

Don't forget -lm: using math.h functions (sqrt, pow, ...) requires the linker flag -lm.

Projectile motion with air drag (Euler integration)

With air drag the trajectory has no closed-form solution, so we use Euler integration: update velocity and position every small time step Δt.

Equations of motion

    ax = -k · v · vx

    ay = -g - k · v · vy

    (v = √(vx² + vy²), k = drag coeff, g = 9.8)

One update step

    vx += ax · Δt

    vy += ay · Δt

    x  += vx · Δt

    y  += vy · Δt

#include <stdio.h>
#include <math.h>

#define G  9.8      // gravity [m/s^2]
#define K  0.01     // drag coefficient
#define DT 0.01     // time step [s]

int main(void) {
    // Initial state: 45 deg, v0=30 m/s
    double x = 0, y = 0;
    double angle = 45.0 * M_PI / 180.0;
    double v0 = 30.0;
    double vx = v0 * cos(angle);
    double vy = v0 * sin(angle);

    printf("t,x,y\n");
    double t = 0;
    while (y >= 0) {
        printf("%.2f,%.3f,%.3f\n", t, x, y);
        double v  = sqrt(vx * vx + vy * vy);
        double ax = -K * v * vx;
        double ay = -G - K * v * vy;
        vx += ax * DT;
        vy += ay * DT;
        x  += vx * DT;
        y  += vy * DT;
        t  += DT;
    }
    printf("# landed: x=%.2fm, t=%.2fs\n", x, t);
    return 0;
}

$ gcc projectile.c -o projectile -lm
$ ./projectile | head -5
t,x,y
0.00,0.000,0.000
0.01,0.212,0.212
0.02,0.424,0.423
0.03,0.636,0.635
$ ./projectile | tail -3
3.86,57.289,0.099
3.87,57.445,0.029
# landed: x=57.60m, t=3.88s

Try K=0: you get the ideal parabola with range v₀²·sin(2θ)/g = 91.8 m. Vary K to see how drag shortens the throw.

Plot the result

Save the output as CSV and plot it with Excel / gnuplot / Python.

$ ./projectile > traj.csv
$ gnuplot -e "set datafile separator ','; plot 'traj.csv' using 2:3 with lines"

Conway's Game of Life

Cells on a 2D grid are either alive or dead. The next generation is computed from the count of live neighbors. Complex patterns (like the glider) emerge from very simple rules.

Rules

Live cell: survives with 2 or 3 live neighbors, dies otherwise.
Dead cell: becomes alive with exactly 3 live neighbors.

#include <stdio.h>
#include <string.h>
#include <unistd.h>     // usleep (POSIX)

#define W 30
#define H 15
#define GENS 40

int grid[H][W];
int next_grid[H][W];

int count_neighbors(int y, int x) {
    int c = 0;
    for (int dy = -1; dy <= 1; dy++) {
        for (int dx = -1; dx <= 1; dx++) {
            if (dy == 0 && dx == 0) continue;
            int ny = y + dy, nx = x + dx;
            if (ny < 0 || ny >= H || nx < 0 || nx >= W) continue;
            c += grid[ny][nx];
        }
    }
    return c;
}

void step(void) {
    for (int y = 0; y < H; y++) {
        for (int x = 0; x < W; x++) {
            int n = count_neighbors(y, x);
            if (grid[y][x]) next_grid[y][x] = (n == 2 || n == 3);
            else           next_grid[y][x] = (n == 3);
        }
    }
    memcpy(grid, next_grid, sizeof(grid));
}

void print_grid(int gen) {
    printf("\033[H\033[2J");  // clear screen
    printf("gen: %d\n", gen);
    for (int y = 0; y < H; y++) {
        for (int x = 0; x < W; x++) {
            putchar(grid[y][x] ? '#' : '.');
        }
        putchar('\n');
    }
}

int main(void) {
    // place a glider
    grid[1][2] = grid[2][3] = 1;
    grid[3][1] = grid[3][2] = grid[3][3] = 1;

    for (int g = 0; g < GENS; g++) {
        print_grid(g);
        usleep(150000);   // 0.15 s pause
        step();
    }
    return 0;
}

gen: 5
..............................
..#...........................
...##.........................
..##..........................
..............................
..............................

Why two grids? The next state of each cell must be computed from the current neighborhood. If you write into the same grid while reading it, already-updated neighbors corrupt the rule. Write to next_grid, then memcpy in one shot.

SIR epidemic model

The simplest mathematical model of disease spread splits the population into three groups:

Susceptible

Not yet infected.

Infected

Currently infectious.

Recovered

Recovered (and immune).

Discrete-time update

    ΔS = -β · S · I / N

    ΔI = +β · S · I / N - γ · I

    ΔR = +γ · I

β: infection rate, γ: recovery rate (1/duration), N: total population.

#include <stdio.h>

#define DAYS 100

int main(void) {
    double N  = 1000.0;      // total population
    double S  = 999.0;       // susceptible
    double I  = 1.0;         // initial infected
    double R  = 0.0;         // recovered
    double beta  = 0.3;      // infection rate
    double gamma = 0.1;      // recovery rate (10-day illness)

    printf("day,S,I,R\n");
    for (int day = 0; day <= DAYS; day++) {
        printf("%d,%.1f,%.1f,%.1f\n", day, S, I, R);
        double dS = -beta * S * I / N;
        double dI =  beta * S * I / N - gamma * I;
        double dR =  gamma * I;
        S += dS;
        I += dI;
        R += dR;
    }
    return 0;
}

$ ./sir
day,S,I,R
0,999.0,1.0,0.0
10,995.6,4.1,0.3
30,880.7,87.0,32.3
50,306.4,311.1,382.5
70,69.5,81.8,848.7
100,56.2,8.7,935.1

Basic reproduction number R₀ = β/γ. Here 0.3/0.1 = 3 — one infected person spreads to three others on average. Lowering β (less contact) or raising γ (faster recovery) flattens the curve.

Find the peak

// track the maximum of I inside the loop
double peak = 0;
int peak_day = 0;
for (int day = 0; day <= DAYS; day++) {
    if (I > peak) { peak = I; peak_day = day; }
    // ... regular update
}
printf("peak on day %d: %.0f people\n", peak_day, peak);

Challenges

Challenge 1: 2D random walk

Pick one of 4 directions (N/S/E/W) each step. Run 10,000 trials and verify the mean distance still scales like √N.

Challenge 2: Optimal launch angle

Run the projectile simulation for angles 5°–85° in 5° steps and find the angle that maximizes range. With drag the optimum is below 45°.

Challenge 3: Random initial grid

Initialize the Life grid randomly with 50% live cells. Report how many cells are alive after 50 generations. See the random-number lesson.

Challenge 4: SIR with vaccination

Add a term that moves v people per day directly from S to R (vaccination). Compare v=10 to v=0 — how does the peak change?

Challenge 5: Lotka-Volterra (predator-prey)

Rabbits (R) and foxes (F):
dR/dt = α R - β R F, dF/dt = δ R F - γ F
Start with R=40, F=9, α=0.1, β=0.02, γ=0.4, δ=0.02 and watch the oscillation.

Next steps

Higher-accuracy integration → try Runge-Kutta (see numerical computing)
Draw graphics → experiment with SDL2 or raylib
Speed up heavy simulations → parallelize with OpenMP