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

Mi, la posedanto de la aŭtorrajto por ĉi tiu verko, ĉi-maniere publikigas tiun laŭ la jena permesilo:

Ĉi tiu dosiero estas disponebla laŭ la Krea Komunaĵo CC0 1.0 Universala Publikaĵiga Dediĉo.

La persono kiu asociis verkon kun ĉi tiu faro dediĉis la verkon kiel publikaĵon forlasante ĉiujn siajn rajtojn al la verko mondvaste sub aŭtorrajta leĝo inkluzivante ĉiujn rilatajn aŭ najbarajn rajtojn permesitajn de leĝo. Oni povas kopii, modifi, distribui kaj ludi la verkon, eĉ por komercaj kialoj, ĉiuj sen peto de permeso. http://creativecommons.org/publicdomain/zero/1.0/deed.enCC0Creative Commons Zero, Public Domain Dedicationfalsefalse

This file, which was originally posted to https://twitter.com/j_bertolotti/status/1044947721696808961, was reviewed on 19 October 2018 by reviewer Ronhjones, who confirmed that it was available there under the stated license on that date.