Блог

Я создаю орбитальную модель солнечной системы, поэтому мне нужно интегрировать второй закон Ньютона для каждой планеты в виде суммы. Я использую ODEINT и получаю эту ошибку значения, но не знаю, где возникает реальная проблема. Это мой первый пост, поэтому есть проблема с отступом под функцией…

 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.
 

Проблема здесь:

 IV = [x0, y0, vx0, vy0] 
 

odeint Вызов определяет ваши выходные размеры как размеры начального значения по размерам времени. Здесь вы установили для своего начального значения значение четвертого размера, даже если 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)