用scipy-odein在python中求解向量常微分方程

Question

我正在尝试解决，ode涉及向量，不能得出一个可行的答案。所以我把它分成6个分量，一个是分量的时间导数，一个是速度分量的时间导数。第一个值似乎是合理的，但后来它跳到了数百万，我不知道为什么。老实说，我不太确定该怎么做，现在就试试看。我似乎在网上找不到任何信息，如果有此类问题的例子，我可以使用一些帮助或链接。任何信息都将是非常感谢如何得到这一点来解决颂歌。在

def dr_dt(y, t):
"""Integration of the governing vector differential equation.
d2r_dt2 = -(mu/R^3)*r with d2r_dt2 and r as vecotrs.
Initial position and velocity are given.
y[0:2] = position components
y[3:] = velocity components"""

    G = 6.672*(10**-11)
    M = 5.972*(10**24)
    mu = G*M
    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2) 

    dy0 = y[3]
    dy1 = y[4]
    dy2 = y[5]
    dy3 = -(mu / (r**3)) * y[0]
    dy4 = -(mu / (r**3)) * y[1]
    dy5 = -(mu / (r**3)) * y[2]
    return [dy0, dy3, dy1, dy4, dy2, dy5]

在这个问题解决后，我想把它画出来。它应该是一个椭圆，但说实话，我也不太清楚该怎么做。我在想，把这个位置的大小，然后用时间来规划它。如果有更好的方法，请随时告诉我。在

谢谢。在

Answer 1

首先，如果您的y[:3]是位置，y[3:]是速度，那么dr_dt函数应该按照这个顺序返回组件。其次，为了绘制轨迹，我们可以使用优秀的matplotlib mplot3d模块，或者省略位置和速度的第z个分量（因此我们的运动在XY平面上），并绘制y与{}。示例代码（返回值的正确顺序）如下所示：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def dr_dt(y, t):
    """Integration of the governing vector differential equation.
    d2r_dt2 = -(mu/R^3)*r with d2r_dt2 and r as vecotrs.
    Initial position and velocity are given.
    y[0:2] = position components
    y[3:] = velocity components"""

    G = 6.672*(10**-11)
    M = 5.972*(10**24)
    mu = G*M
    r = np.sqrt(y[0]**2 + y[1]**2 + y[2]**2).

    dy0 = y[3]
    dy1 = y[4]
    dy2 = y[5]
    dy3 = -(mu / (r**3)) * y[0]
    dy4 = -(mu / (r**3)) * y[1]
    dy5 = -(mu / (r**3)) * y[2]
    return [dy0, dy1, dy2, dy3, dy4, dy5]

t = np.arange(0, 100000, 0.1)
y0 = [7.e6, 0., 0., 0., 1.e3, 0.]
y = odeint(dr_dt, y0, t)
plt.plot(y[:,0], y[:,1])
plt.show()

这就产生了一个漂亮的椭圆：