问题描述
我正在尝试使用integral.solve_ivp绘制木星在太阳周围的轨道(以及在拉格朗日点上的2个恒星簇),但是当我绘制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[2],y[6],y[7],y[8])
Vjg = vectorise(y[3],y[8])
Vst = vectorise(y[0],y[9],y[10],y[11])
Vjt = vectorise(y[3],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,#initial values of y[0-2]; sun position
0.0,5.2,#initial values of y[3-5] jupiter position
-5.2*np.sqrt(3/4),5.2*0.5,#initial values of y[6-8]; greek position
5.2*np.sqrt(3/4),#initial values of y[9-11]; trojan position
0,#initial values of y[12-14]; sun veLocity
2.755,#initial values of y[15-17]; jupiter veLocity in AU per year
2.755*0.5,2.755*np.sqrt(3/4),#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()
解决方法
这些是数值方法错误或方法参数错误的典型症状。
阅读the documentation时,可以使用多种方法。对于默认的"RK45"
,我得到了您所描述的内容。但是,使用
scipy.integrate.solve_ivp(...,method="RK23",...)
我的椭圆很好。