• 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  = 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)の関数としてプロット †

# 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()

宿題 †