N Body Problem in 2 Dimensions using indicies

This notebook demonstrates the use of indicies for the N body problem in 2 dimensions.   Indicies allow the problem to be expressed generally for any number of bodies.The default indicies from Mathematica are not ideal for this purpose. This subscript package by Ted Ersek may be downloaded from MathSource.The package makes it simple to generalize a system of diffrential equations, assign initial values  and clear them. 

William C. Collins
wccollins@starband.net

Packages

Needs["Utilities`SubscriptSymbols`"];

Initial Conditions: M[n], clearM[n] & cv0[n], clearic[n]

Define masses:    M[n] &  ClearM[n]

M[n_] := Table[(Subscript[m, j]), {j, 1, n}];
clearM[n_] := Table[(Subscript[m, j]) =. , {j, 1, n}]

First assign  masses for explicit n : 
            Evaluate[M[8]]  = {1, 2, 3, 4, 5, 6, 7, 
          8} ;

M[8]

Evaluate[M[8]] = Table[i, {i, 1, 8}]

M[8]

To clear use ClearM[
        8]

ClearM[8]

M[8]

Initial cooridinates and velocities: cv0[n]

cv0[n_] := Table[{(Subscript[x0, i]), (Subscript[y0, i]), (Subscript[vx0, i]), (Subscript[vy0, i])}, {i, 1, n}]

Coordinates and velocities at time t  from NDSolve interpolating function: s=-1 reverses velocities at time t: cv[n,t,data,s]

