* Brown運動の解析 [#hd990bc1]
** 粒子位置の分布関数(シミュレーション終了時) [#cc1c1ce8]
# positional distribution of particles at the end of the simulation
fig, ax = plt.subplots(figsize=(7.5,7.5))
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')
plt.show()
** 粒子の平均二乗変位と拡散係数 [#ade886f4]
# mean square displacement vs time
fig, ax = plt.subplots(figsize=(7.5,7.5))
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)
plt.show()
** ランダム力(F)と粒子速度(V)の自己相関関数 [#z279f806]
# auto correlation functions for F(green) and V(blue)
fig, ax = plt.subplots(figsize=(7.5,7.5))
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/dt # corrf(0)*dt=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*time),'r',lw=2)
plt.show()
** 粒子速度(V)のスペクトル密度 [#u6dbde3f]
# power spectrum for V
fig, ax = plt.subplots(figsize=(7.5,7.5))
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
plt.show()
** ランダム力(F)のスペクトル密度 [#oa14e9e1]
# power spectrum for F
fig, ax = plt.subplots(figsize=(7.5,7.5))
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/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/m,'r',lw=2) # -\infty < omega [2\pi*Hz] < \infty
plt.ylim(0,8)
plt.show()
** 宿題 [#ca7cff77]
確率過程(に見えるもの)の実データを解析してみよう.
- 地震(震度)の時系列データ
- 気温変動の時系列データ
- 株価の時系列データ
などなど,種々の実データを取得して,分布関数,自己相関関数,スペクトル密度などを求めてみよう.