使用 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

看起来您正在查看 object 在 Schwarzschild 引力中移动的测地线方程的不变平面中的空间坐标。

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

警告：Y 的初始条件等于 Schwarzschild 半径，因此这些方程可能会失败或需要特殊处理（尤其是方程的时间分量，您未在此处包含！）可能是您必须切换到不同的坐标，以消除偶数 horizon 处的奇点。此外，解决方案可能不是周期性曲线，而是 quasi-periodic，因此它们可能不会很好地关闭。

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

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

关于本征时间tau，然后抵消两边关于r的一阶导数dr / dtau，得到二阶导数为左边的半径 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)

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

使用 Python 耦合微分方程

Couple Differential Equations using Python

python

physics

astronomy

differential-equations