ピタゴラス3体問題 (Burrau’s problem of three bodies) https://taroyabuki.github.io/2009/06/26/burraus-problem-of-three-bodies/

pythagoras3body.m

(*Mass*)
m = {3, 4, 5};
n = Length@m;

(*Coordination*)
p = Array[{Subscript[px, #][t], Subscript[py, #][t]} &, n];
q = Array[{Subscript[ x, #][t], Subscript[ y, #][t]} &, n];

initialCondition = {
    p == {{0, 0}, { 0,  0}, {0,  0}},
    q == {{1, 3}, {-2, -1}, {1, -1}}} /. t -> 0;

(*Total kinetic energy*)
T = Sum[(p[[i]].p[[i]])/2/m[[i]], {i, 1, n}];

(*Total potential energy*)
d[x_] := Sqrt[x.x]
V = Total[With[{i = #[[1]], j = #[[2]]},
      -m[[i]] m[[j]]/d[q[[i]] - q[[j]]]] & /@ Subsets[Range[n], {2}]];

(*Hamiltonian*)
H = T + V;

(*Total momentum*)
P = Total[p];

(*Total angular momentum*)
(*M = Total[Array[Cross[Append[q[[#]], 0], Append[p[[#]], 0]][[3]] &, n]];*)
M = Total[Array[Det[{q[[#]], p[[#]]}] &, n]];

(*Equations of motion*)
eqns = Array[{
     D[p[[#]], t] == -D[H, {q[[#]]}],
     D[q[[#]], t] ==  D[H, {p[[#]]}]} &, n];

tmax = 70;
s = NDSolve[Flatten[{eqns, initialCondition}], Flatten[{p, q}], {t, 0, tmax},
   Method -> {"Projection",
     Method -> "ExplicitRungeKutta",
     "Invariants" -> {H, P[[1]], P[[2]], M}},
   StartingStepSize -> 1/20];

plot = ParametricPlot[Evaluate[q /. s], {t, 0, tmax}]

c = {RGBColor[{94, 129, 181}/255], RGBColor[{215, 156, 36}/255], RGBColor[{143, 177, 49}/255]};

plots = Table[Show[{
     Graphics[Table[{PointSize[0.05], c[[i]], Point[Evaluate[q[[i]] /. s /. t -> time]]}, {i, 1, 3}]], 
     ParametricPlot[Evaluate[q /. s], {t, 0, time}]},
    PlotRange -> {{-5, 5}, {-5, 5}}], {time, 0.01, tmax}];

ListAnimate[plots]

hpm = GraphicsGrid[{Plot[Evaluate[# /. s], {t, 0, tmax}, PlotRange -> All]} & /@ {H, P, M}]