Basically I am generating a grid in space (so a vector where v[[1]]=position 1 in x) and a random number of particles in space (so a vector where p[[1]]= position of particle 1 on space).

Afterwards, I want to see in which spaces of the grid the particles are and do somethings.

The Do is the specific part of the code I want to optimize.

rho = Compile[{{xp, _Real, 1}, {xg, _Real, 
1}, {carga, _Real}, {np, _Real}, {ng, _Real}, {dx, _Real}},

n = ng + 1;
  ρ = Table[carga*np/(dx*ng), n];
  ρ[[ng + 1]] = 0;
  Do[If[Abs[xp[[i]] - xg[[j]]] <= 
     dx, ρ[[j]] = ρ[[j]] - 
      carga*(1 - (Abs[xp[[i]] - xg[[j]]])/dx)/dx; r += 1, 0] , {i, 1, 
    np}, {j, 1, ng + 1}];
  ρ[[1]] = ρ[[1]] + ρ[[ng + 1]];]

Basically I am worried that I am thinking too much in terms of other languages and maybe mathematica has a better way to do this instead of using

If[Abs[xp[[i]] - xg[[j]]] <= 

That manually makes a lot of operations

Here is a part of the code I am using

eps = 8.85*10^(-12);
q = 80;
m = 50;
L = 2;
nparticles = 2000;
ngrid = 200;
vth = 2;
u = 1;

dx = L/ngrid;

x = RandomVariate[UniformDistribution[{0, L}], nparticles];
r = 0;

xgrid = Range[0, L, dx] // N;
rho[x, xgrid, q, nparticles, ngrid, dx]; // AbsoluteTiming
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Commented Jan 22, 2022 at 15:37
  • 1
    $\begingroup$ Minor tips: You should localize your variables inside Compile or there's really no point in compiling the code. If you need to have several return types (which can't be handled strictly within the compiled environment), then set global variables at the very end, after the loop is done. The r += 1 looks like a no-op to me, unless it's counting steps for some hidden purpose. $\endgroup$
    – Michael E2
    Commented Jan 22, 2022 at 17:27
  • $\begingroup$ Please, try to be more precise. Daniel Huber already wasted his time just because you did not describe your problem correctly. With no word do you mention that you are working in 1D. Also your post is lacking example data on which we can run any tests. And your Compile code does not make any sense (e.g., are ng and np supposed to be integers rather then doubles?). And btw., you should probably scope some variables (e.g., ρ,n,r) and initialize r. $\endgroup$ Commented Jan 22, 2022 at 17:27
  • $\begingroup$ The r+=1 is a counter that I used to see how many time the if proved true. As long as I initialize r=0 before I call the rho in the code I can always know if it did the apropiate operations. It will eventually go away but for now it is useful to make sure the code is doing what it should. yeah np and ng are suposed to integers/reals. And my problem with Daniel Huber's answer is that he doesn't create a vector that stores the coordinates. As you can see in the code I need the vectors in order to calculate ρ . And I will add an example of the code. $\endgroup$
    – divica
    Commented Jan 22, 2022 at 19:03

3 Answers 3


This is how I would write your function:

rho2 = Compile[{{p, _Real, 1}, {g, _Real, 
     1}, {c, _Real}, {np, _Integer}, {ng, _Integer}, {dx, _Real}},
   Block[{\[Rho], init, pdiv, pmod, j, a, dxinv},
    init = N[np]/N[ng];
    dxinv = 1./dx;
    a = c dxinv;
    \[Rho] = Table[0., {ng + 1}];
     pdiv = Compile`GetElement[p, i] dxinv;
     j = Ceiling[pdiv];
     pmod = pdiv - j;
     \[Rho][[j]] -= Abs[pmod];
     \[Rho][[j + 1]] -= Abs[pmod + 1];
     {i, 1, np}];
    \[Rho][[1]] += \[Rho][[ng + 1]];
    Do[\[Rho][[j]] = a (Compile`GetElement[\[Rho], j] + init), {j, 1, ng}];
    \[Rho][[-1]] *= a;
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"

Note that I scoped \[Rho] within the function and handed it over as output of the function.

Here is how it compares to the orginal function rho with a substantially finer grid and more particles:

q = 80;
m = 50;
L = 2;
nparticles = 10000;
ngrid = 1000;

dx = L/ngrid;
p = RandomVariate[UniformDistribution[{0, L}], nparticles];
g = Subdivide[0., L, ngrid];

r = 0;
t1 = First@RepeatedTiming[
    rho[p, g, q, nparticles, ngrid, dx];
t2 = First@RepeatedTiming[
    \[Rho]2 = rho2[p, g, q, nparticles, ngrid, dx];
"Relative error" -> Max[Abs[(1 - \[Rho]2/\[Rho])]]
"Speedup" -> t1/t2

"Relative error" -> 4.39937*10^-11

"Speedup" -> 12270.6

Note that the relative error is quite high (I would have expected something of the order of 10^-15). And it increases with ngrid. I think this is due to a serious precision issue in your code: You initialize \[Rho] with values that are orders of magnitude greater than the ones you subtract from it. This is why I initialize \[Rho] by 0., do the subtraction loop, and add the initial values afterwards.

The key idea here is what Daniel Huber suggested: Instead of comparing each particle to all grid points, we can just compute the integers belonging to the two grid points next to the particle. This reduces the original computational complexity from $O(\mathrm{ngrid} \times \mathrm{nparticles})$ to $O(\mathrm{nparticles} + \mathrm{ngrid})$. So the speedup also grows with ngrid.

The remaining changes are just elementary refactoring:

  • Instead of dividing by dx so often, I compute its reciprocal only once and multiply with it. (Floating point division is several times slower than multiplication.)

  • Also, I postpone the multiplication with a until after the contibutions of all particles have been accumulated into \[Rho].

  • $\begingroup$ It is a beautiful piece of code. I have 4 funtions simillar to this in my code and this is one of the most complex and is orders of magnitude faster than everything else. Can you explain in detail the algortihm you used to do this so that I can implement it in my other functions? Also why does you rho go from -1001 to 1001 intead of only 1 to 1001? $\endgroup$
    – divica
    Commented Jan 23, 2022 at 21:02

Instead of creating a space grid, we simply imply integer coordinates from 0 to n. E.g. {0,0,0} indicates the cell at the origin. {1,0,0} the next cell along the x axis etc. Then we can get the cell number for (nearly) "free":

To create some random points:

n = 100;
npts = 10;
pts = RandomReal[{0, n}, {npts, 3}]

(* {{81.7389, 11.142, 78.9526}, {18.7803, 24.1361, 6.57388}, {54.2247, 23.1155, 39.6006}, {70.0474, 21.1826, 74.8657}, {42.2851, 24.7495, 97.7172}, {82.5163, 92.5275, 57.8056}, {29.287, 20.8051, 58.0474}, {12.8821, 30.6427, 71.2012}, {39.0582, 81.9967, 32.5351}, {59.326, 51.8774, 16.9013}} *)

Now, to get the cells in wich the random points are, we only need to calculate the floor of the coordinates:


(* {{81, 11, 78}, {18, 24, 6}, {54, 23, 39}, {70, 21, 74}, {42, 24, 97}, {82, 92, 57}, {29, 20, 58}, {12, 30, 71}, {39, 81, 32}, {59, 51, 16}} *)

  • $\begingroup$ Pretty smart code, but I don't see how I can adapt it to my code. My space is only 1D, but the real problem is that if you look here ``` Do[If[Abs[xp[[i]] - xg[[j]]] <= dx, ρ[[j]] = ρ[[j]] - carga*(1 - (Abs[xp[[i]] - xg[[j]]])/dx)/dx; r += 1, 0] , {i, 1, np}, {j, 1, ng + 1}]; ``` afte knowing where the particle is on the grid I use that position in the grid vector to calculate [rho] in the grid $\endgroup$
    – divica
    Commented Jan 22, 2022 at 16:18
  • 2
    $\begingroup$ For 1 dim, simple write pts = RandomReal[{0, n}, npts]; and then use Round again. The result will be the indices of the cells. $\endgroup$ Commented Jan 23, 2022 at 9:15

A program written using Compile will only run faster if your code does not violate it's limitations. For an explanation of what those limitations are read this.

  • $\begingroup$ Using the compile my program runs 90% faster. my doubt is if there is some other method to make it faster. for example some mathematica way to make the loop Do faster. For example, whenever I try to use parallelDo the code becomes much slower. Or if there is some fuction to compare the elements of xp and xg without using the if for each element... $\endgroup$
    – divica
    Commented Jan 22, 2022 at 21:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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