# Optimizing Code N-Body Simulation

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]


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