使用Integrated.solve_ivp在python中绘制轨道

如何解决使用Integrated.solve_ivp在python中绘制轨道

我正在尝试使用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",...)

我的椭圆很好。