(*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]

我，本作品著作权人，特此采用以下许可协议发表本作品：

	本作品采用知识共享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