Construa Simulações em C

Quatro programas completos: random walk, movimento de projétil, Game of Life e o modelo epidêmico SIR.

O que é simulação?

Uma simulação codifica as regras de um sistema do mundo real e avança no tempo passo a passo, deixando você observar o estado evoluir. É como você estuda comportamentos difíceis de resolver analiticamente.
Estocástico
Random walk
Uma partícula escolhe uma direção aleatória a cada passo. Um modelo clássico de difusão.
Física
Movimento de projétil
Aplica gravidade e arrasto do ar em pequenos passos de tempo (integração de Euler).
Discreto
Game of Life
Um autômato celular 2D com regras simples de nascimento e morte.
População
Epidemia SIR
Acompanha como populações suscetíveis, infectadas e recuperadas mudam com o tempo.

Três passos comuns

  1. Inicializar: monte o estado inicial (posições, velocidades, grade, populações, ...).
  2. Regra de atualização: uma função que avança o estado por um único passo.
  3. Laço: chame a atualização um número fixo de vezes (ou até uma condição ser atingida), imprimindo valores intermediários ao longo do caminho.
Pré-requisitos: for/while, arrays, rand(), math.h, funções.

Random walk

A cada passo, a partícula se move +1 ou -1 com probabilidade igual. É a simulação estocástica mais simples e o modelo 1D clássico de difusão.

Versão 1D

#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;
}
Ponto-chave: rand() % 2 dá 0 ou 1, que mapeamos para +1 / -1. Sem srand, você recebe a mesma trajetória toda execução.

Muitas tentativas e a regra do √N

Para um random walk 1D, a distância média até a origem após N passos escala como √N. Vamos verificar isso com 10.000 tentativas.
#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("distância RMS = %.3f (teoria sqrt(%d) = %.3f)\n",
           rms, STEPS, sqrt((double)STEPS));
    return 0;
}
$ gcc walk_rms.c -o walk_rms -lm $ ./walk_rms distância RMS = 31.585 (teoria sqrt(1000) = 31.623)
Não esqueça do -lm: usar funções de math.h (sqrt, pow, …) exige a flag de linker -lm.

Movimento de projétil com arrasto (integração de Euler)

Com arrasto do ar não há solução fechada, então usamos integração de Euler: atualizamos velocidade e posição a cada pequeno passo de tempo Δt.

Equações do movimento

ax = -k · v · vx
ay = -g - k · v · vy
(v = √(vx² + vy²), k = coef. de arrasto, g = 9,8)

Um passo de atualização

vx += ax · Δt
vy += ay · Δt
x += vx · Δt
y += vy · Δt 
#include <stdio.h>
#include <math.h>

#define G  9.8      // gravidade [m/s^2]
#define K  0.01     // coeficiente de arrasto
#define DT 0.01     // passo de tempo [s]

int main(void) {
    // Estado inicial: 45 graus, 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("# pousou: 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 # pousou: x=57.60m, t=3.88s
Teste K=0: você obtém a parábola ideal com alcance v₀²·sin(2θ)/g ≈ 91,8 m. Varie K para ver como o arrasto encurta o lançamento.

Plotar o resultado

Salve a saída como CSV e plote com Excel, gnuplot ou Python.
$ ./projectile > traj.csv $ gnuplot -e "set datafile separator ','; plot 'traj.csv' using 2:3 with lines"

Game of Life de Conway

Células numa grade 2D estão vivas ou mortas. O próximo estado de cada célula depende da contagem de vizinhos vivos. Padrões complexos (como o glider) emergem de regras surpreendentemente simples.

Regras

#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");  // limpa a tela
    printf("ger: %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) {
    // coloca um 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);   // pausa de 0,15 s
        step();
    }
    return 0;
}
ger: 5 .............................. ..#........................... ...##......................... ..##.......................... .............................. ..............................
Por que duas grades? O próximo estado de cada célula precisa ser calculado a partir da vizinhança atual. Se você escreve na mesma grade que está lendo, vizinhos já atualizados corrompem a regra. Escreva em next_grid e depois faça memcpy de uma só vez.

Modelo epidêmico SIR

O modelo matemático mais simples de propagação de doença divide a população em três compartimentos:
S
Suscetíveis
Ainda não infectados.
I
Infectados
Atualmente infecciosos.
R
Recuperados
Recuperados (e imunes).

Atualização em tempo discreto

ΔS = -β · S · I / N
ΔI = +β · S · I / N - γ · I
ΔR = +γ · I
β: taxa de infecção, γ: taxa de recuperação (1/duração), N: população total.
#include <stdio.h>

#define DAYS 100

int main(void) {
    double N  = 1000.0;      // população total
    double S  = 999.0;       // suscetíveis
    double I  = 1.0;         // infectados iniciais
    double R  = 0.0;         // recuperados
    double beta  = 0.3;      // taxa de infecção
    double gamma = 0.1;      // taxa de recuperação (doença de 10 dias)

    printf("dia,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 dia,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
Número básico de reprodução R₀ = β/γ. Aqui, 0,3/0,1 = 3 — em média, cada infectado transmite para três outros. Reduzir β (menos contato) ou aumentar γ (recuperação mais rápida) achata a curva.

Encontre o pico

// rastreia o máximo de I dentro do laço
double peak = 0;
int peak_day = 0;
for (int day = 0; day <= DAYS; day++) {
    if (I > peak) { peak = I; peak_day = day; }
    // ... atualização regular
}
printf("pico no dia %d: %.0f pessoas\n", peak_day, peak);

Desafios

Desafio 1: Random walk 2D
Escolha uma de quatro direções (N/S/L/O) a cada passo. Rode 10.000 tentativas e confira que a distância média ainda escala como √N.
Desafio 2: Ângulo de lançamento ótimo
Rode a simulação de projétil para ângulos de 5° a 85° em passos de 5°, e encontre o ângulo que maximiza o alcance. Com arrasto, o ótimo fica abaixo de 45°.
Desafio 3: Grade inicial aleatória
Inicialize a grade do Life aleatoriamente com 50% de células vivas. Informe quantas estão vivas após 50 gerações. Veja a aula de números aleatórios.
Desafio 4: SIR com vacinação
Adicione um termo que move v pessoas por dia direto de S para R (vacinação). Compare v=10 com v=0 — como o pico muda?
Desafio 5: Lotka-Volterra (predador-presa)
Coelhos (R) e raposas (F):
dR/dt = α R - β R F, dF/dt = δ R F - γ F
Comece com R=40, F=9, α=0,1, β=0,02, γ=0,4, δ=0,02 e veja a oscilação.

Próximos passos

Quiz de Revisão

Teste seu entendimento desta aula!

Q1. Qual a fórmula básica para atualizar a posição com um passo de tempo dt?

x_{n+1} = x_n + v * dt
x_{n+1} = x_n * v
x_{n+1} = v - dt

Esta é a forma mais simples do Euler method — a posição avança pela velocidade v vezes o passo de tempo dt.

Q2. Qual a maneira padrão de inicializar o gerador para que cada execução produza uma sequência diferente?

srand(time(NULL))
srand(0)
Nenhuma inicialização é necessária

time(NULL) retorna o horário atual em segundos, dando uma semente diferente a cada execução.

Q3. Qual é a ideia central do método de Monte Carlo?

Estimar um valor estatisticamente a partir de muitas tentativas aleatórias
Resolver exatamente uma equação diferencial
Achar o mínimo ordenando

Ele depende de convergência probabilística — estimar π, por exemplo. A precisão melhora conforme o número de tentativas cresce.

Compartilhe este artigo
Compartilhar no X