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

Question

我试图使用 integrate.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[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",...)我得到了很好的省略号。