使用 Python 计算微分方程组

Question

我正在尝试使用 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
Phi(0) = math.pi

轨道方程

我尝试这样做的方式：

def rhs(t, u):
    Y, phi = u
    dY = np.sqrt((E**2 / (m**2 * c**2) - (1 - rs / Y) * (c**2 + h**2 / Y**2)))
    dphi = L / Y**2
    return [dY, dphi]

Y0 = np.array([rs,math.pi])
sol = solve_ivp(rhs, [1, 1000], Y0, method='Radau', dense_output=True)

Answer 1

看起来您正在查看在史瓦西引力中移动的物体的测地方程的不变平面中的空间坐标。

人们可以使用许多不同的方法，尽可能多地保留模型的基础几何结构，例如辛几何积分器或微扰理论。正如 Lutz Lehmann 在评论中指出的那样，“solve_ivp”的默认方法使用默认的 Dormand-Prince (4)5 步进器，该步进器利用外推模式，即 5 阶步骤，步长选择由4阶步骤的误差估计。

警告：您的初始条件Y等于史瓦西半径，因此这些方程可能会失败或需要特殊处理（尤其是方程的时间部分，您没有将其包含在此处！）您可能必须切换到不同的坐标，从而删除偶地平线上的奇点。此外，解可能不是周期曲线，而是准周期曲线，因此它们可能无法很好地闭合。

对于快速而肮脏的处理，但可能是相当准确的处理，我将对第一个方程进行微分

(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求解：vr=r0

(v0)^2 = (E2_mc2 - c2) + (2*GM)/r0 - (h^2)/(r0^2) + (r_schw*h^2)/(r0^3)

也许像这样的某种Python代码可能会起作用：

(dr / dtau)^2 = (E2_mc2 - c2) + (2*GM)/r - (h^2)/(r^2) + (r_schw*h^2)/(r^3)

仔细检查方程式，可能存在错误和不准确。