Goal

We would like to model the motion of a particle in a box of specified dimensions. To do so, we explore random number generation within python to help simulate the motion and we practice its application in other functions.

Overview

Within the numpy module there are functions to generate random numbers. We use these functions for the following simulations:

Die Roll
Random Walk in 1 Dimension
Random Walk in 2 Dimensions
Brownian Motion
Animation

Die Roll

A very simple and intuitive application of random number generation is a six-sided die. First we define the function:

def rolldie():
  '''
  rolls a die
  '''
  roll = np.random.randint(1,7)
  return roll

Then can can simulate specified number of rolls to observe the effects.

rolls = [10,100,1000]
total=[]
plt.figure(figsize=(10,4))
for i in range(0,3):
  for j in range(0,rolls[i]):
    total.append(rolldie())
  plt.subplot(1,3,i+1)
  plt.hist(total,bins=np.arange(1,8))
  plt.xlabel("die roll")
  plt.xticks([1,2,3,4,5,6])
  plt.ylabel("count")
  plt.title(f"# rolls: {rolls[i]}")
plt.subplots_adjust(wspace=0.4)

Resultant plot:

As we can see in the graph, as the number of rolls total increase, the distribution of the number of rolls for each value becomes even. This is expected since there is an equal chance of obtaining any single value on the die.

Random Walk in 1D

To assist in our goal of particle motion simulation, we simulate random movement in one dimension. Essentially, for our randomwalk function, the walker will have equal probability to move either forward of backward. We define our function as follows:

def randomwalker(Nsteps=100,stepsize=1):
  '''
  creates a random position array of 'Nsteps' length with step size of 'stepsize'
  '''
  x=np.zeros(Nsteps)
  for i in range(1,Nsteps):
    direct=np.random.random()
    if direct<0.5:
      x[i]=x[i-1]-stepsize
    else:
      x[i]=x[i-1]+stepsize
  return x

Now, we simulate the walk for two walkers:

p1 = randomwalker()
p2 = randomwalker()

plt.figure()
plt.plot(np.arange(0,100),p1,"bo",label="person 1")
plt.plot(np.arange(0,100),p2,"rs",label="person 2")
plt.xlabel("step")
plt.ylabel("position")
plt.legend()

Resultant plot:

What happens to the mean square distance over time? This is explored below:

Nwalkers = 500
Nsteps = 100
# create an array xsq_sum to store 100 values of the x^2
xsq_sum=np.zeros(Nsteps)
# create a loop to call your function 500 times (or Nwalkers times)
for i in range(0,Nwalkers):
# after each call, square the x values that are returned and add them to xsum
  x = randomwalker(Nsteps)
  xsq_sum+=x**2
# once the loop is over, divide xsum by the number of walkers to get xsq_ave
xsq_ave=xsq_sum/Nwalkers

# create a plot of xsq_ave versus the number of steps.
plt.figure()
plt.plot(np.arange(0,Nsteps),xsq_ave,label="mean square")
plt.xlabel("step")
plt.ylabel("mean square of position")
# fit a polynomial to your data
c=np.polyfit(np.arange(0,Nsteps),xsq_ave,1)
line=np.polyval(c,np.arange(0,Nsteps))
plt.plot(np.arange(0,Nsteps),line,label="best fit",alpha=0.5)
# plot the best-fit polynomial with the data
plt.legend()

slope_bf = (max(line)-min(line))/len(line)
print(slope_bf)
minimum = min(line)
print(minimum)

Resultant plot:

The line represents the square distance the walker is away from the starting point (center) per step.

Random Walk in 2D

As a step up from a one dimension system, this new two-dimensional function has a random walk in a random radial direction. We define our new function as follows:

def randomwalker2d(Nsteps=100,stepsize=1):
  '''
  creates a random position array of 'Nsteps' length with step size of 'stepsize'
  '''
  x=np.zeros(Nsteps)
  y=np.zeros(Nsteps)
  for i in range(1,Nsteps):
    angle=np.random.uniform(0,2*np.pi)
    x[i] = x[i-1] + stepsize*np.cos(angle)
    y[i] = y[i-1] + stepsize*np.sin(angle)
  return x,y

And then we simulate:

x,y=randomwalker2d()
plt.plot(x,y, 'bo-')
plt.xlabel("x")
plt.ylabel("y")

Resultant plot:

What happens to the mean square distance over time? This is explored below:

# determine how <r^2> scales with the number of steps, as we did in Part 1
Nwalkers = 500
Nsteps = 100
# create an array xsq_sum to store 100 values of the x^2
rsq_sum=np.zeros(Nsteps)
# create a loop to call your function 500 times (or Nwalkers times)
for i in range(0,Nwalkers):
# after each call, square the x values that are returned and add them to xsum
  x,y = randomwalker2d(Nsteps)
  rsq_sum+=x**2+y**2
