🇯🇵 日本語 | 🇺🇸 English
広告スペース

C言語でシミュレーションを作る

乱数・物理・セルオートマトン・感染症モデルの4題を完全コード付きで実装する。

シミュレーションとは

現実世界のルールをプログラムで再現し、時間を進めながら状態の変化を観察する手法がシミュレーションです。方程式を手で解くのが難しい現象でも、離散時間で少しずつ状態を更新していけば振る舞いを調べられます。
確率的
ランダムウォーク
乱数で方向を選んで進む粒子。拡散現象のモデル。
物理
放物運動
重力と空気抵抗を微小時間ごとに適用する数値積分。
離散
ライフゲーム
2次元セルの生死をルールで更新するセルオートマトン。
集団
SIR感染症モデル
感受性・感染・回復の人口比率を時間発展させる。

共通する3ステップ

  1. 初期化: 状態(位置、速度、セル配列、人口など)を用意する。
  2. 更新ルール: 1ステップで状態をどう変えるかを1つの関数にまとめる。
  3. ループ: 一定回数(または条件満たすまで)更新を繰り返し、途中経過を出力する。
必要な知識: for/while配列乱数(rand)math.h関数

ランダムウォーク

粒子が毎ステップ ±1 のどちらかにランダムに進む最も単純な確率的シミュレーションです。拡散現象(ブラウン運動)の1次元モデルとしても使われます。

1次元版

#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;
}
ポイント: rand() % 2 で 0 か 1 が出て、それを +1 / -1 に対応付けます。srand を入れないと毎回同じ軌跡が出るので注意。

多数回平均と√N則

1次元ランダムウォークでは、N ステップ後の原点からの距離の平均は √N に比例することが知られています。これを N=1000 回の試行で確認してみます。
#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 距離 = %.3f (理論値 sqrt(%d) = %.3f)\n",
           rms, STEPS, sqrt((double)STEPS));
    return 0;
}
$ gcc walk_rms.c -o walk_rms -lm $ ./walk_rms RMS 距離 = 31.585 (理論値 sqrt(1000) = 31.623)
-lm を忘れずに: math.h の関数(sqrt, pow 等)を使うときはリンカオプション -lm が必要です。

放物運動の数値積分(空気抵抗あり)

空気抵抗がある場合の放物運動は解析的に解きにくいので、微小時間 Δt ごとに速度と位置を更新するオイラー法で近似します。

運動方程式

ax = -k · v · vx
ay = -g - k · v · vy
(v = √(vx² + vy²)、k は空気抵抗係数、g=9.8)

更新ルール(1ステップ)

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

#define G  9.8      // 重力加速度 [m/s^2]
#define K  0.01     // 空気抵抗係数
#define DT 0.01     // 時間刻み [s]

int main(void) {
    // 初期条件: 45度で初速30m/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("# 着地: 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 # 着地: x=57.60m, t=3.88s
K=0 にすると? 空気抵抗なしの理想的な放物運動になり、着地距離は v₀²·sin(2θ)/g = 91.8m になります。K を変えて比較してみましょう。

結果をグラフ化する

出力を CSV としてファイルに保存すれば、Excel / gnuplot / Python などで軌跡を描けます。
$ ./projectile > traj.csv $ gnuplot -e "set datafile separator ','; plot 'traj.csv' using 2:3 with lines"

ライフゲーム(Conway's Game of Life)

2次元グリッドのセルが「生」か「死」の状態を持ち、周囲8マスの状態で次世代が決まるセルオートマトンです。ごく単純なルールから複雑な模様(グライダーなど)が生まれます。

ルール

#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");  // 画面クリア
    printf("世代: %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) {
    // グライダーを配置
    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秒待つ
        step();
    }
    return 0;
}
世代: 5 .............................. ..#........................... ...##......................... ..##.......................... .............................. ..............................
なぜ2つの配列が必要? 次世代のセルは「現世代」の隣接状態から計算します。同じ配列を書き換えながら処理すると、まだ更新していない隣のセルが既に更新されてしまい、ルールが崩れます。next_grid に書いてから一気にコピーするのが定石です。

SIR感染症モデル

感染症の広がりを最も単純化した数理モデルです。人口を次の3群に分けて時間変化を追います。
S
Susceptible
まだ感染していない人(感受性)
I
Infected
現在感染している人
R
Recovered
回復した(または免疫を獲得した)人

更新式(離散時間)

ΔS = -β · S · I / N
ΔI = +β · S · I / N - γ · I
ΔR = +γ · I
β: 感染率(1人の感染者が1日に感染させる平均人数)、γ: 回復率(1/感染期間)、N: 総人口。
#include <stdio.h>

#define DAYS 100

int main(void) {
    double N  = 1000.0;      // 総人口
    double S  = 999.0;       // 感受性
    double I  = 1.0;         // 初期感染者
    double R  = 0.0;         // 回復者
    double beta  = 0.3;      // 感染率
    double gamma = 0.1;      // 回復率 (10日で回復)

    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
基本再生産数 R₀ = β/γ。この例では 0.3/0.1 = 3 で、1人が3人に感染させる強い感染力。β を下げる(接触削減)か γ を上げる(早期治療)と流行曲線は低く緩やかになります。

最大感染者数を求める

// ループの中でピーク値を記録
double peak = 0;
int peak_day = 0;
for (int day = 0; day <= DAYS; day++) {
    if (I > peak) { peak = I; peak_day = day; }
    // 既存の更新処理...
}
printf("ピーク: %d日目に %.0f 人\n", peak_day, peak);

チャレンジ課題

課題1: 2次元ランダムウォーク
東西南北の4方向にランダムに進むプログラムを書け。10,000試行の平均距離が √N に近いか確認する(2次元でも同じスケーリング則が成り立つ)。
課題2: 発射角を最適化
放物運動のコードで発射角を 5°〜85° まで5°刻みで変え、着地距離が最大になる角度を求めよ。空気抵抗ありだと45°より小さい角度が最適になるはず。
課題3: ライフゲームの初期配置
ランダムに50%のセルを生きた状態にして開始せよ。50世代後の生存セル数を表示する。乱数の講座参照。
課題4: ワクチン接種付きSIR
毎日 S のうち v 人がワクチン接種で R に直接移動するようモデルを拡張せよ。v=10 と v=0 で感染ピークがどう変わるか比較する。
課題5: 捕食者-被食者モデル(Lotka-Volterra)
ウサギとキツネの個体数変化を表す古典モデルを実装せよ。
dR/dt = α R - β R F / dF/dt = δ R F - γ F
初期値 R=40, F=9, α=0.1, β=0.02, γ=0.4, δ=0.02 で振動解を観察する。

次のステップ

この記事をシェア
X(Twitter)でシェアはてブ