6
$\begingroup$

I am new in Mathematica, it is my first program. Many codes of diffusion-limited aggregation are available but, they are too difficult to understand. I basically want to create a simulation of the random walk demonstrating the Brownian motion in a visible grid. It is a regular lattice. It starts with a single seed particle at the origin of a lattice. A second particle (random walker), is added at some random site at a long distance from the origin. This particle walks randomly until it reaches a site adjacent to the seed and becomes part of the growing cluster. The particle jumps from the current site to one of their nearest lattice sites at each step until it hits and sticks to the cluster of fixed particles. A third, fourth and so on particles are introduced in the same way and they also join the growing cluster. If the particle reaches the boundaries of the lattice in its random walk, it is killed, and another introduced. The procedure is repeated until a cluster of sufficiently large size is formed. The aggregation process is similar to the DLA model. The aggregation starts with a seed at the site which is nearest to the centre of the square. A particle is added at a random site at large distance from the seed. The particle chooses randomly a neighbour and walks on the geometric graph. When the particle reaches a site that is adjacent to the seed or the cluster, the particle becomes a part of the cluster. I want the end result of this type given in link -> aggregation model

Towards this end, I have created a grid and two particles. But, I am not able to move the particle. Below is my code:

Graphics[{PointSize[Large], Point[{{25/2, 25/2}, 
{RandomInteger [25], RandomInteger [25]}}]},  
GridLines -> {Range[0, 25], Range[0, 25]}, 
PlotRange -> {{0, 25}, {0, 25}}, Axes -> True]
$\endgroup$
  • 1
    $\begingroup$ So, what would be the kinematics that dominate the movement of the particles. You described how this should look like but that is not enough for an implementation. $\endgroup$ – Henrik Schumacher Aug 29 '18 at 7:14
  • 1
    $\begingroup$ Indeed, it's very hard to know what kind of "random walk" - and "sticking" is expected. Random walk on a grid? Or is this supposed to be a physical simulation of some sort, and if so, what sort? Depending on such details this problem can be anything from trivial to essentially impossible to solve. $\endgroup$ – kirma Aug 29 '18 at 7:25
  • $\begingroup$ @BiSarfraz That's a pretty good description. Why don't you add it directly to your question (see the "edit"-link below your post)? $\endgroup$ – Henrik Schumacher Aug 29 '18 at 8:09
  • $\begingroup$ If the walk always arrives to the same endpoint, then it is not a random walk in the usual sense. Some clarification about this would be needed. As for DLA, there's an implementation in this answer: mathematica.stackexchange.com/a/20907/12 $\endgroup$ – Szabolcs Aug 29 '18 at 8:26
8
$\begingroup$

So, here is my reformulations of C.E.'s code.

What has changed:

  • I use reservoirs rand and randpts of large random arrays; generating such arrays is much less expensive than generating every random number on its own.

  • I use RegionMember and Nearest for the collision tests.

Implementation

