Моделирование орбит с использованием законов физики
За последние пару недель я пытался симулировать орбиты в симуляции солнечной системы, которую я выполняю в рамках университетского модуля. Короче говоря, моя симуляция написана на C ++ с использованием механизма рендеринга Ogre3D. Я попытался реализовать орбиты, используя закон всемирного тяготения Ньютона, который заставил мою планету направиться к Солнцу по прямой линии, пройти через Солнце и затем вернуться в исходное положение. Я также попробовал шаги из раздела «Положение как функция времени» этой википедии. статья, но у меня это тоже не сработало.
Я веду симуляцию с помощью простого метода интеграции Эйлера. Если кто-то имеет опыт работы с такого рода симуляциями или просто вообще много знает об этих законах физики, то любая помощь или указание мне в правильном направлении будет принята с благодарностью.
Решение
Вы должны дать планете начальную скорость v = (vx, vy, vz), касательную к желаемой орбите. Если положение солнца — s, а планета — p, то между ними всегда действует сила: та, что на планете, указывает на солнце, вектор t = (s — p) и наоборот. Величина этой силы равна g Ms Mp / (t dot t), где «точка» — это произведение точек, g — стандартное ускорение под действием силы тяжести, а Ms, Mp — соответствующие массы.
Если вы выполняете детальную модель, в которой все тела могут оказывать влияние на все другие тела, то алгоритм состоит в том, чтобы накапливать все попарные силы, чтобы получить один результирующий вектор силы, действующий на каждое тело (планету или солнце). В противном случае вы можете согласиться на приближение, когда только солнце тянет за планеты, а другие силы считаются слишком маленькими, чтобы иметь значение.
Итак, алгоритм такой:
Choose dt, the time step, a small interval.
Set initial positions and velocities
(for planets, velocity is tangent to desired orbit. Sun has velocity zero.)
loop
Accumulate all forces on all bodies with current positions.
For each body position=p, velocity=v, net resultant force=f, mass=m,
update its velocity: v_new = v + f / m dt
update position p_new = p + 0.5 * (v + v_new)
v = v_new; p = p_new
Render
end loop
Как уже упоминалось, Эйлер прост, но требует очень небольшого временного шага, чтобы получить даже разумную точность. Иногда вы можете ввести лишь небольшое сопротивление в системе (умножьте скорость на коэффициент чуть чуть меньше 1), чтобы сохранить стабильность, когда она взорвется.
Другие решения
Проконсультируйтесь с проектом «Движущиеся звезды вокруг«есть старая версия C и модернизированная версия Ruby. (А теперь версия C ++?)
Краткий совет: методы Эйлера вредны для энергосбережения. явный эйлер увеличивает энергию, неявный эйлер уменьшает энергию. Просто проверьте это на изображении фазового пространства гармонического осциллятора y » + y = 0.
Используйте симплектические интеграторы, самым простым и известным является Метод чехарды или верлет что уже Ньютон рассуждал о движении планет.
Вы могли бы попробовать Рунге Кутта 4, но на орбитах все еще будет некоторый «дрейф». Поддержание постоянной энергии системы необходимо для стабилизации системы, но я не уверен, как это сделать. Лучший метод для небесная механика является симплектический интегратор.