我正在尝试使用 python 求解测地线轨道方程组。它们是耦合的普通方程。我尝试了不同的方法，但它们都产生了错误的形状（绘制 r 和 phi 时，形状应该是一些周期函数）。关于如何做到这一点有什么想法吗？这是我的常量
G = 4.30091252525 * (pow(10, -3)) #Gravitational constant in (parsec*km^2)/(Ms*sec^2)
c = 0.0020053761 #speed of light , AU/sec
M = 170000 #mass of the central body, in solar masses
m = 10 #mass of the orbiting body, in solar masses
rs = 2 * G * M / pow(c, 2) #Schwarzschild radius
Lz= 0.000024 #Angular momemntum
h = Lz / m #Just the constant in equation
E= 1.715488e-007 #energy
初始条件为：
Y(0) = rs …