使用Mathematica无法解决双体问题？

Question

编辑:

@auxsvr是正确的,我有力方程式错误,关于-3/2指数.

另一种方法是将它简单地看到2维并考虑从原点起作用的力,与1/r ^ 2成正比,就像重力一样,其中r是距离原点的距离.

在(x,y)处,力作用于方向(-x,-y).然而,这只是方向,而不是幅度.如果我们使用k作为比例常数,则力是(-kx,-ky).

因此,力的大小是Sqrt [( - kx)^ 2 +( - ky)^ 2],或k*Sqrt [x ^ 2 + y ^ 2],或k*Sqrt [r ^ 2]或k*[R

由于力大小也是1/r ^ 2,这给出了k = 1/r ^ 3.

因此力是(-x/r ^ 3,-y/r ^ 3).

由于我最初使用r ^ 2作为我的主要数量,那是(r ^ 2)^( - 3/2),这是3/2来自的地方.

这实际上使我的问题无效,尽管它仍然是一个有趣的理论讨论.

我用正确的方程重试了这个Mathematica,但仍然没有得到答案.正如其他人指出的那样,在某些条件下结果只是一个椭圆(在其他情况下可能是抛物线或双曲线).

另外,虽然最终的轨道是圆锥形截面,但是初始轨道可以螺旋进出,直到达到最终的圆锥截面轨道.

编辑结束在这里

我正在使用Mathematica来解决两体问题:

DSolve[{
    d2[t] == (x1[t]-x0[t])^2 + (y1[t]-y0[t])^2 + (z1[t]-z0[t])^2,
    D[x0[t], t,t] == (x1[t]-x0[t])/d2[t],
    D[y0[t], t,t] == (y1[t]-y0[t])/d2[t],
    D[z0[t], t,t] == (z1[t]-z0[t])/d2[t],
    D[x1[t], t,t] == -(x1[t]-x0[t])/d2[t],
    D[y1[t], t,t] == -(y1[t]-y0[t])/d2[t],
    D[z1[t], t,t] == -(z1[t]-z0[t])/d2[t]
    },
    {x0,y0,z0,x1,y1,x1,d2},
    t
]

但我回来了:

因变量比方程式更少,因此系统超定.

我算7个方程和7个因变量？

事实上,该系统是半未确定的,因为我没有在0时提供位置和速度.

我意识到我的方程本身可能对于两体问题是错误的,但我仍然想知道为什么Mathematica抱怨这个.

Answer 1

怎么样NDSolve？

d2[t_] = (-x0[t] + x1[t])^2 + (-y0[t] + y1[t])^2 + (-z0[t] + 
    z1[t])^2; sol = {x0, y0, z0, x1, y1, z1} /. 
  NDSolve[{x0''[t] == (-x0[t] + x1[t])/d2[t], 
     y0''[t] == (-y0[t] + y1[t])/d2[t], 
     z0''[t] == (-z0[t] + z1[t])/d2[t], x1''[t] == -x0''[t], 
     y1''[t] == -y0''[t], z1''[t] == -z0''[t], x0[0] == 0, y0[0] == 0,
      z0[0] == 0, x1[0] == 1, y1[0] == 0, z1[0] == 0, x0'[0] == -0.5, 
     y0'[0] == 1, z0'[0] == 0.5, x1'[0] == 0.5, y1'[0] == -1, 
     z1'[0] == -0.5}, {x0, y0, z0, x1, y1, z1}, {t, 0, 120}][[1]]

r = 3;
 Animate[
  Graphics3D[
   {
    PointSize -> 0.05,
    Point[{sol[[1]][t], sol[[2]][t], sol[[3]][t]}],
    Point[{sol[[4]][t], sol[[5]][t], sol[[6]][t]}],
    Red,
    Line[Table[{sol[[1]][t1], sol[[2]][t1], sol[[3]][t1]}, {t1, 0, t, 0.1}]],
    Green,
    Line[Table[{sol[[4]][t1], sol[[5]][t1], sol[[6]][t1]}, {t1, 0, t, 0.1}]]
   }, 
   PlotRange -> {{-r, r}, {-r, r}, {-r, r}}
  ], {t, 0, 120}, AnimationRate -> 4
 ]

Answer 2

我很惊讶没有人注意到每个人都错误地写出了运动方程,这从图中可以看出,因为两个物体的引力势中的有界轨道总是闭合的(Bertrand定理).正确的运动方程是

{x0''[t] == (-x0[t] + x1[t])/d2[t]^(3/2), 
 y0''[t] == (-y0[t] + y1[t])/d2[t]^(3/2),
 x1''[t] == -x0''[t], 
 y1''[t] == -y0''[t]}

同

d2[t_]:= (x1[t]-x0[t])^2 + (y1[t]-y0[t])^2

因为对于中心力场,运动是平面的.此外,必须适当地设定初始条件,否则质心移动并且轨道不再是圆锥形截面.