# 3D Random Walk with Periodic Boundary Conditions

I am working on a 3D random walk with periodic boundary conditions and have written a program that will do this for me but it is extremely slow. Does anyone have any suggestions on how I could speed this up?

Timing[
n = 1.;
T = 10000.;
boundary = 5;
initial = boundary;
step = 0.1;
RandomWalk[x_] :=
Accumulate[
Join[{Table[RandomReal[{-initial, initial}], {3}]},
RandomVariate[NormalDistribution[0, step], {x, 3}]]];
i = 1;
(Label[begin];
p[i] = RandomWalk[T + 1];
If[i == n, Goto[end]];
i = i + 1;
Goto[begin];
Label[end];);
t = 0;
(Label[begin1];
i = 1;
Label[begin2];
positionslist = p[i];
If[p[i][[t + 1, 1]] > boundary,
p[i] = Join[Take[positionslist, t + 1],
Accumulate[
Join[{{-boundary, p[i][[t + 1, 2]], p[i][[t + 1, 3]]}},
RandomVariate[
NormalDistribution[0, step], {T - t + 1, 3}]]]]];
If[p[i][[t + 1, 1]] < -boundary,
p[i] = Join[Take[positionslist, t + 1],
Accumulate[
Join[{{boundary, p[i][[t + 1, 2]], p[i][[t + 1, 3]]}},
RandomVariate[
NormalDistribution[0, step], {T - t + 1, 3}]]]]];
If[p[i][[t + 1, 2]] > boundary,
p[i] = Join[Take[positionslist, t + 1],
Accumulate[
Join[{{p[i][[t + 1, 1]], -boundary, p[i][[t + 1, 3]]}},
RandomVariate[
NormalDistribution[0, step], {T - t + 1, 3}]]]]];
If[p[i][[t + 1, 2]] < -boundary,
p[i] = Join[Take[positionslist, t + 1],
Accumulate[
Join[{{p[i][[t + 1, 1]], boundary, p[i][[t + 1, 3]]}},
RandomVariate[
NormalDistribution[0, step], {T - t + 1, 3}]]]]];
If[p[i][[t + 1, 3]] > boundary,
p[i] = Join[Take[positionslist, t + 1],
Accumulate[
Join[{{p[i][[t + 1, 1]], p[i][[t + 1, 2]], -boundary}},
RandomVariate[
NormalDistribution[0, step], {T - t + 1, 3}]]]]];
If[p[i][[t + 1, 3]] < -boundary,
p[i] = Join[Take[positionslist, t + 1],
Accumulate[
Join[{{p[i][[t + 1, 1]], p[i][[t + 1, 2]], boundary}},
RandomVariate[
NormalDistribution[0, step], {T - t + 1, 3}]]]]];
If[i == n, Goto[end2]];
i = i + 1;
Goto[begin2];
Label[end2];
If[t == T, Goto[end1]];
t = t + 1;
Goto[begin1];
Label[end1];);]

{9.44905, Null}


I would like to do this simulation for multiple particles (at least 50) for about 100000 time steps and not have it take all day.

-
Welcome on Mathematica.SE. Have you seen this question here which seems very close to what you ask? –  halirutan Oct 16 '12 at 0:30
"periodic boundary conditions" - you are aware of Mod[]? –  Guess who it is. Oct 16 '12 at 0:43
For a start, (Label[begin]; p[i] = RandomWalk[T + 1]; If[i == n, Goto[end]]; i = i + 1; Goto[begin]; Label[end];); can surely be done more efficiently as a Table command: Table[p[i]=RandomWalk[T+1],{i,n}]. –  Verbeia Oct 16 '12 at 0:46
@halirutan No, I have not seen the question you are referencing. Thanks. I will see if it can help. –  Kane Oct 16 '12 at 0:54
@J.M. No, I am not familiar with Mod[]. How can it help me? –  Kane Oct 16 '12 at 1:03

As I mentioned in the comments, using Mod[] is one good way to enforce your periodic boundary conditions: just generate the random walk as usual, and then apply Mod[] to bring back inside the sections that are outside your box.

Here's an example of what I'm describing:

n = 5*10^3; (* number of steps *)
s = 20; (* cube edge length *)
h = 1/10; (* step bound *)

