題名の通り.擬似乱数を可視化するために,ランジュバン方程式に基づくブラウン運動を計算してみた.*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:これをしているのを昔どこかで見たのだけれど思い出せない