Sample code 7.1. To generate repeatable random numbers, use SeedRandom[X], where X any integer.

For[i=1,i≤5,i++,
Print[i, " ",RandomReal[]]];(*end for i*)

For[i=1,i≤5,i++,
SeedRandom[2];
Print[i, " ",RandomReal[]]];(*end for i*)

Sample code 7.2. Random number is generated using RandomReal[] distributed uniformly in [0,1]. To generated 25 random number distributed according to Gaussian distribution NormalDistribution[μ,σ]
can be generated using RandomVariate[NormalDistribution[μ,σ],25].

data=RandomVariate[NormalDistribution[4,10],25];
Print[data]

Sample code 7.3. We can bin the 10^4 random number generated with RandomVariate[NormalDistribution[1,3],10^4] into 20 bins. This can be acchieved using the command Histogram[data,20,”ProbabilityDenstiy”].

data=RandomVariate[NormalDistribution[4,10],10^4];Histogram[data,100,"ProbabilityDensity"]

Sample code 7.3.1. Plot[PDF[NormalDistribution[4,10],x],{x,-40,40},PlotStyle->Thick] will plot the actual Gaussian Distribution with μ=4,σ=10 for x in [-40,40]
Overlapping the histogram and the Gaussian distribution curve, we can then check that indeed the random numbers are distributed according to the guassian distribution with the corresponding parameters  μ=4,σ=10.

Show[
Histogram[data,100,"ProbabilityDensity"],
Plot[PDF[NormalDistribution[4,10],x],{x,-40,40},PlotStyle->Thick]
](*end show*)

Sample code 7.4. 1-D random walk simulation involving only A SINGLE walker, simulated for a time span of nlast timestep.

(* 1D random walk for 1 walker *)
(*Initialisation*)
nlast=10000;
x[0]=0;sumx=0;sumxsq=0;
(*End of Initialisation*)

For[n=1,n≤nlast,n++,
(*ran=Random[Integer,{-1,1}];*)
ran=Random[Real];
If[
ran<0.5,x[n]=x[n-1]-1,x[n]=x[n-1]+1;
];(*end if *)
sumx=sumx+x[n];
avex[n]=sumx/(n+1);
sumxsq=N[sumxsq+(x[n]^2)];
avexsq[n]=sumxsq/(n+1);
];(*end for n*)
tnx=Table[{n,x[n]},{n,0,nlast}];
tnavex=Table[{n,avex[n]},{n,0,nlast}];
ListPlot[tnx,Joined->True]
ListPlot[tnavex,Joined->True]

Sample code 7.5. 1-D random walk simulation involving jlast walkers, simulated for a time span of nlast timestep.  Basically, code 7.5 and 7.4 are similar except that here the code evaluates jlast walkers instead of only one, as in code 7.4.  The code is obtained by modifying sample code 7.4 by adding an additional loop, indexed by j. We wish to know the following quantities at each time step n:
(i) avex[n] = Sum[x[j],{j,1,jlast}]/jlast, at time n.
(ii) Average displacement squared at timestep n,avexsq[n]=Sum[x[j]^2,{j,1,jlast}]/jlast; 

According to analysis on the 1-D jlast random walker, the variance of position will depend on time via the relation =, where l = 1, the spatial step the walkers are shifted in each simulation step. The variance should be proportional to n with a proportional constant = 1 is a theoretical prediction.

The variance, σ^2, on other other hand can be calculated in our simulation since it is given by σ^2= avexsq[n]. 
We will use simulation to find avexsq[n] as a function of n (=t). We will show that the result of the simulation shows that avexsq[n] = 2Dn, where D =1/2. (this constant is called the diffusion coeffience). n denotes time (in unit of multiple of 1).  This is consistent with what theoretical analysis has predicted.

(*Sample code 7.4.2*)

(* Random walk of a collection of jlast walkers in 1D *)
(*Initialisation*)
nlast=500;jlast=2000;(*total number of walkers*)
L=100;L1=-L;L2=L;
xc=IntegerPart[(L1+L2)/2];

