The Many Body Problem

AP1601
Columbia University

Introduction

The equations of motion for two bodies interacting with a central force can be solved exactly. However, when more than two bodies are interacting, solutions are possible only in spe cial cases or with simplifying assumptions. When we assumed that the mass of the sun was much larger than the planets, we were making one of these simplifying assumptions.

In this notebook, we will explore the solution of several objects interacting through gravitational forces.

Representing Many Objects

How should we represent many objects?

I find it most convenient to represent the objects as arrays indexed by an integer. Arrays can be randomly accessed by using a integer within the range defined for that array.

a[1] = 3; a[4] = 2; a[7][t_] = t^2;

?a

Plot[a[7][t],{t,0,5}]

bList = Table[b[i][t_] = t^i, {i, 1, 10}]

?b

Plot[Evaluate[Table[b[i][t],{i,1,10}]],{t,0,2},PlotStyle->Table[Hue[(i/10)],{i,1,10}]]

Equations of Motion

The Definitions

Imagine a collection of objects of equal mass.

nBodies = 5;

mass = {10.} ~ Join~ Table[1., {i,2,nBodies}]

p = Table[{x[i][t], y[i][t], z[i][t], vx[i][t], vy[i][t],vz[i][t]}, {i, 1, nBodies}] // Flatten;

"p" is a list of all of the unknown dynamical functions of time.

p

The force (per unit mass) on object "i" from object "j" is found from...

distance[i_,j_] := Sqrt[(x[i][t] - x[j][t])^2 + (y[i][t] - y[j][t])^2 + (z[i][t] - z[j][t])^2]

forcex[i_,j_] := (x[j][t] - x[i][t])mass[[j]]/distance[i,j]^3;
forcey[i_,j_] := (y[j][t] - y[i][t])mass[[j]]/distance[i,j]^3;
forcez[i_,j_] := (z[j][t] - z[i][t])mass[[j]]/distance[i,j]^3;

We must sum the forces from all objects, except we must exclude self-forces.

?Complement

