Im trying to plot the orbit of jupiter around the sun (as well as 2 astroid clusters in the Lagrange points) using integrate.solve_ivp, but when i plot a graph of x position and y, I'm getting a spiral, rather than a stable orbit. Can anyone help?
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()
These are typical symptoms of a wrong numerical method or wrong parameters to the method.
Reading the documentation, you can use several methods. For the default "RK45" I got what you described. However, using
scipy.integrate.solve_ivp(...,method="RK23",...)
I got nice ellipses.