#python #orbit
#python #орбита
Вопрос:
Я пытаюсь построить орбиту Юпитера вокруг солнца (а также 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()
Комментарии:
1. Похоже, что проблема точности численного метода..
Ответ №1:
Это типичные симптомы неправильного численного метода или неправильных параметров метода.
Читая документацию, вы можете использовать несколько методов. По умолчанию "RK45"
я получил то, что вы описали. Однако, используя
scipy.integrate.solve_ivp(...,method="RK23",...)
Я получил красивые эллипсы.
Комментарии:
1. Большое вам спасибо! Это намного, намного лучше. Через некоторое время он все еще вращается по спирали, если я запущу его, скажем, на 100000 лет, есть ли метод, который позволил бы ему оставаться стабильным дольше / неопределенно долго?
2. Нет. Теоретически вы можете получить результат численного метода с бесконечно малой ошибкой (например, используя бесконечно малый временной шаг), но 1) это было бы слишком дорого в вычислительном отношении и 2) все равно были бы ошибки компьютерного округления. Если у вас есть два тела, вы можете использовать аналитическое решение. Если у вас более двух тел, вы обречены на не совсем точное численное решение ..
3. Я читал о методах RK, и кажется, что RK45 — более точный метод, вы знаете, почему в этом случае RK23 сработал лучше?
4. Я не математик и знаю только приблизительно предысторию методов. Поскольку на самом деле это не связано с программированием, а является математическим вопросом, может быть, вы можете попробовать задать на math stack exchange или math overflow или аналогичном?