* Brown運動のシミュレーション [#m5769370]

溶媒内にある質量&mathjax{m};の粒子(灰色)は,周囲の溶媒分子(赤色)の衝突に起因する熱揺動力&mathjax{\mathbf{F}(t)};と,粒子の速度に比例した流体抵抗力&mathjax{-\zeta\mathbf{V}(t)};とを受けながら不規則に運動する.これがブラウン運動である.

#ref(Brownian.JPG,center,40%)

- &mathjax{\mathbf{R}(t)};: 粒子の位置座標
- &mathjax{\mathbf{V}(t)};: 粒子の速度
- &mathjax{\mathbf{F}(t)};: 溶媒分子による熱揺動力
- &mathjax{t};: 時間
- &mathjax{m};: 粒子の質量
- &mathjax{\zeta};: 粒子の流体抵抗係数

** Langevin方程式 [#o72ff680]

膨大な数に上る全ての溶媒分子の運動を正確に知ることは非現実的なので,搖動力を確率的な力として定式化する.つまり粒子の運動方程式を,確率微分方程式である以下のLangevin方程式で表すことにする.

#mathjax(m\frac{d\mathbf{V}(t)}{dt}=-\zeta\mathbf{V}(t)+\mathbf{F}(t) \tag{1})

ここで確率的な揺動力&mathjax{\mathbf{F}(t)};は形が決まった関数として書けないので,ランダムなホワイトノイズ(時間相関がない→コンピュータの擬似乱数がそのまま使える)として取り扱う.第2回目の講義で説明したように,以下のようにxyz各成分ともに平均がゼロ,分散が&mathjax{2k_B T\zeta /m};であることが揺動散逸定理により要請される(異なる成分の相関はない).

#mathjax(\langle \mathbf{F}(t)\rangle=\mathbf{0} \tag{2})

#mathjax(\langle \mathbf{F}(t)\mathbf{F}(0)\rangle = \frac{2k_B T\zeta}{m}\mathbf{I}\delta(t)  \tag{3})

** ランダム力の積分 [#xe4b06d8]

時刻&mathjax{t};を小さな間隔&mathjax{\Delta t};で区切り,&mathjax{t_i\equiv i\Delta t \ (i=0,1,2,...)};と離散化して考える.

#mathjax(\mathbf{V}_{i+1}=\mathbf{V}_i-\frac{\zeta}{m}\int_{t_i}^{t_{i+1}} dt\mathbf{V}(t)+\frac{1}{m}\int_{t_i}^{t_{i+1}} dt\mathbf{F}(t)\simeq\left(1-\frac{\zeta}{m}\Delta t\right)\mathbf{V}_i + \frac{1}{m} \Delta \mathbf{W}_i \tag{6})

ここで右辺第1項は被積分関数が滑らかなのでオイラー法を適用できる.しかし激しく振動する第2項には適用できないので,揺動力を&mathjax{\Delta t};の区間で積分した新たな確率変数&mathjax{\Delta \mathbf{W}_i};を導入する.

計算の詳細は省略するが,確率変数&mathjax{\Delta \mathbf{W}_i};のxyz各成分はやはりホワイトノイズであり,ともに平均がゼロ,分散が&mathjax{2k_B T\zeta \Delta t};の正規分布に従う.


** Brown運動する自由粒子のシミュレーション&動画表示 [#oca03ea3]

- 初期時刻&mathjax{t=0};に原点にあるnump個の粒子のBrown運動を同時にシミュレーションし,赤球で表示する.
- n=0の粒子のみ青線で軌跡を表示する.

 % matplotlib nbagg
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.mlab as mlab
 import matplotlib.animation as animation

 dim  = 3
 nump = 1000
 nums = 1024
 box  = 30
 dt   = 0.05
 zeta = 1.0
 m    = 1.0
 kBT  = 1.0
 std  = np.sqrt(2*kBT*zeta*dt)
 np.random.seed(0)
 R = np.zeros([dim,nump])
 V = np.zeros([dim,nump])
 F = np.zeros([dim,nump])
 Rs= np.zeros([dim,nump,nums])
 Vs= np.zeros([dim,nump,nums])
 Fs= np.zeros([dim,nump,nums])
 time = np.zeros([nums])
 omega= np.zeros([nums])

 # 3D Brownian particle simulation with 3D animation
 from mpl_toolkits.mplot3d import Axes3D
 fig = plt.figure(figsize=(10,10))
 ax = fig.add_subplot(111, projection = '3d')
 ax.set_xlim(-box/2,box/2)
 ax.set_ylim(-box/2,box/2)
 ax.set_zlim(-box/2,box/2)
 ax.set_xlabel("x", fontsize=20)
 ax.set_ylabel("y", fontsize=20)
 ax.set_zlabel("z", fontsize=20)
 ax.view_init(elev = 12, azim = 120)
 
 particles, = ax.plot([], [], [], 'ro',ms=8, alpha=0.5)
 title = ax.text(-150., 0.,200.,'', transform = ax.transAxes, va='center')
 line, = ax.plot([], [], [], lw=1, alpha=0.8)
 xdata, ydata, zdata = [], [], []
 n = 0
 
 def init():
     title.set_text('')
     xdata.append(R[0,n])
     ydata.append(R[1,n])
     zdata.append(R[2,n])
     line.set_data(xdata,ydata)
     line.set_3d_properties(zdata)
     particles.set_data(R[0:2, :nump])
     particles.set_3d_properties(R[2, :nump])
     return particles,title,line
 
 def animate(i):
     global R,V,F,Rs,Vs,Fs,time
     F = std*np.random.randn(dim,nump)
     V = V*(1-zeta/m*dt)+F/m
     R = R + V*dt
     title.set_text("t = "+str(i*dt))
     xdata.append(R[0,n])
     ydata.append(R[1,n])
     zdata.append(R[2,n])
     line.set_data(xdata,ydata)
     line.set_3d_properties(zdata)
     particles.set_data(R[0:2, :nump])
     particles.set_3d_properties(R[2, :nump])
     Rs[0:dim,0:nump,i]=R
     Vs[0:dim,0:nump,i]=V
     Fs[0:dim,0:nump,i]=F
     time[i]=i*dt
     return particles,title,line
 
 anim = animation.FuncAnimation(fig, animate, init_func=init,
                                frames=nums, interval=2, blit=True, repeat=False)
 plt.show()

- n=0の粒子のみ軌跡が青線で表示されている.

** 保存された種々のデータ(A(t), B(t), ...)を時刻(t)の関数としてプロット [#ccc9f4ef]

- 1番目(n=0)の粒子の位置(R_x, R_y, R_z)を時刻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()

** 宿題 [#j7077cb3]

- 1番目(n=0)の粒子の速度(V_x, V_y, V_z)を時刻tの関数として表示せよ.グラフが見やすいように表示範囲を調整してみよう.