Given the starting positions, velocities, and masses of three objects interacting via gravity, the **classical three-body problem** involves determining the motions of the three particles throughout time.

What’s cool about the three-body system is that it’s *impossible* to solve for the motions of the objects exactly. That is, we can’t write down an *equation* that describes the system. Instead of finding an exact solution, we solve the system numerically, which amounts to finding an accurate approximation.

The three-body problem is an example of a chaotic system, meaning that even a slight change in the starting conditions *drastically* changes the time-evolution of the system.

The GIF above shows a planar (i.e., two-dimensional) three-body system.

Mathematica code:

G = 1; time = 30;
mA = 1; xA0 = 0; yA0 = 0; vxA0 = 0; vyA0 = 0;
mB = 1; xB0 = 1; yB0 = 0; vxB0 = 0; vyB0 = 0;
mC = 1; xC0 = 0; yC0 = 0.8; vxC0 = 0; vyC0 = 0;
soln1 = NDSolve[
{mA*Derivative[2][xA][t] ==
-((G*mA*mB*(xA[t] - xB[t]))/((xA[t] - xB[t])^2 +
(yA[t] - yB[t])^2)^(3/2)) -
(G*mA*mC*(xA[t] - xC[t]))/((xA[t] - xC[t])^2 +
(yA[t] - yC[t])^2)^(3/2),
mA*Derivative[2][yA][t] ==
-((G*mA*mB*(yA[t] - yB[t]))/((xA[t] - xB[t])^2 +
(yA[t] - yB[t])^2)^(3/2)) -
(G*mA*mC*(yA[t] - yC[t]))/((xA[t] - xC[t])^2 +
(yA[t] - yC[t])^2)^(3/2),
mB*Derivative[2][xB][t] ==
-((G*mB*mC*(xB[t] - xC[t]))/((xB[t] - xC[t])^2 +
(yB[t] - yC[t])^2)^(3/2)) -
(G*mB*mA*(xB[t] - xA[t]))/((xB[t] - xA[t])^2 +
(yB[t] - yA[t])^2)^(3/2),
mB*Derivative[2][yB][t] ==
-((G*mB*mC*(yB[t] - yC[t]))/((xB[t] - xC[t])^2 +
(yB[t] - yC[t])^2)^(3/2)) -
(G*mB*mA*(yB[t] - yA[t]))/((xB[t] - xA[t])^2 +
(yB[t] - yA[t])^2)^(3/2),
mC*Derivative[2][xC][t] ==
-((G*mC*mA*(xC[t] - xA[t]))/((xC[t] - xA[t])^2 +
(yC[t] - yA[t])^2)^(3/2)) -
(G*mC*mB*(xC[t] - xB[t]))/((xC[t] - xB[t])^2 +
(yC[t] - yB[t])^2)^(3/2),
mC*Derivative[2][yC][t] ==
-((G*mC*mA*(yC[t] - yA[t]))/((xC[t] - xA[t])^2 +
(yC[t] - yA[t])^2)^(3/2)) -
(G*mC*mB*(yC[t] - yB[t]))/((xC[t] - xB[t])^2 +
(yC[t] - yB[t])^2)^(3/2),
xA[0] == xA0, yA[0] == yA0,
Derivative[1][xA][0] == vxA0, Derivative[1][yA][0] == vyA0,
xB[0] == xB0, yB[0] == yB0,
Derivative[1][xB][0] == vxB0, Derivative[1][yB][0] == vyB0,
xC[0] == xC0, yC[0] == yC0,
Derivative[1][xC][0] == vxC0, Derivative[1][yC][0] == vyC0
},
{xA, yA, xB, yB, xC, yC},{t, 0, time}, MaxSteps -> 100000]
x1[t_] := Evaluate[xA[t] /. soln1[[1,1]]]
y1[t_] := Evaluate[yA[t] /. soln1[[1,2]]]
x2[t_] := Evaluate[xB[t] /. soln1[[1,3]]]
y2[t_] := Evaluate[yB[t] /. soln1[[1,4]]]
x3[t_] := Evaluate[xC[t] /. soln1[[1,5]]]
y3[t_] := Evaluate[yC[t] /. soln1[[1,6]]]
Manipulate[Show[
{ParametricPlot[
{{x1[t], y1[t]}, {x2[t], y2[t]}, {x3[t], y3[t]}},
{t, tmax - 0.5, tmax},
Axes -> False, PlotRange -> {{-0.55, 1.45},
{-0.55, 1.08}}, PlotStyle -> {Red, Green, Blue},
GridLines -> {Table[0.25*x + 0.07, {x, -100, 100}],
Table[0.25*y + 0.01, {y, -100, 100}]},
GridLinesStyle -> Directive[LightGray]]},
{Graphics[{Opacity[0.7], EdgeForm[Directive[Black]], Red,
Disk[{x1[tmax], y1[tmax]}, 0.03], Green,
Disk[{x2[tmax], y2[tmax]}, 0.03], Blue,
Disk[{x3[tmax], y3[tmax]}, 0.03]}]}, ImageSize -> 600],
{tmax, 6.05, 16.05}]