Optimizing Code N-Body Simulation

Question

below is my code for a gravitational n-body problem. I know Mathematica has a n-body simulation function however I am programming this to learn Mathematica more and programming numerical techniques. Some of my simulations take a while to finish (small dt,large number of particles) and I am trying to speed up the code. The code I provided randomly chooses a position and velocity of each mass then runs under a softened Newtonian potential using the kick-drift-kick method (numerical technique I am currently learning). I have not scaled anything which is why G=1. The code seems to work how it should and you should be able to copy paste the code. I have provided some test cases but the comment marks need to be changed.

I believe the getAcc() can get cleaned up but I am not sure how. I thought of using a double sum with Sum[Sum[]] but this allows for particle particle interaction when i=j. I do not know how to exclude when i=j in the double sum. I believe this is where the majority of the slow down occurs due to the double for loop. Any suggestions are greatly appreciated.

Clear[Mass, Pos, Vel, Acc, NPart];
(*Gravitational Constant, Time Step, Softening, and Number of \
Particles *)
G = 1;(*6.674*10^-11;*)
dt = .1;
Soft = 100;
NPart = 2;

(* Testing Code: Same mass particles at rest falling into each other. \
Change Soft to 10.
Mass = {1,1};
Pos = {{1,0},{-1,0}};
Vel = {{0,0},{0,0}}; *)

(* Testing Code: Large-small mass system. Change Soft to 100 *)

Mass = {1, 1000};
Pos = {{10, 0}, {0, 0}};
Vel = {{.000001, .25}, {0, 0}}; 


(* Randomly generated masses, initial position, and initial \
velocities. Change Soft to 10.
Mass = RandomInteger[{1,1},NPart];
Pos = RandomReal[{-10,10},{NPart,2}];
Vel = RandomReal[{-.2,.2},{NPart,2}];*)

(* Sets all initial velocities to zero for at rest scenario. *)
\
(*Vel= ConstantArray[0,{Length[Mass],2}];*)


(* Finds the total acceleration on each particle for a given set of \
positions, mass, and the softening parameter*)

getAcc[{m_, r_, s_}] :=
 ((* Clears the Acc list and PEtotal value. 
  This is done each time this function is called because the \
acceleration and potential energy is position dependent. 
  Everytime the position is updated the accelerations and potential \
energy needs to be updated. *)
  Acc = {};
  PEtotal = 0; 
  
  (*Calculates acceleration for ith particle *)
  
  For[i = 1, i <= Length[r], i++,
   
   (*Clears {ax,ay} for ith particle*)
   ax = 0; 
   ay = 0;
   
   (* Finds jth particle's acceleration contribution*)
   
   For[j = 1, j <= Length[r], j++,
    
    (*Does not allow for particle to interact with itself*)
    
    If[i != j,
      
      (*Calculates relative distance between particle i and j*)
      
      dx = 0;
      dy = 0;
      dx = r[[j, 1]] - r[[i, 1]];
      dy = r[[j, 2]] - r[[i, 2]];
      
      (* Acceleration denominator*)
      
      denAcc = Sqrt[( dx^2 + dy^2 + s^2)];
      
      (* Adds acceleration contribution to previous acceleration *)
  
          ax = ax + G*(dx/denAcc^3)*m[[j]];
      ay = ay + G*(dy/denAcc^3)*m[[j]];
      
      (* Potential denominator *)
      
      denPot = ( dx^2 + dy^2 + s^2)^-.5;
      
      (* Adds potential energy contribution to previous potential \
energy *)
      PEtotal = PEtotal - G *m[[i]]*m[[j]]*denPot;
      
      , Null
      ];                            
    ];
   
   (* Turns the acceleration components into a vector and adds vector \
to Acc list *)
   a = {ax, ay};
   AppendTo[Acc, a];
   
   ])

(* Finds total kinetic energy and potential energy of the current \
configuration of the system *)

getEnergy[{m_, r_, v_}] :=
 (KEtotal = 0;
  KineticList =  v v;
  For[i = 1, i <= Length[m], i++,
   KEtotal += (1/2) Mass[[i]]*
      Sqrt[KineticList[[i, 1]] + KineticList[[i, 2]]];
   ];
  AppendTo[PotentialData, PEtotal];
  AppendTo[KineticData, KEtotal];
  )

(* Simulation Code *)
Clear[PosSim, VelSim, PosData]
PosSim = Pos; 
VelSim = Vel;
PosData = {};
KineticData = {};
PotentialData = {};

getAcc[{Mass, PosSim, Soft}];
getEnergy[{Mass, PosSim, VelSim}];

For[k = 0, k <= 5000*dt, k += dt,
 VelSim = VelSim + Acc*dt/2;
 PosSim = PosSim + VelSim*dt;
 getAcc[{ Mass, PosSim, Soft}];
 VelSim = VelSim + Acc*dt/2;
 getEnergy[{Mass, PosSim, VelSim}];
 AppendTo[PosData, PosSim];
 ]

ListPlot[Table[PosData[[All, i]], {i, Length[Mass]}]]
ListPlot[{PotentialData, KineticData, KineticData + PotentialData}]
Mean[KineticData + PotentialData]

Answer 1

Here is a simulation to get you going using the gravitation law without changes. I am using a very simple integration scheme, that fails if the bodies come too close. You may improve this if you like.

n = 3; nt = 2000; SeedRandom[1234];
pos = RandomReal[{-1, 1}, {n, 2}];
vel = 1/2 RandomReal[{-1, 1}, {n, 2}];
tvel = Table[0, n, 2];
mass = ConstantArray[1, n];
g = 1;
dt = 0.001;

step[pos0_, vel0_, dt_] := 
 Module[{pos = pos0, vel = vel0, last = 1, sum = 0},
  Do[
   If[i != j,
    If[ i == last
     , sum += 
       mass[[j]] (pos[[j]] - pos[[i]])/Norm[pos[[j]] - pos[[i]]]^3;
     , vel[[last]] += g  sum  dt; 
     pos[[last]] += vel[[last]] dt; sum = 0; last = i;
     ]
    ];
   If[i == n + 1, Break[]]
   , {i, n + 1}, {j, n}];
  {pos, vel}
  ]


res = Table[({pos, vel} = step[pos, vel, dt])[[1]], {nt}];
mima = {MinMax[res[[All, All, 1]]], MinMax[res[[All, All, 2]]]};
Animate[
 Graphics[{Red, PointSize[0.02], Point[res[[i]]]}, PlotRange -> mima, 
  Axes -> True]
 , {i, 1, nt, 1}
 ]