Блог

Я пытаюсь построить орбиту Юпитера вокруг солнца (а также 2 скопления астроидов в точках Лагранжа) с помощью integrate.solve_ivp, но когда я строю график положения x и y, я получаю спираль, а не стабильную орбиту. Кто-нибудь может помочь?

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

#takes the position of two masses and outputs the vector difference, and the modulus
def vectorise(x1,y1,z1,x2,y2,z2):
    v = np.array([x2-x1,y2-y1,z2-z1])
    return v, np.linalg.norm(v)

def derivatives(t, y):
    G =4*np.pi**2
    Mj = 0.001
    Ms = 1

    #Vij gives the vector pointing from i to j (leading to force on j from i)
    Vjs = vectorise(y[3],y[4],y[5],y[0],y[1],y[2])
    Vsg = vectorise(y[0],y[1],y[2],y[6],y[7],y[8])
    Vjg = vectorise(y[3],y[4],y[5],y[6],y[7],y[8])
    Vst = vectorise(y[0],y[1],y[2],y[9],y[10],y[11])
    Vjt = vectorise(y[3],y[4],y[5],y[9],y[10],y[11])

    return [y[12],y[13],y[14],#first differentials of sun position
    y[15],y[16],y[17],#first differentials of Jupiter position
    y[18],y[19],y[20],#first differentials of Greek position
    y[21],y[22],y[23], #first differentials of Trojan position
    -G*Mj*1/(Vjs[1]**3) *Vjs[0][0], #second differentail of y[12] (sun x)
    -G*Mj*1/(Vjs[1]**3) *Vjs[0][1], #second differentail of y[13] (sun y)
    -G*Mj*1/(Vjs[1]**3) *Vjs[0][2], #second differentail of y[14] (sun z)
    G*Ms*1/(Vjs[1]**3) *Vjs[0][0], #second differentail of y[15] (Jupiter x)
    G*Ms*1/(Vjs[1]**3) *Vjs[0][1], #second differentail of y[16] (Jupiter y)
    G*Ms*1/(Vjs[1]**3) *Vjs[0][2], #second differentail of y[17] (Jupiter z)
    -G*(Ms*1/(Vsg[1]**3) * Vsg[0][0]   Mj*1/(Vjg[1]**3) * Vjg[0][0]), #second differentail of y[18] (Greek x)
    -G*(Ms*1/(Vsg[1]**3) * Vsg[0][1]   Mj*1/(Vjg[1]**3) * Vjg[0][1]), #second differentail of y[19] (Greek y)
    -G*(Ms*1/(Vsg[1]**3) * Vsg[0][2]   Mj*1/(Vjg[1]**3) * Vjg[0][2]), #second differentail of y[20] (Greek z)
    -G*(Ms*1/(Vst[1]**3) * Vst[0][0]   Mj*1/(Vjt[1]**3) * Vjt[0][0]), #second differentail of y[21] (Trojan x)
    -G*(Ms*1/(Vst[1]**3) * Vst[0][1]   Mj*1/(Vjt[1]**3) * Vjt[0][1]), #second differentail of y[22] (Trojan y)
    -G*(Ms*1/(Vst[1]**3) * Vst[0][2]   Mj*1/(Vjt[1]**3) * Vjt[0][2])] #second differentail of y[23] (Trojan z)

def solver():
    solution = scipy.integrate.solve_ivp(
        fun = derivatives,
        t_span=(0,118.6),
        # initial velocities are given in sun frame, and are in AU/year, distances in AU (1AU = 1.496 *10^11 m)
        # jupiter has average orbital velocity of 13.07 km/s = 2.755 au/year
        y0 = (0.0, 0.0, 0.0, #initial values of y[0-2]; sun position
        0.0, 5.2, 0.0, #initial values of y[3-5] jupiter position
        -5.2*np.sqrt(3/4), 5.2*0.5, 0,#initial values of y[6-8]; greek position
        5.2*np.sqrt(3/4), 5.2*0.5 , 0,#initial values of y[9-11]; trojan position
        0,0,0,#initial values of y[12-14]; sun velocity
        2.755, 0, 0,#initial values of y[15-17]; jupiter velocity in AU per year
        2.755*0.5, 2.755*np.sqrt(3/4), 0,#initial values of y[18-20]; greek velocity in AU per year
        2.755*0.5, -2.755*np.sqrt(3/4), 0),#initial values of y[21-23]; trojan velocity in AU per year
        t_eval = np.linspace(0,118.6,  1000),
    )
    return solution

plt.plot(solver().y[3],solver().y[4])
plt.show()
  

Это типичные симптомы неправильного численного метода или неправильных параметров метода.

Читая документацию, вы можете использовать несколько методов. По умолчанию "RK45" я получил то, что вы описали. Однако, используя

 scipy.integrate.solve_ivp(...,method="RK23",...)
  

Я получил красивые эллипсы.