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]