* Brown運動のシミュレーション [#m5769370]
** Langevin方程式の数値解法 [#o72ff680]
** ランダム力の積分 [#xe4b06d8]
** Brown運動する自由粒子のシミュレーション＆動画表示 [#oca03ea3]
% 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)の関数としてプロット [#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の関数として表示せよ．グラフが見やすいように表示範囲を調整してみよう．
