我目前遇到的问题是使用黑盒求解器，目前没有给出我应该得到的结果。下面是我用来尝试解决这个问题的代码，开头的大部分内容都以地球为例。
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[0] 是时间，y[1] 是角位置。 r 给出了"合理"的结果，但我希望精度更高，而 y 是可笑的可怕。
另外，这段代码对于我想要的并不理想;f_1(r,phi)我应该有更像dtdr = 1/np.sqrt((2/m)*(E-(-k/r)-(l**2)/(2*m*r**2)))的东西，我甚至不应该在原始代码中有一个np.abs，但是，没有它，黑盒求解器甚至不会运行。任何指向我在程序中有点愚蠢的地方都会有很大帮助。
开始集成系统并不容易，但一旦您拥有合适的模板，它就会变得更容易。您使用的参数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。