3D Random Walk with Periodic Boundary Conditions

Question

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.

Answer 1

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.

Answer 2

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".

Answer 3

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).

Answer 4

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.