捕食・被食関係のシミュレーション

肉食動物と草食動物のような、食う・食われるの関係にある動物の個体数は振動することが知られています。
これをモデル化した有名な例がロトカ＝ヴォルテラの方程式で、これは

というような微分方程式です。

難しく考えなくても、「食われる側の動物が増えると、エサが増えて食う側の動物が増える」「食う側の動物が増えると、食われる側の動物の数は減る」「食われる側の動物の数が減ると、エサがなくなって食う側の動物が減る」「食う側の動物が減ると、食われる側の動物が増える」という関係にありそうだということは直感的に分かると思います。

ここではシミュレーションで振動現象が起こることを確かめてみようと思います。

設定

• 草原に100匹の草食動物Aと20匹の肉食動物Bが生息している。
• 各動物は毎ターン減っていくHPを持っており、HPがゼロになると死ぬ。
• 草原は100x100のタイル状になっていて、全タイルに草が生えている
• 毎ターンごとにAは0∼1タイル、Bは0∼2タイルランダムに移動する。
• 草は毎ターン1ずつ伸びるが、長さ15以上にはならない。
• 草食動物は今いるタイルに長さ10以上の草が生えていたら草食べてHPを回復する。
• 食べられた草は10短くなる
• Bは現在のタイル上にAが存在した場合はAを食べてHPを回復する。
• AもBも毎ターン10%の確率で子供を産む。

その他細かい設定はソースコードを読んでください

ソースコード

画面描画にはOpenGLを使いました。GLUTによる「手抜き」OpenGL入門が参考になりました。
コンパイルは

g++ simulation.cpp -lglut -lGLU -lGL -lm -std=c++11

という感じで行います。

#include <iostream>
#include <random>
#include <vector>
#include <GL/glut.h>

using namespace std;

// 草原のサイズ
const int width = 100;
const int height = 100;

// 草食動物クラス
class herbivore{
public:
    int x;
    int y;
    int HP;
    herbivore(int pos_x, int pos_y);
};
herbivore::herbivore(int pos_x, int pos_y){
    x = pos_x;
    y = pos_y;
    HP = 5;
}

// 肉食動物クラス
class carnivore{
public:
    int x;
    int y;
    int HP;
    carnivore(int pos_x, int pos_y);
};
carnivore::carnivore(int pos_x, int pos_y){
    x = pos_x;
    y = pos_y;
    HP = 15;
}


vector<herbivore> A;
vector<carnivore> B;
int field[width][height];
int frame = 0;

random_device rd;
mt19937 mt(rd());

void init(void) {
    //背景色の設定
    glClearColor(0.0, 0.0, 0.0, 1.0);

    // 草原の設定
    uniform_int_distribution<int> grass(0, 15);
    for(int x=0; x<width; ++x){
        for(int y=0; y<height; ++y){
            field[x][y] = grass(mt); // 草の初期長さはランダム
        }
    }
    // 草食動物の設定
    uniform_int_distribution<int> rand100(0, width);
    for(int i=0; i<100; ++i){
        // ランダムな初期位置
        int x = rand100(mt);
        int y = rand100(mt);
        A.push_back(herbivore(x, y));
    }

    for(int i=0; i<20; ++i){
        int x = rand100(mt);
        int y = rand100(mt);
        B.push_back(carnivore(x, y));
    }
}

