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

Question

我试图使用 integrate.solve_ivp 绘制木星围绕太阳的轨道（以及拉格朗日点中的 2 个星团），但是当我绘制 x 位置和 y 的图形时，我得到的是螺旋线，而不是稳定的轨道。谁能帮忙？import numpy as npimport scipy.integrateimport matplotlib.pyplot as plt#takes the position of two masses and outputs the vector difference, and the modulusdef 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)

Answer 1

这些是数值方法错误或方法参数错误的典型症状。

默认情况下，"RK45"我得到了你所描述的。然而，使用

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

我得到了很好的省略号。