BlockRandom[SeedRandom[6081, Method -> "MersenneTwister"]; (* for reproducibility *)
walk = Accumulate[Join[{RandomReal[{-s, s}/2, 3]},
RandomVariate[NormalDistribution[0, h], {n, 3}]]]];

periodizedWalk = Mod[walk, s, -s/2];
splitPeriodizedWalk = Split[periodizedWalk, EuclideanDistance[#1, #2] < s/2 &];

With[{cube = First[PolyhedronData["Cube"]]},
{Graphics3D[{{Opacity[2/3], Scale[cube, s]}, Line[walk]}, Boxed -> False],
Graphics3D[{{Opacity[2/3], Scale[cube, s]},
Riffle[Line /@ splitPeriodizedWalk, {Blue, Red}, {1, -2, 2}]},
Boxed -> False]}
// GraphicsRow]


Note how the sections that go outside the cube are "folded in".

-

Rapid calculations are afforded by Accumulate to generate the walk and Mod to implement the periodicity. Scaling the entire thing to the unit cube simplifies the code a little. Furthermore, don't generate normally distributed displacements: uniform displacements will do when the increments are small.

With these efficiencies we may generate the coordinates of a random walk of $n$ steps of magnitude at most $e$ extremely quickly and with remarkable brevity of syntax:

points[n_Integer, e_] := Mod[Accumulate /@ RandomReal[{-e, e}, {3, n}] // Transpose, 1];


As an example, let's time the generation of independent random walks of $10^5$ steps for each of $50$ particles. The last one is retained and plotted:

With[{e = 0.01, n = 10^5},
Print[Timing[Table[p = points[n, e], {i, 1, 50}];]];
Graphics3D[{PointSize[0.1/n^(1/3)], Point[p]}]]


{0.359, Null}

It looks like this is around 13,000 times faster than the original. :-)

Starting with this as a base, one can proceed (incrementally) to develop it more fully if desired, such as plotting the values in a prettier fashion, or generating normal rather than uniform increments (in case the time steps are intended to be relatively long), or automatically rescaling the original dimensions to the unit cube, or providing for specified starting points, or inducing correlations among the particles, etc.

-
The use of FractionalPart[] is equivalent to my previous suggestion of using Mod[], since Mod[x, 1] == FractionalPart[x] for x nonnegative... –  Guess who it is. Oct 16 '12 at 17:11
@J.M. Thanks--I have been looking into this and noticed that FractionalPart actually returns negative values for negative numbers! Thus, it is not equivalent to Mod and Mod is the correct choice, not FractionalPart. I will edit the answer. Interestingly, the timing does not noticeably change. –  whuber Oct 16 '12 at 17:21

This is an alternative approach:

boundary = 5;
step = .1;
AbsoluteTiming[
rwtbl = NestList[
Mod[
# + RandomVariate[NormalDistribution[0, step], 3],
boundary] &,
{0, 0, 0},
10000
];
]

Graphics3D@Line[
rwtbl
]


(it's also 2400 times faster than your code on my machine).

-

You are creating the random walk, then imposing the boundary condition and recreating the remaining random steps and accumulating them all over again. This means you are creating and accumulating the steps repeatedly every time you hit the barrier. Surely a simpler way to do this would be to generate the steps and impose the boundary condition at the initial accumulation, using FoldList?

Table[FoldList[ If[Abs[First[#1]+First[#2]]>boundary,
Join[{boundary * Sign[First[#1]+First[#2]]},Rest[#1+#2]] , #1 +#2]&,
someStartingVector, RandomVariate[NormalDistribution[0, step], {x, 3}]], {nbodies}]


You might be able simplify even further using Clip, along the lines of:

Table[FoldList[{Clip[#1[[1]]+#2[[1]],{-boundary,boundary}],
#1[[2]]+#2[[2]],#1[[3]]+#2[[3]]}&,
someStartingVector, RandomVariate[NormalDistribution[0, step], {x, 3}]],{nbodies}]


I have not tested this code but I do think something along the lines of this will work for you, with the added bonus of eliminating all those Labels and Gotos, which are not considered good Mathematica practice.

-