cv[n_, t_, data_, s_:  1] := 
      Table[{(Subscript[x0, i]) = (Subscript[x, i])[t] /. data[[1]], 
          (Subscript[y0, i]) = (Subscript[y, i])[t] /. data[[1]], 
          (Subscript[vx0, i]) = s*(Subscript[x, i])'[t] /. data[[1]], 
          (Subscript[vy0, i]) = s*(Subscript[y, i])'[t] /. data[[1]]}, {i, 1, n}]

Clear Initial Conditions: clearcv0[n]

clearcv0[n_] := 
      Table[{(Subscript[x0, i]) =. , (Subscript[y0, i]) =. , (Subscript[vx0, i]) =. , (Subscript[vy0, i]) =. }, {i, 1, 
          n}]

NDSolve for NBody problem: NBody2D[s,n,G,ms,ag,pg,wp]

n-> number of bodies; G->Gravitational constant; ms-> MaxSteps; ag->AccuracyGoal; pg-> PrecisionGoal; wp->PrecisionGoal

Options[NBody2D] = {n -> 8, G -> 1, ms -> Automatic, 
        ag -> Automatic, pg -> Automatic, 
        wp -> MachinePrecision}; 
    NBody2Dopt[s_, n_, G_, ms_, ag_, pg_, wp_] := 
      NDSolve[Flatten[{Table[{(Subscript[x, i])''[
                    t] == (∑_(j = 1)^(i - 1))((-G (Subscript[m, j]) ((Subscript[x, i])[t] - (Subscript[x, j])[t]))/(((((((Subscript[x, i])[t] - (Subscript[x, j])[t]))^2) + ((((Subscript[y, i])[t] - (Subscript[y, j])[t]))^2)))^(3/2))) + (∑_(j = i + 1)^n)((-G (Subscript[m, j]) ((Subscript[x, i])[t] - (Subscript[x, j])[t]))/(((((((Subscript[x, i])[t] - (Subscript[x, j])[t]))^2) + ((((Subscript[y, i])[t] - (Subscript[y, j])[t]))^2)))^(3/2))), (Subscript[y, i])''[
                    t] == (∑_(j = 1)^(i - 1))((-G (Subscript[m, j]) ((Subscript[y, i])[t] - (Subscript[y, j])[t]))/(((((((Subscript[x, i])[t] - (Subscript[x, j])[t]))^2) + ((((Subscript[y, i])[t] - (Subscript[y, j])[t]))^2)))^(3/2))) + (∑_(j = i + 1)^n)((- G (Subscript[m, j]) ((Subscript[y, i])[t] - (Subscript[y, j])[t]))/(((((((Subscript[x, i])[t] - (Subscript[x, j])[t]))^2) + ((((Subscript[y, i])[t] - (Subscript[y, j])[t]))^2)))^(3/2)))}, {i, 1, n}], 
            Flatten[Table[{(Subscript[x, i])[0] == (Subscript[x0, i]), (Subscript[x, i])'[0] == (Subscript[vx0, i]), 
                  (Subscript[y, i])[0] == (Subscript[y0, i]), (Subscript[y, i])'[0] == (Subscript[vy0, i])}, {i, 1, n}]]}], 
        Flatten[Table[{(Subscript[x, i]), (Subscript[y, i])}, {i, 1, n}]], {t, 0, s}, 
        MaxSteps -> ms, AccuracyGoal -> ag, 
        PrecisionGoal -> pg, WorkingPrecision -> wp]; 
    NBody2D[s_, opts___] := 
      NBody2Dopt[
        s, n /. {opts} /. Options[NBody2D], G /. {opts} /. 
          Options[NBody2D], ms /. {opts} /. 
          Options[NBody2D], ag /. {opts} /. 
          Options[NBody2D], pg /. {opts} /. 
          Options[NBody2D], wp /. {opts} /. Options[NBody2D]];

NBody Plot routine:NBodyPlot2D[t0,ts,data,n,df,mb,pp]

Options: n-> number of bodies, df->$DisplayFunction, mb->MaxBend,  pp->PlotPoints

Options[NBodyPlot2D] = {n -> 8, 
          df -> $DisplayFunction, mb -> 10, 
          pp -> 100};
NBodyPlot2D[t0_, ts_, 
        data_, opts___] := 
      NBodyPlot2Dopt[t0, ts, 
        data, n /. {opts} /. Options[NBodyPlot2D], df /. {opts} /. 
          Options[NBodyPlot2D], mb /. {opts} /. 
          Options[NBodyPlot2D], pp /. {opts} /. 
          Options[NBodyPlot2D]]
NBodyPlot2Dopt[t0_, ts_, data_, n_, df_, mb_, pp_]:=Show[ParametricPlot[
  Evaluate[Table[{{(Subscript[x, i])[t],  (Subscript[y, i])[t]} /.data},{i, 1, n}]],
   {t,  t0,  ts}, AspectRatio  ->  Automatic,Frame -> True,MaxBend -> mb,PlotPoints -> pp,PlotStyle -> 
                  Table[{RGBColor[Abs[Sin[i*Pi/4]], 
                        Abs[Sin[(i + 1)*Pi/4]], 
                        Abs[Sin[(i + 2)*Pi/4]]]}, {i, 1, n}],DisplayFunction -> Identity],	\(Table[
              Graphics[
		{RGBColor[Abs[Sin[i*Pi/4]], 
                    Abs[Sin[\((i + 1)\)*Pi/4]], Abs[Sin[\((i + 2)\)*Pi/4]]], 
                  PointSize[0.04], 
                  Point[Flatten[{x\_i[ts] /. data, y\_i[ts] /. data}]]}], {i, 
                1, n}]\),	DisplayFunction -> df];

Examples: 8 Body circular distributions:

Circular symmetric positions and tangeantal velocities: using r[i].v[i]=0

Symmetric positions:

cirpos[d, n] := 
      Table[{d Sin[k*2  Pi/n], d Cos[k*2  Pi/n]}, {k, 1, n}]

Symmetric tangeantal velocities:

cirvel[n_, v_] := 
      Table[{-v Cos[k*2  Pi/n], v Sin[k*2  Pi/n]}, {k, 1, 
          n}]

Circular symmetric initial conditions:  cir[n,d,v]

cir[n_, d_, v_] := 
      Table[{d Sin[k*2  Pi/n], d Cos[k*2  Pi/n], -v Cos[k*2  Pi/n], 
          v Sin[k*2  Pi/n]}, {k, 1, n}]

CircleA: using default precision

clearM[8]

Evaluate[M[8]] = Table[i/i, {i, 1, 8}]

M[8]

Assign initial conditions for a symmetric circular distribution: cir[n, r, v]

Evaluate[clearcv0[8]]

Evaluate[cv0[8]] = cir[8, 8, 2/5]

run forward to t = 200

circleA = NBody2D[200];

NBodyPlot2D[0, 100, circleA]

NBodyPlot2D[0, 200, circleA]

CircleB: using extended precision

clearM[8]

Evaluate[M[8]] = Table[i/i, {i, 1, 8}]

M[8]

Assign initial conditions

Evaluate[clearcv0[8]]

Evaluate[cv0[8]] = cir[8, 8, 3/5]

Solve & Plot using extended AccuracyGoal and WorkingPrecision

Options[NBody2D]

circleB = NBody2D[800, ag -> 25, wp -> 30];

NBodyPlot2D[0, 200, circleB]

NBodyPlot2D[200, 400, circleB]

NBodyPlot2D[400, 600, circleB]

NBodyPlot2D[600, 800, circleB]

Examples: 8 Body circular distributions with one central body

CentralA: using default precision

clearM[9]

Evaluate[M[9]] = {1, 1, 1, 1, 1, 1, 1, 1, 10}

M[9]

Assign initial conditions for a symmetric circular distribution: cir[n, r, v] with a central heavy body

Evaluate[clearcv0[9]]

cv0[9]

Evaluate[cv0[8]] = cir[8, 12, 3/5]

{(Subscript[x0, 9]), (Subscript[y0, 9]), (Subscript[vx0, 9]), (Subscript[vy0, 9])} = {0, 0, 0, 0}

centralA = NBody2D[200, n -> 9];

NBodyPlot2D[0, 50, centralA, n -> 9]

NBodyPlot2D[0, 200, centralA]