Random Walks

AP 1601y
Columbia University

Introduction

In this notebook, we begin a study of random processes and statistical physics. We explore a simple model for diffusion by considering a simple model for the unbiased, one-dimensional random walk. Diffusion is the process by which a quantity (like particles, heat, or momentum) moves from a region of high concentration to low concretionary by a multi-step random process.  We also briefly introduce the concept of entropy, introduced by Rudolf Clausius in 1865.

Jacob Bernoulli

Jacob Bernoulli was a Swiss mathematician who lived from 1654 to 1705. Online biographies are available from St. Andrews in Scotland and at scienceworld.wolfram.com. Bernoulli's work on combinatorics and probability forms the foundation of today's lecture.

One Dimensional Random Walks

We consider a particle that takes either a right-going step (+1) or a left-going step (-1) along the x-axis. This is the so-called "random walk" that often forms the starting point in a study of statistical physics.

s[]:=2 RandomInteger[{0,1}]-1

Table[s[], {10}]

Taking 500,000 steps is easy…

(w1 = NestList[(# + s[])&, 0, 500000]; )// Timing

ListPlot[w1[[Range[1,500000,1000]]],AxesLabel->{"n/1000","x/Δx"}]

walk1D[n_] := NestList[(# + s[])&, 0, n]

One thousand walkers…

We take 500,000 steps…

nWalkers = 1000;
walkList = Table[walk1D[500],{nWalkers}];

walkList // Shallow (* 1000 walks 500 steps long *)

pos[n_] := walkList[[All, n+1]] (* walker locations after n steps *)

pos[100] // Shallow

Manipulate[ListPlot[pos[n], PlotRange -> {-50, 50}, AxesLabel -> {"walker", "x/Δx, n= " <> ToString[n]}], {n, 0, 500}]

BinCounts

Mathematica has useful tools for statistical analysis. BinCounts[ ] is one such function.

?BinCounts

BinCounts[{0,1,3,1,5,5,4,3,3,2}, {-1/2, 5 + 1/2, 1}]

BinCounts[pos[0], {-50-1/2, +50 + 1/2,1}]

BinCounts[pos[100], {-50-1/2, +50 + 1/2,1}]

grid = Table[BinCounts[pos[n], {-50-1/2, +50 + 1/2,1}], {n, 0, 500}];

ListDensityPlot[grid,Mesh->False,ColorFunction->Hue,PlotRange->{1,100},PerformanceGoal -> "Speed",MaxPlotPoints -> 100,FrameLabel->{"Box(x/Δx), "<>ToString[nWalkers]<>" Walkers","n"}]

p1=ListPlot[grid[[25,All]],PlotRange->All,PlotStyle->{RGBColor[1,0,0],PointSize[0.01`]},DisplayFunction->Identity];
p2=ListPlot[grid[[250,All]],PlotRange->All,PlotStyle->{RGBColor[0,0,1],PointSize[0.01`]},DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction,AxesOrigin->{50,0},AxesLabel->{"Box(x/Δx)","Count"}]

Histogram

?Histogram

p3=Histogram[pos[25],{-50,50,1},PlotLabel->"h[j,n], n= "<>ToString[100],AxesLabel->{"x/Δx",""}]

Manipulate[Histogram[pos[n], {-50,50,1},PlotLabel -> "h[j,n], n= " <> ToString[n], AxesLabel -> {"x/Δx",""}], {n, 0, 500}]

Mean and Standard Deviation

?Mean

Mean /@ Transpose[{{a1, a2, a3, a4}, {b1, b2, b3, b4}, {c1, c2, c3, c4}}]

avgX=Mean/@Transpose[walkList];
ListPlot[avgX,AxesLabel->{"n","<x>/Δx"}]

avgX2=Mean/@Transpose[(walkList^2)];
ListPlot[avgX2,AxesLabel->{"n","<(x^2)>/((Δx)^2)"}]

?Variance

varX=Variance/@Transpose[walkList];
ListPlot[varX,AxesLabel->{"n","Var(x/Δx)"}]

?StandardDeviation

StandardDeviation /@ Transpose[{{a1, a2, a3, a4}, {b1, b2, b3, b4}, {c1, c2, c3, c4}}]

stdDevX=StandardDeviation/@Transpose[walkList];
ListPlot[stdDevX,AxesLabel->{"n","StdDev(x/Δx)"}]

The Binomial Distribution

We can develop a mathematical model for the probability for the random walk problem.

Let's begin by making a "decision tree". If the probability of moving to the right (+1) is equal to the probability of moving to the left (-1), then the probability of reaching a location, x, after n steps is given by P(x,n) = (number of paths arriving at x)/(total number of possible paths). This probability distribution is called the binomial distribution.

The probability of being located at x = 0 after two steps is P(0,2) = 1/2, since there are four possible paths after two steps and two possible ways of reaching x = 0.  (Notice that, for n even, only even positions can be populated, and, for n odd, only odd positions can be populated.)

Now, let's consider an arbitrary number of n steps.

After n steps, the total number of steps to the right, (Subscript[n, R]), and the total number of steps to the left, (Subscript[n, L]),must satisfy, n  = (Subscript[n, R]) + (Subscript[n, L]).In order to reach location j, the difference between right and left directed steps must satisfy, j = (Subscript[n, R]) –  (Subscript[n, L]), or j = 2 (Subscript[n, R]) – n. For example, if j = –3 after 15 steps, we must have (Subscript[n, R]) = 6 and (Subscript[n, L])= 9. It doesn't really matter what order the right and left steps are taken. As long as there are 6 right steps and 9 left steps, j = –3 after 15 steps. The probability for any of these possible paths must be equal to (probability of moving ((right))^(Subscript[n, R]))× (probability of moving ((left))^(Subscript[n, L])) = (((1/2))^(Subscript[n, R]))(((1/2))^(Subscript[n, L])) = (((1/2))^n).

How many paths are there for a given (j, n)?

If there were n distinct elements of the list of steps (i.e. {+1, -1, -1, +1, -1, +1, -1, -1, +1, -1, -1, +1, -1, -1, +1}), then there would be n! = n (n  – 1) (n – 2) … ways of re-arranging these elements, However, the members within the group of all right steps and the group of all left steps are identical. These different re-arrangements don't count. Therefore,
	(number of possible paths) = n! / ((Subscript[n, R])!(Subscript[n, L])!), assuming that n + j  is even and |j| < n.

We can check this model graphically…

prob[n_, j_] := Module[{nR, nL},
nR = (n + j)/2;
nL = n - nR;
If[OddQ[n + j] || Abs[j] > n, 0,
(1/2)^nR (1/2)^nL n!/ nR! / nL!]];

Plot[{prob[5,x],prob[50,x]},{x,-50,50},PlotRange->{0,0.5},AxesLabel->{"x/Δx","prob"},PlotStyle->{{RGBColor[1,0,0]},{RGBColor[0,0,1]}}]

Manipulate[Plot[nWalkers prob[n,x], {x, -n-1, n+1}, PlotRange -> {{-50,50},{0, 500}}, PlotStyle -> Thickness[0.01],DisplayFunction -> Identity], {n, 1, 50}]

The number of possible paths in the decision tree for a one-dimensional random walk is a binomial coefficient.

In Mathematica, the number  n! / ((Subscript[n, R])!(Subscript[n, L])!) =  n! / ((Subscript[n, R])!(n - (Subscript[n, R]))!) is also called "n choose (Subscript[n, R])" or Binomial[n,nR]. Mathematica can compute Binomial[...] many 100's of times faster than  n!/(nR! (n - nR)! especially as n becomes large!

(Binomial[n,nR] == n!/(nR! (n - nR)!)) /. n -> 10 /. nR -> 4

Do[Binomial[n,3], {n,5000}]; // Timing

Do[n!/(3! (n - 3)!), {n,5000}]; // Timing

Averages

We can use the binomial probability distribution function to compute averages of various properties of many random walkers.

The average value of x/Δx is

probX[n_] := Sum[prob[n,x] x, {x,-n, n}]

probX[10]

probX2[n_] := Sum[prob[n,x] x^2, {x,-n, n}]

probX2[10]

ListPlot[Table[{n,probX2[n]},{n,0,100,10}],AxesOrigin->{0,0},AxesLabel->{"n","<(x^2)>/((Δx)^2)"},Joined->True]

Diffusion Equation

If we have a distribution of walkers (or "particles"), how does this distribution evolve in time?

If we consider time and space to change by small steps, Δt and Δx, then we can advance a distribution of particles with the equation
	ρ[x, t + Δt] = 0.5 ρ[x - Δx, t] + 0.5 ρ[x + Δx, t]

If we make a series expansion, this becomes the diffusion equation...

Series[ρ[x,t+Δt], {Δt, 0, 1}]== 0.5 Series[ ρ[x-Δx,t]+ ρ[x+Δx,t], {Δx, 0, 3}] // Simplify

% // Normal

This is the diffusion equation, 
	(∂_t) ρ = D (∂_(x,x))ρ
where D = 0.5 ((Δx)^2)/Δt is the diffusion coefficient.

Entropy

The concept of entropy was developed by the German physicist Rudolf Clausius in 1865. He explained that physical processes always increase the disorder of the system (and of the universe.) The increase in entropy sets the direction of time and defines the second law of thermodynamics. In this section, we will examine and measure the disorder of our collection of random walkers.

Imagine that we start nWalkers at x = 0. What is the probability that they will re-occur back at the origin? Answer: (((P(n, 0)))^nWalkers).

For two walkers and 10 steps, the probability that both will return to 0 after 10 steps is 6%,...

prob[10,0]^2 // N

ListPlot[Table[Log[10,prob[n,0]],{n,2,200,2}],AxesLabel->{"n/2","Log[10,P[n,0]]"}]

The number of times that two particles will meet (at any location) after n steps is

probCollide[n_] := Sum[prob[n,x] prob[n,x], {x, -n, n}]

probCollide[10]

ListPlot[Table[probCollide[n],{n,1,100}],AxesOrigin->{0,0},AxesLabel->{"n","Prob Collide after n"}]

The total number of times a collision will occur (after n steps) is equal to

numCollide[n_] := Sum[probCollide[i], {i, 1, n}]

numCollide[2]

ListPlot[Table[{n,numCollide[n]},{n,1,100,5}],AxesOrigin->{0,0},PlotStyle->PointSize[0.02`],AxesLabel->{"n","Num Collisions after n Steps"}]

This scales like (Sqrt[n]). So two random walkers (in one-dimension) will collide many, many times as they wander!

With 100 walkers and 10 steps, the probability of returning to the original position is

prob[10,0]^100 // N

Requiring many times "longer than the age of the universe."

Entropy

Entropy is a measure of the disorder of the distribution of random walkers. It is defined as S(t) = – ∫P(x,t)  Log[P(x,t)]dx, where P(x,t) is the probability of a walker being at cell x at time t. Initially, all of the walkers are located at x = 0 with probability of 100%, The entropy is zero. As the probability become more uniform in x, then the entropy will increase. According to the ergodic hypothesis, all of the available states of a system in equilibrium will be occupied with equal probability. For our random walkers, x can range from -∞ to +∞ with equal probability as t -> ∞. Therefore, our entropy calculation below is really only valid for a limited time when there are enough walkers within each cell in order to estimate the probability as the ratio, P(x,t) ≈ (# of walkers within grid box x)/(total number of walkers).

(You can show that the diffusion equation implies that entropy must always increase. A  question remains: the diffusion equation is valid only as the number of walkers -> ∞. For a finite number of walkers, fluctuations from the average are large, and sometimes these fluctuations lead to a momentary decrease in entropy.)

We can compute the entropy evolution of our random-walker simulation easily…

{nWalkers, Plus @@ grid[[1,All]],Plus @@ grid[[500,All]]} (* Notice: our "grid" begins to be too small to record all of the walkers. *)

entropyGrid = Table[Plus @@ (If[# > 1,-N[(#/nWalkers) Log[(#/nWalkers)]],0]& /@  grid[[t,All]]), {t, 1, 500}];

ListPlot[entropyGrid,AxesLabel->{"n","Entropy[n], "<>ToString[nWalkers]<>" Walkers"},Joined->True,PlotRange->{0,4}]

Entropy with only 10 Walkers

When there is a large number of walkers, entropy (practically) always increases. But, when the number of walkers decreases (see nWalkers = 10 below), then random fluctuations are large, and the entropy occasionally decreases. When entropy decreases, then the particles tend to return to their originally (ordered) positions.

nWalkers2 = 10;
walk2List = Table[walk1D[500],{nWalkers2}];
grid2 = Table[BinCounts[walk2List[[All, n+1]], {-50-1/2, +50 + 1/2,1}], {n, 0, 500}];
{nWalkers2, Plus @@ grid2[[1,All]],Plus @@ grid2[[500,All]]}

ListDensityPlot[grid2,Mesh->False,ColorFunction->Hue,PlotRange->{0,5},PerformanceGoal -> "Speed",MaxPlotPoints -> 100,FrameLabel->{"Box(x/Δx), "<>ToString[nWalkers2]<>" Walkers","n"}]

entropyGrid2=Table[Plus@@(If[#1>1,-N[((#1 Log[((#1)/nWalkers2)])/nWalkers2)],0]&)/@grid2[[t,All]],{t,1,500}];
ListPlot[entropyGrid2,AxesLabel->{"n","Entropy[n], "<>ToString[nWalkers2]<>" Walkers"},Joined->True,PlotRange->{0,3}]

Summary

We have investigated both one-dimensional and two-dimensional models of random walks. A random walk is a microscopic  description of the macroscopic process of diffusion.

We find the average position (from the starting point) of the walkers increase as the square root of the number of steps,  ⟨x⟩ ~ (Sqrt[n]), or ⟨(x^2)⟩  ~ n . This behavior is observed in all models of random walks and diffusion.

The probability distribution for a random walk process is described by the Binomial Distribution.