ランダム・ウォークのシミュレーション

jupyter notebook で matplotlib を使うときの呪文。(プログラムファイルを編集してターミナルから実行する場合には不要)

%matplotlib inline

ライブラリのimport

import random
import numpy as np
import matplotlib.pyplot as plt

1次元

n = int(input("Enter number of steps:"))

x = 0.0
xs = [x]
for j in range(n):
    x += (random.random() - 0.5) * 2
    xs.append(x)
plt.plot(xs, marker='')

#plt.savefig("rw01.png")
plt.show()

Enter number of steps:300

複数のトラジェクトリを走らせて統計(ヒストグラム)をとる。

m = int(input("Enter number of runs:"))
n = int(input("Enter number of steps:"))

trj = np.zeros((m, n))

for i in range(m):
    x = 0.0
    for j in range(n):
        x += (random.random() - 0.5) * 2
        trj[i,j] = x

for i in range(20):
    plt.plot(trj[i,:])
plt.show()

plt.hist(trj[:,-1], bins=50)
#plt.savefig("rw01hist.png")
plt.show()

Enter number of runs:1000
Enter number of steps:500

平均2乗変位

m = int(input("Enter number of runs:"))
n = int(input("Enter number of steps:"))

nplot = 20
pintv = m / nplot
xsq = np.zeros(n)

for i in range(m):
    x = 0.0
    xs = [x]
    for j in range(n):
        x += (random.random() - 0.5) * 2
        xsq[j] += x**2
        xs.append(x)
    if i % pintv == 0: plt.plot(xs) 

plt.show()

plt.plot(xsq/m)
#plt.savefig("rw01msd.png")
plt.show()

Enter number of runs:1000
Enter number of steps:500

2次元

n = int(input("Enter number of steps:"))

x = 0.0
y = 0.0
xlist = [x]
ylist = [y]

for i in range(n):
    x += (random.random() - 0.5) * 2
    y += (random.random() - 0.5) * 2
    xlist.append(x)
    ylist.append(y)

plt.plot(xlist, ylist, marker='+', color = 'r')
#plt.savefig("rw02.png")
plt.show()

Enter number of steps:200

3次元

from mpl_toolkits.mplot3d import Axes3D

n = int(input("Enter number of steps:"))

fig = plt.figure()
ims = []

x = 0.0
y = 0.0
z = 0.0
xlist = [x]
ylist = [y]
zlist = [z]

#ax = plt.subplot(1,1,1,projection='3d')
ax = fig.gca(projection='3d')

for i in range(n):
    x += (random.random() - 0.5) * 2
    y += (random.random() - 0.5) * 2
    z += (random.random() - 0.5) * 2
    xlist.append(x)
    ylist.append(y)
    zlist.append(z)

ax.plot(xlist, ylist, zlist, marker='+', color = 'r')
plt.show()

Enter number of steps:200

上のコードをエディタに貼り付け、ターミナルから起動してもよい。その場合は、グラフィック・ウィンドウ上でドラッグすると図を回転できる。

アニメーション

jupyter notebook 内ではアニメーションの表示が出来ないので、以下のコードをエディタに貼り付け、ターミナルから実行するとよい。

1次元

import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

n = int(input("Enter number of steps:"))

fig = plt.figure()
ims = []

x = 0.0
xlist = [x]

for i in range(n):
    x += (random.random() - 0.5) * 2
    xlist.append(x)
    im = plt.plot(xlist, marker='+', color = 'r')
    ims.append(im)

ani = animation.ArtistAnimation(fig,ims,interval=100)
#ani.save('rw01anim.gif', writer='imagemagick')
plt.show()

2次元

import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

n = int(input("Enter number of steps:"))

fig = plt.figure()
ims = []

x = 0.0
y = 0.0
xlist = [x]
ylist = [y]

for i in range(n):
    x += (random.random() - 0.5) * 2
    y += (random.random() - 0.5) * 2
    xlist.append(x)
    ylist.append(y)
    im = plt.plot(xlist, ylist, marker='+', color = 'r')
    ims.append(im)

ani = animation.ArtistAnimation(fig,ims,interval=100)
#ani.save('rw02anim.gif', writer='imagemagick')
plt.show()

3次元

import random
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

n = int(input("Enter number of steps:"))

fig = plt.figure()
ims = []

x = 0.0
y = 0.0
z = 0.0
xlist = [x]
ylist = [y]
zlist = [z]

ax = fig.gca(projection='3d')

for i in range(n):
    x += (random.random() - 0.5) * 2
    y += (random.random() - 0.5) * 2
    z += (random.random() - 0.5) * 2
    xlist.append(x)
    ylist.append(y)
    zlist.append(z)
    im = ax.plot(xlist, ylist, zlist, marker='+', color = 'r')
    ims.append(im)

ani = animation.ArtistAnimation(fig,ims,interval=100)
#ani.save('rw03anim.gif', writer='imagemagick')
plt.show()