🙇‍♂️

【トイモデルで遊びたかった】Sandpile Model【1回】

に公開
数理モデル
sandpile
トイモデル
冪分布
tech

はじめに

みなさんトイモデル、やってますか？
私は学部生のころEhrenfest's Urn Modelを触ってからというもの、そこそこご無沙汰です。

トイモデルとは、数理モデルの一種で、現実の現象を完全に再現することよりも、モデルをできるだけ簡素化することによって、現象の本質を捉えようとするモデルと私は思っています。

今回は雪崩や連鎖的な崩壊をみることができるSandpile Modelを触ってみようと思います。

Sandpile Modelって？

手順

トイモデルだけあってとても単純です。

  1. 一辺L四方のマス目付きボードを用意します。
  2. ランダムに1つマスを選び、そこに一つの石を配置することを繰り返す。
  3. ある閾値Z以上の石が同じマスに積まれた時、積まれている石を上下左右4マスに配ります。
  4. 石をもらった上下左右のマスでも閾値Zを超えた場合は、以降連鎖させていく。

基本これだけです。
今回はボードの外に出た石は消えるようになっていますが、周期境界にするパターンや、ある一定方向に重点的に配る石を増やすパターンなど、いくつか派生したものがあります。

Z=4つ乗っかると崩れる

実装

import numpy as np
from collections import deque
import matplotlib.pyplot as plt

class Sandpile:
    def __init__(self, L=64, mode="btw", seed=0):
        self.L = L
        self.mode = mode
        self.rng = np.random.default_rng(seed)
        self.h = np.zeros((L, L), dtype=np.int32)

        if mode == "btw":
            # 四方に散らす
            self.neigh = [(1,0,1), (-1,0,1), (0,1,1), (0,-1,1)]
        else:  
            # 下マスに2つ散らす
            self.neigh = [(1,0,2), (0,1,1), (0,-1,1)]
        self.zc = sum(w for _,_,w in self.neigh)

    def _relax(self):
        """連鎖が終わるまで緩和させる。総転倒回数 s と、転倒した固有サイト数 a を返す。"""
        L, h = self.L, self.h
        q = deque(map(tuple, np.argwhere(h >= self.zc)))
        s = 0
        visited = set()
        while q:
            i, j = q.popleft()
            if h[i, j] < self.zc:
                continue
            h[i, j] -= self.zc
            s += 1
            visited.add((i, j))
            for di, dj, w in self.neigh:
                ii, jj = i + di, j + dj
                if 0 <= ii < L and 0 <= jj < L:
                    before = h[ii, jj]
                    h[ii, jj] += w
                    # 崩壊したらキュー追加
                    if before < self.zc and h[ii, jj] >= self.zc:
                        q.append((ii, jj))
            # 境界外にいったら消える
        return s, len(visited)

    def drop_and_relax(self):
        i = self.rng.integers(0, self.L)
        j = self.rng.integers(0, self.L)
        self.h[i, j] += 1
        return self._relax()

def simulate(L=64, steps=100000, mode="directed", seed=0):
    sp = Sandpile(L=L, mode=mode, seed=seed)
    sizes, areas = [], []
    for _ in range(steps):
        s, a = sp.drop_and_relax()
        if s > 0:
            sizes.append(s)
            areas.append(a)
    return np.array(sizes), np.array(areas), sp.h, sp

sizes, areas, H, sp = simulate(L=64, steps=150, mode="btw", seed=0)

時間発展させてみる
100回石を落とした 1000回石を落とした
2300回石、崩れそう 2500回、２連鎖が一度起きたところ

以下は10000回までの間、10ステップごとに保存したgifになっています。
8000回程度のところから連鎖の大きな雪崩が見えると思います。

雪崩の大きさの分布

横軸は雪崩の連鎖回数、縦軸は発生確率。この2つを両対数で描くと、直線に近い領域が現れ、べき分布っぽいものが見えるかと思います。

べき分布は以下のようなものです。

P(s)∝s^{−\tau}

シンプルに指数 \tau だけで形が決まるのが特徴です。典型的なスケールを持たず、場合によっては平均や分散が発散することすらあります（例えば \tau \leq 2 なら平均が、\tau \leq 3 なら分散も発散します）。これが「スケールを持たない分布」と呼ばれる理由です。
では、なぜ Sandpile でべき分布が現れるのでしょうか。

このモデルは石を一つずつ加え、閾値を超えたら崩壊が連鎖し、境界から外に出た石は散逸して消えます。その結果、系は平均的な石の高さが閾値より少し下の定常値に落ち着きます。外から特別な調整をしなくても自然に臨界に張り付く、この現象を自己組織化臨界 (Self-Organized Criticality, SOC) と呼びます。

臨界状態にあると、相関が系全体に広がります。つまり、どこから始まった崩壊でも全体に及ぶ可能性があり、小さな崩壊も巨大な崩壊も、確率的には同じメカニズムで出現します。有限サイズの系では、最大雪崩が制限されるため分布の右端が折れます。

結局のところ、べき分布は「臨界に張り付いた系のサイン」であり、Sandpile のようなトイモデルで自然に現れるのはとても面白い現象だと思います。

おわりに

今回見た Sandpile Model は、シンプルなルールにもかかわらず、自己組織化臨界性という現象を示すことを見てきました。やっぱり現象の本質を考える小さな実験台としてトイモデルはとても面白いと思います。

第一回と書いてある通り、次回はもっと詳しくこのモデルを掘り下げたり、べき指数の推定などにも触れてみようと思います。

Discussion