• List of Backups
• View the diff.
• View the diff current.
• View the source.
• Go to ry/Ilas2_05.
  • 1 (2016-06-22 (Wed) 05:04:07)
  • 2 (2016-06-22 (Wed) 05:27:34)
  • 3 (2016-06-23 (Thu) 06:21:24)
  • 4 (2016-06-25 (Sat) 06:00:04)
  • 5 (2016-06-25 (Sat) 09:00:06)
  • 6 (2016-06-25 (Sat) 11:17:07)
  • 7 (2016-06-25 (Sat) 15:16:50)

---

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

Langevin方程式の数値解法 †

ランダム力の積分 †

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の関数として表示せよ．グラフが見やすいように表示範囲を調整してみよう．