Simulations of Brownian motion

PHYS-436 STATISTICAL PHYSICS IV
Vivishek Sudhir, Tobias Kippenberg (LPQM)

(* generates n Brownian motion paths in the interval [0,1], with npoints time points *)
BrownianWalks[npoints_,n_]:=Map[
Accumulate,Join[
Transpose[{Table[0,{n}]}],
RandomReal[NormalDistribution[0,Sqrt[1.0/(npoints-1)]],{n,npoints-1}]
,2]
];

(* calculates ensemble average of function F applied to Brownian walk in [0,1] with npoints, using nensemble realizations *)
BrownianEnsembleAverage[F_,npoints_,nensemble_]:=Module[{walks,avg},
	walks=BrownianWalks[npoints,nensemble];
	avg=(Plus@@(F/@walks))/nensemble;
	Return[avg]
];

(* generates an Ornstein/Uhlenbeck trajectory in the interval [0,1] with npoints time points *)
(* Ito equation is: dx = -γ x dt + α dw, where w is a Wiener process so that dw = b (Sqrt[dt]) where b ∈ N(0,1)  *)
OUWalk[γ_,α_,npoints_]:=Module[{b,x,dt,i},
	x=ConstantArray[0,npoints];
	dt=(((1/2))^(npoints-2));
	b=Table[RandomReal[NormalDistribution[0,1]],{npoints}];
	(* apply update *)
	(* this is the brute-force approach, for OU process, there exists an exact update equation which is quicker *)
	For[i=2,i≤Length[x],i++,
		x[[i]]=x[[i-1]]-γ dt x[[i-1]]+α b[[i]] (Sqrt[dt]);];
	Return[x];
];

(* generates n OU trajectories *)
OUWalks[γ_,α_,npoints_,n_]:= Module[{walks},
	DistributeDefinitions[OUWalk]; (* distribute function definition to all processor cores *)
	walks = ParallelTable[OUWalk[γ,α,npoints],{n}]; (* evaluate in parallel to speed-up *)
	Return[walks]
];

(* calculates ensemble average of F over OU trajectories *)
OUEnsembleAverage[F_,γ_,α_,npoints_,nensemble_]:=Module[{walks,avg},
	walks=OUWalks[γ,α,npoints,nensemble];
	avg=(Plus@@(F/@walks))/nensemble;
	Return[avg]
];

(* settings for plotting *)
SetOptions[{ListLinePlot,ListPlot,Histogram},
Frame->True,FrameStyle->Directive[Gray,11],
LabelStyle->Directive[Black,FontFamily->"Helvetica"]];

Some sample Brownian paths

ListLinePlot[BrownianWalks[500,10],
	FrameLabel->{"time steps","x(t)"}]

Variance of Brownian walk does satisfy (Sqrt[t]) law

f[x_]:=x^2;
ListLinePlot[Sqrt/@BrownianEnsembleAverage[f,1000,1000],PlotStyle->Blue,
	FrameLabel->{"time steps","(Sqrt[<(x^2)>])"}]

Probability distribution of the Brownian walk

walks=BrownianWalks[10,500];
Animate[
	Histogram[
		Select[Flatten[walks[[1;;n]]],(#≠0)&],
		{0.15},"Probability",
		ChartStyle->Pink,FrameLabel->{"x","P(x)"},PlotLabel->"PDF of x(t)",PlotRange->{{-3,3},{0,0.15}}],
	{n,1,500,1},AnimationRate->30,AnimationRepetitions->1
]

Divergence in the x-correlation function

f[x_,n_]:=x[[2]]x[[2+n]];(* starting from 2nd element since first is always zero *)
ListLinePlot[
	Table[{n,BrownianEnsembleAverage[f[#,n]&,1000,100]},{n,0,990}],FrameLabel->{"n","<x(0)x(n)>"}
]

Velocity of the Brownian process

npoints=1000;
v[x_]:=Differences[x]/(npoints-1); (* velocity of Brownian walk *)
vcorr[x_,n_]:=v[x][[2]]*v[x][[2+n]];
GraphicsRow[
{ListLinePlot[
	v/@BrownianWalks[npoints,1],
	FrameLabel->{"time steps","v(t)"},PlotLabel->"Sample velocity trajectory"
],
ListLinePlot[
	Table[{n,BrownianEnsembleAverage[vcorr[#,n]&,npoints,1000]},{n,0,50}],
	PlotRange->Full,PlotStyle->Blue,
	FrameLabel->{"n","<v(0)v(n)>"},PlotLabel->"Decay of velocity correlations"
]}
]

General OU process with drift and damping

γ=10;α=1;
npoints=1000;
ListLinePlot[OUWalks[γ,α,npoints,5],
	FrameLabel->{"time steps","x(t)"},PlotLabel->"Sample OU trajectories starting at x=0"]

Brownian walk in a double-well potential

(* double well *)
V[x_,a_,h_]:=h ((((x+a)(x-a)))^2); (* zeros at x={a,-a}, and V(x=0)=(ha^4) *)
F[x_,a_,h_]:=-D[V[y,a,h],y]/.y->x;

(* Smoluchowski equation and its boundary/initial conditions *)
smolPDE[F_,Dc_]:=D[P[x,t],t]==Dc D[P[x,t],{x,2}]-D[F P[x,t],x]; (* works with arbitrary force field F *)
bcs[xmin_,xmax_]:={ (* BC for closed system - probability flux at end points is zero *)
Derivative[1,0][P][xmin,t]==0,
Derivative[1,0][P][xmax,t]==0
};
ic[x0_,dx0_]:={(* IC - probability is localized at some point x=x0, with width dx0 *)
P[x,0]==Exp[((-1)/2)(((((x-x0)/dx0)))^2)]/(Sqrt[2π (dx0^2)])
};  

(* specific example *)
a=2;h=1;
sol=NDSolve[{
	smolPDE[F[x,a,h],3],
	bcs[-5,5],
	ic[-a,0.1]
},P[x,t],{x,-5,5},{t,0,50}]//Flatten;
Animate[
	ListLinePlot[{
		Table[{x,V[x,a,h]},{x,-5,5,0.05}],
		Table[{x,(P[x,t]/.sol)/.t->T},{x,-5,5,0.02}]},
		PlotRange->{{-5,5},{0,4}},PlotStyle->{Red,Blue},
		FrameLabel->{"x","P(x,t)"}],
{T,0,50,0.2},AnimationRate->10,AnimationRepetitions->1]