python - 使用 Python 的耦合微分方程
问题描述
我正在尝试使用 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)
解决方案
似乎您正在查看在施瓦西重力中移动的物体的测地线方程的不变平面中的空间坐标。
人们可以使用许多不同的方法,尽可能多地保留模型的基本几何结构,例如辛几何积分器或微扰理论。正如 Lutz Lehmann 在评论中指出的那样,“solve_ivp”的默认方法默认使用 Dormand-Prince (4)5 步进器,该步进器利用外推模式,即 5 阶步进,步长选择由阶数 4 的误差估计。
警告:您的初始条件Y
等于 Schwarzschild 的半径,因此这些方程可能会失败或需要特殊处理(尤其是方程的时间分量,您没有包括在此处!)这可能是您必须切换到不同的坐标,从而删除偶数地平线上的奇点。此外,解可能不是周期曲线,而是准周期曲线,因此它们可能不会很好地接近。
对于一个快速而肮脏的处理,但可能是一个相当准确的处理,我会区分第一个方程
(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)
仔细检查方程,错误和不准确是可能的。
推荐阅读
- angular - 带有 Yarn 2 (dlx) 的 Angular CLI
- java - 用equals方法比较组合对象的问题
- javascript - JavaScript + discord.js:如何使用一个对象的属性在另一个对象中同时存储在同一个对象中?
- prisma - Prisma2:在创建函数中使用相关表中的用户 ID
- python - 如何在 django 中将满足特定条件的行从一个表复制到另一个表?
- python - 删除外键值
- asp.net-core - .Net Core 3.1 中是否有任何内置函数/方法可以返回,IActionResult/ActionResult 而不是 HttpResponseMessage
- node.js - 获取 ejs 文件中未定义的错误消息(在 express 中使用 flash 时)
- regex - 查找没有指定字符串的行并删除空行
- javascript - Jquery 转换为 Javascript 导致标题未定义?如何正确转换它?