Блог

Я пытаюсь построить орбиту Луны вокруг Юпитера, используя гравитационное ускорение. Кажется, я не могу определить, как правильно использовать функцию solve_ivp. Что-то просто не щелкает… Я создал ODE для луны, связанной с Юпитером в начале координат.

 year = np.arange(0,31536000,3600)
G = 6.67408e-11
jupiter_xpos = 0
jupiter_ypos = 0
jupiter_vel = (0,0)
jupiter_mass = 1.89819e27
Io_orbit = 421700000
Io_xpos = -421700000
Io_ypos = 0
Io_xvel = 0
Io_yvel =  -1773400
Io_mass = 8.9319e22
Io = [Io_xpos,Io_xvel,Io_ypos,Io_yvel]


def ode(Moon,t_max):
    #Moon[0,1,2,3]=[x,v_x,y,v_y]
    output=[0,0,0,0]
    R = ((Moon[0]**2   Moon[2]**2)**(3/2))
    output[0]=Moon[1]
    output[1]= -G*jupiter_mass*Moon[0]/R
    output[2]=Moon[3]
    output[3]= -G*jupiter_mass*Moon[2]/R
    return output

#This is where the problem is
sol= solve_ivp(ode,Io,year)


plt.plot(sol[:,0],sol[:,2],0,0,'ro')
plt.axis('equal')
plt.grid(True)
plt.show()
  

Я надеюсь получить такой 2D-орбитальный график…

а также отслеживать и отображать каждое изменение положения и скорости x и y во времени.

Документация для solve_ivp показывает, что параметры

 sol = solve_ivp(ode, [t0,tf], u0, tspan=year, atol = 100, rtol=1e-9)
  

где year=np.arange(t0,tf,hour) . Затем вы находите значения решения sol.t для повторяющихся раз и sol.y для значений.