格子モデルシミュレーションにおける拡散
シミュレーション
このポストでは、格子空間でのある成分の拡散現象を模倣しようとする。これは同時にSI病気の拡散モデルのシミュレーションでもあり、空間が限られているという点でSIRモデルとも見ることができる。
変数
- $t$: 現在のターンを意味する。
- $I(t) \in \mathbb{N}$: $t$ターンで拡散している成分の量を表す。
- $S(t) \in \mathbb{N}$: $t$ターンで成分がない空き空間の量を表す。
パラメータ
- $\beta \in [0,1]$: 拡散率で、成分がどのように拡散していくかを示す。
アクション
(すべてのセルで、毎ターン)
- 計算: 各セルを基準に上下左右に成分が何個あるか数える。$i$行$j$列で計算された数を$n_{ij}$としよう。
- 拡散: 各セルは$1-(1-\beta)^{n_{ij}}$の確率で成分を持つことになる。
コードレビュー
ステップ1. 初期化と格子空間の生成
julia> row\_size = 10
10
julia> column\_size = 10
10
julia> β = 0.5 # 확산률
0.5
julia> Random.seed!(0);
julia> stage\_lattice = rand(['S'], row\_size, column\_size)
10×10 Array{Char,2}:
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
julia> stage\_lattice[rand(1:row\_size), rand(1:column\_size)] = 'I'; stage\_lattice
10×10 Array{Char,2}:
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'I' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S' 'S'
格子ベースのシミュレーションチュートリアルと似ていて、格子空間を生成してランダムな位置に成分を入れる。上述のケースでは、$6$行$6$列に生成されたことが確認できる。これをヒートマップで視覚化すると次のようになる。
ステップ2. 計算
julia> I\_lattice = (stage\_lattice .== 'I')
10×10 BitArray{2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
julia> count\_lattice =
vcat(I\_lattice[2:end, :], zeros(Int64, 1, column\_size)) +
vcat(zeros(Int64, 1, column\_size), I\_lattice[1:end-1, :]) +
hcat(I\_lattice[:, 2:end], zeros(Int64, row\_size, 1)) +
hcat(zeros(Int64, column\_size, 1), I\_lattice[:, 1:end-1])
10×10 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 1 0 1 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
周りに成分がどれだけあるかを数える行列$(n_{ij})$を作る。例ではベクター演算を使って計算したが、必ずしもこのように計算する必要はない。
ステップ3. 拡散
julia> probability\_lattice = 1 .- (1-β).^count\_lattice
10×10 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> hit\_lattice = probability\_lattice .> rand(row\_size, column\_size)
10×10 BitArray{2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 1 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
各セル毎に拡散を受ける確率を計算して行列$\left( 1-(1-\beta)^{n_{ij}} \right)$を作る。この計算理由は、周りに成分が多いほど拡散率が高くなることを反映するためだ。空き空間が成分を受ける確率は$\beta$、受けない確率は$(1-\beta)$だ。空き空間$i$行$j$列が周りの$n_{ij}$個の成分から拡散を受けない確率は$(1-\beta)^{n_{ij}}$で、結局、その中のどれかから拡散を受けた確率がちょうど$\left( 1-(1-\beta)^{n_{ij}} \right)$である。
julia> probability\_lattice = 1 .- (1-β).^count\_lattice
10×10 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.5 0.0 0.5 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
julia> hit\_lattice = probability\_lattice .> rand(row\_size, column\_size)
10×10 BitArray{2}:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0
0 0 0 0 1 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
今、一様分布$U(0,1)$から同じサイズの行列を作って、拡散が起きたかどうかを判断する。上記のケースでは、上、左、右へ成分が広がるのが見られる。実際にヒートマップで描いてみると次のようになる。
このプロセスを繰り返すと、本文や上部のようなウェブ画像が得られる。