📂動力学

格子モデルシミュレーションにおける拡散

シミュレーション

このポストでは、格子空間でのある成分の拡散現象を模倣しようとする。これは同時に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)$から同じサイズの行列を作って、拡散が起きたかどうかを判断する。上記のケースでは、上、左、右へ成分が広がるのが見られる。実際にヒートマップで描いてみると次のようになる。

このプロセスを繰り返すと、本文や上部のようなウェブ画像が得られる。