问题描述
我正在尝试使用 Runge-Kutta 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)
plt.plot(sol.t,sol.y[2],color='g')
plt.subplot(2,3)
plt.plot(sol.t,sol.y[4],color='r')
plt.show()
解决方法
使用给定的 ODE 函数，您正在求解系统的第一个组件
xdot = x
ydot = y
具有众所周知的指数解决方案。由于指数因子与两个解的长度相同，xy-plot 将沿一条穿过原点的线移动。
解决方案当然是用 zdot[0:2] 填充 dx,dy，用 zdot[2:4] 或 ax,ay 填充 ddx,ddy 或者你想命名加速度。那么初始状态也只有 4 个组件。或者您需要将位置和速度视为 3 维。
您需要将单位放在常量中，并注意所有使用相同的单位。引用的 G 在 m^3/kg/s^2 中，因此您定义的任何 M 都在 kg 中，任何长度都在 m 中，任何速度在 m/s 中}}。在这种情况下，您的常量可能会显得小得可笑。
不管代码中的注释怎么说，都不会有神奇的转换。您需要使用实际的转换计算来获得真实的数字。例如使用数字
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
这给出了一个非常薄的椭圆，这并不令人惊讶，因为初始速度几乎直接指向重心。
| 标记了半天步长的轨道变体，长度单位为 AU |
|---|
 |
如果交换速度分量，就会得到
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
这样比较平衡。