Вселенная на моем столе Universe On My Table Программа для решения  знаменитой задачи небесной механики, известной как  "задача N тел" Program for PC to solve the famous problem of celestial mechanics known as "N-body problem" Download: Newton-x64.zip (64 bit version)   Введение. Постановка задачи 1. Как это было. От теории к практике 2. Несерьезные эксперименты 3. Новые горизонты! Три измерения 4. Традиционый тест. Солнечная система 5. Язык описания сцен Заключение

Введение. Постановка задачи

Закон всемирного тяготения (или 4-й закон Ньютона) гласит: сила притяжения двух материальных тел прямо пропорциональна произведению их масс и обратно пропорциональна квадрату расстояния между ними. 

Ньютон установил, что два тела под действием взаимного притяжения могут двигаться только по эллипсу, параболе или гиперболе около общего центра масс. Задача о движении трех и более тел оказалась намного более трудной. Доказано, что она не может быть решена путем получения соответствующих математических формул, описывающих движение тел в каждый момент времени, кроме отдельных частных случаев. Таковы, например, "треугольник Лагранжа" и, найденное в 2000 г. решение - "невероятная восьмерка".

С развитием вычислительной техники стал возможен "численный эксперимент". Зная взаимное расположение тел, их массы и скорости в начальный момент, легко рассчитать силу воздействия на каждое тело со стороны всех остальных. Полагая, что в короткий отрезок времени движение под этим воздействием происходит по прямой, определим новое положение и скорость каждого тела. Затем цикл расчета повторяется. Если шаг по времени выбран малым, то мы вправе ожидать, что получим достаточно точные траектории движения тел. Проверка правильности расчета состоит в его повторении для вдвое меньшего и вдвое большего временных интервалов. Если результат не меняется, значит, все в порядке.

Июньской короткой ночью над городом скрежетала гроза. С очередной зарницей я оставил напрасный счет слонов и позволил мыслям парить свободно.  И тут...  осенило.  Это можно сделать!  Когда алгоритм, короткий и ясный, окончательно сложился, я не заметил, как уснул.

"Изделие" вчерне было закончено за 2 дня, оставалось только "зачистить напильником". Почему так споро? Я - не гений, просто для альфа-версии проекта использовал язык программирования Python и готовый компонент FloatCanvas для графической библиотеки wxPython.

Впоследствии программа была переписана на Delphi, что сократило размеры дистрибутива в 10 раз, и увеличило быстродействие раз в 15-20. Python очень красивый и необычайно мощный язык, но, к сожалению, медленный. Поэтому Python часто используют для предварительного программирования сложных проектов, когда надо вовремя предъявить готовый и работающий продукт. Потом критичные блоки перекодируют на Си или Дельфи.

Практика - критерий истины, кто это сказал? Если программа будет правильно отображать движение небесных тел для классических, еще Кеплером и Ньютоном рассмотренных случаев, то заложенные в алгоритм посылки, скорее всего, верны.

А они вот такие. Сперва я рассматривал плоскую задачу многих тел, исключительно потому, что мониторы на компьютерах, сами понимаете, тоже плоские. Это не слишком большое ограничение - планеты Солнечной системы лежат практически в одной плоскости, системы спутников планет обладают тем же свойством, да и галактики при взгляде издалека в определенных ракурсах напоминают блин. Так что, на первую попытку, хватит нам икса с игреком. Ось X направлена горизонтально вправо, Y - вверх. Но, с самого начала, алгоритм был разработан для трех измерений, просто для первых экспериментальных расчетов, координата Z полагалась равной 0.

Итак, каждое небесное тело характеризуется положением (X,Y,Z), радиусом, массой и цветом, чтобы отличать от других. Если цвет не задан, он назначается программой произвольно. Все тела рисуются на экране. Затем рассчитывается смещение каждого тела под воздействием остальных - предполагается, что в краткий отрезок времени оно происходит по прямой. Возможные коллизии обрабатываются программой, как абсолютно неупругие. Шаг по времени (dt)чисто условно полагается равным 1, это упрощает запись формул.

Последовательность шагов такова (метод Эйлера-Кромера):

Если p_t, v_t, F(p_t) - векторы координат, скоростей и ускорений тел в момент t, то:

v_t+1 = v_t+F(p_t); p_t+1 = p_t+v_t+1

А можно так (метод Рунге-Кутта 4-го порядка):

k₁ = F(p_t); k₂ = F(p_t+k₁/2); k₃ = F(p_t+k₂/2); k₄ = F(p_t+k₃); v_t+1 = v_t+(k₁+2k₂+2k₃+k₄)/6; p_t+1 = p_t+v_t+1

 Или вот так (метод Верле в скоростной формулировке):

p_t+1 = p_t+v_t+F(p_t)/2; v_t+1 = v_t+F(p_t)/2+F(p_t+1)/2

Экран очищается и рисуются цветные кружочки в их новых позициях. И так много раз в секунду. Мультик готов. Уменьшая величину шага по времени, можно достичь впечатляющей точности - здесь все указанные методы равноценны.

