Exercises (6)

Assigned: Feb 21, 2012
Due: March 6, 2012

Rolling Dice (in-class)

What is the probability of throwing 8 or less with three dice, each with six sides numbered 1 through 6? Write a Mathematica simulation (with at least 1000 trial dice rolls) to numerically estimate the probability. (Can you estimate the probability analytically?)

Solution

diceRoll[]:=RandomInteger[{1,6}]

Table[diceRoll[], {10}]

eightOrLess[] := If[diceRoll[] + diceRoll[] + diceRoll[] ≤ 8, 1, 0]

rolls = Table[eightOrLess[], {1000}];
(Plus @@ rolls)/1000 // N

This probability can be estimated by counting all of the ways to roll 8 or less and all of the ways to roll three dice.

The number of ways to roll eight with three dice:
	1, 1, 6		2, 1, 5		3, 1, 4		4, 1, 3		5, 1, 2		6, 1, 1
	1, 2, 5		2, 2, 4		3, 2, 3		4, 2, 2		5, 2, 1
	1, 3, 4		2, 3, 3		3, 3, 2		4, 3, 1
	1, 4, 3		2, 4, 2		3, 4, 1
	1, 5, 2		2, 5, 1
	1, 6, 1
is 6 + 5 + 4 +3 + 2 +1 = 6 (6 + 1)/2 = 21
The number of ways to roll seven with three dice:
	1, 1, 5		2, 1, 4		3, 1, 3		4, 1, 2		5, 1, 1
	1, 2, 4		2, 2, 3		3, 2, 2		4, 2, 1
	1, 3, 3		2, 3, 2		3, 3, 1
	1, 4, 2		2, 4, 1
	1, 5, 1
is 5 + 4 + 3 + 2 + 1 = 5 (5 + 1)/2 = 15
The number of ways to roll six is 4(4 + 1)/2 = 10.
The number of ways to roll five is 3(3 + 1)/2 = 6.
The number of ways to roll four is 2(2 + 1)/2 = 3.
The number of ways to roll three is 1.
Therefore, the number of ways to roll eight or less is 21 +15 + 10 + 6 + 3 + 1 = 56.

The total number of ways to roll three dice is (6^3)= 216.

56/216 // N

Sum[n(n+1)/2, {n, 1, 6}]/(6^3)

% // N

Langevin Equations for Brownian Motion (in-class)

The Langevin equation for Brownian motion can be written as a sequence of random steps. Each time step, random motion of molecules cause particles to change their velocities. Friction with the surroundings cause the velocity to slow back towards zero. This is a model for the random motion of fine dust particles on the surface of liquid, for example.

The Langevin equation describes the trajectory of an object with a random force, s, and a frictional damping, γ. In differential form, the Langevin equation is
	(∂_t) v = - γ v + s
	(∂_t) x = v
The Langevin equation can be written in finite-difference form as a two-dimensional map:
	(Subscript[v, n]) = (1 – γ Δt) (Subscript[v, n-1]) + (Subscript[s, n])
	(Subscript[x, n]) = (Subscript[x, n-1]) + Δt (Subscript[v, n-1])
This map has a similar form to the deterministic equations for planetary orbits and pendulums. The big difference here is that the change in velocity cause by (Subscript[s, n])is random!

A trajectory following the Langevin equations can be found from the map below, where γ = 0.1, Δt = 1, and where (Subscript[s, n])= ± 1.

langevin[{x_,v_}]:={x+v,0.9 v+(2 RandomInteger[]-1)}

langevin[langevin[{0, 10}]]

Part a:
Construct a simulation using langevin[…] and NestList[...] beginning with {(Subscript[x, 0]), (Subscript[v, 0])} = {0, 10} for about 50 steps and for about 1000 "walkers".

Part b:
What is the evolution of the average position and average velocity?

Part c:
Plot the mean-squared position and velocity of the walkers. Does the velocity and position diffuse? If so, can you estimate the diffusion coefficient?

Solution

Part a

nWalkers = 1000;

walkList = Table[NestList[langevin, {0, 10}, 50], {nWalkers}];

walkList // Dimensions

Part b

avgX = Table[Mean[walkList[[All, n, 1]]], {n, 1, 51}];

ListPlot[avgX,AxesLabel->{"n","⟨x⟩"}]

avgV = Table[Mean[walkList[[All, n, 2]]], {n, 1, 51}];

ListPlot[avgV,AxesLabel->{"n","⟨v⟩"}]