getCluster[startcluster_, nparticles_, L_] := Module[{
    getStartpoint, getDirection,
    directions, cluster, nf, box, inboxQ, nrand, rand, nrandpts, 
    randpts, iparticle, irand, irandpts, particle, inbox
    },

   getStartpoint[] := Block[{p},
     If[irandpts < nrandpts,
      irandpts++,
      irandpts = 1;
      randpts = RandomInteger[{-L, L}, {nrandpts, 2}];
      ];
     p = randpts[[irandpts]];
     While[
      Length[nf[p, {1, 0}]] > 0,
      If[irandpts < nrandpts,
       irandpts++,
       irandpts = 1;
       randpts = RandomInteger[{-L, L}, {nrandpts, 2}];
       ];
      p = randpts[[irandpts]];
      ];
     p
     ];

   getDirection[] := (
     If[irand < nrand,
      irand++,
      irand = 1;
      rand = RandomInteger[{1, 4}, nrand];
      ];
     directions[[rand[[irand]]]]
     );

   directions = Developer`ToPackedArray[{{0, 1}, {1, 0}, {0, -1}, {-1, 0}}];

   (* set up a sufficiently large array to store the cluster in order to avoid the costly Append.*)
   cluster = Join[startcluster, ConstantArray[{0, 0}, {nparticles}]];
   nf = Nearest[cluster -> Automatic];

   (* define the box and a collision checker *)

   box = Rectangle[{-L, -L}, {L, L}];
   inboxQ = RegionMember[box];

   (* counters for the reservoirs for random numbers and points *)

   irand = nrand = 1000000;
   irandpts = nrandpts = 10000;

   iparticle = Length[startcluster];

   (* share some progress info with the user *)       
   PrintTemporary[
    Dynamic[
     Grid@
      Transpose[{{"iparticle", "irand", "irandpts"}, {iparticle, 
         irand, irandpts}}]
     ]
    ];
   While[iparticle < Length[cluster],
    Check[particle = getStartpoint[];, Print["!"]];
    inbox = True;
    While[inbox,
     Check[particle += getDirection[];, Print["?"]];
     inbox = inboxQ[particle];
     If[inbox,
      If[Length[nf[particle, {1, 1}]] > 0,
       Check[cluster[[iparticle]] = particle;, Print["."]];
       nf = Nearest[cluster -> Automatic];
       iparticle++;
       inbox = False
       ]
      ]
     ]
    ];
   Association[
    "Cluster" -> cluster,
    "Box" -> box,
    "N" -> nparticles
    ]
   ];

And here a function to plot the particles with color according to their age:

showCluster[a_Association] := With[{L = a[["Box", 2, 2]]},
   ArrayPlot[
    Transpose@SparseArray[
      (a[["Cluster"]] + (L + 1)) -> 
       Rescale[Range[Length[a[["Cluster"]]]]],
      {2 L + 1, 2 L + 1},
      2.
      ],
    ColorFunction -> "DeepSeaColors",
    PlotRange -> {0, 1},
    ClippingStyle -> White,
    DataReversed -> {True, False}
    ]
   ];

Usage example

A test run with 2000 particles, a single condensation core at {0,0}, and a box of edgelength 101 = 2 * 50 + 1 (yes, this takes a bit):

data = getCluster[{{0, 0}}, 10000, 100]; // AbsoluteTiming // First

184.823

And here a plot of the result:

showCluster[data]

enter image description here

You can continue the simulation with

data2 = getCluster[data[["Cluster"]], 10000, 100]; // AbsoluteTiming // First
showCluster[data2]

28.3286

enter image description here

Surprisingly complex patterns!

$\endgroup$
  • $\begingroup$ It's very nice, +1 $\endgroup$ – C. E. Aug 29 '18 at 16:39
  • $\begingroup$ @C.E. Thanks! I really appreciate that. $\endgroup$ – Henrik Schumacher Aug 29 '18 at 16:42
  • $\begingroup$ @HenrikSchumacher Thank you so very much for your kind help and time. $\endgroup$ – BiSarfraz Aug 29 '18 at 17:30
  • $\begingroup$ You're welcome! $\endgroup$ – Henrik Schumacher Aug 29 '18 at 17:30
  • $\begingroup$ @BiSarfraz Oh, sorry. I made a mistake when copying the code to the site. If fixed it. Please try again. $\endgroup$ – Henrik Schumacher Aug 30 '18 at 5:57
3
$\begingroup$

This is like what you described, except that I also made sure that the particle didn't wander outside of the region $-100 \leq x, y \leq 100$ to ensure that convergence wouldn't take too long.

particles = {{0, 0}};
nf = Nearest[particles -> Automatic];
nextStep := RandomChoice[{{0, 1}, {1, 0}, {0, -1}, {-1, 0}}];
nParticles = 100;
l = 100;

inBounds[{x_, y_}, l_] := -l < x < l && -l < y < l

Do[
  particle = RandomInteger[{-l, l}, 2];
  While[
   Length[nf[particle, {1, 1}]] == 0,
   particle = RandomInteger[{-l + 1, l - 1}, 2];
   ];
  step = nextStep;
  While[
   Length[nf[particle, {1, 1}]] == 0,
   particle += step;
   step = nextStep;
   While[
    ! inBounds[particle + step, l],
    step = nextStep;
    ]
   ];
  AppendTo[particles, particle];
  nf = Nearest[particles -> Automatic],
  {nParticles}
  ];

Graphics[{
  Point[particles]
  }]

Mathematica graphics

$\endgroup$
  • $\begingroup$ I'd suggest to do the collision test with Nearest: Define a nearest function nf = Nearest[particles -> Automatic]; and do the actual collision test with Length[nf[particle, {1, 1}]] == 0. $\endgroup$ – Henrik Schumacher Aug 29 '18 at 8:31
  • $\begingroup$ I meant to use Length[nf[particle, {1, 1}]] == 0 as replacement of ! MemberQ[particles, particle + step]. Moreover, nf = Nearest[particles -> Automatic]; has to run whenever particles changes, so right in the beginning and right after AppendTo[particles, particle]. $\endgroup$ – Henrik Schumacher Aug 29 '18 at 14:32
  • $\begingroup$ @HenrikSchumacher Thank you for pointing this out, I incorporated it into the answer. $\endgroup$ – C. E. Aug 29 '18 at 16:35
  • $\begingroup$ @BiSarfraz I made an update to the answer to address performance issues. $\endgroup$ – C. E. Aug 29 '18 at 16:40
  • $\begingroup$ @C.E. Thank you so much, I really appreciate your help and concern. $\endgroup$ – BiSarfraz Aug 29 '18 at 17:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.