For[j=1,j≤jlast,j++,
x[j]=xc;
];(*end of j*)
avex[0]=Sum[x[j],{j,1,jlast}]/jlast;
avexsq[0]=Sum[x[j]^2,{j,1,jlast}]/jlast;
tnx[0]=Table[{x[j],0},{j,1,jlast}];
tx[0]=Table[x[j],{j,1,jlast}];
(*End of Initialisation*)

For[n=1,n≤nlast,n++,
For[j=1,j≤jlast,j++,
ran=Random[Real];
If[ran≤0.5,x[j]=x[j]-1
,(*else*)
x[j]=x[j]+1
](*end if *)
];(*end for j*)

tnx[n]=Table[{x[j],0},{j,1,jlast}];
tx[n]=Table[x[j],{j,1,jlast}];

(*evaluate the average of x[j] and xsq[j] at fixed n*)
avex[n]=Sum[x[j],{j,1,jlast}]/jlast;
avexsq[n]=Sum[x[j]^2,{j,1,jlast}]/jlast;
(*evaluate the average of x[j] and xsq[j] at fixed n*)
];(*end for n *)
tnavex=Table[{n,avex[n]},{n,0,nlast}];
tnavexsq=Table[{n,avexsq[n]},{n,0,nlast}];
lptnavexsq=ListPlot[tnavexsq,PlotLabel->"avexsq, vs n"];

yth[n_]:=2*Dx*n;Dx=0.5;
tny=Table[{n,yth[n]},{n,0,nlast}];
lptny=ListPlot[tny,PlotLabel->"yth vs n",PlotStyle->Red];
Show[{lptny,lptnavexsq}]
(*Manipulate[
ListPlot[tnx[n],PlotRange->{{L1,L2},{-1,1}}],
{n,0,nlast,1}]*)
Manipulate[Show[Histogram[tx[n],50]],{n,0,nlast,1}];

Manipulate[
ListPlot[N[BinCounts[tx[n],{L1,L2,2}]/jlast],PlotRange->{0,1},Filling->Axis]
,{n,0,nlast,1}]

Sample code 7.6. 2-D random walk simulation involving jlast walkers, simulated for a time span of nlast timestep. We would like to visual the walkers on a 2D plan, using the Manipulate function of mathematica. We assume all walkers are located at r0={x0,y0} at n=0.

After the 2D random walk is done, we have collected all the various value of r[j] for all n. Based on these values we can caculate easily <r^2>, which is stored in tnrsqave. The graph of <r^2> is plotted against n to show a linear relation with n. It is easily found that the slop for the 2D random walk in <r^2> vs. n is given by <r^2> = 2*D*n, where D = 1 in this case.

(*2D random walk*)

(*Initialisation*)
nlast=500;ilast=jlast=500;(*total number of walkers*)
L=100;L1=-L;L2=L;
xc=IntegerPart[(L1+L2)/2];

For[j=1,j≤jlast,j++,
y[j]=xc;
For[i=1,i≤ilast,i++,
x[i]=xc;
];(*end of j*)
];(*end of j*)

avex[0]=Sum[x[j],{j,1,jlast}]/jlast;
avexsq[0]=Sum[x[j]^2,{j,1,jlast}]/jlast;

avey[0]=Sum[y[j],{j,1,jlast}]/jlast;
aveysq[0]=Sum[y[j]^2,{j,1,jlast}]/jlast;

tnx[0]=Table[{x[j],0},{j,1,jlast}];
tx[0]=Table[x[j],{j,1,jlast}];
tny[0]=Table[{y[j],0},{j,1,jlast}];
ty[0]=Table[y[j],{j,1,jlast}];

tnr[0]=Table[Sqrt[x[j]^2+y[j]^2],{j,1,jlast}];

tnr2[0]=Table[{x[j],y[j]},{j,1,jlast}];
(*End of Initialisation*)

