数学或python中的轨道

Question

我有 4 个差异。方程（代表植物的轨道方程）

x'[t] == px[t] + y[t]
y'[t] == py[t] - x[t]
px'[t] == py[t] - dVx[t]
py'[t] == -px[t] - dVy[t]

我想在任何时间 t 求解 x[t] 和 y[t]。给定的变量是

x[0]==0
y[0]==0
px[0]==0
py[0]==2.0
\[Epsilon]==0.2

-dVx[t] == x[t] - (1 - \[Epsilon])*(x[t] + \[Epsilon])/((x[t] + \[Epsilon])^2 + 
           y[t]^2)^(3/2) - \[Epsilon] (x[t] + \[Epsilon] - 
           1)/(((x[t] + \[Epsilon] - 1)^2 + y[t]^2)^(3/2))

-dVy[t] == y[t]*(1 - (1 - \[Epsilon])/((x[t] + \[Epsilon])^2 + y[t]^2)^(3/
           2) - \[Epsilon]/((x[t] + \[Epsilon] - 1)^2 + y[t]^2)^(3/2))

我怎样才能得到 x,y 任何时间并在 x,y 平面上绘制图。我用 NDSolve 尝试过，但失败了。我的代码是

In[49]:= -dVx[t] == x[t] - (1 - \[Epsilon])*(x[t] + \[Epsilon])/((x[t] + \ 
        [Epsilon])^2 + y[t]^2)^(3/2) - \[Epsilon] (x[t] + \[Epsilon] - 
         1)/(((x[t] + \[Epsilon] - 1)^2 + y[t]^2)^(3/2))

Out[49]= -dVx[t] == -(0.16/(0.04 + y[t]^2)^(3/2)) + 
          0.16/(0.64 + y[t]^2)^(3/2)

In[50]:= -dVy[t] == 
          y[t]*(1 - (1 - \[Epsilon])/((x[t] + \[Epsilon])^2 + y[t]^2)^(3/
          2) - \[Epsilon]/((x[t] + \[Epsilon] - 1)^2 + y[t]^2)^(3/2))

Out[50]= -dVy[t] == 
         y[t] (1 - 0.8/(0.04 + y[t]^2)^(3/2) - 0.2/(0.64 + y[t]^2)^(3/2))

In[56]:= DSolve[{x'[t] == px[t] + y[t], y'[t] == py[t] - x[t], 
         px'[t] == py[t] - dVx[t], py'[t] == -px[t] - dVy[t], px[0] == 0, 
         y[0] == 0, py[0] == 2.0, x[0] == 0, \[Epsilon] == 0.2}, {x[t], 
         y[t]}, t]

During evaluation of In[56]:= DSolve::dsfun: 0 cannot be used as a function.

Out[56]= DSolve[{Derivative[1][x][t] == px[t] + y[t], 
         Derivative[1][y][t] == 2., Derivative[1][px][t] == 2. - dVx[t], 
         Derivative[1][py][t] == -dVy[t] - px[t], px[0] == 0, y[0] == 0, 
         py[0] == 2., True, True}, {0, y[t]}, t]

我是 mathematica 的新手，很高兴能得到任何帮助。如果这样更容易，我可以使用 python

Answer 1

许多语法错误。试试这个：

\[Epsilon] = 0.2;
dVx = -(x[
      t] - (1 - \[Epsilon])*(x[
          t] + \[Epsilon])/((x[t] + \[Epsilon])^2 + y[t]^2)^(3/
          2) - \[Epsilon] (x[t] + \[Epsilon] - 
         1)/(((x[t] + \[Epsilon] - 1)^2 + y[t]^2)^(3/2)));
dVy = -(y[
      t]*(1 - (1 - \[Epsilon])/((x[t] + \[Epsilon])^2 + y[t]^2)^(3/
           2) - \[Epsilon]/((x[t] + \[Epsilon] - 1)^2 + y[t]^2)^(3/
           2)));
NDSolve[{
  x'[t] == px[t] + y[t],
  y'[t] == py[t] - x[t],
  px'[t] == py[t] - dVx,
  py'[t] == -px[t] - dVy,
  px[0] == 0,
  y[0] == 0,
  py[0] == 2,
  x[0] == 0
  }, {x[t], y[t], px[t], py[t]}, {t, 0, 1}]