今回は,数学定数の各桁を外力としてみる.
数学定数のn桁目をdとしたときに,時間n-1からnまでに粒子に与える外力を
(F_x, F_y) = (cos(2 * pi * d / 10), sin(2 * pi * d / 10)
とする.
試した数学定数は
である.参考までにメルセンヌツイスターから得られた擬似乱数も試す.
とりあえず,100万桁計算した
結果の図
凡例が小さくて見づらいが,薄い水色が円周率,緑が自然対数の底,青が2の平方根,紫が3の平方根,赤がメルセンヌツイスターである.
だから?って感じだが,円周率はメルセンヌツイスターと同程度に偏りが少ない感じに見える.
前回とほとんど同じだが,参考までにコードは以下:
import math import random import sys 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): # 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) if len(sys.argv) == 1: random.seed(0) # random integer in range [0, 9] def nextint(): return random.randint(0, 9) else: # read a file including a sequence of number in range [0,9] f = open(sys.argv[1], 'r') def nextint(): while(True): i_char = f.read(1) if i_char == '': raise EOFError try: i = int(i_char) return i except ValueError: pass # 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))
データは以下から得た.
http://apod.nasa.gov/htmltest/gifcity/e.1mil
http://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil
http://apod.nasa.gov/htmltest/gifcity/sqrt3.1mil
このサイトにはもっと多くの桁と種類があるので,次に試してみる.