ODEINT错误-ValueError:用序列设置数组元素(太阳系的Python模型)

Question

我正在创建一个太阳系的轨道模型，所以我需要将每颗行星的牛顿第二定律积分为一个和。我正在使用ODEINT并得到这个ValueError，但不知道实际问题发生在哪里。这是我的第一篇文章，所以函数下的缩进有问题。。。

def func(Q,t):
G = const.G.value #[m3/(kg s2)]
masses = np.array([1.9891*10**6, 0.33011, 4.8675, 5.9724, 0.64171])*10**24 #[kg]
(x_sun, y_sun, vx_sun, vy_sun,
x_merc, y_merc, vx_merc, vy_merc, 
x_ven, y_ven, vx_ven, vy_ven,
x_earth, y_earth, vx_earth, vy_earth,
x_mars, y_mars, vx_mars, vy_mars) = Q
x = np.array(Q[0::4])
y = np.array(Q[1::4])
vx = np.array(Q[2::4])
vy = np.array(Q[3::4])
accx = np.zeros(len(masses))
accy = np.zeros(len(masses))
for p in range(len(masses)):
for m in masses:
i = list(masses).index(m)
r = np.sqrt((x[2]-x[i])**2 + (y[2]-y[i])**2)
if r == 0:
accx[p] = 0
accy[p] = 0
else:
accx[p] += -G*m * (x[p]-x[i])/r**3
accy[p] += -G*m * (y[p]-y[i])/r**3
dQdt = np.zeros(len(Q))
dQdt[0::4] = vx
dQdt[1::4] = vy     
dQdt[2::4] = accx
dQdt[3::4] = accy
return dQdt 
'''set initial conditions'''
names = np.array(['Sun','Mercury', 'Venus', 'Earth', 'Mars'])
G = const.G.value                                                    #[m3/(kg s2)]
semimajor = np.array([0, 57.909, 108.209, 149.596, 227.923])*10**9      
#[m]
masses = np.array([1.9891*10**6, 0.33011, 4.8675, 5.9724, 0.64171])*10**24  #[kg]
x0 = np.array([0, 46.002, 107.476, 147.092, 206.617])*10**9             #[m]
y0 = 0                                     
vx0=0
vy0 = 1000*(G*const.M_sun.value*(1/1e6)*((2/x0)-(1/semimajor)))**0.5 #[m/s]
IV = [x0, y0, vx0, vy0] 
'''set time parameters'''
t_i = 0
year = 31557600 
t_f = 8*year
N = 100
t = np.linspace(t_i, t_f, N)
'''find solution with built-in integrator'''
soln = spi.odeint(func, IV, t)
...

ValueError                                Traceback (most recent 
call last)
<ipython-input-209-681766e01086> in <module>
62 
63 '''find solution with built-in integrator'''
---> 64 soln = spi.odeint(func, IV, t)
65 solnT = np.transpose(soln)
66 x = solnT[0::4]/AU
~/anaconda3/lib/python3.7/site-packages/scipy/integrate/odepack.py 
in odeint(func, y0, t, args, Dfun, col_deriv, full_output, ml, mu, r
tol, atol, tcrit, h0, hmax, hmin, ixpr, mxstep, mxhnil, mxordn, mxords, printmessg, tfirst)
231                              full_output, rtol, atol, tcrit, h0, hmax, hmin,
232                              ixpr, mxstep, mxhnil, mxordn, mxords,
--> 233                              int(bool(tfirst)))
234     if output[-1] < 0:
235         warning_msg = _msgs[output[-1]] + " Run with full_output = 1 to get quantitative information."
ValueError: setting an array element with a sequence.

Answer 1

问题就在这里：

IV = [x0, y0, vx0, vy0] 

odeint调用将输出维度确定为初始值的维度乘以时间的维度。在这里，您已经将初始值设置为大小为4，即使x0, y0, vx0, vy0本身就是数组，IV也不会连接数组，它只是将它们放在长度为4的列表中。还要注意，您的IV需要与您的函数预期的形式相同，即按照的顺序和形状

(x_sun, y_sun, vx_sun, vy_sun,
x_merc, y_merc, vx_merc, vy_merc, 
x_ven, y_ven, vx_ven, vy_ven,
x_earth, y_earth, vx_earth, vy_earth,
x_mars, y_mars, vx_mars, vy_mars)

改变你的IV和变量，使其具有正确的形状，并将它们连接起来，而不是将它们放在列表中，你的问题就会消失。

x0 = np.array([0, 46.002, 107.476, 147.092, 206.617])*10**9             #[m]
y0 = np.full(x0.shape, 0.)                               
vx0= np.full(x0.shape, 0.)      
vy0 = 1000*(G*const.M_sun.value*(1/1e6)*((2/x0)-(1/semimajor)))**0.5 #[m/s]
IV = np.vstack([x0, y0, vx0, vy0]).reshape(-1)