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 s 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*)

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 a 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 for a 1-D random walk. It simulates the spreading of an initial distribution of walkers concentrated in the center of the simulation box of size L. This walkers are then undergo random walks and spread out in the simulation box, and 0<= x<= L. L has to be an integer in our case here since the stepping size of the walker is ±1.
 
A histogram of N bins is constructed to display the number of walkers in each stem step n. The width of the bins in the histogram is simply dx=L/N. 

To build the histrogram, we use the command: Histogram[data,dx,”Probability”], where data is a list containing the information of the position of the walkers, data={2,1,5,8,9,0,1,3,...}, dx the width of bin.

The following code is modified from the Sample code 7.5. Note that we will use the new command Export and Import for data manipulation. Export saves the numerical data in the hard disk in form of *.dat. To use the same set of data, simply Import the *.dat file.

(* Sample code 7.4.3 *)
(*Begin of initialisation*)
SetDirectory["C:\\Users\\Yoon\\Documents\\teaching\\ZCE111\\1011SEM2\\notes\\lecture7"];
x0=0.0;nlast=200;(*nlast = total time = 100 by default*)
jlast=6000;(*jlast = total number of walker = 5000 by default*)
(*Size of simulation box*)
L=100;
dx=1;(*width of bin*)
(*End of initialisation*)
Export["nlast.dat",nlast];
Export["L.dat",L];
Export["dx.dat",dx];


(*Setting the initial conditions for the walker at n=0*)
For[j=1,j≤jlast,j++,
(*The walkers are initially concentrated in the central point of simulation box*)
x[j]=IntegerPart[L/2];
];(*end for j*)
tx[0]=Table[x[j],{j,1,jlast}];

(*End of setting the initial conditions for the walker at n=0*)


For[n=1,n≤nlast,n++,
tx[n]={0};
For[j=1,j≤jlast,j++,
q=Random[Integer];
If[q==0,x[j]=x[j]+1,x[j]=x[j]-1];(*end if*)(* Shift the walker randomly*)
tx[n]=Append[tx[n],x[j]];
];(*end for j*)
BC[n]=BinCounts[tx[n],{0,L,dx}];(*Count the number of walkers in each bin in the 
histrogram, where the range of x is in [0,L], and the width of the bin is dx=1*)
(*BC[n] is a list containing L elements*)
BC[n]=BC[n]/jlast;(*normalise the number of walker to make it a fractional valaue*)
];(*end of n*)

For[n=1,n≤nlast,n++,
tnBC[n]={0};
For[k=1,k≤L/dx,k++,
nBC=N[BC[n][[k]]];(*BC[n][[k]] = # of walkers in bin indexed by k, where k = 1,2,...,L*)
tnBC[n]=Append[tnBC[n],nBC];(*tnBC[n] is the histrogram to be ploted*)
];(*end for k*)
name="tnBC"<>ToString[n]<>".dat";(*export as "tnBC1.dat", "tnBC2.dat"m ...*)
Export[name,tnBC[n]];
];(*end for n*)

Manipulate[
ListPlot[tnBC[n],PlotRange->{0,0.3},(*Mesh->Automatic,Joined->True,*)Filling->Axis],{n,1,nlast,1}]

Now we “smoothen” the histogram:

Sample code 7.4.3a

(* Sample code 7.4.3a *)
SetDirectory["C:\\Users\\Yoon\\Documents\\teaching\\ZCE111\\1011SEM2\\notes\\lecture7"];
L=Import["L.dat","List"][[1]];
nlast=Import["nlast.dat","List"][[1]];
dx=Import["dx.dat","List"][[1]];

For[n=1,n≤nlast,n++,
name2="tnBC"<>ToString[n]<>".dat";
datt[n]=Import[name2,"List"];
];(*end for n*)

For[n=1,n≤nlast,n++,
tnBCsm[n]={0.0};
For[k=1,k≤IntegerPart[(L/dx)]-1,k++,
(*nBCsm=N[(BC[n][[k]]+BC[n][[k+1]])/2];*)
nBCsm=N[(datt[n][[k]]+datt[n][[k+1]])/2];
tnBCsm[n]=Append[tnBCsm[n],nBCsm];(*tnBC[n] is the histrogram to be ploted*)
(*Print[BCtx[n]][[k]]*)
];(*end for k*)
name2="tnBCsm"<>ToString[n]<>".dat";
Export[name2,tnBCsm[n]]; (*export as "tnBCsm1.dat", "tnBCsm2.dat"m ...*)
];(*end for n*)

Manipulate[
ListPlot[tnBCsm[n],PlotRange->{0,0.08}(*,Mesh->Automatic,Joined->True*),Filling->Axis],{n,1,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 match the σ (t) to show that it is proportional Sqrt[t], as expected from the diffusion equation. For this purpose, you may need to use the command FindFit. The following sample code illustrate how we use FindFit to fit the the walker' s data (imported from tnBC*.dat) against the Gaussian parameterss. 

In the following code, we fit the smoothened data as generated by Sample code 7.4.3a against the function to fit,
y(x) = Aexp[((-(((x-μ))^2))/(2(σ^2)))], where μ=L/2 is at the center of the simulation box.

(*Sample code 7.4.4 on Findfit*)

SetDirectory["C:\\Users\\Yoon\\Documents\\teaching\\ZCE111\\1011SEM2\\notes\\lecture7"];

L=Import["L.dat","List"][[1]];
nlast=Import["nlast.dat","List"][[1]];
dx=Import["dx.dat","List"][[1]];

For[n=1,n≤nlast,n++,
name3="tnBCsm"<>ToString[n]<>".dat";
datt[n]=Import[name3,"List"];
];(*end for n*)

sol=FindFit[datt[nlast], A*Exp[-(x - 0.5*L)^2/(2wid^2)], 
    {A, wid}, x];
Print[sol];

y[x1_]:=A1*Exp[-(x1 - 0.5*L)^2/(2wid1^2)]/.{A1->A/.sol,wid1->wid/.sol};
Print[y[x]];
ListPlot[{Table[{k,y[k]},{k,0,L/dx}],datt[nlast]},Joined->True,Mesh->Automatic,PlotStyle->{Red,Black},PlotRange->All]