砂山崩壊モデル(Sandpile CA)

code: https://colab.research.google.com/drive/1lY2-h-m74U7ODkFLVCwAsf5xkMEChjni

cf. 『Mathematica複雑系のシミュレーション』R.J.ゲイロード・P.R.ウェリン著

自己組織臨界性(self-organized critical) SOAにあるものは、わずかな系の変化（たとえば一粒の砂）が加わるだけで、激変に導びく連鎖反応を引き起こす。

例。地学（地震、火山爆発）、天文学（パルサーの自然周期の突然増加）、経済学（株式相場のうらぎ）、生態系（進化）、流体力学（乱流）

条件

  • 内部サイトは１から８までの乱数整数値をとる。
  • ５以上の値は不安定とする。
  • 境界サイトの値は０。
  • ただし、最初の初期値は３か４。
  • そして、どこか１カ所の値が５以上になるまで１ずつ増やし、その後、ルールにしたがって計算を開始する。

ルールspread:

  • 自分自身安定してる場合(r<5)の時は、自分自身の値(r)と最近接サイトの不安定サイトの個数(t)を合計したものと置き換える。
  • 自身、不安定サイトの場合、r+tの合計から４を引く。これは崩壊の可能性を示唆するため。
  • それ以外で、自分自身が安定（０）の時は、そのまま。
In [ ]:
Sandpile[snd_, m_] := 
     Module[{absorbBC, landscape, spread, rule, foreststeps},
         absorbBC = 
            (Prepend[Append[Map[Prepend[Append[#, 0], 0] &, #],
                               Table[0, {Length[#] + 2}]], 
                       Table[0, {Length[#] + 2}]]) &;
         landscape = absorbBC[RandomInteger[{3, 4}, {snd, snd}]];
   
         While[Max[landscape]  <  5,
              randx = RandomInteger[ {2, snd + 1}];
              randy = RandomInteger[ {2, snd + 1}];
              landscape[[randx, randy]]++
              ];
   
   spread[r_, neighbors_] := Module[{t},
     If[r == 0, 0,
      t = Total[Floor[neighbors/5.0]];
      If[r < 5, r + t, r - 4 + t]]
     ];
   rule = {
     {{_, n_, _}, {w_, c_, e_}, {_, s_, _}} :> 
      spread[c, {n, w, e, s}]};
   foreststeps =
    CellularAutomaton[rule, landscape, m]];
In [ ]:
In [ ]:
sandsteps = Sandpile[50, 100];
In [ ]:
ListDensityPlot/@Take[sandsteps, {10,20,2}]
Out[ ]:
Output
In [ ]:
ListDensityPlot/@Take[sandsteps, {40,60,4}]
Out[ ]:
Output
In [ ]:
ListDensityPlot/@Take[sandsteps, {90,100,2}]
Out[ ]:
Output
In [ ]: