{x10, y10, z10} = {2, 2, 0}; {vx10, vy10, vz10} = {1, 2.9, 0.1};
{x20, y20, z20} = {-3, -3, 0}; {vx20, vy20, vz20} = {-1, 2, -0.2};
{x30, y30, z30} = {5, -2, 0}; {vx30, vy30, vz30} = {1.7, -0.1, 0.2};
m1 = 1; m2 = 1; m3 = 1;
G = 15;
TMAX = 200;
DeltaTime = TMAX/200;
MyLTBR = 100;
vcenter = (m1*{vx10, vy10, vz10} + m2*{vx20, vy20, vz20} +
m3*{vx30, vy30, vz30})/3;{vx10, vy10, vz10} -= vcenter;
{vx20, vy20, vz20} -= vcenter;
{vx30, vy30, vz30} -= vcenter;
F21x[t_] =
G*m1*m2*(x2[t] -
x1[t])/((y2[t] - y1[t])^2 + (x2[t] - x1[t])^2 + (z2[t] -
z1[t])^2)^1.5;
F21y[t_] :=
G*m1*m2*(y2[t] -
y1[t])/((y2[t] - y1[t])^2 + (x2[t] - x1[t])^2 + (z2[t] -
z1[t])^2)^1.5;
F21z[t_] :=
G*m1*m2*(z2[t] -
z1[t])/((y2[t] - y1[t])^2 + (x2[t] - x1[t])^2 + (z2[t] -
z1[t])^2)^1.5;
F31x[t_] :=
G*m1*m3*(x3[t] -
x1[t])/((y3[t] - y1[t])^2 + (x3[t] - x1[t])^2 + (z3[t] -
z1[t])^2)^1.5;
F31y[t_] :=
G*m1*m3*(y3[t] -
y1[t])/((y3[t] - y1[t])^2 + (x3[t] - x1[t])^2 + (z3[t] -
z1[t])^2)^1.5;
F31z[t_] :=
G*m1*m3*(z3[t] -
z1[t])/((y3[t] - y1[t])^2 + (x3[t] - x1[t])^2 + (z3[t] -
z1[t])^2)^1.5;
F32x[t_] :=
G*m2*m3*(x3[t] -
x2[t])/((y3[t] - y2[t])^2 + (x3[t] - x2[t])^2 + (z3[t] -
z2[t])^2)^1.5;
F32y[t_] :=
G*m2*m3*(y3[t] -
y2[t])/((y3[t] - y2[t])^2 + (x3[t] - x2[t])^2 + (z3[t] -
z2[t])^2)^1.5;
F32z[t_] :=
G*m2*m3*(z3[t] -
z2[t])/((y3[t] - y2[t])^2 + (x3[t] - x2[t])^2 + (z3[t] -
z2[t])^2)^1.5;
p = NDSolve[
{m1*x1''[t] == F21x[t] + F31x[t],
m1*y1''[t] == F21y[t] + F31y[t],
m1*z1''[t] == F21z[t] + F31z[t],
m2*x2''[t] == -F21x[t] + F32x[t],
m2*y2''[t] == -F21y[t] + F32y[t],
m2*z2''[t] == -F21z[t] + F32z[t],
m3*x3''[t] == -F31x[t] - F32x[t],
m3*y3''[t] == -F31y[t] - F32y[t],
m3*z3''[t] == -F31z[t] - F32z[t],
x1[0] == x10, x1'[0] == vx10,
y1[0] == y10, y1'[0] == vy10,
z1[0] == z10, z1'[0] == vz10,
x2[0] == x20, x2'[0] == vx20,
y2[0] == y20, y2'[0] == vy20,
z2[0] == z20, z2'[0] == vz20,
x3[0] == x30, x3'[0] == vx30,
y3[0] == y30, y3'[0] == vy30,
z3[0] == z30, z3'[0] == vz30}, {x1, x2, x3, y1, y2, y3, z1, z2,
z3}, {t, 0, TMAX}];
{r1, r2, r3} = {m1^(1/3), m2^(1/3), m3^(1/3)};
{{px1, px2, px3, py1, py2, py3, pz1, pz2, pz3}} = {x1, x2, x3, y1, y2,
y3, z1, z2, z3} /. p;
GraphList = Table[
Show[{ParametricPlot3D[{{px1[t], py1[t], pz1[t]}, {px2[t], py2[t],
pz2[t]}, {px3[t], py3[t], pz3[t]}}, {t, 0, tmax},
PlotStyle -> {Darker[Red], Darker[Green], Darker[Blue]},
PlotRange -> All],
Graphics3D[
{Sphere[{px1[tmax], py1[tmax], pz1[tmax]}, r1],
Sphere[{px2[tmax], py2[tmax], pz2[tmax]}, r2],
Sphere[{px3[tmax], py3[tmax], pz3[tmax]}, r3]}]}], {tmax,
DeltaTime, TMAX, DeltaTime}]
两个小球重合到一起了
而且最后的运动也不符合力学原理啊，都在走远？？？求大神

绝对刚体？..

有没有想过，这种三体问题，涉及混沌现象的可能有微扰放大效应。

与天体运行问题相关的微分方程求解中的一个常见问题是，变量单位设置不当以致出现“黑洞”。同学你那个等于15的万有引力常量和个位数的间距真的没问题吗？

完全没有问题。
相撞/逃逸是几乎必然的现象，除非是极其少数的精心设计的数据。
逃逸也是完全符合力学原理的；二体运动高中应该学过，其中有天体无限相互远离的情况，三体运动也类似。