#

# 1D-randomwalk: A walker makes several steps,

# with a given number of walks pr. trial

#

import numpy, sys

from  matplotlib import pyplot as plt

import numpy as np

​

​

def mc_trial(number_walks,move_probability,walk_cum,walk2_cum):

    """

    Do a MonteCarlo trial, that is,

    random-walk one particle.

​

    Input:

    - number_walks:     Number of steps to walk the particle

    - move_probability: Probability that the particle

                        will step right when doing a step

    - walk_cum:         Numpy-array of length number_walks + 1,

                        containing the sum of the position

                        of the particles as a function of time

                        (usefull to calculate mean pos. as a function

                        of time)

    - walk2_cum:        Same as walk_cum, but with the sum of the

                        positions squared

​

    Output: As walk_cum and walk2_cum are numpy arrays, they are altered.

    """

    #Initial pos. As walk_cum[0]=walk2_cum[0] = 0.0

    #by initialization, it is uneccessary to add this step to

    #the arrays...

    pos = 0;

    for walk in range(number_walks+1):

        if numpy.random.random() <= move_probability:

            pos += 1

        else:

            pos -= 1

        walk_cum[walk]  += pos

        walk2_cum[walk] += pos**2

​

def mc_sample(trials, number_walks, move_probability):

    """

    Run as many trials as asked for by input variable trials.

    Wrapper to mc_trial, split out for easier paralellization

​

    Output: NumPy arrays walk_cum and walk2_cum, length number_walks + 1

    """

​

    walk_cum  = numpy.zeros(number_walks+1)

    walk2_cum = numpy.zeros(number_walks+1)

    for trial in range(trials):

        mc_trial(number_walks,move_probability,walk_cum,walk2_cum)

​

    return (walk_cum,walk2_cum)

​

#Main program

​

# initialize, values can easily be changed

trials           =   10000

number_walks     =   100

move_probability = 0.5

​

#Do the MC

(walk_cum,walk2_cum) = mc_sample(trials,number_walks, move_probability);

Dim = len(walk_cum)

x = np.zeros(Dim)

xaverage = np.zeros(Dim)

variance = np.zeros(Dim)

#Output

for i in range(Dim):

    x[i] = i

    #Normalize to number of trials (= number of walkers)

    xaverage[i]  = walk_cum[i]/float(trials)

    x2average = walk2_cum[i]/float(trials)

    variance[i] = x2average - xaverage[i]**2

plt.figure(1)

plt.subplot(211)

plt.xlabel(r'Steps')

plt.ylabel(r'Average displacement  $\Delta x$')

plt.plot(x, xaverage, 'b-')

plt.subplot(212)

plt.plot(x, variance, 'r-')

plt.ylabel(r'Variance  $\langle\Delta x)^2\rangle-\langle x\rangle^2$')

plt.savefig('rw.pdf')

plt.show()