Speeding up a loop?

Questions
,
1

Given that Mathematica has a wide range of functions I would like to know if there is a better way to make this loop using Mathematica.

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

Basically I am comparing the position i of vector xp with the position j of vector xg. If the condition is verified then I use the iterator j to calculate the function [Rho] of j.
My question is if Mathematica doesn’t have a better way to do this operation that would take less time to run. Like some built in function or smarter method to do operations between vectors?

Here is a snippet of my code if anyone need better context to understand.

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

dx = L/ngrid;

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]];]

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

xgrid = Range[0, L, dx] // N;
rho[x, xgrid, q, nparticles, ngrid, dx]; // AbsoluteTiming
r
2

Just a clarification… Is it expected that each position of \[Rho] can get updated multiple times during each inner loop?

3

Yes. ρ is calculated at each point xg[[j]]. The number of times ρ is updated depends on the number of particles that has a position between xg[[j]] + or - dx.

4

Diogo,

What you seem to be doing is finding where each point lies on a grid. There is a much simpler way to do this rather than to compare the value to every point on the grid. I would first compute the integer part of the x values to find where on the grid you start:

pos = Floor[x/dx] + 1;

I did not use Ceil because I think you always round up even if the error is zero (but if I’m wrong you can change this)

Next get the fractional part:

frac = FractionalPart[x/dx];

or get both at the same time with MixedFractionParts[].

With this information you can set your rho values with

 MapThread[rhoFormula,{pos,frac}]

in which you index into rho using pos and do your equivalent computation of the Do loop using frac.This works because you do all this looping to find the position in the grid just below and just above your x values. pos is the position in the grid below your value and pos+1 is the position in the grid above your value. frac and (1 - frac) gives you the expression Abs[xp[[i]] - xg[[j]]])/dx which is used in your calculations.

I hope this helps.

Regards,

Neil

5

I am confused.
so you are saying my code should be something like this?

Do[pos=Floor[x[[i]]/dx];rho[[pos]]=q/dx * Abs[FractionalPart[x[[i]]/dx], {i,nparticles}];
 Do[pos=Ceiling[x[[i]]/dx];\[rho][[pos]]=q/dx * Abs[FractionalPart[x[[i]]/dx], {i,nparticles}];

Because the particles influence the point of the grid before and after them.

6

Related question on MSE.

7

Diogo,

I am confused. so you are saying my code should be something like this?

Almost. However, Do[] loops are generally not recommended in Mathematica. Also, the FractionalPart is always positive so the Abs is not necessary. I would use Map or MapThread instead of the loop.

For example, your posted code takes 0.2 seconds to run compiled. This code gets the same result but takes only 0.007 seconds on my M1 machine without compilation:

AbsoluteTiming[rr = Table[q*nparticles/(dx*ngrid), n];
 rr[[ngrid + 1]] = 0;
 MapThread[(rr[[#1]] = rr[[#1]] - q*(1 - (#2))/dx) &, {pos, frac}];
 MapThread[(rr[[#1]] = rr[[#1]] - q*((#2))/dx) &, {pos + 1, frac}];
 rr[[1]] = rr[[1]] + rr[[ngrid + 1]];]

Regards,

Neil

8

Another thing you might try is Nearest.

Nearest[deltaSequence -> {"Index", "Distance"}, particleValues, {All, deltaParam}]

I renamed some variables because I was having trouble keeping track, but I think it would be this with your original variables:

Nearest[xg -> {"Index", "Distance"}, xp, {All, dx}]

I included my renamed version so that if I messed it up, you can hopefully correct it.

Anyway, what this will give you is something like the following:

{{{6, 0.0122529}, {7, 0.0877471}}, {{18, 0.00370914}, {17, 
   0.0962909}}, {{18, 0.0247811}, {17, 0.0752189}}, {{4, 
   0.0364263}, {5, 0.0635737}}, {{8, 0.0122352}, {7, 0.0877648}}, <...etc...>}

You should be able to write a function that uses this data to do your update to your rho table. I haven’t done a performance comparison.

9

Okay, I went ahead and looked at performance a bit. First, here’s an alternate to your rho function:

rho2[xp_, xg_, carga_, np_, ng_, dx_] :=
 Module[
  {rhoTable = Table[carga*np/(dx*ng), 1 + ng]},
  rhoTable[[1 + ng]] = 0;
  With[
   {neighbors = Nearest[xg -> {"Index", "Distance"}, xp, {All, dx}]},
   Map[(rhoTable[[#[[1]]]] -= (carga*(1 - (#[[2]])/dx)/dx)) &, 
    neighbors, {-2}];
   rhoTable[[1]] += rhoTable[[ng + 1]];
   {Length[Flatten@neighbors]/2, rhoTable}]]

I used local variables rather than update global ones. It looked like you were interested in number of updates, so I returned both that and the rho table, although I don’t know if number of updates would ever be anything other than 2*nparticles.
Here’s the timing I got on my machine for your version:

(result1 = rho[x, xgrid, q, nparticles, ngrid, dx]; {r, \[Rho]}) // AbsoluteTiming // Short
(*{0.11498, {4000, {-12152.9, 8155.76, -4859.24, <<195>>, -11347.9, 468.88, -60255.7}}}*)

And here’s the timing with the version using Nearest:

(result2 = rho2[x, xgrid, q, nparticles, ngrid, dx]) // AbsoluteTiming // Short
(*{0.02478, {4000, {-12152.9, 8155.76, -4859.24, <<195>>, -11347.9, 468.88, -60255.7}}}*)

And a quick check:

result1[[2]] == result2[[2]]
(*True*)