N-Body Simulator!

In this code, we calculate the dynamics of a system of N-point masses (currently the code is set up to handle two point masses) experiencing a mutual attraction via a gravitational 1/r potential, but confined to move in 2D.

We discretize time in order to approximately solve the differential equation, separating it into equal small intervals:

TSTEP=1/10;NSTEP=10000;

We specify the initial position (x,y) coordinates of the two masses at time t=0:

POS={{{5.,5.},{-5.,-5.}}};

We specify the initial velocity of the two masses

VEL={{{.1,-.1},{-.1,.1}}}; (* Made sum zero, so CM is not moving if masses are equal *)

We now calculate the new positions of the two objects, using the initial velocity to get the position at time t=TSTEP:

AppendTo[POS,POS[[1]]+VEL[[1]] * TSTEP]

Here are the masses, to be used in the code:

MASS={1,1};

The next part is the meat of the code, which calculates first the force on each object (the force here is an array of two vectors, the total force on object 1, followed by the total force on object 2

For[i=1,i<NSTEP,i++,
(*  At each value of i, we are at some time t  *)
(*  LASTPOS is an array containing only the last two elements of the trajectory array, the positions at time t, and time t-STEP  *)
LASTPOS=Take[POS,-2];
(*  The FORCE is obtained by taking the gradient of the gravitational potential (G_Newton has been set to one) - it sums over all the pairs of masses, although in this case, there is only one pair - check that you agree with the physics!  I.e. F = -Grad V = r/|r|^3 (F and r are vectors) *)
FORCE=Table[Sum[If[j==k,{0,0},MASS[[j]]*MASS[[k]]/Dot[LASTPOS[[2]][[j]]-LASTPOS[[2]][[k]],LASTPOS[[2]][[j]]-LASTPOS[[2]][[k]]]^(3/2)*(LASTPOS[[2]][[j]]-LASTPOS[[2]][[k]])],{j,1,2}],{k,1,2}];
(*  Now we apply Newton's 2nd law to get the acceleration of each mass  *)
ACC=Table[FORCE[[j]]/MASS[[j]],{j,1,2}];
(*  Now we use the fact that ACC = (POS(t+TSTEP)-2 POS(t)+POS(t-TSTEP))/TSTEP^2  is a good approximation to the acceleration in the small TSTEP limit, and solve for POS(t+TSTEP), which we call NEWPOS *)
NEWPOS=ACC*TSTEP^2+2LASTPOS[[2]]-LASTPOS[[1]];
(*  Finally we add the new position to the array, and iterate to the next time-step  *)
AppendTo[POS,NEWPOS];
]

Here we make a big table of plots, where each plot has one mass shown red, and another mass shown blue, and each plot corresponds to one value of the time

plots=Table[ListPlot[{Style[POS[[i]][[1]],Red],Style[POS[[i]][[2]],Blue]},PlotRange->{{-10,10},{-10,10}},AspectRatio->1,PlotStyle->PointSize[Medium]],{i,1,Dimensions[POS][[1]],50}];

Here, we create a new table of plots that shows not just the positions at one time, but also the positions for all times before.

frames=FoldList[Show,plots];

ListAnimate[frames]

Here is the same animated plots, but not showing the previous positions of the masses:

Animate[ListPlot[{Style[POS[[i]][[1]],Red],Style[POS[[i]][[2]],Blue]},PlotRange->{{-10,10},{-10,10}},AspectRatio->1,PlotStyle->PointSize[Medium]],{i,1,Dimensions[POS][[1]],1}]

.