For[n=1,n≤nlast,n++,

For[j=1,j≤jlast,j++,
ran=Random[Real];
If[ran≤0.25,x[j]=x[j]-1];(*end if *)
If[0.25<=ran≤0.5,x[j]=x[j]+1];(*end if *)
If[0.5<=ran≤0.75,y[j]=y[j]+1];(*end if *)
If[0.75<=ran≤1.0,y[j]=y[j]-1];(*end if *)
r[j]=Sqrt[x[j]^2+y[j]^2]
];(*end for j*)

tnr[n]=Table[r[j],{j,1,jlast}];
tnr2[n]=Table[{x[j],y[j]},{j,1,jlast}];

(*evaluate the average of r[j] and rsq[j] at fixed n*)
aver[n]=Sum[r[j],{j,1,jlast}]/jlast;
aversq[n]=Sum[r[j]^2,{j,1,jlast}]/jlast;
(*evaluate the average of r[j] and rsq[j] at fixed n*)
];(*end for n *)
tnaver=Table[{n,aver[n]},{n,0,nlast}];
tnaversq=Table[{n,aversq[n]},{n,0,nlast}];

Manipulate[
ListPlot[tnr2[n],AspectRatio->1,PlotRange->{{L1,L2},{L1,L2}}],
{n,0,nlast,1}]

lptnaversq=ListPlot[tnaversq,PlotLabel->"aversq, vs n"];
yth[n_]:=2*Dx*n;Dx=0.5;
tny=Table[{n,yth[n]},{n,0,nlast}];
lptny=ListPlot[tny,PlotLabel->"yth vs n",PlotStyle->Red];
Show[{lptny,lptnaversq}]

Section 7.4a: Diffusion

Consider a 3-D cubic lattice, with a lattice constant 1. 
Each lattice site (i,j,k) has 6 nearest neighbours. If probability a neighbour, say (i-1,j,k) is occupied by a walker at time n-1 is P(i-1,j,k, n-1), the probability this walker to hope to site (i,j,k) at time n is 1/6. Since there are 6 neighbours to the site (i,j,k), the probability that it is occupied at time n is simply
P(i,j,k,n) = (1/6)[P(i-1,j,k, n-1)+P(i+1,j,k, n-1)+P(i,j-1,k, n-1)+P(i,j+1,k, n-1)+P(i,j,k-1, n-1)+P(i,j,k+1, n-1)]. 
Now, consider the difference, P(i,j,k,n) -P(i,j,k,n-1) =  P(i,j,k,n) + (1/6)[-2P(i,j,k,n-1) - 2P(i,j,k,n-1) - 2P(i,j,k,n-1)].



Convince yourself that indeed Eq. (7.16) is the discretised difference equation for (assumed ((Δx)^3)=((Δz)^3)=((Δy)^3)).
, where .

Say there is a total N walkers, N = (∑_(all i,j,k))m(i,j,k,n), where m(i,j,k,n)=N·P(i,j,k,n) is the number of walkers at site (i,j,k) at time n. 
m(i,j,k,n)/(((Δx)^3))=N·P(i,j,k,n)/(((Δx)^3))≡ρ(i,j,k,n) is the walker density = number of walkers at site (i,j,k) per unit volume at time n. We see that ρ(i,j,k,n) ∝ P(i,j,k,n). The density,  ρ(i,j,k,n) obeys the same equation as for the probability, .

ρ tell us how the walkers are distributed in space at any given time n, and also how is this distribution evolves in time.

Diffusion equaiton in 1D

In 1-D, . ,  is
Given an initial distribution ρ(i,n=0) is known, we can evaluate  ρ(i,n) for all n in the future.  
For numerical programming, numerical stability requires that the time step and spatial step must obey . 
A time step that is larger than  has the meaning that information of the density distribution propagates faster than theoretically allowed by the interpretation of the equation . You can varify that if time step is choosen larger than allowed, your simulation indeed become not stable. 

