Моделирование орбиты Земли с использованием симплекси c Эйлера в Python

Question

Я пытаюсь смоделировать движение Земли вокруг Солнца. (это задача, которую я пытаюсь выполнить) Это то, что я придумал до сих пор

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 является совершенно правильной ? С этим кодом я думаю, что орбита всегда будет круговой. Любые советы, чтобы улучшить это будет высоко ценится

Answer 1

Ваш текст немного неправильный. Symplecti c Эйлер - это метод порядка 1, а Штормер-Верле - порядок 2. Можно было бы получить правильный метод с помощью схемы скачкообразного перегиба, где скорости измеряются в средних точках временных интервалов. На практике это означает, что начальная скорость должна быть смещена к скорости в момент времени 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

Ваша реализация силы неверна, поскольку она использует постоянный радиус и не фактический радиус PL anet в настоящее время. Изменить на

#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