• The added line is THIS COLOR.
• The deleted line is THIS COLOR.
• Go to ry/Ilas2_06.
• Deleting diff of ry/Ilas2_06

* Brown運動の解析 [#hd990bc1]

前回の「Brown運動のシミュレーション」のつづき．「Brown運動する自由粒子のシミュレーション＆動画表示」までを予め実行してください．そこで作成されたデータを今回の解析で使用します．

** 粒子位置の分布関数[#cc1c1ce8]

- 動画で軌跡が青線で示されているn=0の粒子について，位置&mathjax{R_x(t),R_y(t),R_z(t) };を横軸，時刻&mathjax{t};を縦軸にとってプロットする．

 # particle positions vs time
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("$R_{x,y,z}(t)$", fontsize=20)
 ax.set_ylabel("$t$", fontsize=20)
 n = 0
 ax.plot(Rs[0,n,0:nums],time,'r')
 ax.plot(Rs[1,n,0:nums],time,'b')
 ax.plot(Rs[2,n,0:nums],time,'g')
 plt.xlim(-40, 40)
 ax.legend(['x','y','z'], fontsize=14)
 plt.show()

- nump個の粒子の位置&mathjax{R_x(t),R_y(t),R_z(t)};を同時にプロットし，時刻&mathjax{t};の経過による粒子位置の分布の変化をみる．

 # particle positions vs time
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("$R_{x,y,z}(t)$", fontsize=20)
 ax.set_ylabel("$t$", fontsize=20)
 for n in range(nump):
     ax.plot(Rs[0,n,0:nums],time,'r', alpha=0.1)
     ax.plot(Rs[1,n,0:nums],time,'b', alpha=0.1)
     ax.plot(Rs[2,n,0:nums],time,'g', alpha=0.1)
 plt.xlim(-40, 40)
 plt.show()

- シミュレーション終了時&mathjax{t=t_{\rm end}};における粒子位置の分布関数&mathjax{P(R_x),P(R_y),P(R_z)};を，平均&mathjax{<P>=0};，分散&mathjax{\sigma^2=\frac{2k_BTt_{\rm end}}{\zeta}};の正規分布&mathjax{P(R)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{(P-<P>)^2}{2\sigma^2}\right]};と比較する．

 # positional distribution of particles at the end of the simulation
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("$R_{x,y,z}(t=t_\mathrm{end})$", fontsize=20)
 ax.set_ylabel("$P(R_{x,y,z})$", fontsize=20)
 plt.hist(Rs[0,0:nump,nums-1], bins=50,normed=True,facecolor='r',alpha=0.5)
 plt.hist(Rs[1,0:nump,nums-1], bins=50,normed=True,facecolor='b',alpha=0.5)
 plt.hist(Rs[2,0:nump,nums-1], bins=50,normed=True,facecolor='g',alpha=0.5)
 sig2=2*kBT/zeta*dt*nums
 ave=0.0
 x = np.arange(-40,40,1)
 y = np.exp(-(x-ave)**2/2/sig2)/np.sqrt(2*np.pi*sig2)
 plt.plot(x,y,lw=4,color='k')
 ax.legend(['Gauss $(\sigma^2=2k_BT t/\zeta)$','x','y','z'], fontsize=14)
 plt.xlim(-40,40)
 plt.show()

** 粒子速度の分布関数 (Maxwell-Boltzmann分布) [#v8e807ad]

- nump個の粒子の速度&mathjax{V_x(t),V_y(t),V_z(t)};を時刻&mathjax{t};の関数として同時にプロットし，粒子速度の分布の時間変化をみる．

 # particle velocities vs time
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("$V_{x,y,z}(t)$", fontsize=20)
 ax.set_ylabel("$t$", fontsize=20)
 for n in range(nump):
     ax.plot(Vs[0,n,0:nums],time,'r', alpha=0.1)
     ax.plot(Vs[1,n,0:nums],time,'b', alpha=0.1)
     ax.plot(Vs[2,n,0:nums],time,'g', alpha=0.1)
 plt.xlim(-6, 6)
 plt.show()

- シミュレーション終了時&mathjax{t=t_{\rm end}};における粒子速度の分布関数&mathjax{P(V_x),P(V_y),P(V_z)};を，平均&mathjax{<V>=0};，分散&mathjax{\sigma^2=\frac{2k_BT}{m}};の正規分布&mathjax{P(V)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{(V-<V>)^2}{2\sigma^2}\right]};と比較する（Maxwell-Boltzmann分布）．

 # velocity distribution of particles at the end of the simulation
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("$V_{x,y,z}(t=t_\mathrm{end})$", fontsize=20)
 ax.set_ylabel("$P(V_{x,y,z})$", fontsize=20)
 plt.hist(Vs[0,0:nump,nums-1], bins=30,normed=True,facecolor='r',alpha=0.5)
 plt.hist(Vs[1,0:nump,nums-1], bins=30,normed=True,facecolor='b',alpha=0.5)
 plt.hist(Vs[2,0:nump,nums-1], bins=30,normed=True,facecolor='g',alpha=0.5)
 sig2=kBT/m
 ave=0.0
 x = np.arange(-10,10,0.1)
 y = np.exp(-(x-ave)**2/2/sig2)/np.sqrt(2*np.pi*sig2)
 plt.plot(x,y,lw=4,color='k')
 ax.legend(['Baxwell-Boltzmann $(\sigma^2=k_BT/m)$','x','y','z'], fontsize=14)
 plt.xlim(-6, 6)
 plt.ylim(0, 0.7)
 plt.show()

** 粒子の平均二乗変位と拡散係数 [#ade886f4]

- 粒子の平均二乗変位&mathjax{\langle[\mathbf{R}(t)-\mathbf{R}(0)]^2\rangle};を時間&mathjax{t};の関数としてプロットし，Einsteinの式&mathjax{D=\frac{k_BT}{\zeta}};を用いた理論値&mathjax{6Dt};と比較する．

 # mean square displacement vs time
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("$t$", fontsize=20)
 ax.set_ylabel("mean square displacement", fontsize=16)
 msd = np.zeros([nums])
 for i in range(nums):
     for n in range(nump):
         for d in range(dim):
             msd[i]=msd[i]+Rs[d,n,i]**2
     msd[i] = msd[i]/nump
 ax.plot(time,msd)
 ax.plot(time,6*kBT/zeta*time)
 ax.legend(['$< R^2(t)>$','$6Dt={6k_BT t}/{\zeta}$'], fontsize=16)
 plt.show()

** 自己相関関数 [#z279f806]

- 粒子速度の自己相関関数&mathjax{C_V(t)=\langle\mathbf{V}(t)\cdot \mathbf{V}(0)\rangle};を求め，理論値&mathjax{C_V(t)=\frac{3k_BT}{m}\exp\left[-\frac{\zeta}{m}t\right]};と比較する．
- 揺動力による力積の自己相関関数&mathjax{C_W(t)=\langle\Delta\mathbf{W}(t)\cdot \Delta\mathbf{W}(0)\rangle};を求め，ホワイトノイズであることを確認する．

 # auto correlation functions for dW(green) and V(blue)
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("$t$", fontsize=20)
 ax.set_ylabel("auto correlation", fontsize=16)
 def auto_correlate(dt):
     cor = np.correlate(dt,dt,mode="full")
 #    return cor[cor.size/2:]      # This is original but creates a warning
     return cor[nums-1:]          # cor[0:2*nums-1] is an even function centered at nums-1
 X     = np.zeros([nums])
 corrv = np.zeros([nums])
 corrf = np.zeros([nums])
 for n in range(nump):
     for d in range(dim):
         X = Fs[d,n,0:nums]
         corrf=corrf+auto_correlate(X)/nums
         X = Vs[d,n,0:nums]
         corrv=corrv+auto_correlate(X)/nums
 corrf=corrf/nump                 # corrf(0)*dt/m=2*dim*kBT*zeta/m
 corrv=corrv/nump                 # corrv(0)=dim*kBT/m
 plt.xlim(-1, 10)
 #plt.ylim(-1, 8)
 ax.plot(time,corrf,'g',lw=4)
 ax.plot(time,corrv,'b',lw=4)
 ax.plot(time,dim*kBT/m*np.exp(-zeta/m*time),'r',lw=2)
 ax.legend(['$C_W(t)=<\Delta\mathbf{W}(t)\cdot \Delta\mathbf{W}(0)>$','$C_V(t)=<\mathbf{V}(t)\cdot \mathbf{V}(0)>$','$(3k_BT/m)\exp(-\zeta t/m)$'], fontsize=16)
 plt.show()


** スペクトル密度 [#u6dbde3f]

- 粒子速度のスペクトル密度&mathjax{S_V(\omega)=\langle |\mathbf{V}(\omega)|^2\rangle};を求め，理論値&mathjax{S_V(\omega)=\frac{6k_BT}{m^2\omega^2+\zeta^2}};と比較する．


 # power spectrum for V
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("angular frequency, $\omega$", fontsize=16)
 ax.set_ylabel("spectrum density", fontsize=16)
 X = np.zeros([nums])
 Y = np.zeros([nums])
 Z = np.zeros([nums])
 omega=np.zeros(nums)
 for n in range(nump):
     for d in range(dim):
         X = Vs[d,n,0:nums]
         Y,omega=mlab.psd(X,NFFT=1024, Fs=1/dt,noverlap=1024/4)
         Z[:len(Y)] = Z[:len(Y)] + Y
 Z = Z/nump
 ax.plot(omega*2*np.pi,Z[0:len(Y)]/2,'b',lw=4)           # 0 < omega [Hz] < \infty, Z[0]/2=2*dim*kBT/zeta**2
 plt.xlim(0, 10)
 plt.ylim(0, 8)
 ax.plot(omega,6*kBT/(m**2*omega**2+zeta**2),'r',lw=2)  # -\infty < omega [2\pi*Hz] < \infty
 ax.legend(['$S_V(\omega)=<|\~\mathbf{V}(\omega)|^2>$','$6k_BT/(m^2\omega^2+\zeta^2)$'], fontsize=16)
 plt.show()

- 揺動力による力積スペクトル密度&mathjax{S_W(\omega)=\langle |\mathbf{\Delta W}(\omega)|^2\rangle};を求め，理論値&mathjax{S_W(\omega)=6k_BT\zeta\Delta t;と比較する．
- 揺動力による力積スペクトル密度&mathjax{S_W(\omega)=\langle |\mathbf{\Delta W}(\omega)|^2\rangle};を求め，理論値&mathjax{S_W(\omega)=6k_BT\zeta\Delta t};と比較する．


 # power spectrum for dW
 fig, ax = plt.subplots(figsize=(7.5,7.5))
 ax.set_xlabel("angular frequency, $\omega$", fontsize=16)
 ax.set_ylabel("spectrum density", fontsize=16)
 X = np.zeros([nums])
 Y = np.zeros([nums])
 Z = np.zeros([nums])
 omega=np.zeros(nums)
 for n in range(nump):
     for d in range(dim):
         X = Fs[d,n,0:nums]
         Y,omega=mlab.psd(X,NFFT=1024, Fs=1/dt,noverlap=1024/4)
         Z[:len(Y)] = Z[:len(Y)] + Y
 Z = Z/nump/dt
 ax.plot(omega*2*np.pi,Z[0:len(Y)]/2,'b',lw=4)               # 0 < omega [Hz] < \infty, Z=2*dim*kBT*zeta/m
 ax.plot(omega*2*np.pi,0*omega+2*dim*kBT*zeta*dt,'r',lw=2)   # -\infty < omega [2\pi*Hz] < \infty
 plt.ylim(0,1)
 ax.legend(['$S_W(\omega)=<|\Delta\~\mathbf{W}(\omega)|^2>$','$6k_BT\zeta\Delta t$'], fontsize=16)
 plt.show()

** 宿題 [#ca7cff77]

確率過程（に見えるもの）の実データを解析してみよう．

- 地震（震度）の時系列データ

- 気温変動の時系列データ

- 株価の時系列データ

などなど，種々の実データを取得して，分布関数，自己相関関数，スペクトル密度などを求めてみよう．