16

6

I am trying to speed up this code for many particles to take a random walk. I'm not sure why it is so slow for such a simple task.

I got a few hints from colleagues to reduce the precision of the calculations, compile some of the functions, parallelize the code, or implement C code. I tried to implement each of these, but I most likely don't know what I'm doing because I couldn't get it to work. The particles can and should be simulated in parallel to take advantage of the two cores.

I have a bunch (10,000) particles that each take 1000 fixed steps in random directions. I have a bunch of lists with the length of number of particles, and I operate on those lists. This is faster than doing a nested loop, where I first loop over every particle and then loop over each time step. Two random angles are generated for each particle at each time step. Each particle's position is then incremented. I also want to know how many times each particle has crossed a spherical shell. The history of each particle is important: eventually, I will distinguish particles by the number of times they cross the shell so that they take bigger steps the more they have crossed.

I tried timing each step, and it seems like the trig functions slow things down. Is there any way to speed this up? Is there a better way to write this code to make it much faster? Eventually, I would like to add some more complications to this code and run it for many more time steps and particles. So, any general speed tips would help.

For comparison, the absolute timing of the Do loop is 3 seconds on my laptop with Core i5 M 540 @ 2.53 GHz, with 4 GB RAM.

```
numparticles = 10^4;
numsteps = 10^3;
radius = 1.;(*radius of shell*)
particlesx = Table[1.001, {numparticles}];(*x-coordinate for each particle*)
particlesy = Table[0., {numparticles}];(*y-coordinate for each particle*)
particlesz = Table[0., {numparticles}];
numcrossings = Table[0., {numparticles}];(*counts the number of crossings for each particle*)
Do[
\[Theta]rand = RandomReal[{0., Pi}, numparticles];(*random polar angle for each particle*)
\[Phi]rand = RandomReal[{0, 2.*Pi}, numparticles];(*random azimuthal angle for each particle*)
rold = particlesx^2 + particlesy^2 + particlesz^2;(*original distance from origin for each particle*)
particlesx = particlesx + 0.01*Sin[\[Theta]rand]*Cos[\[Phi]rand];(*update each particle's position*)
particlesy = particlesy + 0.01*Sin[\[Theta]rand]*Sin[\[Phi]rand];
particlesz = particlesz + 0.01*Cos[\[Theta]rand];
rnew = particlesx^2 + particlesy^2 + particlesz^2;(*new distance from origin*)
isitoutside = (rold - radius)*(rnew - radius);(*yields negative number if the particle crossed the sphere, else positive*)
switchsides = Map[If[# < 0, 1, 0] &, isitoutside];(*if the particle crossed, it increments numcrossings*)
numcrossings = numcrossings + switchsides;
, {step, 1, numsteps}]; // AbsoluteTiming
```

Not an answer to the question, but perhaps it's good to mention that a simple and fast way to generate an

`n`

-dimensional (or if you wish:`n`

-particle) random walk with`m`

steps is`Accumulate@RandomReal[{-1,1}, {m, n}]`

. – Szabolcs – 2012-05-12T10:48:59.403