Исходные данные для расчета выглядят, прямо-таки, как шпионский шифр. Либо... как тарабарский язык описания сцен. (Скажу по секрету, что это - популярный среди разработчиков компьютерных игр скриптовый язык Lua). Итак, начали...

Это не рецепт спагетти, а описание начальных условий хоровода из 8 звезд:

Устойчивость метода можно оценить по тому, как он моделирует классические конфигурации. Так, после 500 витков (свыше миллиона итераций) "невероятная восьмерка" ни на йоту не изменила формы. Преисполнившись энтузиазма, продолжим...             

Гипотеза Канта. Вначале был Хаос. Из него, путем объединения частиц, возникло Солнце и планеты. Круговое движение последних есть результат осреднения скоростей.

Сказано - сделано. 150 тел с экспоненциальным распределением масс хаотично движутся в согласии с учением великого философа. Старт... 3 минуты, 16000 итераций... финиш.

Комментарий. Из ничего и не будет ничего. Суммарный вращательный момент системы если и отличен от 0, то очень ненамного, за счет статистических отклонений. Все массы сливаются в одно Солнце и только несколько мелких, случайно уцелевших ошметков, болтаются вокруг по сильно вытянутым орбитам. Повторы эксперимента ничего более интересного не дали.

Гипотеза О.Ю.Шмидта. Известно, что вращательный момент планет составляет 98% вращательного момента всей солнечной системы. Предположим, что околосолнечное протопланетное облако обладало им с самого начала. Как?! Было захвачено пролетавшим сбоку Солнцем и заверчено вокруг него, как ловит партнершу фигурист на льду? Законы небесной механики этого не допускают. Сыграла роль "помощь" другой звезды? Маловероятно. В общем, примем это, как данность.

Сказано - сделано. 200 хаотически расположенных тел вращаются в одну сторону вокруг звезды. Начальные скорости  генерируются случайным образом в интервале 0.6..1.0 круговой. Распределение масс - экспоненциальное. Общая масса облака  составляет 0.5% от массы звезды. Размеры протопланетных  тел сильно преувеличены, чтобы а) видеть их на экране, б) обеспечить достаточную вероятность слипания при сближении. (Или, можно полагать, что звезда и облако очень уж карликовые. Для точного соблюдения масштабов надо брать в расчет многие тысячи протопланетных тел и иметь доступ к суперкомпьютеру ).

Вот как в нашем примере задаются начальные условия на "языке описания сцен":

h = {}
h.METHOD = 0; h.SHOWFRAME = 4; h.ASPECT = 1; h.STEPLIMIT = 40000
h.MESS = 'Otto J. Shmidt hypothesis'

h.GAMMA = 5e-7

R, M, X, Y = 2, 2e6, 50, 50

body = {}
body[1] = {R,M,X,Y,0,0,'YELLOW'}

math.randomseed(os.time())
n = 200
for i = 2, n+1 do
m = randexp(50); r = 0.1*m^0.3; x = 100*math.random(); y = 100*math.random()
vx, vy = Vcirc(h.GAMMA,M,X,Y,x,y)
vx = (1-0.4*math.random())*vx; vy = (1-0.4*math.random())*vy
body[i] = {r,m,x,y,vx,vy,'WHITE'}
end

Scene(h,body)

Расчет продлевался дважды по истечении очередных 40000 итераций. Через несколько минут,  в районе 100000-й, я решил: хватит.

Окончательное сотворение планетной системы требовало до 200000 итераций. По мере слипания протопланет число "действующих лиц" уменьшалось и первоначальное легкое торможение при смене кадров исчезало. Когда от вращения остатков облака начинало мельтешить в глазах, я использовал такой прием. Еще больше ускорял демонстрацию нажатием клавиш "Ctrl+]" до тех пор, пока не получал статическую картинку, меняющуюся не чаще раза в 1-2 сек. Таким способом 40000 итераций могут пролететь за несколько секунд. Когда становилось ясно, что процесс завершился, то возвращал нормальную смену кадров неоднократным нажатием на "Ctrl+[". Затем F10, чтобы включить режим эскиза. При остановке сцены (Ctrl+Q) через пару-тройку секунд на экране возникали планетные орбиты во всей своей красе.

Формировалось, как правило 5-6 планет, 1-2 самые крупные всегда находились в середине системы. Обычным явлением была небольшая внешняя планета с заметно вытянутой орбитой. Для больших полуосей орбит в 75% случаев удовлетворительно выполнялся закон планетных расстояний: L_n = akⁿ, где a > 0 и k > 1 - константы. Еще в 15% случаев из гармонического ряда выпадала лишь последняя планета.

А вот оставшаяся десятая доля в ЗПР не вписывалась. Все такие системы имели особенность: пару планет-сестер, обращающихся по очень близким симметричным орбитам.  Вместе с солнцем они образовывали фигуру, похожую на треугольник Лагранжа, с той разницей, что угол между радиус-векторами в некоторых случаях периодически менялся. (От 30 до 120 градусов).

