如您所显示的，您可以将其编写为一个包含六个一阶ode的系统：
x' = x2y' = y2z' = z2x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * xy2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * yz2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z
您可以将其另存为矢量：
u = (x, y, z, x2, y2, z2)
并因此创建一个返回其派生函数：
def deriv(u, t): n = -mu / np.sqrt(u[0]**2 + u[1]**2 + u[2]**2) return [u[3], # u[0]' = u[3] u[4], # u[1]' = u[4] u[5], # u[2]' = u[5] u[0] * n, # u[3]' = u[0] * n u[1] * n, # u[4]' = u[1] * n u[2] * n] # u[5]' = u[2] * n
给定一个初始状态
u0 = (x0, y0, z0, x20, y20,z20)和一个time变量
t，可以这样输入
scipy.integrate.odeint：
u = odeint(deriv, u0, t)
u上面的列表在哪里。或者你也可以解包
u从一开始，并忽略这些值
x2，
y2和
z2（你必须先转输出
.T）
x, y, z, _, _, _ = odeint(deriv, u0, t).T
