正如您所展示的,您可以将其写为六个一阶颂歌的系统:
x' = x2
y' = y2
z' = z2
x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x
y2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y
z2' = -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)
,以及时间变量t
,这可以输入scipy.integrate.odeint
像这样:
u = odeint(deriv, u0, t)
where u
将是上面的列表。或者你可以解压u
从一开始就忽略这些值x2
, y2
, and z2
(您必须先将输出转置为.T
)
x, y, z, _, _, _ = odeint(deriv, u0, t).T