dp = Table[{vx[i][t], vy[i][t], vz[i][t], Plus @@ (forcex[i,#]& /@ Complement[Range[nBodies], {i}]),
Plus @@ (forcey[i,#]& /@ Complement[Range[nBodies], {i}]),
Plus @@ (forcez[i,#]& /@ Complement[Range[nBodies], {i}])}, {i, 1, nBodies}] // Flatten;

A particular initial condition...

p0=Flatten[Table[{x[i][0]==RandomReal[{-1.`,1.`}],y[i][0]==RandomReal[{-1.`,1.`}],z[i][0]==RandomReal[{-1.`,1.`}],vx[i][0]==RandomReal[{-1.`,1.`}],vy[i][0]==RandomReal[{-1.`,1.`}],vz[i][0]==RandomReal[{-1.`,1.`}]},{i,1,nBodies}]];

eqs = Thread[D[p,t] == dp]~ Join ~ p0;

"eqs" are the equations of motion!

Example Solution

Print["nBodies = ",nBodies];(sol1 = NDSolve[eqs , p, {t, 0, 5}, MaxSteps -> 80000] // Flatten;) // Timing

How long does Mathematica require to solve these equations as the number of planetary objects change? Try this with nBodies = 5, 10, 15. (Be warned: the solution time will change greatly depending upon the initial condition! If two objects are too close, then the solution requires very many time steps to maintain accuracy.)

Visualizing the Solution

To view the solution, we need 3D graphics.

threeRange = {{-3.0, 3.0}, {-3.0, 3.0}, {-3.0, 3.0}};

objects[tt_] := Graphics3D[{PointSize[0.04],Red, Point[{x[1][t], y[1][t], z[1][t]}], PointSize[0.02],Blue} ~ Join ~Table[Point[{x[i][t], y[i][t], z[i][t]}] , {i, 2, nBodies}] /. sol1 /. t -> tt, PlotRange->threeRange]

Show[objects[1.`]]

Animation

Animate[objects[t], {t, 0.0, 5.0, 0.1},AnimationRunning-> False]

We made a oversight. In order to follow the moving objects, we must change our frame of reference so that we are moving and located at the "center of gravity" of the many bodies.

Averages

We can find the average radius and the average position.

avgRadius[tt_] = (Plus @@Table[Sqrt[x[i][t]^2 + y[i][t]^2 + z[i][t]^2] /. sol1/. t -> tt, {i, 1, nBodies}] )/nBodies;

The Apply[...] function is used in its "infix" form...

Plus @@ {a1, a2, a3, a4}

avgRadius[0.2]

Plot[avgRadius[t],{t,0,5},PlotLabel->"Average Radius"]

The average position is a vector...

avgPosition[tt_] = (Plus @@ Table[{x[i][t], y[i][t],z[i][t]}/. sol1 /. t -> tt, {i, 1, nBodies}])/nBodies;

Plus @@ {{x1,y1,z1},{x2,y2,z2},{x3,y3,z3}}

avgPosition[0.2]

sp1=ListPointPlot3D[Table[avgPosition[t],{t,0,5,0.1`}],PlotStyle->{Blue,PointSize[0.02`]},PlotLabel->"Average Position"]

The center of mass (COM) is the density-weighted position...

avgCOM[tt_] = (Plus @@(mass Table[{x[i][t], y[i][t],z[i][t]}/. sol1 /. t -> tt, {i, 1, nBodies}]))/(Plus @@ mass);

sp2=ListPointPlot3D[Table[avgCOM[t],{t,0,5,0.1`}],PlotStyle->{Red,PointSize[0.02`]},PlotLabel->"COM"]

Show[sp1,sp2,PlotLabel->"Averages"]

Removing the COM

Since the mass of each object is identical, we can easily sum (and average) the initial velocities and positions of the objects in order to find the center of gravity and the center of mass velocity.

Print["Number of planetary objects is ",nBodies]

Defining the COM

v0=Table[RandomReal[{-1.,1.},3],{i,1,nBodies}]

vcom = MapThread[Plus, v0 mass]/(Plus @@ mass)

r0=Table[RandomReal[{-1.,1.},3],{i,1,nBodies}];

rcom = MapThread[Plus, r0 mass]/(Plus @@ mass)

p0 = Table[{x[i][0] == r0[[i,1]] - rcom[[1]],
y[i][0] == r0[[i,2]] - rcom[[2]],
z[i][0] == r0[[i,3]] - rcom[[3]], vx[i][0] == v0[[i,1]] - vcom[[1]], vy[i][0] == v0[[i,2]] - vcom[[2]],vz[i][0] == v0[[i,3]]- vcom[[3]]}, {i, 1, nBodies}] // Flatten;

eqs2 = Thread[D[p,t] == dp]~ Join ~ p0;

Example Solution

Print["nBodies = ",nBodies];(sol2 = NDSolve[eqs2 , p, {t, 0, 5}, MaxSteps -> 50000] // Flatten;) // Timing

Animation

objects2[tt_] := Graphics3D[{PointSize[0.04],Red, Point[{x[1][t], y[1][t], z[1][t]}], PointSize[0.02],Blue} ~ Join ~Table[Point[{x[i][t], y[i][t], z[i][t]}] , {i, 2, nBodies}] /. sol2 /. t -> tt, PlotRange->threeRange]

Animate[objects2[t], {t, 0.0, 5.0, 0.1},AnimationRunning-> False]

ParametricPlot3D[Evaluate[{x[i][t],y[i][t],z[i][t]}/.i->1/.sol2],{t,0,5},PlotPoints->200,PlotLabel->"Heavy Object"]

ParametricPlot3D[Evaluate[{x[i][t],y[i][t],z[i][t]}/.i->3/.sol2],{t,0,5},PlotPoints->200,PlotLabel->"Light Object"]

Now, what happens when nBodies -> 2?

Averages

We can find the average radius and the average position.

avgRadius[tt_] = (Plus @@Table[Sqrt[x[i][t]^2 + y[i][t]^2 + z[i][t]^2] /. sol2/. t -> tt, {i, 1, nBodies}] )/nBodies;

Plot[avgRadius[t],{t,0,5},PlotLabel->"Average Radius"]

avgPosition[tt_] = (Plus @@ Table[{x[i][t], y[i][t],z[i][t]}/. sol2 /. t -> tt, {i, 1, nBodies}])/nBodies;

sp1=ListPointPlot3D[Table[avgPosition[t],{t,0,5,0.1`}],PlotStyle->{Blue,PointSize[0.02`]},PlotLabel->"Average Position"]

avgCOM[tt_] = (Plus @@(mass Table[{x[i][t], y[i][t],z[i][t]}/. sol2 /. t -> tt, {i, 1, nBodies}]))/(Plus @@ mass);

sp2=ListPointPlot3D[Table[avgCOM[t],{t,0,5,0.1}],PlotStyle->{Red,PointSize[0.02]},PlotLabel->"COM"]

The COM remains "zero" within a small numerical error.

Show[sp1,sp2,PlotLabel->"Averages"]

Nearby Orbits

Let's choose two nearby initial conditions: the second differing only by a slightly different velocity.

tMax = 5.0;

v0a=Table[RandomReal[{-1.`,1.`},3],{i,1,nBodies}]

v0b = Drop[v0a,-1] ~ Join~ {Last[v0a] {1.001,1,1}}

vcoma = MapThread[Plus, v0a mass]/(Plus @@ mass);
vcomb = MapThread[Plus, v0b mass]/(Plus @@ mass);

r0=Table[RandomReal[{-1.`,1.`},3],{i,1,nBodies}];

rcom = MapThread[Plus, r0 mass]/(Plus @@ mass);

p0a = Table[{x[i][0] == r0[[i,1]] - rcom[[1]],
y[i][0] == r0[[i,2]] - rcom[[2]],
z[i][0] == r0[[i,3]] - rcom[[3]], vx[i][0] == v0a[[i,1]] - vcoma[[1]], vy[i][0] == v0a[[i,2]] - vcoma[[2]],vz[i][0] == v0a[[i,3]]- vcoma[[3]]}, {i, 1, nBodies}] // Flatten;
p0b = Table[{x[i][0] == r0[[i,1]] - rcom[[1]],
y[i][0] == r0[[i,2]] - rcom[[2]],
z[i][0] == r0[[i,3]] - rcom[[3]], vx[i][0] == v0b[[i,1]] - vcomb[[1]], vy[i][0] == v0b[[i,2]] - vcomb[[2]],vz[i][0] == v0b[[i,3]]- vcomb[[3]]}, {i, 1, nBodies}] // Flatten;

eqsA = Thread[D[p,t] == dp]~ Join ~ p0a;
eqsB = Thread[D[p,t] == dp]~ Join ~ p0b;

solA = NDSolve[eqsA , p, {t, 0, tMax}, MaxSteps -> 50000] // Flatten;
solB = NDSolve[eqsB , p, {t, 0, tMax}, MaxSteps -> 50000] // Flatten;

Orbit Separation

dr[t_] = ({x[i][t], y[i][t], z[i][t]} /. i->nBodies /. solB) - ({x[i][t], y[i][t], z[i][t]} /. i->nBodies /. solA)

ParametricPlot3D[dr[t],{t,0,tMax},PlotPoints->200,PlotLabel->"Orbital Difference"]

Plot[Evaluate[Log[10,(Sqrt[dr[t].dr[t]])]],{t,0,tMax},PlotPoints->200,PlotLabel->"Log Orbital Difference",Frame->True,PlotRange->{-6,0}]

Two objects that begin nearby deviate exponentially fast!! (This is a measure of chaos.)

Summary

We have shown how to use Mathematica to solve the nonlinear, coupled differential equations for many interacting objects in three dimensions. This is a very complicated problem. For n interacting objects, there are 6 × n equations to solve. Additionally, computation of the force acting on any object grows as (n - 1). This implies that computation of the interacting forces grows as n × (n -1) ~ (n^2). Finally, these equations can not be solved in closed form, and they exhibit sensitivity to very small changes in the initial conditions.

Extra notes:

Schemes for many body problems with computations scaling like N log N are possible using a recursive "cell averaging" procedure. 
Say we divide N bodies into M cells with about N/M bodies per cell. We need to sort these bodies so that all of the cells are separated while the bodies within a cell are relatively close. Then, the acceleration of a body within a cell equals to the accelerations due to the other (M - 1) cells plus the accelerations due to other bodies within the cell. The number of operations is roughly of order (M^2)+ (N^2)/M. When M ~ (N^(2/3)), then the number of operations are minimized, saving considerable computational time.