期末レポート課題

期末レポート課題の内容については、こちらを参照してください。

内容

• 課題 (1) のヒント
• 課題 (2) のヒント

解説動画

課題 (1) のヒント

まず第12回の課題１で取り組んだ「相互作用を考慮した2次元ランダムウォークモデル（周期境界条件）」の プログラム（一部省略）から考えてみましょう。

# 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.99
# 領域の大きさ
XMAX = 20
YMAX = 20
# 粒子数　=　格子の総数　*　密度（ただし、粒子数は整数）
P_N = int((XMAX+1)*(YMAX+1)*d)
print("粒子数：",P_N)


# 乱数のシードを指定（結果の再現性を得るためには必要）
rd.seed(1)

# 0~M-1の整数をランダムにT*P_N個生成して、リスト r に格納
M = 10
r = rd.randint(0,M,(T,P_N))



# 各粒子の位置 px, py の初期化（粒子数が要素数である1次元配列）
px = np.zeros(P_N)
py = np.zeros(P_N)

# それぞれの位置に粒子が存在するかどうかを記憶する二次元配列の初期化
# １：粒子が存在する　0：粒子が存在しない
pos = np.zeros((XMAX+1, YMAX+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




# ここで配列 w と変数 wt の初期化




# グラフの描画領域を作成
fig = plt.figure(figsize=[6,6])

# 描画領域内での描画場所を指定（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, pos

  
  # 描画が重ならないようにするための処理
  if len(p) > 0:
    p.pop().remove()
  # 現在の位置をplot
  p.append(ax.scatter(px,py,c="red",s=100))
  

  # P_N 個の粒子のランダムウォークについての繰り返し処理
  for j in range(P_N):

    # 相互作用を考慮した２次元ランダムウォーク（周期境界条件）の処理
    # この部分は省略
    
    # 右端（x座標=XMAX）から左端（x座標=0）への移動（ワープ）処理の中に、ワープをしたチェックを追加
    


  
  
  # 全ての粒子のランダムウォークの処理（1 step 分）が終わったら、
  # 各粒子がワープしたかを判定して、全ての粒子がワープ済みなら現在時刻を記憶
  



# 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 は縦１列なので、YMAX + 1 個であるのは明らかです。
次に P_N 個の粒子の座標は「中央の縦１列」であることを考えれば、x 座標は int(XMAX/2)、 y 座標は i (0 $\le$ i $\le$ YMAX) とすれば良さそうです。

以上まとめると、次のようなアルゴリズムになります。

• for などを用いた繰り返し処理で、 以下のように i 番目 (0 $\le$ i $\le$ YMAX) の粒子の初期配置をおこなう
1. x 座標を px[i] に代入
2. y 座標を py[i] に代入
3. 粒子の位置を記憶するための二次元配列 pos 内の上記の x, y 座標に相当する要素に、1 を代入
   （この位置に粒子が存在するという意味）

全ての粒子がワープするまでにかかるステップ数の計算

課題のヒントにもあるように、各粒子が１回でもワープしたかどうかを記憶するために w という配列（要素数 P_N）を用意します。
初期値では配列の全ての値を 0 に設定しておいて、j 番目の粒子がワープした場合には w[j] = 1 を代入します。
その上で、各ステップですべての粒子の移動判定を終えた後に配列 w の全ての値が 1 かどうかをチェックすれば、 求めたいステップ数が計算できそうです。

しかし、この方法には少し問題点があります。
例えば、43 ステップ目に w の全ての値が 1 になったとすると、 それ以降では常にこの w の値に関する条件はクリアできてしまいます。
つまり、一番最初にこの条件を達成したステップ数を記憶するためには、何か別の変数を別途用意する必要があります。
ここでは、wt という変数にステップ数を記憶することにしましょう。

以上を踏まえて、全ての粒子が少なくとも１回ワープするまでにかかるステップ数を表示するアルゴリズムは以下のようになります。

• 配列 w (要素数 P_N)と変数 wt に 0 を代入して初期化
  青字で書かれた「# ここで配列 w と変数 wt の初期化」あたりに追加すればOK
• j 番目の粒子のランダムウォーク処理中で、粒子が右端から左端にワープしたかのチェック
  右端（x座標=XMAX）から左端（x座標=0）への移動（ワープ）処理の中に、ワープをしたチェックを追加
  具体的には、上記の移動処理に加えて w[j] = 1 （ワープ済みという意味）を加えればOK
• 各粒子がワープしたかの判定
  
  • wt = 0 （初期値）ならば、1 step前までは全粒子がワープ済みではない状態なので、以下の判定処理を実行
  1. 各粒子がワープしているかの判定を行うループ（ここでは後のために while を使った方が良い）
     もし１つでもワープしていない粒子が存在すれば、途中でループを抜け出す
     具体的には、while ループの中で k 番目の粒子について（0 $\le$ k < P_N）以下の処理を実行
     • もしワープしていないなら（if w[k] == 0）、ループを抜け出す（break）
     • もしワープ済みなら（else）、次の k に移動（k += 1）
  2. 1 のループを途中で抜け出さなければ、全ての粒子がワープ済みなのでステップ数を取得
     この場合は k = P_N となっているはずなので、これを判定条件（if k == P_N）にしてステップ数を記憶
     ただし、ここでワープのチェックしているのは移動処理を行った後の粒子についてなので、
     wt = i + 1 とすることに注意（i は現在の step 数）

（注）配列 w と変数 wt は関数内でも利用しているので、px, py, pos と同様に関数内で global 変数に設定する。

おまけ：より高速に多くのシミュレーションを行いたい場合のヒント

上で紹介したアルゴリズムで求めたい時間を計算することは可能ですが、全ての粒子がワープした後も ランダムウォークの処理は実行され続けます。 今知りたいのは全ての粒子がワープするのにかかるステップ数なので、それを求めたら処理を終えてしまえば 計算時間は短くなります。
また、描画処理も計算時間がかかるので、時間を求めるだけであればアニメーションを描画させない方が 圧倒的に計算時間は短くなります。

以上を踏まえて、アニメーションを描画させるプログラムとは別に、ステップ数を計算するためだけのプログラムを作成した方が 条件を変更した場合のシミュレーションを効率的に行うことができます。 描画に関する処理は、上のサンプルコードでは緑色の文字で書かれたブロックです。
また、関数 FuncAnimation では１ステップごとに関数 PlotRW_2D を呼び出していますが、 この処理自体にかなりの時間がかかります。
そこで、FuncAnimationを使わずに関数 PlotRW_2D で１ステップごとに行ってる処理を for 文などで繰り返し処理を行えば、計算時間はかなり短縮されます。

さらに、wt の値を取得したらループを抜け出す処理を書き加えれば、 その時点で全ての処理を終了させることができるので、より計算時間を短縮することができます。

課題 (2) のヒント

課題（2）は発展的な内容の課題として設定しているので、ヒントは少なめにします。

２種類の粒子の設定

課題のヒントにもあるように、２種類の粒子を扱うために一番簡単な拡張は粒子に関する様々な変数を ２つ用意すれば良いです。
位置に関する変数 px, py のほか、ワープしたかを記憶する変数 w や ランダムウォークのための乱数を保持する変数 r も２種類用意するのが分かりやすいでしょう。
またランダムウォークについても、粒子1と粒子2それぞれに関して分けて処理を書くのが分かりやすいでしょう。

２種類の粒子の描画

２種類の粒子を区別するためには、種類ごとに粒子の描画色を変えるなどの工夫が必要になります。
ここでは、以下のようにして粒子1を赤、粒子2を青にして描画してみます。

p.append(ax.scatter(px1,py1,c="red",s=100))
p.append(ax.scatter(px2,py2,c="blue",s=100))

粒子が１種類の場合には、p.append(ax.scatter(px,py,c="red",s=100)) などと書いていた場所を 上の２行に書き換えればOKです。

また、PlotRW_2D 関数内の

# 描画が重ならないようにするための処理
if len(p) > 0:
  p.pop().remove()

という所については、以下のように p.pop().remove() を２つ書く必要があります。

# 描画が重ならないようにするための処理
if len(p) > 0:
  p.pop().remove()
  p.pop().remove()

理由としては、p.append の方で２つの描画を追加しているので、remove も２つする必要があるためです。

戻る