1

我正在尝试模拟地球围绕太阳的运动。（这是我正在尝试做的任务） 这是我到目前为止想出的

import numpy as np
import matplotlib.pyplot as plt
#Set parameters:
N = 365                     # Earth days in a year
dt = 1.00/N                 # Time Step: Fractions of a year - 1 Earth day (i.e. 1/365)
mu = 4 * np.pi**2           # mu=4pi^2 is the Gravitational Parameter: mu = GM where G=6.67e-11 is the Universal Gravitational Constant and M is the mass of the body
rEar = 1

#Create an array, for all variables, of size N with all entries equal to zero:
xEar = np.zeros((N,))
yEar = np.zeros((N,))
vxEar = np.zeros((N,))
vyEar = np.zeros((N,))

# Initial Conditions:
xEar[0] = rEar                   # (x0 = r, y0 = 0) in AU
vyEar[0] = np.sqrt(mu/rEar)      # (vx0 = 0, v_y0 = sqrt(mu/r)) AU/yr

#Implement Verlet Algorithm:
for k in range(0,N-1):
    vxEar[k+1] = vxEar[k] - (mu*xEar[k]) / (rEar**3)*dt
    xEar [k+1] = xEar[k] + vxEar[k+1]*dt
    vyEar[k+1] = vyEar[k] - (mu*yEar[k]) / (rEar**3)*dt
    yEar [k+1] = yEar[k] + vyEar[k+1]*dt

#Plot:
plt.plot(xEar, yEar, 'go')
plt.title ('Circular Orbit of Earth')
plt.xlabel ('x')
plt.ylabel ('y')
plt.axis('equal')
plt.show()

但我认为我的 verlet 实现不太正确？使用此代码，我认为轨道将始终是圆形的。任何改善这一点的提示将不胜感激

4

1 回答 1

1

你的文字有点错误。Symplectic Euler 是 1 阶方法，而 Stormer-Verlet 是 2 阶方法。人们可以通过实施跨越式 Verlet 方案来弥合方法之间的差距，其中速度取自时间间隔的中点。在实践中，这意味着必须将初始速度转移到 time t0-dt/2，其他一切都保持不变。

# Initial Conditions:
xEar[0] = rEar                   # (x0 = r, y0 = 0) in AU
vyEar[0] = np.sqrt(mu/rEar)      # (vx0 = 0, v_y0 = sqrt(mu/r)) AU/yr,
vxEar[0] = 0 - mu / (rEar**2)*(-0.5*dt) # corrected for time -dt/2

您对力的实施是错误的，因为它使用恒定半径而不是当前行星的实际半径。改成

#Implement Verlet Algorithm:
for k in range(0,N-1):
    rEar = (xEar[k]**2+yEar[k]**2)**0.5
    vxEar[k+1] = vxEar[k] - (mu*xEar[k]) / (rEar**3)*dt
    xEar [k+1] = xEar[k] + vxEar[k+1]*dt
    vyEar[k+1] = vyEar[k] - (mu*yEar[k]) / (rEar**3)*dt
    yEar [k+1] = yEar[k] + vyEar[k+1]*dt
于 2020-04-02T09:15:29.513 回答