使用rk4将行星轨道显示为线性图

Question

我试图用龙格库塔4方法模拟行星绕恒星的轨道。在与导师交谈后，我的代码应该是正确的。然而，我并没有生成预期的2D轨道图，而是一个线性图。这是我第一次用solve_ivp来解二阶微分。有人能解释为什么我的阴谋是错误的吗

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# %% Define derivative function
def f(t, z):
    
    x = z[0] # Position x
    y = z[1] # Position y
    dx = z[2] # Velocity x
    dy = z[3] # Velocity y
    
    G = 6.674 * 10**-11 # Gravitational constant
    M = 2 # Mass of binary stars in solar masses
    c = 2*G*M 
    r = np.sqrt(y**2 + x**2) # Distance of planet from stars
    
    zdot = np.empty(6) # Array for integration solutions
    zdot[0] = x
    zdot[1] = y
    zdot[2] = dx # Velocity x
    zdot[3] = dy #Velocity y
    zdot[4] = (-c/(r**3))*(x) # Acceleration x
    zdot[5] = (-c/(r**3))*(y) # Acceleration y
    
    return zdot
# %% Define time spans, initial values, and constants
tspan = np.linspace(0., 10000., 100000000)
xy0 = [0.03, -0.2, 0.008, 0.046, 0.2, 0.3] # Initial positions x,y in R and velocities
# %% Solve differential equation
sol = solve_ivp(lambda t, z: f(t, z), [tspan[0], tspan[-1]], xy0, t_eval=tspan)
# %% Plot 
#plot
plt.grid()
plt.subplot(2, 2, 1)
plt.plot(sol.y[0],sol.y[1], color='b')
plt.subplot(2, 2, 2)
plt.plot(sol.t,sol.y[2], color='g')
plt.subplot(2, 2, 3)
plt.plot(sol.t,sol.y[4], color='r')
plt.show()

Answer 1

使用给定的ODE函数，您将在系统的第一个组件中求解

xdot = x
ydot = y

它有著名的指数解。由于两种解决方案的指数因子相同，xy图将沿穿过原点的直线移动

当然，解决方案是用dx,dy填充zdot[0:2]，用ax,ay或ddx,ddy填充zdot[2:4]，或者你想如何命名加速度的组成部分。那么初始状态也只有4个分量。或者你需要把位置和速度做成三维的

---

你需要把单位放在常数里，注意所有的单位都是一样的G引用的是{}，因此您定义的任何{}都是{}，任何长度都是{}，任何速度都是{}。在这种情况下，您的常量可能看起来小得可笑

不管代码中的注释是什么，都不会有神奇的转换。您需要使用实际的转换计算来获得真实的数字。例如使用数字

G       = 6.67408e-11 # m^3 s^-2 kg^-1
AU      = 149.597e9 # m
Msun    = 1.988435e30 # kg
hour = 60*60 # seconds in an hour
day = hour * 24 # seconds in one day
year = 365.25*day # seconds in a year (not very astronomical)

人们可以猜测，对于一个由两颗质量相等的恒星组成的敏感双星系统，一个

M = 2*Msun # now actually 2 sun masses
x0 = 0.03*AU
y0 = -0.2*AU
vx0 = 0.008*AU/day
vy0 = 0.046*AU/day

由于位置仅以AU为单位，因此速度也可以为AU/hour。通过https://math.stackexchange.com/questions/4033996/developing-keplers-first-law和Cannot get RK4 to solve for position of orbiting body in Python，半径为R=0.2AU的圆形轨道围绕2*M组合质量的速度为

sqrt(2*M*G/R)=sqrt(4*Msun*G/(0.2*AU)) = 0.00320 * AU/hour = 0.07693 AU/day

这是。。。如果给定的速度实际上在AU/day范围内，则不太合理。调用https://math.stackexchange.com/questions/4050575/application-of-the-principle-of-conservation中的计算来计算开普勒椭圆是否合理

r0 = (x0**2+y0**2)**0.5
dotr0 = (x0*vx0+y0*vy0)/r0
L = x0*vy0-y0*vx0           # r^2*dotphi = L constant, L^2 = G*M_center*R
dotphi0 = L/r0**2

R = L**2/(G*2*M)
wx = R/r0-1; wy = -dotr0*(R/(G*2*M))**0.5
E = (wx*wx+wy*wy)**0.5; psi = m.atan2(wy,wx)
print(f"half-axis R={R/AU} AU, eccentr. E={E}, init. angle psi={psi}")
print(f"min. rad. = {R/(1+E)/AU} AU, max. rad. = {R/(1-E)/AU} AU")

返回

half-axis R=0.00750258 AU, eccentr. E=0.96934113, init. angle psi=3.02626659
min. rad. = 0.00380969 AU, max. rad. = 0.24471174 AU

这给出了一个非常薄的椭圆，这并不令人惊讶，因为初始速度几乎直接指向重心

^{tb1}$

如果交换速度分量，则会得到

half-axis R=0.07528741 AU, eccentr. E=0.62778767, init. angle psi=3.12777251
min. rad. = 0.04625137 AU, max. rad. = 0.20227006 AU

这是比较平衡的