void update(void){
    frame++;

    // 草が伸びる(最長15)
    for(int x=0; x<width; ++x){
        for(int y=0; y<height; ++y){
            if(field[x][y] < 15) field[x][y] += 1;
        }
    }
    uniform_int_distribution<int> move(-1, 1);
    for(int i=0; i<A.size(); ++i){
           A[i].HP -= 2; // HP消耗
           // 移動する 
           A[i].x += move(mt);
           A[i].x = max(0, A[i].x);
           A[i].x = min(width-1, A[i].x);

           A[i].y += move(mt);
           A[i].y = max(0, A[i].y);
           A[i].y = min(height-1, A[i].y);
    }

    // 草を食う
    for(int i=0; i<A.size(); ++i){
        // HP満タンではなく、草が10以上伸びているとき
        if(A[i].HP <= 10 && field[A[i].x][A[i].y] >= 10){
           A[i].HP += 3;                // HP回復
           field[A[i].x][A[i].y] -= 10; // 草が短くなる
        }
    }

    uniform_int_distribution<int> move2(-2, 2);
    for(int i=0; i<B.size(); ++i){
           B[i].HP -= 2; // HP消耗
           // 移動する 
           B[i].x += move2(mt);
           B[i].x = max(0, B[i].x);
           B[i].x = min(width-1, B[i].x);

           B[i].y += move2(mt);
           B[i].y = max(0, B[i].y);
           B[i].y = min(height-1, B[i].y);
    }

    // 肉食動物が草食動物を食う
    for(int i=0; i<B.size(); ++i){
        for(int j=0; j<A.size(); ++j){
            if(A[j].x==B[i].x && A[j].y == B[i].y){
                A.erase(A.begin()+j);
                B[i].HP += 12; // HP回復
                break;
            }
        }
    }

    // 子供を産む
    uniform_int_distribution<int> dice(0, 100);
    int size = A.size();
    for(int i=0; i<size; ++i){
        if(dice(mt) < 10) {             // 確率10%
            A.push_back(herbivore(A[i].x, A[i].y));
        }
    }
    size = B.size();
    for(int i=0; i<size; ++i){
        if(dice(mt) < 10){
            B.push_back(carnivore(B[i].x, B[i].y));
        }
    }

    // 生死判定
    for(vector<herbivore>::iterator it=A.begin(); it!=A.end();){
        if ((*it).HP <= 0) it = A.erase(it);
        else ++it;
    }
    for(vector<carnivore>::iterator it=B.begin(); it!=B.end();){
        if((*it).HP <= 0) it = B.erase(it);
        else ++it;
    }
}

void display(void) {
    update();
    glClear(GL_COLOR_BUFFER_BIT);
    cout << frame << "," << A.size() << "," << B.size() << endl;
    for(int i=0; i<A.size(); ++i){
        double x = (double) (A[i].x*2.0/width-1.0);
        double y = (double) (A[i].y*2.0/width-1.0);
        glBegin(GL_QUADS);
        glColor3d(0, 1.0, 0);
        glVertex2d(x, y);
        glVertex2d(x, y-2.0/height);
        glVertex2d(x+2.0/width, y-2.0/height);
        glVertex2d(x+2.0/width, y);
        glEnd(); 
    }
    for(int i=0; i<B.size(); ++i){
        double x = (double) (B[i].x*2.0/width-1.0);
        double y = (double) (B[i].y*2.0/width-1.0);
        glBegin(GL_QUADS);
        glColor3d(1.0, 0, 0);
        glVertex2d(x, y);
        glVertex2d(x, y-2.0/height);
        glVertex2d(x+2.0/width, y-2.0/height);
        glVertex2d(x+2.0/width, y);
        glEnd(); 
    }
    glutSwapBuffers();
}

void timer(int value){
    glutPostRedisplay();
    glutTimerFunc(30, timer, 0);
}

int main(int argc, char *argv[]) {
    // 700x700の画面
    glutInitWindowSize(700, 700);
    glutInit(&argc, argv);
    // RGBA表示 ダブルバッファリング
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
    glutCreateWindow("OpenGL");
    glutDisplayFunc(display);
    init();
    glutTimerFunc(1, timer , 0);
    glutMainLoop();
}

こんな感じで動きます。赤が肉食動物、緑が草食動物のつもりです。

結果

個体数の推移をグラフにしたらこんな感じでした。

草食動物の個体数の変化に遅れて肉食動物の個体数が変化していることが分かると思います。

追記：
Roombaさんのコメントを参考に相空間上での変化を見てみました。
アニメーションも作りました。

この振動する結果に至るまでにはいろいろとパラメータの調整を行う必要がありました。
うまいパラメータを設定しないと動物が増えすぎて処理が追いつかなくなったり、どちらかが全滅して振動しなかったりといった問題が起こりました。
個体群密度の影響など自然にある制約を導入していないのが複雑なパラメータ調整が必要になった原因かもしれません。

草食動物と肉食動物の個体数の比が20:1になるなどいろいろと面白いことが知られていますが、この実験では振動現象くらいしか再現することができませんでした。パラメータ次第ではそのような現象も再現できるかもしれませんが、残念ながら今のところできていません。

OpenGLで適当に描画すると現在の状態が確認できて面白いので、今後も機会があれば遊んでみようと思います。
今回の記事は以上です。