* 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 = 10
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])
# First set up the figure, the axis, and the plot element we want to animate
fig, ax = plt.subplots(figsize=(7.5,7.5))
ax = plt.axes(xlim=(-box/2,box/2),ylim=(-box/2,box/2))
particles, = ax.plot([], [], 'ro', ms=10)
title = ax.text(0.7, 1.05, '', transform = ax.transAxes, va='center')
def init():
particles.set_data([], [])
title.set_text('')
return particles,title
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
particles.set_data(R[0,0:nump], R[1,0:nump])
title.set_text("t = "+str(i*dt))
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
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=nums, interval=5, blit=False, repeat=False)
** 保存された種々のデータ(A(t), B(t), ...)を時刻(t)の関数としてプロット [#ccc9f4ef]
# particle positions vs time
fig, ax = plt.subplots(figsize=(7.5,7.5))
for n in range(nump):
ax.plot(Rs[0,n,0:nums],time,'r', alpha=0.05)
ax.plot(Rs[1,n,0:nums],time,'b', alpha=0.05)
ax.plot(Rs[2,n,0:nums],time,'g', alpha=0.05)
plt.show()
** 宿題 [#j7077cb3]