# once the loop is over, divide xsum by the number of walkers to get xsq_ave
rsq_ave=rsq_sum/Nwalkers

# create a plot of xsq_ave versus the number of steps.
plt.figure()
plt.plot(np.arange(0,Nsteps),rsq_ave,label="mean square")
plt.xlabel("step")
plt.ylabel("mean square of position")
# fit a polynomial to your data
c=np.polyfit(np.arange(0,Nsteps),rsq_ave,1)
line=np.polyval(c,np.arange(0,Nsteps))
plt.plot(np.arange(0,Nsteps),line,label="best fit",alpha=0.5)
# plot the best-fit polynomial with the data
plt.legend()

slope_bf = (max(line)-min(line))/len(line)
print(slope_bf)
minimum = min(line)
print(minimum)

Resultant plot:

Brownian Motion

This is similar to the random walk in 2D, but restricted within a boundary. So we adjust our code as appropriate:

def brownian(Nsteps=100,stepsize=1,length=10):
  x = np.zeros(Nsteps)
  y = np.zeros(Nsteps)
  for i in range(1,Nsteps):
    angle=np.random.uniform(0,2*np.pi)
    x[i] = (x[i-1] + stepsize*np.cos(angle))
    y[i] = (y[i-1] + stepsize*np.sin(angle))
    while (x[i]>length/2) or (x[i]<-(length/2)) or (y[i]>length/2) or (y[i]<-(length/2)):
      angle=np.random.uniform(0,2*np.pi)
      x[i] = (x[i-1] + stepsize*np.cos(angle))
      y[i] = (y[i-1] + stepsize*np.sin(angle))

  return x,y

Running the simulation:

x,y = brownian()
plt.plot(x,y,"bo-")
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.xticks(np.arange(-5,6))
plt.yticks(np.arange(-5,6))
plt.xlabel("x")
plt.ylabel("y")
plt.grid()

Resultant plot:

It is hard to see the confinement of the particle with this number of steps, so we increase it:

x,y = brownian(1000)
plt.plot(x,y,"bo-")
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.xticks(np.arange(-5,6))
plt.yticks(np.arange(-5,6))
plt.xlabel("x")
plt.ylabel("y")
plt.grid()

Resultant plot:

Animation

We can also animate the motion of our particle. The code to do this is as follows:

from matplotlib.animation import FuncAnimation
from matplotlib import rc

# ANIMATION CODE
# number of steps of random walker
nsteps = 10
print(f"Our animation will have {nsteps} timesteps")
L = 10 # length of the box that the particle is in


# call your random walker function to get the
# x and y positions
# pass in the length of the box L and the number of steps nsteps
xrand, yrand = brownian(nsteps,stepsize=1,length=L)

fig = plt.figure()
# set the range of our axes to the size of the box
ax = plt.axes(xlim=(-(L/2), L/2), ylim=(-(L/2),L/2))


# We also need to create an empty point that we'll draw into each frame
# This looks weird but it's just a regular plot...with no data in it!
point, = ax.plot([], [], 'ro', markersize=20)

# We need an initialization function that empties the point at the start
# of the animation
#
# These initialization routines are fairly common in other programs and
# routines.
#
# Note that all this is doing is emptying out the data!
def init():
    point.set_data([], [])
    return point,

# Now we need an animate function that takes in the step number
# and sets the values for the points
def animate(i):
    xpos = [xrand[i]]
    ypos = [yrand[i]]
    point.set_data(xpos, ypos)
    return point,

# Note that the number of frames is the number of steps here, though it doesn't have to be
# interval is the time in milliseconds between frames
# blit is whether or not to *only* draw the new stuff, see what happens when it's False
anim_walker = FuncAnimation(fig, animate, init_func=init,
                               frames=nsteps, interval=100, blit=True)

# We need this so it doesn't draw an extra figure
plt.close()

rc('animation', html='jshtml')
anim_walker

This will play an animation, like below, but moving.

Lab 10 ‐ Random Processes - MohamedS-Siena/computational_physics GitHub Wiki

Goal

Overview

Die Roll

Random Walk in 1D

Random Walk in 2D

Brownian Motion

Animation

⚠️ GitHub.com Fallback ⚠️

Lab 10 ‐ Random Processes - MohamedS-Siena/computational_physics GitHub Wiki

Goal

Overview

Die Roll

Random Walk in 1D

Random Walk in 2D

Brownian Motion

Animation

⚠️ **GitHub.com Fallback** ⚠️

⚠️ GitHub.com Fallback ⚠️