Python 开普勒定律绘图

Question

我将绘制围绕太阳的地球。因此，该任务分为 2 个子任务。 在第一个任务中，我将近似认为运动是一个圆圈。

我使用以下代码来获得解决方案，但不知何故程序会编辑一个点而不是几个点。你能帮我解决我的算法吗？

所以我的代码：

npoints = 360    
x= np.zeros((npoints,1))
y= np.zeros((npoints,1))
v_x=np.zeros((npoints,1))
v_y=np.zeros((npoints,1))
r=1
dt=1
x[0]=1.
y[0]=0.
v_x[0]=-1.
v_y[90]=-1.
v_x[180]=1.
v_y[270]=1.
for step in range(0,npoints-1):
    v_x[step+1]=v_x[step]-4*pi**2*x[step]/(r**3)*dt
    x[step+1]=x[step]+v_x[step+1]*dt
    v_y[step+1]=v_y[step]-4*pi**2*y[step]/(r**3)*dt
    y[step+1]=y[step]+v_y[step+1]*dt


plt.plot(x, y)
plt.axis([-100, 100, -100, 100])
plt.ylabel('y-axis')
plt.xlabel('x-axis')
plt.show()  

感谢您的帮助:)

Answer 1

我认为你的代码有两个错误：

1. 您没有使用正确的度量单位（或者，您没有始终如一地使用它们。）
   据我从您发布的代码中可以看出，与太阳的距离应以Astronomical Units 为单位，时间以一年的几分之一为单位，因此r == 1. 表示距离为 1AU（~1.49e8 公里），@ 987654324@ 是一年。 （dt == 1.太大，可以除以npoints：dt = 1./npoints。如果npoints == 360，时间步长为一天。）
   此外，您还需要将重力参数mu 表达为一致的测量单位。使用轨道周期表达式T = 2*pi*sqrt(r**3 / mu) 并加上T=1. 和r=1.，我们得到mu = 4 * pi**2。
2. 您对速度施加了错误的初始条件。
   让我们假设（如您所做的那样）地球的初始位置具有坐标（x=1 AU，y=0 AU）。速度与轨道相切，因此在所选参考系中它只有一个垂直分量 (v_y)，其模块由 Circular Velocity 的等式给出。所以你强加（v_x=0 AU/yr，v_y=sqrt(mu/r) AU/yr）。 请注意，如果您施加这组初始条件，则不必施加任何其他条件，因为问题已经很好地提出了。 （此外，v_y[90]=-1. 之类的条件会在 for 循环中被覆盖，根本不会影响您的计算。）

完整代码如下：

import numpy as np              # please next time include the relevant
import matplotlib.pyplot as plt # `import` statements and variable
pi = np.pi                      # definitions

npoints = 360
r = 1.          # AU
dt = 1./npoints # fractions of a year
mu = 4 * pi**2  # 
x = np.zeros(npoints)
y = np.zeros(npoints)
v_x = np.zeros(npoints)
v_y = np.zeros(npoints)

# Initial Conditions
x[0] = r               # (x0 = r, y0 = 0) AU
v_y[0] = np.sqrt(mu/r) # (v_x0 = 0, v_y0 = sqrt(mu/r)) AU/yr

for step in range(0,npoints-1):
    v_x[step+1]=v_x[step]-4*pi**2*x[step]/(r**3)*dt
    x[step+1]=x[step]+v_x[step+1]*dt
    v_y[step+1]=v_y[step]-4*pi**2*y[step]/(r**3)*dt
    y[step+1]=y[step]+v_y[step+1]*dt

plt.plot(x, y, 'bo')
plt.axis('equal')
plt.show()