3 Body Problem 输出尖球而不是轨道路径

Question

我正在尝试用 solve_ivp 及其 runge kutta sim 解决 3 体问题，但它输出的不是一个漂亮的轨道路径，而是一个尖刺的死亡球。我试过各种改变步长和步长，我不知道为什么图表如此尖锐，这对我来说毫无意义。

我现在已经按照建议实现了速度，但我可能做错了

我做错了什么？

更新代码：

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

R = 150000000 #radius from centre of mass to stars orbit
#G = 1/(4*np.pi*np.pi) #Gravitational constant in AU^3/solar mass * years^2
G = 6.67e-11
M = 5e30 #mass of the stars assumed equal mass in solar mass
Omega = np.sqrt(G*M/R**3.0) #inverse of the orbital period of the stars
t = np.arange(0, 1000, 1)
x = 200000000
y = 200000000
vx0 = -0.0003
vy0 = 0.0003

X1 = R*np.cos(Omega*t)
X2 = -R*np.cos(Omega*t)
Y1 = R*np.sin(Omega*t)
Y2 = -R*np.sin(Omega*t) #cartesian coordinates of both stars 1 and 2
r1 = np.sqrt((x-X1)**2.0+(y-Y1)**2.0) #distance from planet to star 1 or 2
r2 = np.sqrt((x-X2)**2.0+(y-Y2)**2.0)
xacc = -G*M*((1/r1**2.0)*((x-X1)/r1)+(1/r2**2.0)*((x-X2)/r2))
yacc = -G*M*((1/r1**2.0)*((y-Y1)/r1)+(1/r2**2.0)*((y-Y2)/r2)) #x double dot and y double dot equations of motions

#when t = 0 we get the initial contditions

r1_0 = np.sqrt((x-R)**2.0+(y-0)**2.0)
r2_0 = np.sqrt((x+R)**2.0+(y+0)**2.0)
xacc0 = -G*M*((1/r1_0**2.0)*((x-R)/r1_0)+(1/r2_0**2.0)*((x+R)/r2_0))
yacc0 = -G*M*((1/r1_0**2.0)*((y-0)/r1_0)+(1/r2_0**2.0)*((y+0)/r2_0))

#inputs for runge-kutta algorithm
tp = Omega*t
r1p = r1/R
r2p = r2/R
xp = x/R
yp = y/R
X1p = X1/R
X2p = X2/R
Y1p = Y1/R
Y2p = Y2/R

#4 1st ode 

#vx = dx/dt
#vy = dy/dt

#dvxp/dtp = -(((xp-X1p)/r1p**3.0)+((xp-X2p)/r2p**3.0))
#dvyp/dtp = -(((yp-Y1p)/r1p**3.0)+((yp-Y2p)/r2p**3.0))

epsilon = x*np.cos(Omega*t)+y*np.sin(Omega*t)
nave = -x*np.sin(Omega*t)+y*np.cos(Omega*t)


# =============================================================================
# def dxdt(x, t):
#     return vx
# 
# def dydt(y, t):
#     return vy
# =============================================================================

def dvdt(t, state):
    xp, yp = state
    X1p = np.cos(Omega*t)
    X2p = -np.cos(Omega*t)
    Y1p = np.sin(Omega*t)
    Y2p = -np.sin(Omega*t)
    r1p = np.sqrt((xp-X1p)**2.0+(yp-Y1p)**2.0)
    r2p = np.sqrt((xp-X2p)**2.0+(yp-Y2p)**2.0)
    return (-(((xp-X1p)/(r1p**3.0))+((xp-X2p)/(r2p**3.0))),-(((yp-Y1p)/(r1p**3.0))+((yp-Y2p)/(r2p**3.0)))) 
    
def vel(t, state): 
    xp, yp, xv, yv = state
    return (np.concatenate([[xv, yv], dvdt(t, [xp, yp]) ]))

p = (R, G, M, Omega)
initial_state = [xp, yp, vx0, vy0]

t_span = (0.0, 1000) #1000 years

result_solve_ivp_dvdt = solve_ivp(vel, t_span, initial_state, atol=0.1) #Runge Kutta
fig = plt.figure()
plt.plot(result_solve_ivp_dvdt.y[0,:], result_solve_ivp_dvdt.y[1,:])
plt.plot(X1p, Y1p)
plt.plot(X2p, Y2p)

输出：绿色是星星图，蓝色是速度 Km and seconds Years, AU and Solar Masses

Answer 1

你已经给出了等式

dv/dt = a(x)

但是后来你用加速度，速度的导数，作为位置的导数。这在物理上是错误的。

你应该传递函数

lambda t, xv: np.concantenate([xv[2:], dvdt(xv[:2]) ])

解算器，除了位置分量外，还包含速度分量的合适初始状态。

在固定轨道的2星系统中，恒星之间的距离为1。这个距离，而不是到中心的0.5距离，应该进入angular速度的计算。

z_1 = 0.5*exp(2*pi*i*t),  z_2 = -z_1  ==>  z_1-z_2=2*z_1, abs(z_1-z_2)=1

z_1'' = -GM * (z_1-z_2)/abs(z_1-z_2)^3

-0.5*4*pi^2 = -GM  or  GM = 2*pi^2

现在将一颗卫星插入半径为 R 的圆半径中，就好像只有一个中心质量 2M 静止在原点

z_3 = R*exp(i*w*t)

z_3'' = -2GM * z_3/abs(z_3)^3

R^3*w^2=2GM

position (R,0), velocity (0,w*R)=(0,sqrt(2GM/R))

在python代码中

GM = 2*np.pi**2
R = 1.9

def kepler(t,u):
    z1 = 0.5*np.exp(2j*np.pi*t)
    z3 = u[0]+1j*u[1]
    a = -GM*((z3-z1)/abs(z3-z1)**3+(z3+z1)/abs(z3+z1)**3)
    return [u[2],u[3],a.real,a.imag]

res = solve_ivp(kepler,(0,17),[R,0,0,2*np.pi*(1/R)**0.5], atol=1e-8, rtol=1e-11)
print(res.message)

这给出了

的轨迹图

双星系统对卫星的影响是一系列连续的 swing-by 机动，加速 angular 速度直到达到逃逸速度。在 R=1.5 或更小的情况下，这发生在卫星与最近的恒星第一次近距离接触时，因此卫星会立即从系统中弹出。

Never-the-less，仍然可以获得“spiky-ball”个轨道。在上面的代码中设置 R=1.6，具有更严格的误差容限并集成到 t=27 给出轨迹

3 Body Problem 输出尖球而不是轨道路径

3 Body Problem Outputs a spikey ball rather than an orbital path

python

plot

physics

runge-kutta