Моя главная проблема на данный момент: что орбитальный спутник вращается вокруг целевой планеты, хотя он должен находиться на стабильной круговой орбите длиной 320 км.
Я создал четыре разные функции для четырех разных интеграций. Эйлер, Чехарда, Верлет и RK4. Я ожидаю, что Эйлер и чехарда будут иметь некоторую неточность, но не RK4 или Верлет. Мои познания в области математического анализа ограничены, поэтому мне нужны дополнительные глаза для проверки.
Они следующие:
Код: Выделить всё
def leapfrog_integration(satellite, planet, dt): #most accurate under 5 orbits
# Update velocity by half-step
satellite.velocity += 0.5 * satellite.acceleration * dt
# Update position
satellite.position += (satellite.velocity / 1000)
# Calculate new acceleration
satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
# Update velocity by another half-step
satellite.velocity += 0.5 * satellite.acceleration * dt
def euler_integration(satellite, planet, dt):
# satellite.accleration = satellite.calculate_gravity(planet)
satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
satellite.velocity += satellite.acceleration * dt
satellite.position += (satellite.velocity / 1000) #convert into kilometers
def verlet_integration(satellite, planet, dt):
acc_c = (satellite.acceleration_due_to_gravity(planet) / 1000)#convert to km/s
satellite.velocity = (satellite.position - satellite.previous_position)
new_pos = 2 * satellite.position - satellite.previous_position + (acc_c * dt)
satellite.previous_position = satellite.position #km
satellite.position = new_pos #km
satellite.velocity = (satellite.position - satellite.previous_position)
def rk4_integration(satellite, planet, dt):# need to resolve the conversion to km for position. If i remove the DT from the kx_r then it's the excat same as Verlet and Euler
def get_acceleration(position, velocity):
temp_mass = PointMass(position,satellite.mass,satellite.radius,satellite.colour,(velocity), np.zeros_like(satellite.acceleration))
return temp_mass.acceleration_due_to_gravity(planet)
k1_v = dt * get_acceleration(satellite.position, (satellite.velocity))
k1_r = (satellite.velocity / 1000)
k2_v = dt * get_acceleration(satellite.position + 0.5 * k1_r, (satellite.velocity) + 0.5 * k1_v)
k2_r = (satellite.velocity + 0.5 * k1_v) / 1000
k3_v = dt * get_acceleration(satellite.position + 0.5 * k2_r, (satellite.velocity) + 0.5 * k2_v)
k3_r = (satellite.velocity + 0.5 * k2_v) / 1000
k4_v = dt * get_acceleration(satellite.position + 0.5 * k3_r, (satellite.velocity) + 0.5 * k3_v)
k4_r = (satellite.velocity + 0.5 * k3_v) / 1000
satellite.position +=(k1_r + 2*k2_r + 2*k3_r + k4_r) / 6
satellite.velocity +=(k1_v + 2*k2_v + 2*k3_v + k4_v) / 6
Код: Выделить всё
class PointMass:
def __init__(self, position, mass, radius, colour, velocity, acceleration):
self.position = np.array(position) #in KM
self.mass = mass #in Kilograms
self.radius = radius #in meters
self.colour = colour
self.velocity = np.array(velocity) #in m/s
self.acceleration = np.array(acceleration) #in m/s per second
self.gForce = None #This is in Newtons
self.previous_position = self.position - (self.velocity / 1000) # Initialize previous position for Verlet integration
def acceleration_due_to_gravity(self,other):
dVector = self.position - other.position # distance vector from self to the other point mass in pixels(km)
distance_km = np.linalg.norm(dVector) #Compute the Euclidean distance in km
distance_m = distance_km * 1000 + other.radius #the distance including the radius to the centre in meters
unit_vector = (dVector / distance_km) #the unit vector for the direction of the force
acceleration_magnitude = -Constants.mu_earth / distance_m**2
return acceleration_magnitude * (unit_vector * 1000) #Return the acceleration vector by multiplying the magnitude with the unit vector(converted to meters)
#the returned acceleration vector is in m/s
Код: Выделить всё
planet = PointMass(
position=[600.0,400.0,0.0],
mass=5.9722e24,
radius=6.371e6,
velocity=[0.0,0.0,0.0], #in m/s
acceleration=[0.0,0.0,0.0], #in m/s^2
colour=[125,100,100]
)
satellite = PointMass(
position=[280.0,400.0,0.0],
mass=100,
radius=1,
velocity=[0.0,7718.0,0.0], #need an intial velocity or else it'll jsut fall to the central mass
acceleration=[1.0,0.0,0.0],#added an non-zero acceleration jsut to make sure there's no issues with the integrations.
colour=[255,255,255]
)
Теперь мне не 100 Проблема связана с интеграцией, но я думал, что начну с этого и перейду к этому. Это может быть моя функция силы тяжести, неверное применение дельты времени, проблема с преобразованием единиц измерения или тем, как она размещается на экране.
Ссылка на GithubВыше приведена ссылка на репозиторий, чтобы все было в порядке, поскольку это может быть что угодно. Снимки экрана тестов и т. д.
Теперь можно было бы просто ожидать поведения от Python с потерей точности с числами с плавающей запятой, но я бы нажал X по этому поводу, поскольку это, скорее всего, мои плохие знания в области математического анализа.
После тщательного тестирования я получил DT 0,021 как лучший результат для разницы времени. Чехарда кажется наиболее точной как для малого количества витков, так и для большого количества витков.
Эйлер и Верле, похоже, даже меньше 100 витков, где Эйлер немного отстает, а RK4 кажется нестабильным, поскольку он медленный. замедляется, и, таким образом, орбиты становятся все больше и больше.
У меня есть преобразование метров в км в разных местах для каждой интеграции, поскольку они не будут работать правильно в зависимости от того, где я его поместил.
Мне пришлось удалить DT из некоторых интеграционных частей, так как, если бы я оставил их в объекте, объект просто слегка закрутился бы и упал бы на центральную массу.>
Подробнее здесь: https://stackoverflow.com/questions/785 ... -in-python