One of the possible analytical solution to the diffusion equation is given by . Note that σ is time dependent. It means the size of the ‘drop of milk’ (as represented by a collection of walkers) gets spread out.

Sample code 7.4.1 The code below numerically solve the diffusion equation in 1-D by assuming that all walkers are initially concentrated at a point.

(* Sample code 7.4.1 *)
nlast=100;
size=10.0;(*in arbitrary unit*)
Deltax=0.1;(*in arbitrary unit*)
Dcoef=1.0;(*dimensionless*)
Deltat=Deltax^2/(2*Dcoef);
rho0=50;
iini=0;ilast=IntegerPart[size/Deltax];
ic=IntegerPart[(ilast+iini)/2];

(*Setting initial profile*)
For[i=iini-1,i≤ilast+1,i++,
rho[i,0]=0;
];(*end for i*)
rho[ic,0]=rho0;(*end for i*)
tnrho[0]=Table[{i*Deltax,rho[i,0]},{i,iini,ilast}];
(*End of setting initial profile*)


dmy=Dcoef*Deltat/Deltax^2;
For[n=0,n≤nlast,n++,
	For[i=iini,i≤ilast,i++,
		rho[i,n+1]=rho[i,n]+dmy*(rho[i+1,n]+rho[i-1,n]-2rho[i,n]);
	];(*end for i*)

(*boundary condition*)
	rho[iini-1,n+1]=0;rho[ilast+1,n+1]=0;
(*boundary condition*)
tnrho[n]=Table[{i*Deltax,rho[i,n]},{i,iini,ilast}];
];(*end for n*)
Print[
Manipulate[
ListPlot[tnrho[n],
Mesh->Automatic,Joined->True,
PlotRange->{{iini*Deltax,ilast*Deltax},{0,rho0}}]
,{n,0,nlast,1}](*end Manipulate*)
](*end print*)

However, the output of the Sample code 7.4.1 are ‘spiky’, as there are many values of rho that are zero. We can ‘smoothen’ out the output by  the following code (Sample code 7.4.2).

Sample code 7.4.2 The code is modifed from sample code 7.4.1 by averaing the results in code 7.4.1 over adjacent grid elements. This is done by “post-processing” the output rho[i,n]: Make the value of rho[i,n] to be defined as the average of that from the both sides of i.

(* Sample code 7.4.2 *)
nlast=100;
size=10.0;(*in arbitrary unit*)
Deltax=0.1;(*in arbitrary unit*)
Dcoef=1.0;(*dimensionless*)
Deltat=Deltax^2/(2*Dcoef);
rho0=50;
iini=0;ilast=IntegerPart[size/Deltax];
ic=IntegerPart[(ilast+iini)/2];

(*Setting initial profile*)
For[i=iini-1,i≤ilast+1,i++,
rho[i,0]=0;
];(*end for i*)
rho[ic,0]=rho0;(*end for i*)
tnrho[0]=Table[{i*Deltax,rho[i,0]},{i,iini,ilast}];
(*End of setting initial profile*)


dmy=Dcoef*Deltat/Deltax^2;
For[n=0,n≤nlast,n++,
	For[i=iini,i≤ilast,i++,
		rho[i,n+1]=rho[i,n]+dmy*(rho[i+1,n]+rho[i-1,n]-2rho[i,n]);
	];(*end for i*)

(*boundary condition*)
	rho[iini-1,n+1]=0;rho[ilast+1,n+1]=0;
(*boundary condition*)


(* average rho for adjacent spatial sites*)
	For[i=iini,i≤ilast,i++,
		rho[i,n+1]=(rho[i-1,n+1]+rho[i+1,n+1])/2;
	];(*end for i*)
(* average rho for adjacent spatial sites*)
tnrho[n]=Table[{i*Deltax,rho[i,n]},{i,iini,ilast}];
];(*end for n*)
Print[
Manipulate[
ListPlot[tnrho[n],
Mesh->Automatic,Joined->True,
PlotRange->{{iini*Deltax,ilast*Deltax},{0,rho0}}]
,{n,0,nlast,1}](*end Manipulate*)
](*end print*)

