4
$\begingroup$

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]
$\endgroup$

1 Answer 1

1
$\begingroup$

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

enter image description here

$\endgroup$

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.