Brown運動のシミュレーション †

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

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

Langevin方程式 †

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

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

ここで確率的な揺動力&mathjax{\mathbf{F}(t)};は形が決まった関数として書けないので，ランダムなホワイトノイズ（時間相関がない→コンピュータの擬似乱数がそのまま使える）として取り扱う．第２回目の講義で説明したように，以下のように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})

↑

ランダム力の積分 †

時刻&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})

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

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

↑

Brown運動する自由粒子のシミュレーション＆動画表示 †

% 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()
↑

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

  • 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()
↑

宿題 †

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