- List of Backups
- View the diff.
- View the diff current.
- View the source.
- Go to ry/Ilas2_06.
Brown運動の解析 †
粒子位置の分布関数 †
# 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()
# 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','x','y','z'], fontsize=14) plt.xlim(-40,40) plt.show()
粒子速度の分布関数 (Baxwell-Boltzmann分布) †
# 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()
# 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','x','y','z'], fontsize=14) plt.xlim(-6, 6) plt.show()
粒子の平均二乗変位と拡散係数 †
# 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$'], fontsize=16) plt.show()
自己相関関数 †
# auto correlation functions for F(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/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/m*time),'r',lw=2) ax.legend(['$C_F(t)=<\mathbf{F}(t)\cdot \mathbf{F}(0)>$','$C_V(t)=<\mathbf{V}(t)\cdot \mathbf{V}(0)>$','$(3k_BT/m)\exp(-\zeta t)$'], fontsize=16) plt.show()
スペクトル密度 †
# 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()
# power spectrum for F 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/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) ax.legend(['$S_F(\omega)=<|\~\mathbf{F}(\omega)|^2>$','$6k_BT\zeta/m$'], fontsize=16) plt.show()
宿題 †
確率過程(に見えるもの)の実データを解析してみよう.
- 地震(震度)の時系列データ
- 気温変動の時系列データ
- 株価の時系列データ
などなど,種々の実データを取得して,分布関数,自己相関関数,スペクトル密度などを求めてみよう.