轨道力学的scipy.integrate.solve_ivp

Question

我目前遇到的问题是使用黑盒解算器，目前它没有给出我应该得到的结果。下面是我用来尝试和解决这个问题的代码，开头的大部分时间都是以地球为例。

import numpy as np
from scipy import integrate
from matplotlib import pyplot

m = 5.972*10e24
M = 1.989*10e30
G = 6.67430*10e-11
k = G*M*m
y = np.array([1.495979*10e12,20*2*np.pi/360])
z = np.array([1.07*10e8, 2*10e-7])
r_0, phi_0 = y[0], y[1]
r_dot_0, phi_dot_0 = z[0], z[1]
l = m*r_0**2*phi_dot_0
V=-k/r_0
E=1/2*m*(r_dot_0**2+r_0**2*phi_dot_0**2)+(1/2)*(l**2)/(m*r_0**2)+V
assert(E>0)
r_min = (-k/E-np.sqrt(k**2/E**2+2*l**2/(m*E)))/2
r_max = (-k/E+np.sqrt(k**2/E**2+2*l**2/(m*E)))/2
r= np.linspace(r_min,r_max,10000)
def f_1(r,phi):
    drdt=np.sqrt((2/m)*np.abs((E-(-k/r)-(l**2)/(2*m*r**2))))
    d0dt = phi_0 + l/(m*r**2)
    return drdt, d0dt
def theta(f_1,r_min,r_max,r_0,phi_0):
    return integrate.solve_ivp(f_1,[r_min,r_max],[r_0,phi_0])

我得到的结果是

message: 'The solver successfully reached the end of the integration interval.'
     nfev: 188
     njev: 0
      nlu: 0
      sol: None
   status: 0
  success: True
        t: array([-4.17982328e+11, -4.17982328e+11, -4.17982328e+11, -4.17982328e+11,
       -4.17982328e+11, -4.17982327e+11, -4.17982317e+11, -4.17982215e+11,
       -4.17981721e+11, -4.17979977e+11, -4.17974974e+11, -4.17961411e+11,
       -4.17924396e+11, -4.17822864e+11, -4.17543973e+11, -4.16777612e+11,
       -4.14671109e+11, -4.08877180e+11, -3.92913104e+11, -3.48712226e+11,
       -2.24166605e+11, -8.84161983e+10, -3.66420882e+10,  4.01903359e+10,
        1.28704958e+11,  2.83822588e+11,  4.17982305e+11])
 t_events: None
        y: array([[1.49597900e+13, 1.49597900e+13, 1.49597900e+13, 1.49597901e+13,
        1.49597920e+13, 1.49598540e+13, 1.49618143e+13, 1.50193097e+13,
        1.56994667e+13, 2.05914847e+13, 4.61047292e+13, 1.64356184e+14,
        7.03558012e+14, 3.16004985e+15, 1.43535511e+16, 6.53961403e+16,
        2.98620094e+17, 1.37024180e+18, 6.37281115e+18, 3.08691143e+19,
        1.74995579e+20, 5.63096052e+20, 9.53269689e+20, 2.83709416e+21,
        3.34912831e+21, 3.65873923e+21, 3.74798929e+21],
       [3.49065850e-01, 3.86709595e-01, 7.63147037e-01, 4.52752146e+00,
        4.21712657e+01, 4.18608708e+02, 4.18298313e+03, 3.98437512e+04,
        2.13698548e+05, 8.26969054e+05, 2.58606034e+06, 7.35526027e+06,
        2.03705358e+07, 5.60724135e+07, 1.54139230e+08, 4.23620702e+08,
        1.16438445e+09, 3.20214248e+09, 8.81913334e+09, 2.43925440e+10,
        6.85803126e+10, 1.19038660e+11, 1.44275525e+11, 2.85105393e+11,
        3.23736253e+11, 3.79785425e+11, 4.27122185e+11]])

其中t表示径向位置，y表示时间，y1表示角位置。R给出了“合理的”结果，但我希望精度要大得多，而y是可怕得可笑的。

此外，这段代码并不是我想要的理想代码；在f_1(r,phi)中，我应该有更多类似于dtdr = 1/np.sqrt((2/m)*(E-(-k/r)-(l**2)/(2*m*r**2)))的东西，并且在原始代码中甚至不应该有np.abs，然而，黑盒求解器甚至不能在没有它的情况下运行。任何关于我在程序中的愚蠢之处的提示都将是一个很大的帮助。

Answer 1

开始集成系统并不容易，但是一旦你有了合适的模板，它就会变得容易起来。您在f_1中使用的参数不正确。第一个条目是时间，第二个条目是要集成的数组：

def f_1(t, arg):
    r, phi = arg
    drdt=np.sqrt((2/m)*(E - (l**2 / (2*m*r**2) - k/r)))
    d0dt = l/(m*r**2)
    return [drdt, d0dt]

要获得更多的可定制性，您可以使用ode模块。使用此模块进行集成的适当方式可能是：

r = integrate.ode(f_1)
r.set_integrator('dopri5')
r.set_initial_value([r_0, phi_0], 0)
t1 = 356 * 24 * 60 * 60 # one year
dt = 24 * 60 * 60 # steps of one day
Y = [r.y]
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    Y.append(r.y)
Y = np.array(Y)

使用修改后的f_1版本，您的集成方式也很好。请重新检查E和d0dt的定义。