Part c

avgΔX2 = Table[Variance[walkList[[All, n, 1]]], {n, 1, 51}];

ListPlot[avgΔX2,AxesLabel->{"n","⟨(x - ⟨x⟩\!\(\*SuperscriptBox[\()\), \(2\)]\)⟩"}]

The particle trajectories begin very little spatial diffusion, and later diffuse strongly.

avgΔV2 = Table[Variance[walkList[[All, n, 2]]], {n, 1, 51}];

ListPlot[avgΔV2,AxesLabel->{"n","⟨(v - ⟨v⟩\!\(\*SuperscriptBox[\()\), \(2\)]\)⟩"},PlotRange->All]

The velocity initially diffuses, and then finds an equilibrium where the average velocity fluctuations are constant in time.

Part c (Estimating the diffusion coefficient)

gridX = Table[BinCounts[walkList[[All, n, 1]], {-0-1/2, +100 + 1/2,1}], {n, 1, 51}];

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

gridV = Table[BinCounts[walkList[[All, n, 2]], {-20-1/2, +20 + 1/2,1}], {n, 1, 51}];

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

An estimate of the spatial diffusion coefficient is found from the linear part of the ⟨((Δx)^2)⟩ curve, where ⟨((Δx)^2)⟩ ∝ 2 D t.

Fit[avgΔX2[[Range[30,50]]], {1, t}, t]

Therefore, D ~ 50.

Notice also that ⟨((Δv)^2)⟩ initially increases with time. The velocity diffusion coefficient is found from

Fit[avgΔV2[[Range[1,4]]], {1, t}, t]

Therefore, (Subscript[D, v]) ~ 0.4, which is about D (γ^2). The equilibrium mean-square velocity is ⟨(v^2)⟩ ~ γ D ~ (Subscript[D, v])/γ. This relationship between the diffusion and dissipation was first described by Albert Einstein.

The Game of Craps

Create a simulation and estimate the probability that you will win in a "game of craps."

The game of craps is a game played by rolling two dice (numbered 1 to 6). If you roll 2, 3, or 12 on the first try, you lose. If you roll 7 or 11 on the first try, you win. If you roll anything else (i.e. 4, 5, 6, 8, 9, 10), to win you must roll that number again before you roll 7. If you roll 7 first, then you lose. You can roll as many times as you need until you either win or lose.

Solution

crapRoll[] := diceRoll[] + diceRoll[];

crapWin[] := Module[{roll, nextRoll},
roll = crapRoll[];
If[(roll == 2) || (roll == 3) || (roll == 12), Return[0]];
If[(roll == 7) || (roll == 11), Return[1]];
While[nextRoll = crapRoll[];(nextRoll ≠ roll) && (nextRoll ≠ 7)];
If[nextRoll == roll, 1, 0]]

Table[crapWin[], {10}]

(Plus @@ Table[crapWin[], {1000}])/1000 // N

Two-Dimensional Random Walk

For the two-dimensional random walk, particles make unit steps in one of four orthogonal directions: up, down, right, or left. A simple 2D walker function is

step2D[]:={{0,1},{1,0},{0,-1},{-1,0}}[[RandomInteger[{1,4}]]]

Table[step2D[], {5}]

If a walker starts at the origin {x, y} == {0, 0}, then the following command generates a trajectory with n steps.

walk2D[n_] := NestList[(# + step2D[])&, {0,0}, n]

walk2D[5]

Question: what is the mean distance of a two-dimensional random walk as a function of n? (Hint, define a function radius[{x_, y_}] that returns the radius from the origin, and Map[...] this function to every coordinate-pair of each walker's trajectory.) You should find that the mean radius increases in proportion to (Sqrt[n]).

Solution

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

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

Average position ("center of mass")

avgXY[walkers_List] := N[(Plus @@ walkers)/Length[walkers]]

avgXYList = avgXY /@ Transpose[walkList];

avgXYList // Shallow

ListPlot[avgXYList,PlotRange->{{-1,1},{-1,1}},AspectRatio->Automatic,Joined->True,PlotLabel->"Average Position of All Walkers vs Time"]

The average position of all walkers remains near the origin (as expected).

Average radius

radius[{x_, y_}] := Sqrt[N[x^2 + y^2]]

rList = Map[radius, walkList,{2}];

rList//Shallow

avgR=Mean/@Transpose[rList];
ListPlot[avgR]

Fit[avgR, {(Sqrt[n])}, n]