Похоже, что вы смотрите на пространственные координаты в инвариантной плоскости уравнений геодезии c объекта, движущегося в гравитации Шварцшильда.
Можно использовать много разных методов, которые сохраняют как можно большую часть базовой геометрии c структуры модели, например, симплект c геометрии c интеграторов или теории возмущений. Тем не менее, я подозреваю, что метод default_ivp () по умолчанию - это Runge-Kutta 4-го порядка, который является довольно хорошим, прямым и простым методом.
Warnng: ваше начальное условие для Y равно радиусу Шварцшильда, поэтому эти уравнения могут не работать или требовать особой обработки (особенно временная составляющая уравнений, которую вы здесь не включили!) Возможно, вы приходится переключаться на разные координаты, которые убирают сингулярность на ровном горизонте. Более того, решения могут быть не периодическими c кривыми, а квазипериодическими c, поэтому они могут не очень хорошо закрываться.
Для быстрой и грязной обработки, но, возможно, довольно точной, я бы дифференцировал первое уравнение
(dr / dtau)^2 = (E2_mc2 - c2) + (2*GM)/r - (h^2)/(r^2) + (r_schw*h^2)/(r^3)
относительно подходящего времени tau, затем отменил бы первое производная dr / dtau по r с обеих сторон, и заканчивается уравнение со второй производной для радиуса r слева. Затем превратите это второе уравнение в производную в пару уравнений первой производной для r и скорости его изменения v, то есть
dphi / dtau = h / (r^2)
dr / dtau = v
dv / dtau = - GM / (r^2) + h^2 / (r^3) - 3*r_schw*(h^2) / (2*(r^4))
, и рассчитайте по исходному уравнению для r и его первой производной dr / dtau начальное значение для скорости изменения v = dr / dtau, то есть я бы решил для v уравнения с r=r0:
(v0)^2 = (E2_mc2 - c2) + (2*GM)/r0 - (h^2)/(r0^2) + (r_schw*h^2)/(r0^3)
Может быть, какой-то код python, подобный этому, может работа:
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
#from ode_helpers import state_plotter
# u = [phi, Y, V, t] or if time is excluded
# u = [phi, Y, V]
def f(tau, u, param):
E2_mc2, c2, GM, h, r_schw = param
Y = u[1]
f_phi = h / (Y**2)
f_Y = u[2] # this is the dr / dt auxiliary equation
f_V = - GM / (Y**2) + h**2 / (Y**3) - 3*r_schw*(h**2) / (2*(Y**4))
#f_time = (E2_mc2 * Y) / (Y - r_schw) # this is the equation of the time coordinate
return [f_phi, f_Y, f_V] # or [f_phi, f_Y, f_V, f_time]
# from the initial value for r = Y0 and given energy E,
# calculate the initial rate of change dr / dtau = V0
def ivp(Y0, param, sign):
E2_mc2, c2, GM, h, r_schw = param
V0 = math.sqrt((E2_mc2 - c2) + (2*GM)/Y0 - (h**2)/(Y0**2) + (r_schw*h**2)/(Y0**3))
return sign*V0
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
Lz= 0.000024 #Angular momemntum
h = Lz / m #Just the constant in equation
E= 1.715488e-007 #energy
c2 = c**2
E2_mc2 = (E**2) / (c2*m**2)
GM = G*M
r_schw = 2*GM / c2
param = [E2_mc2, c2, GM, h, r_schw]
Y0 = r_schw
sign = 1 # or -1
V0 = ivp(Y0, param, sign)
tau_span = np.linspace(1, 1000, num=1000)
u0 = [math.pi, Y0, V0]
sol = solve_ivp(lambda tau, u: f(tau, u, param), [1, 1000], u0, t_eval=tau_span)
Перепроверьте уравнения, возможны ошибки и неточности.