#!/usr/bin/env python
# coding: utf-8
#
#
#
# # A brief introduction to Brownian motion and diffusion
#
# ### Examples - Statistical Mechanics
#
# By Andreas S. Krogen, Jonas Tjemsland, Håkon W. Ånes and Jon Andreas Støvneng.
#
# Last edited: March 22nd 2018
# ___
# ## Introduction
# In the year of 1827, a Scottish botanist named Robert Brown made a series of observations of pollen grains under his microscope [1]. The pollen grains were immersed in water and he discovered that their motion was of a rather erratic character, however he was not able to come up with a theory that could explain these observations. Today we know that this random motion is caused by a vast number of collisions between the grains and the water molecules, and the phenomenon is called *Brownian motion* in honour of Robert Brown. The first person to come up with a proper theoretical treatment of Brownian motion was Albert Einstein, who in 1905 published a paper providing a quantitative theory of the phenomenon (*Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen*).
# ## Brownian motion in one dimension
# We consider a very simple probability model of the random motion of a Brownian particle proposed by M. Smoluchowski [2]. Nevertheless, this model gives rise to some interesting physics, as we are about to discover soon. The time evolution of the position of the particle takes the form of discrete steps, and in each step there is an equal probability of the particle moving a distance $h$ to the left or a distance $h$ to the right. The total number of steps taken is denoted by $N = N_l + N_r$, where $N_l$ is the number of steps to the left and $N_r$ the number of steps to the right. The probability that the particle moves $n$ steps to the right is thus given by the binomial distribution
#
# \begin{align*}
# \mathrm{P}_N(N_r = n) &= \binom{N}{n} \left(\frac{1}{2}\right)^n \left(1 - \frac{1}{2}\right)^{N-n} \\
# &= \binom{N}{n} \left(\frac{1}{2}\right)^N \quad.
# \end{align*}
#
# The position of the particle after $N$ steps is equal to the distance travelled to the right minus the distance travelled to the left, assuming that the particle is initially located at $x = 0$
#
# \begin{align*}
# x &= (N_r - N_l)\,h \\
# &= (2N_r - N)\,h \quad,
# \end{align*}
#
# where we have used that $N_l = N - N_r$. Let us plot the probability distribution for a few different values of $N$.
#
# In[1]:
get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import binom
# In[2]:
# Set common figure parameters
newparams = {'font.size': 14, 'lines.linewidth': 1.0, 'figure.figsize': (16, 10),
'legend.frameon': True}
plt.rcParams.update(newparams)
# In[3]:
def calc_pos_distribution(N=10, p=0.5, h=1.0):
""" Creates an array containing the probabilities associated with each position
of the particle after a number N steps.
Parameters:
N: int (non-negative). Number of steps in the time evolution.
p: float (0 <= p <= 1). The probability that the particle moves one step to the right.
h: float. The distance the particle moves during each step.
Returns: (xs, ps)
xs: float array. The possible positions of the particle.
ps: float array. The probabilities of each possible position of the particle.
"""
ns = np.arange(0, N+1)
ps = binom(N, ns) * p**ns * (1 - p)**(N - ns)
xs = 2*ns - N
return (xs, ps)
# In[4]:
fig = plt.figure()
fig.suptitle("Probability distributions of the position of a Brownian particle after $N$ steps with $h = 1$.")
# Handles and data arrays for four subplots
axes = [None] * 4
pos_arrays = [None] * 4
prob_arrays = [None] * 4
# Iterate over each subplot
for i in range(0, 4):
N = 10 * i
pos_arrays[i], prob_arrays[i] = calc_pos_distribution(N)
axes[i] = fig.add_subplot(2, 2, i+1)
axes[i].bar(pos_arrays[i], prob_arrays[i], width=2.0, align='center')
axes[i].set_ylim([0, 1])
axes[i].set_xlim([-20, 20])
axes[i].set_xlabel('$x$')
axes[i].set_ylabel('P($x$)')
axes[i].legend(['$N$ = {0}'.format(N)])
# The expected value of the position of the particle remains at $x = 0$ for all $N$, but as time evolves ($N$ grows bigger), the probability distribution spreads out more and more. This closely resembles the macroscopic phenomenon of *diffusion*, where particles move from a region of high concentration to a region of lower concentration.
# To see this mathematically, we investigate the limit of the binomial distribution where it approaches a continuous probability distribution: Let $N$ increase and let the step size $h$ decrease in such a way that $hN$ is kept constant.
# The number of steps $N$ is related to time $t$ as follows: $N = t/\tau$, where $\tau$ is the average time between collisions of the Brownian particle and a water molecule.
#
# The result is a normal distribution when $N$ gets large (a derivation is given in the appendix)
#
# \begin{equation*}
# \mathrm {P}(x) = \frac{1}{\sqrt{2\pi h^2 t/\tau}} \exp\left({-\frac{1}{2} \frac{x^2}{h^2 t/\tau}}\right) \quad .
# \end{equation*}
#
# The expected value of the position is $x = 0$ at all times, and we have a time-dependent variance $\sigma^2 = h^2 t/\tau$.
# If we introduce the constant $D \equiv h^2/\tau$, the normal distribution can be written in the following form
#
# \begin{equation*}
# \mathrm{P}(x) = \frac{1}{\sqrt{4\pi Dt}} \exp\left({-\frac{1}{2} \frac{x^2}{2Dt}}\right) \quad .
# \end{equation*}
#
# Suppose we have a macroscopic system of $\mathcal{N}$ independent particles.
# The concentration of the particles can then be written as
#
# \begin{equation*}
# n(x, t) = \mathcal{N} \, \mathrm{P}(x) \quad .
# \end{equation*}
#
# A straight-forward calculation shows that $n(x, t)$ satisfies the one-dimensional diffusion equation
#
# \begin{equation*}
# \frac{\partial}{\partial t} n(x, t) = D \frac{\partial^2}{\partial x^2} n(x, t) \quad ,
# \end{equation*}
#
# hence our simple microscopic model gives rise to the macroscopic phenomenon of diffusion!
# ## Numerical simulation of Brownian motion
#
# In order to further acknowledge the connection between Brownian motion and diffusion, it would be helpful to visualise the time development of a system of independent Brownian particles. By independent we mean that the motion of a Brownian particle does not influence the motion of the others. Suppose we have a collection of $\mathcal{N}$ Brownian particles that can move in two dimensions as described in the model above. We are going to plot the positions of the particles after each time step and indeed observe that they tend to move from regions of high concentration to regions of lower concentration.
#
# During the simulation we are also going to compute and display the arithmetic mean and the sample variance of the positions
#
# \begin{align*}
# \overline{x} &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} x_k \\
# s_x^2 &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} (x_k - \overline{x})^2 \quad .
# \end{align*}
#
# The expected values of these quantities are
#
# \begin{align*}
# \langle\overline{x}\rangle &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \langle x_k \rangle \\
# &= 0 \\
# \langle s_x^2 \rangle &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \langle(x_k - \overline{x})^2 \rangle \\
# &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \left(\langle x_k^2 \rangle - \langle \overline{x}\rangle^2\right) \\
# &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} \langle x_k^2 \rangle \\
# &= \frac{1}{\mathcal{N}} \sum_{k = 1}^\mathcal{N} h^2 N \\
# &= h^2 N \quad .
# \end{align*}
# In[5]:
from matplotlib import animation
from IPython.display import HTML
plt.rcParams.update({'animation.html':'html5', 'savefig.dpi': 60})
def populate_particles(nParticles, R=0):
""" Creates two arrays containing randomised positions of the particles.
The positions are randomised within a circle of radius R about the origin.
Parameters:
nParticles: int (positive). The number of Brownian particles.
R: float (positive). The radius within which the randomised positions are confined.
Returns: (xs, ys)
xs: float array. The x-coordinates of the particles.
ys: float array. The y-coordinates of the particles.
"""
angles = 2*np.pi*np.random.rand(nParticles)
radii = R*np.random.rand(nParticles)
xs = radii*np.cos(angles)
ys = radii*np.sin(angles)
return (xs, ys)
def evolve_particles(xs, ys, h=0.1, pr=0.5, pu=0.5, vx=0, vy=0):
""" Evolves the positions of the particles by one time step.
Parameters:
xs: float array. The x-coordinates of the particles.
ys: float array. The y-coordinates of the particles.
h: float. The distance a particle moves during a time step.
pr: unsigned int (0 <= pr <= 1). The probability that a particle moves to the right along the x-axis.
pu: unsigned int (0 <= pr <= 1). The probability that a particle moves upwards along the y-axis.
vx: float. The x-component of the average drift velocity (length h per unit time).
vy: float. The y-component of the average drift velocity (length h per unit time).
Returns: None. The updated positions are stored in xs and ys.
"""
if xs.size != ys.size:
print("Error: The sizes of xs and ys do not match.")
return -1
randnumbers = np.random.random(xs.size)
for i, step_right in enumerate(randnumbers < pr):
if step_right:
xs[i] += h
else:
xs[i] -= h
randnumbers = np.random.random(ys.size)
for i, step_up in enumerate(randnumbers < pu):
if step_up:
ys[i] += h
else:
ys[i] -= h
# Particle drift
xs += vx*h
ys += vy*h
def init_anim():
""" Initialises the animation.
"""
global ax, line, textbox_sim, textbox_theory
line, = ax.plot([], [], 'bo')
ax.set_xlim([-20, 20])
ax.set_ylim([-20, 20])
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Numerical simulation of Brownian motion')
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
# A text box that will display the average and the variance of the positions of the simulated particles.
textbox_sim = ax.text(0.05, 0.95, '', transform=ax.transAxes, fontsize=12,
verticalalignment='top', bbox=props)
# A text box that will display the average and the variance of the theoretical normal distribution of the positions.
textbox_theory = ax.text(0.775, 0.95, '', transform=ax.transAxes, fontsize=12,
verticalalignment='top', bbox=props)
return line, textbox_sim, textbox_theory
def animate(N, xs, ys, h, nParticles, drift_vx=0, drift_vy=0):
""" Runs the simulation and updates the plot.
Parameters:
N: int. The number of time steps taken so far.
xs: float array. The x-coordinates of the particles.
ys: float array. The y-coordinates of the particles.
h: float. The distance a particle moves during a time step.
nParticles: int (positive). The number of Brownian particles.
Returns: Handles to the objects that are going to be updated in the plot.
"""
global ax, line, textbox_sim, textbox_theory
# Evolve the system by one time step.
evolve_particles(xs, ys, h, vx=drift_vx, vy=drift_vy)
N += 1
# Update the plot of the positions.
line.set_data(xs, ys)
# Compute and display quantities associated with the simulated particles.
x_avg, y_avg = np.average(xs), np.average(ys)
x_var = np.sum((xs - x_avg)**2) / nParticles
y_var = np.sum((ys - y_avg)**2) / nParticles
x_avg_str = '$\\langle x \\rangle = %+.2f$' % (x_avg)
y_avg_str = '$\\langle y \\rangle = %+.2f$' % (y_avg)
x_var_str = '$\mathrm{Var}(x) = %.2f$' % (x_var)
y_var_str = '$\mathrm{Var}(y) = %.2f$' % (y_var)
textbox_sim.set_text('Simulation\n' + x_avg_str + '\n' + y_avg_str + '\n' + x_var_str + '\n' + y_var_str)
# Compute and display the theoretical quantities.
x_var = h**2 * N
y_var = h**2 * N
x_avg_str = '$\\langle x \\rangle = %.2f$' % (drift_vx * h * N)
y_avg_str = '$\\langle y \\rangle = %.2f$' % (drift_vy * h * N)
x_var_str = '$\mathrm{Var}(x) = %.2f$' % (x_var)
y_var_str = '$\mathrm{Var}(y) = %.2f$' % (y_var)
textbox_theory.set_text('Theoretical\n' + x_avg_str + '\n' + y_avg_str + '\n' + x_var_str + '\n' + y_var_str)
return line, textbox_sim, textbox_theory
# In[6]:
# The number of Brownian particles.
nParticles = 500
# The total number of steps taken during the simulation.
nSteps = 1000
# Create arrays containing the x and y positions of the particles.
xs, ys = populate_particles(nParticles)
# The distance that the particles move during each time step.
h = 0.1
# Run the simulation and visualise the system as an animation.
fig, ax = plt.subplots()
h_anim = animation.FuncAnimation(fig, animate, init_func=init_anim, frames=nSteps, interval=50, blit=True,
fargs=(xs, ys, h, nParticles))
plt.close(h_anim._fig)
HTML(h_anim.to_html5_video())
# ## Brownian motion with drift
# So far we have looked at a system of independent Brownian particles randomly colliding with the surrounding water molecules. The probability model above thus accounts for the thermal motion of the particles.
# We may however encounter situations where it is necessary to include other mechanisms in the treatment of our system. For instance, the Brownian particles could be ions in a solution subject to an external electric field, making the particles drift in the field direction. Another example would be when we take into account the collisions between the Brownian particles themselves and there is a uniform concentration gradient across the system, resulting in a constant drift velocity.
#
# First, the derivation of a more general expression for the time derivative of the particle concentration is given. This equation can be thought of as an extension of the diffusion equation that also includes the possibility of drift of the particles. Then we shift our focus back on the microphysics and add a particle drift term to the probability model. As we will see, the microscopic model is once again consistent with the macroscopic model – this time described by the drift-diffusion equation.
#
# #### Deriving the drift-diffusion equation
# Let $\rho(\mathbf{r}, t)$ denote the particle density. We assume that the total number of particles in the system is conserved, hence for an arbitrary volume, the change in the number of particles within the volume must equal the number of particles leaving the volume. Mathematically, this can be expressed as
#
# \begin{equation*}
# \int_\mathrm{V} \frac{\partial}{\partial t}\rho(\mathbf{r}, t) \, \mathrm{d}V = - \oint_\mathrm{S} \mathbf{j} \cdot \hat{n}\, \mathrm{d}A \quad ,
# \end{equation*}
#
# where $\mathbf{j}$ denotes the particle flux and the right-hand side of the equation is an integral taken over the closed surface enclosing the volume $V$. Using the divergence theorem, the surface integral can be rewritten as
#
# \begin{equation*}
# \int_\mathrm{V} \frac{\partial}{\partial t}\rho(\mathbf{r}, t) \, \mathrm{d}V = - \int_\mathrm{V} \nabla\cdot {\mathbf{j}} \, \mathrm{d}V \quad .
# \end{equation*}
#
# Since the volume was assumed to be arbitrary, the integrands must be equal, which yields the *continuity equation*
#
# \begin{equation*}
# \frac{\partial}{\partial t}\rho(\mathbf{r}, t) + \nabla\cdot \mathbf{j} = 0 \quad .
# \end{equation*}
#
# This is the starting point of the derivation. We consider two mechanisms that contribute to the particle flux:
#
# 1. The diffusion due to the thermal motion of the particles.
# 2. The drift of the particles.
#
# Fick's first law of diffusion states that $\, \mathbf{j_\mathrm{diff}} = -D \nabla\rho$,
# where the factor $D$ is the diffusion coefficient.
# An intuitive way of thinking about this law is: In regions of higher particle concentration, there will be a higher collision rate than in the regions of lower concentration. It is thus more likely that particles are scattered from regions of higher concentration into regions of lower concentration rather than the opposite direction, resulting in a net particle flux proportional to $-\nabla\rho$.
#
# The drift term is given by $\, \mathbf{j_\mathrm{drift}} = \rho \mathbf{v}$, where $\mathbf{v}$ is the drift velocity of the particles. One can think about this as the average velocity of a group of particles.
# The total particle flux can thus be written as
#
# \begin{align*}
# \mathbf{j} &= \mathbf{j_\mathrm{diff}} + \mathbf{j_\mathrm{drift}} \\
# &= -D \nabla\rho + \rho \mathbf{v} \quad .
# \end{align*}
#
# Inserting the expression for the particle flux into the continuity equation yields
#
# \begin{align*}
# \frac{\partial}{\partial t}\rho(\mathbf{r}, t) &= - \nabla \cdot \mathbf{j} \\
# &= - \nabla\cdot (-D\nabla\rho + \rho\mathbf{v}) \\
# &= \nabla\cdot (D\nabla\rho) - \mathbf{v}\cdot \nabla\rho - \rho \nabla\cdot \mathbf{v} \\
# &= D\nabla^2\rho - \mathbf{v}\cdot \nabla\rho \quad .
# \end{align*}
#
# In the last step we assumed that $D$ is constant and that the system is incompressible ($\nabla\cdot \mathbf{v} = 0$).
# The drift-diffusion equation in one dimension thus reads
#
# \begin{equation*}
# \frac{\partial}{\partial t}\rho(x, t) = D\frac{\partial^2}{\partial x^2}\rho(x, t) - v_x \frac{\partial}{\partial x}\rho(x, t) \quad .
# \end{equation*}
#
# #### Adding a drift term to the probability model
# The position of the Brownian particle with drift is given by $x = (2N_r - N)\,h + v_x t$,
# where $v_x$ is the drift velocity and $t$ is the time parameter.
# It is shown in the appendix that for large $N$, the probability density of $x$ approaches the following normal distribution
#
# \begin{equation*}
# \mathrm{p}(x) \approx \frac{1}{\sqrt{4\pi Dt}} \exp\left({-\frac{1}{2} \frac{(x-v_x t)^2}{2Dt}}\right) \quad ,
# \end{equation*}
#
# which satisfies the drift-diffusion equation!
#
#
#
# ## Numerical simulation of Brownian motion with drift
# In[7]:
# The number of Brownian particles.
nParticles = 500
# The total number of steps taken during the simulation.
nSteps = 250
# Create arrays containing the x and y positions of the particles.
xs, ys = populate_particles(nParticles)
# The distance that the particles move during each time step.
h = 0.1
# The components of the average drift velocity (length h per unit time)
drift_vx = 0.5
drift_vy = 0
# Run the simulation and visualise the system as an animation.
fig, ax = plt.subplots()
h_anim = animation.FuncAnimation(fig, animate, init_func=init_anim, frames=nSteps, interval=50, blit=True,
fargs=(xs, ys, h, nParticles, drift_vx, drift_vy))
plt.close(h_anim._fig)
HTML(h_anim.to_html5_video())
# ## Conclusion
# We have demonstrated how a simple probability model of the motion of Brownian particles can be used to explain why particles tend to move towards regions of lower concentration. The microscopic probability model with an added drift term is consistent with the macroscopic drift-diffusion equation.
# ### Further work and improvements
# * Create a more sophisticated simulation model of this phenomenon. For instance, instead of just assigning probabilities that the particles move in one direction or another, model the water molecules as spherical particles of finite size with a given velocity distribution and simulate elastic collisions between the water molecules and the Brownian particles. Using this approach one may also allow for collisions between the Brownian particles themselves. Will the results be the same?
# ## Appendix
#
# ### The binomial distribution approaches a normal distribution when $N$ grows large and $hN$ is held constant.
#
# The very essence of this derviation is to look at the moment-generating function of the binomially distributed position of the particle, then let $N$ grow large and see that we end up with the moment-generating function of a normal distribution. We shall first state a few general results.
#
# The moment-generating function of a random variable $X$ is defined by
#
# \begin{equation*}
# \mathrm{M}_X(k) = \langle e^{kX} \rangle \quad .
# \end{equation*}
#
# For a linear transformation $aX + b$, where $a$ and $b$ are constants, we have that
#
# \begin{equation*}
# {\mathrm M}_{aX+b}(k) = e^{-bk} \, {\mathrm M}_X(ak) \quad .
# \end{equation*}
#
# #### A theorem of moment-generating functions
# Suppose that $X$ and $Y$ are random variables for which ${\mathrm M}_X(k)$ = ${\mathrm M}_Y(k)$
# for some interval containing $k = 0$. Then the probability distributions of $X$ and $Y$ are equal. [3]
#
# #### The moment-generating functions of the binomial distribution and the normal distribution [4]
# For a binomial distribution of $n$
#
# \begin{equation*}
# \mathrm{P}_N(n) = \binom{N}{n} p^n (1 - p)^{N-n} \quad ,
# \end{equation*}
#
# the moment-generating function reads
#
# \begin{equation*}
# \mathrm{M}_n(k) = (p e^k + 1-p)^N \quad .
# \end{equation*}
#
# For a normal distribution of $x$
#
# \begin{equation*}
# \mathrm{p}(x) = \frac{1}{\sqrt{2\pi} \sigma} \exp\left(-\frac{1}{2} \frac{(x-\mu)^2}{\sigma^2}\right) \quad ,
# \end{equation*}
#
# the moment-generating function reads
#
# \begin{equation*}
# \mathrm{M}_x(k) = e^{\mu k + \sigma^2 k^2 / 2} \quad .
# \end{equation*}
#
# #### From binomial distribution to normal distribution
# Recall that the position of the Brownian particle was given by $x = (2N_r - N)\,h$.
# With $p = 1/2$ we have that
#
# \begin{align*}
# \mathrm{M}_x(k) &= e^{-hNk} \, {\mathrm M}_{2hN_r}(k) \\
# &= e^{-hNk} \left( \frac{1}{2} e^{2hk} + \frac{1}{2} \right)^N \\
# &= e^{-hNk} \left( \frac{1}{2} e^{2hNk/N} + \frac{1}{2} \right)^N \\
# &= e^{-yk} \left( \frac{1}{2} e^{2yk/N} + \frac{1}{2} \right)^N \quad ,
# \end{align*}
#
# where $y \equiv hN$ is held constant. As $N$ grows large, we do the following approximation
#
# \begin{equation*}
# e^{2yk/N} \approx 1 + \frac{2yk}{N} + \left(\frac{yk}{N}\right)^2 \quad ,
# \end{equation*}
#
# which implies that
#
# \begin{equation*}
# \mathrm{M}_x(k) \approx e^{-yk} \left( 1 + \frac{yk + y^2 k^2 / 2N}{N} \right)^N \quad .
# \end{equation*}
#
# We also make use of another approximation
#
# \begin{equation*}
# \left( 1 + \frac{u}{N} \right)^N \approx e^u \quad ,
# \end{equation*}
#
# which is valid when $\frac{u}{N} \ll 1$, hence for large $N$
#
# \begin{align*}
# \mathrm{M}_x(k) &\approx e^{-yk} e^{yk + y^2 k^2 / 2N} \\
# &= e^{(y^2/N) k^2 / 2} \quad .
# \end{align*}
#
# This we immediately recognise as the moment-generating function of $x$ for a normal distribution with $\mu = 0$ and
# $\sigma^2 = y^2/N$. The domain of this moment-generating function contains $k = 0$, hence as $N$ grows large and $hN$ is held constant
#
# \begin{align*}
# \mathrm{p}(x) &\approx \frac{1}{\sqrt{2\pi y^2/N}} \exp\left(-\frac{1}{2} \frac{x^2}{y^2/N}\right) \\
# &= \frac{1}{\sqrt{2\pi h^2 t/\tau}} \exp\left(-\frac{1}{2} \frac{x^2}{h^2 t/\tau}\right) \\
# &= \frac{1}{\sqrt{4\pi Dt}} \exp\left({-\frac{1}{2} \frac{x^2}{2Dt}}\right) \quad .
# \end{align*}
#
# An alternative and more direct approach using Stirling's approximation can be found in [5].
#
# #### Adding a drift term
# We now consider Brownian motion with drift.
# The position of the Brownian particle is then given by $\tilde{x} = (2N_r - N)\,h + v_x t$,
# where $v_x$ is the drift velocity and $t = N\tau$ is the time parameter.
# The term $v_x t$ is a constant (not a random variable), hence the moment-generating function of $\tilde{x}$ becomes
#
# \begin{equation*}
# \mathrm{M}_\tilde{x}(k) = e^{-v_x tk} \, {\mathrm M}_{(2N_r-N)h}(k) \quad.
# \end{equation*}
#
# Using our previous result for ${\mathrm M}_{(2N_r-N)h}(k)$ for large $N$ yields
#
# \begin{equation*}
# \mathrm{M}_\tilde{x}(k) \approx e^{-v_x tk + (y^2/N) k^2 / 2} \quad ,
# \end{equation*}
#
# which we recognise as the moment-generating function of $x$ for a normal distribution with $\mu = v_x t$ and $\sigma^2 = y^2/N$.
#
# The probability density of the position of the Brownian particle with drift is thus
#
# \begin{equation*}
# \mathrm{p}(x) \approx \frac{1}{\sqrt{4\pi Dt}} \exp\left({-\frac{1}{2} \frac{(x-v_x t)^2}{2Dt}}\right) \quad .
# \end{equation*}
# ## References
# [1] R. Brown, "A brief account of microscopical observations made in the months of June, July and August, 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies", *The miscellaneous botanical works of Robert Brown*, (The Ray Society, London, 1866), pp. 463. Available: [The miscellaneous botanical works of Robert Brown](https://archive.org/details/miscellaneousbot01brow) \[11.09.2016\].
#
# [2] M. Smoluchowski, "Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen", *Annalen der Physik* **21** (14), pp. 756 (1906).
#
# [3] R. J. Larsen & M. L. Marx, *An introduction to Mathematical Statistics and Its Applications*, 5th Edition, 2012, Pearson, p. 214.
#
# [4] Institutt for matematiske fag NTNU, *Tabeller og formler i statistikk*, 6. opplag, 2011, Fagbokforlaget, s. 25.
#
# [5] R. K. Pathria & P. D. Beale, *Statistical Mechanics*, 3rd Edition, 2011, Butterworth–Heinemann, p. 588.