前回扱った２次元ランダムウォークのシミュレーションは粒子が1つでしたが、 今回は複数の粒子が同時にランダムウォークするシミュレーションを作成してみましょう。 複数の粒子が同時に動くランダムウォークのシミュレーションは、気体分子のシミュレーションや人の集団の動きのシミュレーションなど、 様々な集団運動のモデルの基礎となります。
ここでは、同じ格子上に複数の粒子が存在できる場合を考えてみましょう。 これは粒子同士の衝突や相互作用を考えない場合、あるいは粒子の大きさが存在しない場合と言い換えることもできます （いわゆる理想気体の状態）。 同じ格子上に複数の粒子が存在できるということは、 それぞれの粒子の移動を考える際に移動先の格子上にすでに粒子が存在しているかどうかを考慮する必要がない、ということです。
この場合には、前回の１粒子のランダムウォークモデルを基にして、比較的簡単に実装することができます。 以下では、複数の粒子を扱うために前回のプログラムから拡張するために必要となる主な要素を考えてみることにしましょう。
上記のうち、２番目の「それぞれの格子上に粒子が存在しているかどうかを記憶する２次元配列 pos」は、 「粒子が重なることができる」場合には導入する必要はないのですが、 この後の課題で取り組む「粒子が重なることができない」場合に必要となるため、あらかじめ導入することにします。
これを踏まえて、粒子が重なることができる場合のプログラムは以下のようになります。
# 2次元ランダムウォーク（複数粒子、相互作用なし）
import numpy as np
import numpy.random as rd
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML
# 繰り返しの回数
T = 200
# 粒子の密度
d = 0.1
# 領域の大きさ
XMAX = 20
YMAX = 20
# 粒子数 = 格子の総数 * 密度（ただし、粒子数は整数）
P_N = int((XMAX+1)*(YMAX+1)*d)
print("粒子数：",P_N)
# 乱数のシードを指定（結果の再現性を得るためには必要）
rd.seed(1)
# 0~4の整数をランダムにT個生成して、リスト r に格納
r = rd.randint(0,5,(T,P_N))
# 各粒子の位置 px, py の初期化（粒子数が要素数である1次元配列）
px = np.zeros(P_N)
py = np.zeros(P_N)
# それぞれの格子に粒子が存在するかどうかを記憶する二次元配列の初期化
# １：粒子が存在する 0：粒子が存在しない
pos = np.zeros((XMAX+1, YMAX+1))
# 時刻 t を格納した配列の生成（0-Tまでの整数を順に格納）
t = np.arange(T+1)
# 各粒子をランダムに初期配置（重ならないように配置）
i = 0
while i < P_N:
# 粒子のx座標とy座標を0-XMAX&0-YMAXの乱数として取得
rx = rd.randint(0,XMAX+1)
ry = rd.randint(0,YMAX+1)
# 取得した座標に粒子が存在しなければ、粒子を配置
if pos[rx][ry] == 0:
# 取得した座標を各粒子の位置に代入
px[i] = rx
py[i] = ry
# 粒子の存在する位置を記憶
pos[rx][ry] = 1
i += 1
# グラフの描画領域を作成
fig = plt.figure(figsize=[5,5])
# 描画領域内での描画場所を指定（1行1列で配置）
ax = fig.add_subplot(1,1,1)
# 各軸のラベル
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("y", fontsize=15)
# グラフの表示域の設定
ax.set_xlim([-1,XMAX+1])
ax.set_ylim([-1,YMAX+1])
# グラフ描画のために、空のリスト作成
p=[]
# 初期化のための関数（初期位置描画のために必要）
def PlotInit_2D():
# 現在の位置をplot
p.append(ax.scatter(px,py,c="red",s=100))
# ランダムウォークの実行するための関数
def PlotRW_2D(i):
global px, py
# 描画が重ならないようにするための処理
if len(p) > 0:
p.pop().remove()
# 現在の位置をplot
p.append(ax.scatter(px,py,c="red",s=100))
for j in range(P_N):
# 0なら、右に移動
if r[i][j] == 0:
px[j] += 1
# 1なら、左に移動
elif r[i][j] == 1:
px[j] -= 1
# 2なら、上に移動
elif r[i][j] == 2:
py[j] += 1
# 3なら、下に移動
elif r[i][j] == 3:
py[j] -= 1
# 境界条件の処理
# 左にはみ出た場合
if px[j] < 0:
px[j] = XMAX
# 右にはみ出た場合
elif px[j] > XMAX:
px[j] = 0
# 下にはみ出た場合
if py[j] < 0:
py[j] = YMAX
# 上にはみ出た場合
elif py[j] > YMAX:
py[j] = 0
# FuncAnimationとしてアニメーションを作成
anim = animation.FuncAnimation(fig,PlotRW_2D,init_func=PlotInit_2D,frames=range(T),interval=50,repeat=False)
# Notebook 上で表示するためのおまじない
plt.close()
HTML(anim.to_jshtml())
前回作成した１粒子のプログラムに比べて非常に長くなっているように感じるかもしれませんが、 ランダムウォークの処理の本体部分は「粒子の総数 P_N」の分だけ繰り返すという処理を加えているだけで、 非常にシンプルな拡張になっています。
なお、「粒子の総数 P_N」の分だけ繰り返すという処理については「粒子の番号」の順に処理を行っていますが、
「粒子が重なることができない」場合には「粒子の番号」の順に処理を行うことで、
意図せず粒子の移動についての優先順位がついてしまう可能性があります。
これを防ぐためには、「粒子の番号」の順に繰り返しの処理を行うのではなく、
移動の処理を行う「粒子の番号」を乱数で決定してランダムな順に動かして繰り返し処理を行う方法（ランダムサンプリング）があります。
例えば、以下のように書き換えることでランダムサンプリングによる処理を行うことができます。
# ランダムウォークの実行するための関数
def PlotRW_2D(i):
global px, py
# 描画が重ならないようにするための処理
if len(p) > 0:
p.pop().remove()
# 現在の位置をplot
p.append(ax.scatter(px,py,c="red",s=100))
for j in range(P_N):
k = rd.randint(0,P_N)
# 0なら、右に移動
if r[i][j] == 0:
px[k] += 1
# 1なら、左に移動
elif r[i][j] == 1:
px[k] -= 1
# 2なら、上に移動
elif r[i][j] == 2:
py[k] += 1
# 3なら、下に移動
elif r[i][j] == 3:
py[k] -= 1
# 境界条件の処理
# 左にはみ出た場合
if px[k] < 0:
px[k] = XMAX
# 右にはみ出た場合
elif px[k] > XMAX:
px[k] = 0
# 下にはみ出た場合
if py[k] < 0:
py[k] = YMAX
# 上にはみ出た場合
elif py[k] > YMAX:
py[k] = 0
ただし、ランダムサンプリングを用いた場合には 1 step の間に全ての粒子が1回ずつ動くことが保証されている訳ではありません。 2回以上動く粒子もあれば、1回も動かない粒子も出てくることがほとんどです。 これは、各粒子の持つ速度が異なっていると考えることもできます。 気体分子のシミュレーションの場合であればこれで問題ありませんが、人の動きのシミュレーション、 特に避難のシミュレーションなどを考える場合には動かない人が出てくるのは妥当なシミュレーションと言えないかもしれません。
そのような場合には、全ての人が1回だけ動き、かつランダムな順に移動の処理を行うような方法を考える必要が出てきます。 今回の授業ではこのような方法については扱いませんが、興味のある人はどのようなプログラムを組めば実装できるか考えてみるのも良いでしょう。
次に、同じ格子上に複数の粒子が存在できない（粒子間の相互作用がある）場合を考えてみましょう。 これは粒子に大きさが存在する場合には、自然なモデルだと言えます。 例えば、実在気体分子集団や人の動きのシミュレーションを行う場合にはこちらが適切なモデルとなるでしょう。 （参考：ファン・デル・ワールスの状態方程式）
粒子が重なることができない場合には、粒子が移動する際に移動先の格子にすでに粒子が存在しているかどうかの情報が必要になります。 もし移動先に粒子が存在していれば、そこに移動することはできないので動かそうとした粒子はその場に留まったままになります。 この処理を行うために、あらかじめ導入していた「それぞれの格子上に粒子が存在しているかどうかを記憶する２次元配列 pos」を用います。
以下では、粒子の相互作用がある場合のランダムウォークのアルゴリズムを考えてみましょう。
上記のアルゴリズムについては、移動の処理を判定する際に境界条件の処理を行っているので、 「相互作用のない」場合のアルゴリズムのように移動の処理を行った後に別途で境界条件の処理を行う必要はないことに注意してください。
複数の粒子が同時にランダムウォークするシミュレーションを行うプログラムを作成せよ。 ただし、同じ格子上で粒子は重なることはできないとし、境界条件は反射境界条件とする。 課題が早く終わった場合には、周期境界条件の場合のプログラムも作成せよ。 （最終課題で周期境界条件におけるシミュレーションを実装する際に使えます）
繰り返し回数や粒子の総数は、粒子が重なることができる場合のサンプルプログラムと同じで良い。 （ただし、自分で密度などのパラメータを変えて色々と実験してみるのも良いだろう）
# Notebook 上で表示するためのおまじない plt.close() HTML(anim.to_jshtml())
「課題１」で作成したプログラムを提出してください。
ただし、コメント文として
学生番号、名前、プログラムの簡単な説明
を書く事。
ファイル名は、11-rep.ipynb として、 提出したファイル単体で実行できるようにすること。
提出締め切り： 2024年12月16日19時
相互作用を考慮した２次元ランダムウォークのモデルを改良して、人の歩行シミュレーションのプログラムを作成せよ。
ここでは、以下のような人の歩行シミュレーションのモデルを考える。
人が左から右に動いていく動きを再現するために、例えば以下のような移動条件を設定する。
例えば、0 ~ 9 の整数が出る乱数を用意して、以下のルールで上下左右に人を移動させる。
このようなモデルを用いて、以下の課題に取り組んで シミュレーションを実装したプログラムファイル(.ipynb ファイル)および 様々なシミュレーションを試した結果についての考察レポート(pdf ファイル)を提出せよ。
プログラム中には適宜コメント文を挿入して、どのような処理を行っているかについて分かりやすく説明すること。
プログラムのファイル名は、repFinalProg.ipynb として、
考察レポートのファイル名はrepFinal.pdf とせよ。
考察レポートについては、必ず pdf ファイルで提出すること。
word ファイルなど他の形式のファイルで提出された場合には、大幅に減点をするので注意！
(1) XMAX = 20, YMAX = 20 で周期境界条件とした上で、中央の縦１列に配置した21個の粒子全てが少なくとも1回
右端から左端の境界に移動（ワープ）するまでにかかるステップ数を計算せよ。
（シミュレーションにかかる実際の時間ではなく、シミュレーション内でのステップ数である事に注意）
また、右に移動する確率を変化させたときに上記で計算したステップ数がどのように変化していくかについても調べて、
変化のグラフ等を作成した上でその理由を考察せよ。
余力のある人は、粒子数を変化させるなど様々な条件を変化させて、
そのときの上記で計算したステップ数の変化について調べた結果についても考察してみよ。
ヒント
それぞれの粒子が1回でもワープしたかどうかを記録するために、例えば w という配列(要素数は粒子の総数)を用意して、 初期値では配列の全ての値を 0 に設定した上で、j 番目の粒子がワープした場合には w[j] = 1 として、 全ての粒子が少なくとも1回ワープしたかどうかをチェックすれば良い。
(2) XMAX = 20, YMAX = 20 で反射境界条件とした上で、右側に移動しやすい粒子と左側に移動しやすい粒子の２種類を用意する。
右に移動する21個の粒子は左端の縦１列に配置し、左に移動する21個の粒子は右端の縦１列に配置する。
その上で、右に移動する全ての粒子が少なくとも1回右側の境界に到達し、かつ左に移動する全ての粒子が少なくとも1回左側の境界に到達する
までにかかるステップ数を計算せよ。
また、左右に移動する確率を変化させたときに上記で計算したステップ数がどのように変化していくかについても調べて、
変化のグラフ等を作成した上でその理由を考察せよ。
さらに、粒子数を変化させるなど様々な条件を変化させて、
そのときの上記で計算したステップ数の変化について調べた結果についても考察してみよ。
ヒント
２種類の粒子を用意する場合には、例えば位置に関する配列を px1, py1, px2, py2 として増やすほか、 格子上の配置を記憶する２次元配列 pos を使って、粒子１が存在する格子は pos[x][y] = 1、 粒子２が存在する格子は pos[x][y] = 2 として、粒子の種類を区別すれば良い。
これに加えて、これら人の歩行の問題に関して自分で自由に問題設定をし、シミュレーション実験・研究を行い、 結果や考察等を記載するとなお良く、採点においては評価・加点する。
注意：
提出するプログラムファイルは、JupyterLab の Notebook 上で実行できるものに限る。
実行できない場合は、大幅に評価は下がる。
必ず実行ができる事を確認すること。
当然のことながら、レポートのプログラムは自分で作成したものを提出すること。 友達や教員、TAの方々に分からないことを聞いたり、相談することはもちろん構わないが、 他人の制作物の丸写し（あるいは、それと同等と見なせるもの）は厳禁！ レポートにおける他人の丸写しやそれに類する行為は、 テストで言えばカンニング行為に当たるので、 見つけた場合にはカンニングと同等の厳しい処分を行う。 （実質上、単位取得は不可と思って良い） 疑わしい場合には、作成したプログラムやレポートについて追加で尋ねることもある。 それを防ぐためにも、自分の言葉で十分に説明を記載しておくこと。
提出上の注意提出期限：2025年1月22日（水）23時
締め切り後の提出については、いかなる理由があっても認めませんので、必ず余裕をもって提出するようにしてください。
より詳しいヒントは、こちらにあります。 （動画は以前に撮ったものですが、内容としては問題ないはず）
最終レポートについて、採点基準は以下の通りです。