# coding: utf-8
# # Brownian Motion: A Numerical Experiment
# Suppose we want to describe the motion of a particle in a gas:
#
# ![Brownian motion](brownian-motion.gif)
# * We could model all those little particles and their collisions.
#
# * Model the motion of the big particle as random.
#
# * *Assume* the particle gets a 'bump' in X and Y that is uniformly distributed between $-0.1$ and $0.1$.
#
# Draw 10,000 random numbers modeling the bumps, print them, store them in the array `bumps`.
# In[2]:
from random import uniform
bumps = []
for i in range(10000):
b = uniform(-0.1, 0.1)
bumps.append(b)
# Plot a histogram to confirm the distribution of the bumps:
# In[3]:
import matplotlib.pyplot as pt
pt.hist(bumps)
# **Question 1:** Where does the particle end up, on average?
#
# **Question 2:** How does the average distance from the origin depend on the time?
#
# **Approach:** Sum the bumps.
# In[4]:
sum(bumps)
# How did it get there?
#
# Plot the path.
# In[5]:
pos = 0
path = []
for bump in bumps:
path.append(pos)
pos += bump
path.append(pos)
pt.plot(path)
# Individual realizations: not that helpful.
#
# Do 1000 of these. Store the endpoints in `endpoints`.
# In[9]:
from random import normalvariate
endpoints = []
for realization in range(1000):
bumps = []
for i in range(10000):
b = normalvariate(0, 0.1)
bumps.append(b)
endpoint = sum(bumps)
endpoints.append(endpoint)
# Plot a histogram of the endpoints:
# In[7]:
pt.hist(endpoints)
# ## Answering Question 1
#
# To answer our question 1: What is the *average* endpoint?
# In[35]:
sum(endpoints)/len(endpoints)
# But we got even more information:
#
# * It looks like the endpoints follow a [Normal distribution](https://en.wikipedia.org/wiki/Normal_distribution).
# * Is this result robust with respect to distribution choice? (We somewhat arbitrarily chose 'uniform'.)
# ## Answering Question 2
#
# Let's compute the average distance from the origin based on the length of the simulation.
# In[24]:
from random import normalvariate
avg_dists = []
for nsteps in range(0, 1500, 100):
print(nsteps, "steps")
endpoints = []
for realization in range(1000):
bumps = []
for i in range(nsteps):
b = normalvariate(0, 0.1)
bumps.append(b)
endpoint = sum(bumps)
endpoints.append(endpoint)
avg_dists.append(sum(ep**2 for ep in endpoints) / len(endpoints))
# In[25]:
pt.plot(avg_dists)
# Got a hypothesis?