(*Staring positions in a triangle*)
x10 = -1;
y10 = -1;
x20 = 1;
y20 = -1;
x30 = 1;
y30 = 1;
(*Initial total momentum is zero, so the center of mass does not \
drift away*)
vx10 = 0.2;
vy10 = 0;
vx20 = -0.1;
vy20 = 0;
vx30 = 0;
vy30 = -0.1;
(*max time the system evolves (in arbitrary units)*)
T = 40;
(*All three bodies have the same mass*)
m1 = 1;
m2 = 1;
m3 = 1;
(*Setting up of the equations copied from \
http://demonstrations.wolfram.com/PlanarThreeBodyProblem/
There are more elegant and compact ways of doing this, but I wasn't \
interested in optimizing the code.*)
nds = NDSolve[
   {x1'[t] == vx1[t], y1'[t] == vy1[t],
    x2'[t] == vx2[t], y2'[t] == vy2[t],
    x3'[t] == vx3[t], y3'[t] == vy3[t],
    m1 vx1'[t] == -((
       m1 m2 (x1[t] - 
          x2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2)) - (
      m1 m3 (x1[t] - x3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(
      3/2), m1 vy1'[t] == -((
       m1 m2 (y1[t] - 
          y2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(3/2)) - (
      m1 m3 (y1[t] - y3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(
      3/2), m2 vx2'[t] == (
      m1 m2 (x1[t] - x2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(
      3/2) - (m2 m3 (x2[t] - 
         x3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2), 
    m2 vy2'[t] == (
      m1 m2 (y1[t] - y2[t]))/((x1[t] - x2[t])^2 + (y1[t] - y2[t])^2)^(
      3/2) - (
      m2 m3 (y2[t] - y3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(
      3/2), m3 vx3'[t] == (
      m1 m3 (x1[t] - x3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(
      3/2) + (m2 m3 (x2[t] - 
         x3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2), 
    m3 vy3'[t] == (
      m1 m3 (y1[t] - y3[t]))/((x1[t] - x3[t])^2 + (y1[t] - y3[t])^2)^(
      3/2) + (m2 m3 (y2[t] - 
         y3[t]))/((x2[t] - x3[t])^2 + (y2[t] - y3[t])^2)^(3/2),
    x1[0] == x10, y1[0] == y10, x2[0] == x20, y2[0] == y20, 
    x3[0] == x30, y3[0] == y30,
    vx1[0] == vx10, vy1[0] == vy10, vx2[0] == vx20, vy2[0] == vy20, 
    vx3[0] == vx30, vy3[0] == vy30},
   {x1, x2, x3, y1, y2, y3, vx1, vx2, vx3, vy1, vy2, vy3}, {t, 0, 
    T}];
funsToPlot = {{x1[t], y1[t]}, {x2[t], y2[t]}, {x3[t], y3[t]}} /. 
   nds[[1]];
evo = Table[funsToPlot /. {t -> j}, {t, 0, T, 0.01}];
dim = Dimensions[evo][[1]];
(*For the elastic force case I used a Verlet integration, as this \
case is numerically very stable.*)
np = 3;
k0 = 1;
(*Same initial condition as the gravitational case*)

pos = {{x10, y10}, {x20, y20}, {x30, y30}};
v0 = {{vx10, vy10}, {vx20, vy20}, {vx30, vy30}};
acc = Table[
   Sum[If[j == k, 0, -k0 (pos[[j]] - pos[[k]])], {k, 1, np}], {j, 1, 
    np}];
dt = 0.005;
posold = pos;
pos = posold + v0 dt + acc/2 dt^2;
range = 5;

evoe = Reap[Do[
      acc = 
       Table[Sum[
         If[j == k, 0, -k0 (pos[[j]] - pos[[k]])], {k, 1, np}], {j, 1,
          np}];
      posoldold = posold;
      posold = pos;
      pos = 2 posold - posoldold + acc dt^2;
      Sow[pos];
      , dim];][[2, 1]];
plots = Table[
   GraphicsRow[{
     Show[
      ListPlot[{evo[[All, 1]][[1 ;; j]], evo[[All, 2]][[1 ;; j]], 
        evo[[All, 3]][[1 ;; j]]}, PlotStyle -> {Purple, Orange, Cyan},
        PlotRange -> {{-range, range}, {-range, range}}, 
       Joined -> True, Axes -> False, 
       PlotLabel -> "Gravitational 3-body problem", 
       LabelStyle -> {Bold, Black}],
      Graphics[{PointSize[0.03], Purple, Point[evo[[All, 1]][[j]]], 
        Orange, Point[evo[[All, 2]][[j]]], Cyan, 
        Point[evo[[All, 3]][[j]]]} , 
       PlotRange -> {{-range, range}, {-range, range}}], 
      ImageSize -> Medium
      ]
     ,
     Show[
      ListPlot[{evoe[[All, 1]][[1 ;; j]], evoe[[All, 2]][[1 ;; j]], 
        evoe[[All, 3]][[1 ;; j]]}, 
       PlotStyle -> {Purple, Orange, Cyan}, 
       PlotRange -> {{-range, range}, {-range, range}}, 
       Joined -> True, Axes -> False, 
       PlotLabel -> "Elastic 3-body problem", 
       LabelStyle -> {Bold, Black}],
      Graphics[{PointSize[0.03], Purple, Point[evoe[[All, 1]][[j]]], 
        Orange, Point[evoe[[All, 2]][[j]]], Cyan, 
        Point[evoe[[All, 3]][[j]]]} , 
       PlotRange -> {{-range, range}, {-range, range}}], 
      ImageSize -> Medium
      ]
     }], {j, 1, dim, 20}];
ListAnimate[plots]

我，本作品的著作權持有者，決定用以下授權條款發佈本作品：

	此檔案於創用 CC CC0 1.0 通用公有領域貢獻宣告下分發。
	在此宣告之下分發本作品者，已依據各國著作權法，在全世界放棄其對本作品所擁有的著作權及所有相關相似的法律權利，從而將本作品貢獻至公有領域。您可以複製、修改、分發和演示該作品，用於任何商業用途，所有這些都不需要請求授權。 http://creativecommons.org/publicdomain/zero/1.0/deed.enCC0Creative Commons Zero, Public Domain Dedicationfalsefalse

圖片審查員Ronhjones，確認本圖片於2018年10月19日可在下列站點找到並符合所選許可證： https://twitter.com/j_bertolotti/status/1044947721696808961