ブラウン運動をシミュレーションしてみる

題名の通り．擬似乱数を可視化するために，ランジュバン方程式に基づくブラウン運動を計算してみた．*1

結果の絵．

コードは以下．

import math
import random

class Particle(object):
    '''
    A particle moving in 2D space.
    '''
    position_x = 0.0
    position_y = 0.0
    momentum_x = 0.0
    momentum_y = 0.0

    mass = 1.0
    viscous = 1.0

    def evolute(self, force_x, force_y, dt):
        # Solve Langevin equation with leap-flog method
        self.position_x = (self.position_x +
                           self.momentum_x / self.mass * dt)
        self.position_y = (self.position_y +
                           self.momentum_y / self.mass * dt)
        self.momentum_x = (self.momentum_x + 
                           (force_x -
                            self.viscous * self.momentum_x / self.mass) * dt)
        self.momentum_y = (self.momentum_y + 
                           (force_y -
                            self.viscous * self.momentum_y / self.mass) * dt)

random.seed(0)

# random integer in range [0, 9]
def nextint():
    return random.randint(0, 9)

# number of steps per a unit time
idt = 10
dt = 1.0 / idt

print 'position_x    position_y'

particle = Particle()

for t in xrange(1000000):
    i = nextint()
    
    # Get an external force vector (force_x, force_y) from int in range [0, 9].
    # The length is unity.
    force_x = math.cos(2.0*math.pi/10.0*i)
    force_y = math.sin(2.0*math.pi/10.0*i)

    for t_ in range(idt):
        particle.evolute(force_x, force_y, dt)
        
    print ('{0:e}    {1:e}'
            .format(particle.position_x, particle.position_y))

擬似乱数はとりあえずPythonの組み込みのもの（Mersenne Twister）である．比較がないので，なんともいえないが，なんかランダムっぽい動きをしている．

円周率とか自然対数の底とかの数学定数を外力に使ってみて，どんな感じになるのか遊んでみる予定．

*1:これをしているのを昔どこかで見たのだけれど思い出せない