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

Tôi, người giữ bản quyền tác phẩm này, từ đây phát hành nó theo giấy phép sau:

Tập tin này được phân phối theo Creative Commons Hiến tặng vào Phạm vi Công cộng Toàn thế giới CC0.

Người nào gán tài liệu này với tác phẩm nghĩa là đã hiến tác phẩm cho phạm vi công cộng bằng cách từ bỏ mọi quyền lợi của người đó đối với tác phẩm theo quy định của luật bản quyền, có hiệu lực trên toàn thế giới và các quyền lợi pháp lý phụ mà người đó có được trong tác phẩm, đến mức độ mà luật pháp cho phép. Bạn được tự do sao chép, phân phối, và biểu diễn tác phẩm này, tất cả đều không bắt buộc ghi công. http://creativecommons.org/publicdomain/zero/1.0/deed.enCC0Creative Commons Zero, Public Domain Dedicationfalsefalse

Tập tin này, vốn được đăng tải tại https://twitter.com/j_bertolotti/status/1044947721696808961, đã được người duyệt hình Ronhjones kiểm tra vào ngày 19 tháng 10 năm 2018 và xác nhận rằng nó đã được phát hành dưới giấy phép tương ứng trong ngày hôm đó.