読者です 読者をやめる 読者になる 読者になる

数学定数でブラウン運動

計算機 手慰み Python

擬似乱数を可視化するためにブラウン運動を計算している.

今回は,数学定数の各桁を外力としてみる.

数学定数のn桁目をdとしたときに,時間n-1からnまでに粒子に与える外力を

(F_x, F_y) = (cos(2 * pi * d / 10), sin(2 * pi * d / 10)

とする.

試した数学定数は

である.参考までにメルセンヌツイスターから得られた擬似乱数も試す.

とりあえず,100万桁計算した

結果の図

f:id:fjkz:20141207123818p:plain

凡例が小さくて見づらいが,薄い水色が円周率,緑が自然対数の底,青が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))

データは以下から得た.

A million digits of Pi

http://apod.nasa.gov/htmltest/gifcity/e.1mil

http://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil

http://apod.nasa.gov/htmltest/gifcity/sqrt3.1mil

このサイトにはもっと多くの桁と種類があるので,次に試してみる.

Mathematical Constants - Millions of Digits