1

我将绘制围绕太阳的地球。因此，该任务分为 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()

谢谢你的帮助 ：）

4

1 回答 1

1

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

  1. 您没有使用正确的测量单位（或者，您没有始终如一地使用它们。）
    据我从您发布的代码中可以看出，与太阳的距离应该以天文单位和时间以一年的分数来测量，所以这r == 1.意味着距离为 1AU（~1.49e8 公里）并且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方程给出。所以你强加（v_x=0 AU/yr，v_y=sqrt(mu/r) AU/yr）。请注意，如果您施加这组初始条件，则不必施加任何其他条件，因为问题已经很好地提出了。（此外，类似的条件在 -loop 中被覆盖，根本不会影响您的计算。）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()
于 2016-02-09T22:58:51.893 回答