Descriere fișier

Mathematica 11.0 code

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

Licențiere

Eu, deținătorul drepturilor de autor ale acestei opere, prin prezenta îmi public lucrarea sub următoarea licență:

Acest fișier a fost eliberat sub licența Creative Commons CC0 1.0 Universal Public Domain Dedication.

Persoana care a asociat o operă cu acest document o oferă domeniului public, renunțând la toate drepturile asupra operei, în toată lumea, atât în ce privește drepturile de autor cât și orice alte drepturi juridice conexe pe care le avea asupra operei, în măsura permisă de lege. Puteți copia, modifica sau distribui opera, inclusiv în scopuri comerciale, fără a fi necesară permisiunea autorului. 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.