Комментарий. Не будем бомбардировать телеграммами Академию наук! Вспомним, что наши доморощенные эксперименты - это всего лишь завораживающая игра. Так к ним и нужно относиться. И, вообще, задача N тел может иметь "хаотический" характер, когда ничтожное отклонение в начальных условиях совершено меняет вид траекторий. Даже изменение шага расчета в этих примерах приводило к построению другой планетной системы. Хотя качественно характер решения сохранялся.

<<

3. Новые горизонты! Три измерения

 

Попробуем теперь задавать координаты и скорости небесных тел в 3-х измерениях.  Переключим режим показа на "Изометрию" и посмотрим. Первое впечатление: решения задачи N тел в пространстве намного богаче, чем на плоскости!

Еще примеры. Четыре тела после сумасшедшего вальса с обменом партнерами, разбиваются, в итоге, на 2 пары (красно-синяя, желто-зеленая), обращающиеся вокруг друг друга по вытянутым орбитам. После чего обе пары сливаются в экстазе. Вот такая вселенская любовь.

Два одинаковых тела вращаются по эллиптическим орбитам вокруг общего центра масс. Третье (ничтожно малой массы)болтается взад-вперед,  наподобие маятника, на прямой, перпендикулярной плоскости вращения и проходящей через центр. Что-то я раньше не слышал о таком очевидном решении. Еще более впечатляющая картина получается, когда третье тело имеет массу, сравнимую с первыми двумя. Здесь массы тел равны: $0|4|(-6,-1,-11),(6,10,11),(0,0,0),2|0|0|1|68740|1| TERRA|40.0 0.3|1.0|-5|0.0|0.0|0.0|0.0|-0.012|YELLOW 0.3|1.0|5|0.0|0.0|0.0|0.0|0.012|MAGENTA 0.3|1|0.0|9.0|0.0|0.0|0.0|0.0|GREEN "Естественного" развала системы от накопления вычислительных ошибок не происходит даже после миллиона итераций, что в реальном времени эквивалентно примерно году. Говорит ли это в пользу устойчивости подобного решения? А теперь... вполне традиционный тест. << 4. Традиционный тест. Солнечная система      "Вселенная на моем столе" - это не программа-планетарий, в которой все орбиты заранее заданы. Но если взять координаты и векторы скоростей планет в некоторый момент времени, и запустить расчет...  Должно получиться именно то, что мы изучаем в школе. У меня были необходимые данные на момент 1.01.07 00:00, хотя не слишком точные. Впрочем, получилось очень похоже. Для удобства внутренние планеты с "большой четверкой" астероидов и внешние с планетоидом Плутон - считались отдельно. Уж больно велика разница масштабов. Меркурий, Венера, Земля, Марс. Астероиды: Церера, Паллада, Юнона, Веста. Шаг расчета - 0.1 суток = 2.4 часа. Время расчета - 1 мин. Гиганты: Юпитер, Сатурн, Уран, Нептун и планетоид Плутон, недавно исключенный астрономами из списка планет. Запомните: планет в Солнечной системе - восемь! Для сравнения масштабов показана и орбита Цереры. Здесь шаг расчета равный одним суткам дает такую же хорошую точность. Время расчета - 1 мин.      Вид на солнечную систему "сбоку" доказывает, что Плутон - явный чужак, хулиган, затесавшийся в приличное общество. В заключение. "Полнокомплектная" Солнечная система со всеми планетами, астероидом Церера и бродягой Плутоном считалась 15 минут при шаге в 0.1 суток. (Хотелось сохранить точность расчета орбит внутренних планет). За это время промчались 250 виртуальных лет... Рисунок не привожу. Мысленно уменьшите в 13 раз первый и вставьте во второй, внутрь орбиты Цереры. Много увидите? То-то. Но в ходе расчета и по его окончании картинку в окне можно увеличивать, уменьшать, переключать проекцию и двигать, как хочешь. А потом сохранять выбранный фрагмент в буфер обмена по Alt-PrintScreen для создания скриншотов. << 5. Язык описания сцен Это набор команд, автоматическое выполнение которых приводит к созданию соответствующего nwt-файла. В основном это нужно, когда в сцене участвует множество небесных тел и nwt-файл трудно создать "вручную". Как уже упоминалось, в качестве языка описания сцен используется язык программирования Lua, предназначенный для встраивания в другие приложения, чтобы дать пользователям возможность писать свои высокоуровневые сценарии. Интерпретатор Lua является неотемлемой частью программы "Вселенная на моем столе". <<

Заключение

Автор должен признаться, что ему нечего сказать в качестве последнего слова. Он сделал то, что мог, другие пусть сделают лучше. Закон всемирного тяготения никто не отменял и в ближайшие миллиарды лет этого не предвидится. Поэтому задача расчета движения многих тел в поле взаимного тяготения всегда будет актуальной. <<