Section 7.4 b : Randam walks and Diffusion

Random walks and diffusion are closedly connected. We can generate the diffusive behavior of an initially concentrated walkers (particles) at certain region using random walk. By constructing the histogram of the number of walkers on the simulation sample at a fixed time n, we can obtain the distribution profile at a fixed time n. The histogram can then be compared with the density profile as obtained from the numerical distribution from diffusion equation.

Cream in your coffee: random walk approach

Let a lot of walkers initially concentrated at a poin, and then allowed to random walk. We wish to obtain the density profile of the walker at time t = n*Deltat. That is, we



The random walk results exhibit the Gaussian spreadinng of the particle distribution that we have come to expect.

Sample code 7.4.3 is just a simplified version of Sample code 7.4.2.  which simulates the spreading of an initial distribution of walkers concentrated in the center of the simulation box of size 2L. These walkers are then undergo random walks and spread out in the simulation box, and L1=-L<= x<= L2=L. L has to be an integer in our case here since the stepping size of the walker is ±1.
 
A histogram of totalbin bins is constructed to display the distribution of walkers in these bins at each stem step n. The width of the bins in the histogram is simply bw=(L2-L1)/totalbin. The position of the walkers are contained in the list tx[n]. BinCounts[tx[n],{L1,L2,bw}]tells Mathematica to output the numbers of walkers in each bin. If you Listplot BinCounts[tx[n],{L1,L2,bw}]you’ll get a dotted curve showing the number of walkers in the y-axis and index of the bin in the x-axis. The index of the bin range from 1 to totalbin.  In Sample code 7.4.3, we list plot the output of  BinCounts[tx[n],{L1,L2,bw}]/jlast instead of just  BinCounts[tx[n],{L1,L2,bw}].

(*Sample code 7.4.3*)

(*Initialisation*)
nlast=100;jlast=5000;(*total number of walkers.*)
L=50;L1=-L;L2=L; 
totalbin=100;(*number of bins*)
bw=(L2-L1)/totalbin;(*bin width*)
xc=IntegerPart[(L1+L2)/2];

For[j=1,j≤jlast,j++,
x[j]=xc;
];(*end of j*)
tx[0]=Table[x[j],{j,1,jlast}];
(*End of Initialisation*)

For[n=1,n≤nlast,n++,
For[j=1,j≤jlast,j++,
ran=Random[Real];
If[ran≤0.5,x[j]=x[j]-1
,(*else*)
x[j]=x[j]+1
](*end if *)
];(*end for j*)

(*tx[n]=Table[x[j],{j,1,jlast}];*)
tx[n]=Table[{(x[j-1]+x[j+1])/2},{j,2,jlast-1}];
];(*end for n *)
Manipulate[
ListPlot[N[BinCounts[tx[n],{L1,L2,bw}]/jlast],PlotRange->{0,1},Filling->Axis,PlotLabel->"(# walkers/total walker) in a bin vs bin index",AxesLabel->{"bin index","# of walkers/total walker"}]
,{n,0,nlast,1}]

For some choice of parameters (jlast, bw,L) the distribution curve could be spiky. In such a case  we may need to “smoothen” the histogram by taking replacing x[i] by the average evaluated using adjacent values of (x[i-1],x[i+1])/2, i.e., x[i] -> (x[i-1],x[i+1])/2. This was done in the line tx[n]=Table[{(x[j-1]+x[j+1])/2},{j,2,jlast-1}] in Sample code 7.4.3.

However, what we really want is a graph with : 
y - axis : (number of walkers/Total walkers) in each bin,
x - axis : the corresponding ranges ofxfor each bin (instead of the index of the bins).

