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 km)，而 dt == 1. 表示一年。 (由于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()