Sample code 7.4 .4 below takes in the data tx[n] generated from Sample code 7.4 .3 where the set of the positions of all walkers are contained in tx[n]. 

Basically we simply wish to convert the bin index j to the corressponding ranges of position, xbin[j], so that we can get the desired curve as mentioned above. 

As a output of Sample code 7.4.4, the horizontal axis is expressed in values of xbin[j]. The vertical axis is the normalised number of walker in bin xbin *)

(* Sample code 7.4.4 *)
(** totalbin = Total number of bins you decide to use; bw is the width of the bin; both are defined in Sample code 7.4.3*)

xbin[j_]:=N[L1+j*bw];(*1 <= j ≤ totalbin)*)
NormalisedBC[j_,n_]:=BinCounts[tx[n],{L1,L2,bw}][[j]]/jlast;(* 1<=j≤totalbin *)

aveNormalisedBC[j_,n_]:=
(NormalisedBC[j-1,n]+NormalisedBC[j+1,n])/2.0;(*(NormalisedBC[j-1,n]+NormalisedBC[j+1,n]+NormalisedBC[j,n])/3.0;*)

(* 2<=j≤totalbin-1 *)

(*aveNormalisedBC[j_,n_]:=N[NormalisedBC[j,n]];*)

txbinnwinbinj[n_]:=Table[{xbin[j],aveNormalisedBC[j,n]},{j,2,totalbin-1}];

Manipulate[ListPlot[txbinnwinbinj[n],Joined->True,Filling->Axis,Mesh->All,PlotRange->{{L1,L2},{0,1}}],{n,0,nlast,1}]

As a homework assignment, you will be required to write a program to investigate how the width of the density distribution of the random walkers widens as a function of time σ = σ (t), and then show that σ (t) is proportional Sqrt[t], as expected from the diffusion equation. For this purpose, you may need to use the command FindFit. 
Take for example the data set at n = n0 ,txbinnwinbinj[100]. 
This is a set of paired data points of rho[j] vs xbin[j], where rho[j] = aveNormalisedBC[j] .We wish to fit this set of data  against a Gaussian  of the form rho(x)=Aexp[((-(((x-μ))^2))/(2(σ^2)))]. We want Mathematica to tell us what is the best-fit values of A and 
 σ. μ is the value of x where the curve is divided symmetrically. In our case μ is ideally μ=xc = (L2 - L1)/2. But due to the finite width of wd (due to the binning procedure), we got to correct μ to μ -> μ + wd/2. 

 Sample code 7.4.5   illustrates how we use FindFit to fit the the walker' s data against the Gaussian parameterss. We fit the  smoothened data for n = 100  generated by Sample code 7.4.4, txbinnwinbinj[100], against the function to fit,rho(x)=Aexp[((-(((x-μ+bw/2))^2))/(2(σ^2)))].

(* Sample code 7.4.5  *)
(*Use FindFit to find the parameters A (normalisation constant), and width (=sigma) of the curve*)

lptxbinnwinbinj=ListPlot[txbinnwinbinj[100],Mesh->All,PlotStyle->Blue]
sol=FindFit[txbinnwinbinj[100], A*Exp[-(x -(xc+bw/2))^2/(2wid^2)], 
    {A, wid}, x];
Print[sol];

(* Plot function *)
A1=A/.sol;wid1->wid/.sol;
y[x1_]:=A1*Exp[-(x1 -(xc+bw/2))^2/(2wid1^2)]/.{A1->A/.sol,wid1->wid/.sol};
Print["y[x]= ",y[x]];
txy=Table[{xbin[ind],y[xbin[ind]]},{ind,1,totalbin}];
lptxy=ListPlot[txy,Joined->True,Mesh->All,PlotStyle->Red]
Show[{lptxbinnwinbinj,lptxy}]
(*ListPlot[{Table[{k,y[x]},{x,L1,L2,bw}],txbinnwinbinj[100]},Joined->True,Mesh->Automatic,PlotStyle->{Red,Black},PlotRange->All]*)