Моделирование динамических систем: Как движется Луна?

  • Tutorial
Светлой памяти моего учителя — первого декана физико-математического факультета Новочеркасского политехнического института, заведующего кафедрой «Теоретическая механика» Кабелькова Александра Николаевича

Введение


Август, лето подходит к концу. Народ яростно рванул на моря, да оно и неудивительно — самый сезон. А на Хабре, тем временем, буйным цветом распускается и пахнет лженаука. Если говорить о теме данного выпуска «Моделирования...», то в нем мы совместим приятное с полезным — продолжим обещанный цикл и совсем чуть-чуть поборемся с этой самой лженаукой за пытливые умы современной молодежи.


А вопрос ведь действительной не праздный — со школьных лет мы привыкли считать, что наш ближайший спутник в космическом пространстве — Луна движется вокруг Земли с периодом 29,5 суток, особенно не вдаваясь в сопутствующие подробности. На самом же деле наша соседка своеобразный и в какой-то степени уникальный астрономический объект, с движением которого вокруг Земли не всё так просто, как, возможно хотелось бы некоторым моим коллегам из ближайшего зарубежья.

Итак, оставив полемику в стороне, попытаемся с разных сторон, в меру своей компетенции, рассмотреть эту безусловно красивую, интересную и очень показательную задачу.

1. Закон всемирного тяготения и какие выводы мы можем из него сделать


Открытый ещё во второй половине 17 века, сэром Исааком Ньютоном, закон всемирного тяготения говорит о том, что Луна притягивается к Земле (и Земля к Луне!) с силой, направленной вдоль прямой, соединяющей центры рассматриваемых небесных тел, и равной по модулю

$F_{1,2} = G \, \frac{m_1 \, m_2}{r_{1,2}^2} $


где m1, m2 — массы, соответственно Луны и Земли; G = 6,67e-11 м3/(кг * с2) — гравитационная постоянная; r1,2 — расстояние между центрами Луны и Земли. Если принимать во внимание только эту силу, то, решив задачу о движении Луны как спутника Земли и научившись рассчитывать положение Луны на небе на фоне звезд, мы довольно скоро убедимся, путем прямых измерений экваториальных координат Луны, что в нашей консерватории не всё так гладко как хотелось бы. И дело здесь не в законе всемирного тяготения (а на ранних этапах развития небесной механики такие мысли высказывались весьма нередко), а в неучтенном возмущении движения Луны со стороны других тел. Каких? Смотрим на небо и наш взгляд сразу упирается в здоровенный, массой аж 1,99e30 килограмм плазменный шар прямо у нас под носом — Солнце. Луна притягивается к Солнцу? Ещё как, с силой, равной по модулю

$F_{1,3} = G \, \frac{m_1 \, m_3}{r_{1,3}^2} $


где m3 — масса Солнца; r1,3 — расстояние от Луны до Солнца. Сравним эту силу с предыдущей

$\frac{F_{1,3}}{F_{1,2}} = \frac{G \, \frac{m_1 \, m_3}{r_{1,3}^2}}{G \, \frac{m_1 \, m_2}{r_{1,2}^2}} = \frac{m_3}{m_2} \, \left(\frac{r_{1,2}}{r_{1,3}}\right)^2$


Возьмем положение тел, в котором притяжение Луны к Солнцу будет минимальным: все три тела на одной прямой и Земля располагается между Луной и Солнцем. В этом случае наша формула примет вид:

$\frac{F_{1,3}}{F_{1,2}} = \frac{m_3}{m_2} \, \left(\frac{\rho}{a + \rho}\right)^2$


где $\rho = 3,844 \cdot 10^{8}$, м — среднее расстояние от Земли до Луны; $a = 1,496\cdot10^{11}$, м — среднее расстояние от Земли до Солнца. Подставим в эту формулу реальные параметры

$\frac{F_{1,3}}{F_{1,2}} = \frac{1.99 \cdot 10^{30}}{5.98\cdot10^{24}} \, \left( \frac{3.844\cdot10^{8}}{1.496\cdot10^{11} + 3.844\cdot10^{8}}\right)^2 = 2.19$


Вот это номер! Получается Луна притягивается к Солнцу силой, более чем в два раза превышающей силу её притяжения к Земле.

Подобное возмущение уже нельзя не учитывать и оно определенно повлияет на конечную траекторию движения Луны. Пойдем дальше, принимая во внимание допущение о том, что орбита Земли круговая с радиусом a, найдем геометрическое место точек вокруг Земли, где сила притяжения любого объекта к Земле равна силе его притяжения к Солнцу. Это будет сфера, с радиусом

$R = \frac{a \, \sqrt{\gamma}}{1 - \gamma}$


смещенная вдоль прямой, соединяющей Землю и Солнце в сторону противоположенную направлению на Солнце на расстояние

$l = R \, \sqrt{\gamma}$


где $\gamma = m_2 / m_3 $ — отношение массы Земли к массе Солнца. Подставив численные значения параметров получим фактические размеры данной области: R = 259300 километров, и l = 450 километров. Эта сфера носит название сферы тяготения Земли относительно Солнца.

Известная нам орбита Луны лежит вне этой области. То есть в любой точке траектории Луна испытывает со стороны Солнца существенно большее притяжение, чем со стороны Земли.

2. Спутник или планета? Гравитационная сфера действия


Эта информация, часто порождает споры, о том, что Луна не спутник Земли, а самостоятельная планета Солнечной системы, орбита которой возмущена притяжением близкой Земли.

Оценим возмущение, вносимое Солнцем в траекторию Луны относительно Земли, а так же возмущение, вносимое Землей в траекторию Луны относительно Солнца, воспользовавшись критерием, предложенным П. Лапласом. Рассмотрим три тела: Солнце (S), Землю (E) и Луну (M).
Примем допущение, что орбиты Земли относительно Солнца и Луны относительно Земли являются круговыми.


Рассмотрим движение Луны в геоцентрической инерциальной системе отсчета. Абсолютное ускорение Луны в гелиоцентрической системе отсчета определяется действующими на неё силами тяготения и равно:

$\vec a_1 = \vec a_1^{(3)} + \vec a_1^{(2)} = \frac{1}{m_1} \, \vec F_{1,3} + \frac{1}{m_1} \, \vec F_{1,2}$


С другой стороны, в соответствии с теоремой Кориолиса, абсолютное ускорение Луны

$\vec a_1 = \vec a_2 + \vec a_{1,2}$


где $\vec a_2$ — переносное ускорение, равное ускорению Земли относительно Солнца; $\vec a_{1,2}$ — ускорение Луны относительно Земли. Ускорения Кориолиса здесь не будет — выбранная нами система координат движется поступательно. Отсюда получаем ускорение Луны относительно Земли

$\vec a_{1,2} = \frac{1}{m_1} \, \vec F_{1,3} + \frac{1}{m_1} \, \vec F_{1,2} - \vec a_2$


Часть этого ускорения, равная $ \vec a_1^{(2)} = \frac{1}{m_1} \, \vec F_{1,2}$ обусловлена притяжением Луны к Земле и характеризует её невозмущенное геоцентрическое движение. Оставшаяся часть

$\Delta \vec a_{1,3} = \frac{1}{m_1} \, \vec F_{1,3} - \vec a_2$


ускорение Луны, вызванное возмущением со стороны Солнца.

Если рассматривать движение Луны в гелиоцентрической инерциальной системе отсчета, то всё намного проще, ускорение $\vec a_1^{(3)} = \frac{1}{m_1} \, \vec F_{1,3} $ характеризует невозмущенное гелиоцентрическое движение Луны, а ускорение $\Delta \vec a_{1,2} = \frac{1}{m_1} \, \vec F_{1,2} $ — возмущение этого движения со стороны Земли.

При существующих в текущую эпоху параметрах орбит Земли и Луны, в каждой точке траектории Луны справедливо неравенство

$\frac{|\Delta \vec a_{1,3}|}{|\vec a_{1}^{(2)}|} < \frac{|\Delta \vec a_{1,2}|}{|\vec a_{1}^{(3)}|}\quad\quad(1) $


что можно проверить и непосредственным вычислением, но я сошлюсь на источник, дабы излишне не загромождать статью.

Что означает неравенство (1)? Да то, что в относительном выражении эффект от возмущения Луны Солнцем (причем очень существенно) меньше эффекта от притяжения Луны к Земле. И наоборот, возмущение Землей геолиоцентрической траектории Луны оказывает решающее влияние на характер её движения. Влияние земной гравитации в данном случае более существенно, а значит Луна «принадлежит» Земле по праву и является её спутником.

Интересным является другое — превратив неравенство (1) в уравнение можно найти геометрическое место точек, где эффекты возмущения Луны (да и любого другого тела) Землей и Солнцем одинаковы. К сожалению это у же не так просто, как в случае со сферой тяготения. Расчеты показывают, что данная поверхность описывается уравнением сумасшедшего порядка, но близка к эллипсоиду вращения. Всё что мы может сделать без лишних заморочек, это оценить общие габариты этой поверхности относительно центра Земли. Решая численно уравнение

$\frac{|\Delta \vec a_{1,3}|}{|\vec a_{1}^{(2)}|} = \frac{|\Delta \vec a_{1,2}|}{|\vec a_{1}^{(3)}|}\quad\quad(2) $


относительно расстояния от центра Земли до искомой поверхности на достаточном количестве точек, получаем сечение искомой поверхности плоскостью эклиптики


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

Для практических расчетов данную поверхность удобно аппроксимировать сферой с центром в центра Земли и радиусом равным

$r = a \, \left(\frac{m}{M} \right)^{\frac{2}{5}}\quad\quad(3)$


где m — масса меньшего небесного тела; M — масса большего тела, в поле тяготения которого движется меньшее тело; a — расстояние между центрами тел. В нашем случае

$r = a \, \left(\frac{m_2}{m_3} \right)^{\frac{2}{5}} = 1.496\cdot10^{11} \, \left(\frac{5.98\cdot10^{24}}{1.99\cdot10^{30}} \right)^{\frac{2}{5}} = 925000,\, км$



Вот этот недоделанный миллион километров и есть тот теоретический предел, за который власть старушки Земли не распространяется — её влияние на траектории астрономических объектов настолько мало, что им можно пренебречь. А значит, запустить Луну по круговой орбите на расстоянии 38,4 млн. километров от Земли (как делают некоторые лингвисты) не получится, это физически невозможно.

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

Например, теперь понятно, для того чтобы иметь теоретическую возможность совершить маневры для выхода на окололунную орбиту, космический аппарат должен попасть внутрь сферы действия Луны относительно Земли. Её радиус легко рассчитать по формуле (3) и он равен 66 тысяч километров.

Таким образом, Луна справедливо может считаться спутником Земли. Однако, ввиду существенно влияния гравитационного поля Солнца она движется не в центральном гравитационном поле, а значит её траектория не является коническим сечением.

3. Задача трех тел в классической постановке


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


Тела считаем материальными точками. Положение тел будем отсчитывать в произвольном базисе, с которым связана инерциальная система отсчета Oxyz. Положение каждого из тел задается радиус-вектором соответственно $\vec r_1$, $\vec r_2$ и $\vec r_3$. На каждое тело действует сила гравитационного притяжения со стороны двух других тел, причем в соответствии с третьей аксиомой динамики точки (3-й закон Ньютона)

$\vec F_{i,j} = -\vec F_{j,i}\quad\quad(4)$



Запишем дифференциальные уравнения движения каждой точки в векторной форме

$\begin{align} & m_1 \, \frac{d^2 \vec r_1}{dt^2} = \vec F_{1,2} + \vec F_{1,3} \\ & m_2 \, \frac{d^2 \vec r_2}{dt^2} = \vec F_{2,1} + \vec F_{2,3} \\ & m_3 \, \frac{d^2 \vec r_3}{dt^2} = \vec F_{3,1} + \vec F_{3,2} \end{align}$



или, с учетом (4)

$\begin{align} & m_1 \, \frac{d^2 \vec r_1}{dt^2} = \vec F_{1,2} + \vec F_{1,3} \\ & m_2 \, \frac{d^2 \vec r_2}{dt^2} = -\vec F_{1,2} + \vec F_{2,3} \\ & m_3 \, \frac{d^2 \vec r_3}{dt^2} = -\vec F_{1,3} - \vec F_{2,3} \end{align}$


В соответствии с законом всемирного тяготения, силы взаимодействия направлены вдоль векторов

$\begin{align} & \vec r_{1,2} = \vec r_2 - \vec r_1 \\ & \vec r_{1,3} = \vec r_3 - \vec r_1 \\ & \vec r_{2,3} = \vec r_3 - \vec r_2 \\ \end{align}$


Вдоль каждого из этих векторов выпустим соответствующий орт

$\vec e_{i,j} = \frac{1}{r_{i,j}} \, \vec r_{i,j}$


тогда каждая из гравитационных сил рассчитывается по формуле

$\vec F_{i,j} = G\,\frac{m_i \, m_j}{r_{i,j}^2}\,\vec e_{i,j}$


С учетом всего этого система уравнений движения принимает вид

$\begin{align} & \frac{d^2 \vec r_1}{dt^2} = \frac{G\,m_2}{r_{1,2}^3} \, \vec r_{1,2} + \frac{G\,m_3}{r_{1,3}^3} \, \vec r_{1,3} \\ & \frac{d^2 \vec r_2}{dt^2} = -\frac{G\,m_1}{r_{1,2}^3} \, \vec r_{1,2} + \frac{G\,m_3}{r_{2,3}^3} \, \vec r_{2,3} \\ & \frac{d^2 \vec r_3}{dt^2} = -\frac{G\,m_1}{r_{1,3}^3} \, \vec r_{1,3} - \frac{G\,m_2}{r_{2,3}^3} \, \vec r_{2,3} \end{align}$


Введем обозначение, принятое в небесной механике

$\mu_i = G\,m_i$


— гравитационный параметр притягивающего центра. Тогда уравнения движения примут окончательный векторный вид

$\begin{align} & \frac{d^2 \vec r_1}{dt^2} = \frac{\mu_2}{r_{1,2}^3} \, \vec r_{1,2} + \frac{\mu_3}{r_{1,3}^3} \, \vec r_{1,3} \\ & \frac{d^2 \vec r_2}{dt^2} = -\frac{\mu_1}{r_{1,2}^3} \, \vec r_{1,2} + \frac{\mu_3}{r_{2,3}^3} \, \vec r_{2,3} \\ & \frac{d^2 \vec r_3}{dt^2} = -\frac{\mu_1}{r_{1,3}^3} \, \vec r_{1,3} - \frac{\mu_2}{r_{2,3}^3} \, \vec r_{2,3} \end{align}$



4. Нормирование уравнений к безразмерным переменным


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

Итак, введем некое абстрактное небесное тело с гравитационным параметром $\mu$, такое, что период обращения спутника по эллиптической орбите с большой полуосью $a$ вокруг него равен $T$. Все эти величины, в силу законов механики, связаны соотношением

$T = 2\,\pi\,\left(\frac{a^3}{\mu}\right)^{\frac{1}{2}}$


Введем замену параметров. Для положения точек нашей системы

$\vec r_i = a \, \vec\xi_i$


где $\vec\xi_i$ — безразмерный радиус-вектор i-й точки;
для гравитационных параметров тел

$\mu_i = \varkappa_i \, \mu$


где $\varkappa_i$ — безразмерный гравитационный параметр i-й точки;
для времени

$t = T \, \tau$


где $\tau$ — безразмерное время.

Теперь пересчитаем ускорения точек системы через эти безразмерные параметры. Применим прямое двукратное дифференцирование по времени. Для скоростей

$\vec v_i = \frac{d \vec r_i}{dt} = a \, \frac{d\vec\xi_i}{dt}=\frac{a}{T} \, \frac{d\vec\xi_i}{d\tau}=\frac{1}{2\,\pi} \, \sqrt{\frac{\mu}{a}}\,\frac{d\vec\xi_i}{d\tau}.$


Для ускорений

$\vec a_i = \frac{d\vec v_i}{dt} = \frac{1}{2\,\pi} \, \sqrt{\frac{\mu}{a}}\,\frac{1}{dt} \left(\frac{d\vec\xi_i}{d\tau}\right) = \frac{1}{4\,\pi^2} \, \frac{\mu}{a^2}\,\frac{d^2\vec \xi_i}{d\tau^2}$



При подстановке полученных соотношений в уравнения движения всё элегантно схлопывается в красивые уравнения:

$\begin{align} &\frac{d^2 \vec\xi_1}{d\tau^2} = 4\,\pi^2 \, \varkappa_2\,\frac{\vec \xi_2 - \vec \xi_1}{|\vec \xi_2 - \vec \xi_1|^3} + 4\,\pi^2 \, \varkappa_3\,\frac{\vec \xi_3 - \vec \xi_1}{|\vec \xi_3 - \vec \xi_1|^3}\\ & \frac{d^2 \vec\xi_2}{d\tau^2} = -4\,\pi^2 \, \varkappa_1\,\frac{\vec \xi_2 - \vec \xi_1}{|\vec \xi_2 - \vec \xi_1|^3} + 4\,\pi^2 \, \varkappa_3\,\frac{\vec \xi_3 - \vec \xi_2}{|\vec \xi_3 - \vec \xi_2|^3}\quad\quad(5) \\ & \frac{d^2 \vec\xi_3}{d\tau^2} = -4\,\pi^2 \, \varkappa_1\,\frac{\vec \xi_3 - \vec \xi_1}{|\vec \xi_3 - \vec \xi_1|^3} - 4\,\pi^2 \, \varkappa_2\,\frac{\vec \xi_3 - \vec \xi_2}{|\vec \xi_3 - \vec \xi_2|^3} \end{align}$



Данная система уравнений до сих пор считается не интегрируемой в аналитических функциях. Почему считается а не является? Потому что успехи теории функции комплексного переменного привели к тому, что общее решение задачи трех тел таки появилось в 1912 году — Карлом Зундманом был найден алгоритм отыскания коэффициентов для бесконечных рядов относительно комплексного параметра, теоретически являющихся общим решением задачи трех тел. Но… для применения рядов Зундмана в практических расчетах с требуемой для них точностью требует получения такого числа членов этих рядов, что эта задача во много превосходит возможности вычислительных машин даже на сегодняшний день.

Поэтому численное интегрирование — единственный способ анализа решения уравнения (5)

5. Расчет начальных условий: добываем исходные данные


Как я уже писал ранее, прежде чем начинать численное интегрирование, следует озаботится расчетом начальных условий для решаемой задачи. В рассматриваемой задаче поиск начальных условий превращается в самостоятельную подзадачу, так как система (5) дает нам девять скалярных уравнений второго порядка, что при переходе к нормальной форме Коши повышает порядок системы ещё в 2 раза. То есть нам необходимо рассчитать целых 18 параметров — начальные положения и компоненты начальной скорости всех точек системы. Где мы возьмем данные о положении интересующих нас небесных тел? Мы живем в мире, где человек ходил по Луне — естественно человечество должно обладать информацией, как эта самая Луна движется и где она находится.

То есть, скажете вы, ты, чувак, предлагаешь нам взять с полок толстые астрономические справочники, сдуть с них пыль… Не угадали! Я предлагаю сходить за этими данными к тем, кто собственно ходил по Луне, к NASA, а именно в Лабораторию реактивного движения, Пасадена, штат Калифорния. Вот сюда — JPL Horizonts web interface.

Здесь, потратив немного времени на изучение интерфейса, мы добудем все необходимые нам данные. Выберем дату, например, да нам всё равно, но пусть это будет 27 июля 2018 года UT 20:21. Как раз в этот момент наблюдалась полная фаза лунного затмения. Программа выдаст нам огромную портянку

Полный вывод для эфемерид Луны на 27.07.2018 20:21 (начало координат в центре Земли)
*******************************************************************************
 Revised: Jul 31, 2013             Moon / (Earth)                           301
 
 GEOPHYSICAL DATA (updated 2018-Aug-13):
  Vol. Mean Radius, km  = 1737.53+-0.03    Mass, x10^22 kg       =    7.349
  Radius (gravity), km  = 1738.0           Surface emissivity    =    0.92
  Radius (IAU), km      = 1737.4           GM, km^3/s^2          = 4902.800066
  Density, g/cm^3       =    3.3437        GM 1-sigma, km^3/s^2  =  +-0.0001  
  V(1,0)                =   +0.21          Surface accel., m/s^2 =    1.62
  Earth/Moon mass ratio = 81.3005690769    Farside crust. thick. = ~80 - 90 km
  Mean crustal density  = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km 
  Heat flow, Apollo 15  = 3.1+-.6 mW/m^2   k2                    = 0.024059
  Heat flow, Apollo 17  = 2.2+-.5 mW/m^2   Rot. Rate, rad/s      = 0.0000026617
  Geometric Albedo      =    0.12

  Mean angular diameter = 31'05.2"         Orbit period          = 27.321582 d
  Obliquity to orbit    = 6.67 deg         Eccentricity          = 0.05490
  Semi-major axis, a    = 384400 km        Inclination           = 5.145 deg
  Mean motion, rad/s    = 2.6616995x10^-6  Nodal period          = 6798.38 d
  Apsidal period        = 3231.50 d        Mom. of inertia C/MR^2= 0.393142
  beta (C-A/B), x10^-4  = 6.310213         gamma (B-A/C), x10^-4 = 2.277317
 
                                 Perihelion  Aphelion    Mean
  Solar Constant (W/m^2)         1414+-7     1323+-7     1368+-7
  Maximum Planetary IR (W/m^2)   1314        1226        1268
  Minimum Planetary IR (W/m^2)      5.2         5.2         5.2
 *******************************************************************************
 
 
*******************************************************************************
Ephemeris / WWW_USER Wed Aug 15 20:45:05 2018 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Moon (301)                      {source: DE431mx}
Center body name: Earth (399)                     {source: DE431mx}
Center-site name: BODY CENTER
*******************************************************************************
Start time      : A.D. 2018-Jul-27 20:21:00.0003 TDB
Stop  time      : A.D. 2018-Jul-28 20:21:00.0003 TDB
Step-size       : 0 steps
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii    : 6378.1 x 6378.1 x 6356.8 km     {Equator, meridian, pole}    
Output units    : AU-D                                                         
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0                                                 
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch                 
*******************************************************************************
JDTDB
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
*******************************************************************************
$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB 
 X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06
 VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05
 LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06
$$EOE
*******************************************************************************
Coordinate system description:

  Ecliptic and Mean Equinox of Reference Epoch

    Reference epoch: J2000.0
    XY-plane: plane of the Earth's orbit at the reference epoch
              Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76)
    X-axis  : out along ascending node of instantaneous plane of the Earth's
              orbit and the Earth's mean equator at the reference epoch
    Z-axis  : perpendicular to the xy-plane in the directional (+ or -) sense
              of Earth's north pole at the reference epoch.

  Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]:

    JDTDB    Julian Day Number, Barycentric Dynamical Time
      X      X-component of position vector (au)                               
      Y      Y-component of position vector (au)                               
      Z      Z-component of position vector (au)                               
      VX     X-component of velocity vector (au/day)                           
      VY     Y-component of velocity vector (au/day)                           
      VZ     Z-component of velocity vector (au/day)                           
      LT     One-way down-leg Newtonian light-time (day)                       
      RG     Range; distance from coordinate center (au)                       
      RR     Range-rate; radial velocity wrt coord. center (au/day)            

Geometric states/elements have no aberrations applied.

 Computations by ...
     Solar System Dynamics Group, Horizons On-Line Ephemeris System
     4800 Oak Grove Drive, Jet Propulsion Laboratory
     Pasadena, CA  91109   USA
     Information: http://ssd.jpl.nasa.gov/
     Connect    : telnet://ssd.jpl.nasa.gov:6775  (via browser)
                  http://ssd.jpl.nasa.gov/?horizons
                  telnet ssd.jpl.nasa.gov 6775    (via command-line)
     Author     : Jon.D.Giorgini@jpl.nasa.gov
*******************************************************************************


Бр-р-р, что это? Без паники, для того, кто хорошо учил в школе астрономию, механику и математику тут боятся нечего. Итак, самое главное конечное искомые координаты и компоненты скорости Луны.

$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB 
X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06
VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05
LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06
$$EOE

Да-да-да, они декартовы! Если внимательно прочесть всю портянку, то мы узнаем, что начало этой системы координат совпадает с центром Земли. Плоскость XY лежит в плоскости земной орбиты (плоскости эклиптики) на эпоху J2000. Ось X направлена вдоль линии пересечения плоскости экватора Земли и эклиптики в точку весеннего равноденствия. Ось Z смотрит в направлении северного полюса Земли перпендикулярно плоскости эклиптики. Ну а ось Y дополняет всё это счастье до правой тройки векторов. По-умолчанию единицы измерения координат: астрономические единицы (умнички из NASA приводят и величину автрономической единицы в километрах). Единицы измерения скорости: астрономические единицы в день, день принимается равным 86400 секундам. Полный фарш!

Аналогичную информацию мы можем получить и для Земли

Полный вывод эфемерид Земли на 27.07.2018 20:21 (начало координат в центре масс Солнечной системы)
*******************************************************************************
 Revised: Jul 31, 2013                   Earth                              399
 
 GEOPHYSICAL PROPERTIES (revised Aug 13, 2018):
  Vol. Mean Radius (km)    = 6371.01+-0.02   Mass x10^24 (kg)= 5.97219+-0.0006
  Equ. radius, km          = 6378.137        Mass layers:
  Polar axis, km           = 6356.752          Atmos         = 5.1   x 10^18 kg
  Flattening               = 1/298.257223563   oceans        = 1.4   x 10^21 kg
  Density, g/cm^3          = 5.51              crust         = 2.6   x 10^22 kg
  J2 (IERS 2010)           = 0.00108262545     mantle        = 4.043 x 10^24 kg
  g_p, m/s^2  (polar)      = 9.8321863685      outer core    = 1.835 x 10^24 kg
  g_e, m/s^2  (equatorial) = 9.7803267715      inner core    = 9.675 x 10^22 kg
  g_o, m/s^2               = 9.82022         Fluid core rad  = 3480 km
  GM, km^3/s^2             = 398600.435436   Inner core rad  = 1215 km
  GM 1-sigma, km^3/s^2     =      0.0014     Escape velocity = 11.186 km/s
  Rot. Rate (rad/s)        = 0.00007292115   Surface Area:
  Mean sidereal day, hr    = 23.9344695944     land          = 1.48 x 10^8 km
  Mean solar day 2000.0, s = 86400.002         sea           = 3.62 x 10^8 km
  Mean solar day 1820.0, s = 86400.0
  Moment of inertia        = 0.3308          Love no., k2    = 0.299
  Mean Temperature, K      = 270             Atm. pressure   = 1.0 bar
  Vis. mag. V(1,0)         = -3.86           Volume, km^3    = 1.08321 x 10^12
  Geometric Albedo         = 0.367           Magnetic moment = 0.61 gauss Rp^3
  Solar Constant (W/m^2)   = 1367.6 (mean), 1414 (perihelion), 1322 (aphelion)

 ORBIT CHARACTERISTICS:
  Obliquity to orbit, deg  = 23.4392911  Sidereal orb period  = 1.0000174 y
  Orbital speed, km/s      = 29.79       Sidereal orb period  = 365.25636 d
  Mean daily motion, deg/d = 0.9856474   Hill's sphere radius = 234.9       
*******************************************************************************
 
 
*******************************************************************************
Ephemeris / WWW_USER Wed Aug 15 21:16:21 2018 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Earth (399)                     {source: DE431mx}
Center body name: Solar System Barycenter (0)     {source: DE431mx}
Center-site name: BODY CENTER
*******************************************************************************
Start time      : A.D. 2018-Jul-27 20:21:00.0003 TDB
Stop  time      : A.D. 2018-Jul-28 20:21:00.0003 TDB
Step-size       : 0 steps
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii    : (undefined)                                                  
Output units    : AU-D                                                         
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0                                                 
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch                 
*******************************************************************************
JDTDB
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
*******************************************************************************
$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB 
 X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05
 VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07
 LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05
$$EOE
*******************************************************************************
Coordinate system description:

  Ecliptic and Mean Equinox of Reference Epoch

    Reference epoch: J2000.0
    XY-plane: plane of the Earth's orbit at the reference epoch
              Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76)
    X-axis  : out along ascending node of instantaneous plane of the Earth's
              orbit and the Earth's mean equator at the reference epoch
    Z-axis  : perpendicular to the xy-plane in the directional (+ or -) sense
              of Earth's north pole at the reference epoch.

  Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]:

    JDTDB    Julian Day Number, Barycentric Dynamical Time
      X      X-component of position vector (au)                               
      Y      Y-component of position vector (au)                               
      Z      Z-component of position vector (au)                               
      VX     X-component of velocity vector (au/day)                           
      VY     Y-component of velocity vector (au/day)                           
      VZ     Z-component of velocity vector (au/day)                           
      LT     One-way down-leg Newtonian light-time (day)                       
      RG     Range; distance from coordinate center (au)                       
      RR     Range-rate; radial velocity wrt coord. center (au/day)            

Geometric states/elements have no aberrations applied.

 Computations by ...
     Solar System Dynamics Group, Horizons On-Line Ephemeris System
     4800 Oak Grove Drive, Jet Propulsion Laboratory
     Pasadena, CA  91109   USA
     Information: http://ssd.jpl.nasa.gov/
     Connect    : telnet://ssd.jpl.nasa.gov:6775  (via browser)
                  http://ssd.jpl.nasa.gov/?horizons
                  telnet ssd.jpl.nasa.gov 6775    (via command-line)
     Author     : Jon.D.Giorgini@jpl.nasa.gov
*******************************************************************************


Здесь в качестве начала координат выбран барицентр (центр масс) Солнечной системы. Интересующие нас данные

$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB 
 X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05
 VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07
 LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05
$$EOE

Для Луны нам понадобятся координаты и скорость относительно барицентра Солнечной системы, мы можем их посчитать, а можем попросит NASA дать нам такие данные

Полный вывод эфемерид Луны на 27.07.2018 20:21 (начало координат в центре масс Солнечной системы)
*******************************************************************************
 Revised: Jul 31, 2013             Moon / (Earth)                           301
 
 GEOPHYSICAL DATA (updated 2018-Aug-13):
  Vol. Mean Radius, km  = 1737.53+-0.03    Mass, x10^22 kg       =    7.349
  Radius (gravity), km  = 1738.0           Surface emissivity    =    0.92
  Radius (IAU), km      = 1737.4           GM, km^3/s^2          = 4902.800066
  Density, g/cm^3       =    3.3437        GM 1-sigma, km^3/s^2  =  +-0.0001  
  V(1,0)                =   +0.21          Surface accel., m/s^2 =    1.62
  Earth/Moon mass ratio = 81.3005690769    Farside crust. thick. = ~80 - 90 km
  Mean crustal density  = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km 
  Heat flow, Apollo 15  = 3.1+-.6 mW/m^2   k2                    = 0.024059
  Heat flow, Apollo 17  = 2.2+-.5 mW/m^2   Rot. Rate, rad/s      = 0.0000026617
  Geometric Albedo      =    0.12

  Mean angular diameter = 31'05.2"         Orbit period          = 27.321582 d
  Obliquity to orbit    = 6.67 deg         Eccentricity          = 0.05490
  Semi-major axis, a    = 384400 km        Inclination           = 5.145 deg
  Mean motion, rad/s    = 2.6616995x10^-6  Nodal period          = 6798.38 d
  Apsidal period        = 3231.50 d        Mom. of inertia C/MR^2= 0.393142
  beta (C-A/B), x10^-4  = 6.310213         gamma (B-A/C), x10^-4 = 2.277317
 
                                 Perihelion  Aphelion    Mean
  Solar Constant (W/m^2)         1414+-7     1323+-7     1368+-7
  Maximum Planetary IR (W/m^2)   1314        1226        1268
  Minimum Planetary IR (W/m^2)      5.2         5.2         5.2
 *******************************************************************************
 
 
*******************************************************************************
Ephemeris / WWW_USER Wed Aug 15 21:19:24 2018 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Moon (301)                      {source: DE431mx}
Center body name: Solar System Barycenter (0)     {source: DE431mx}
Center-site name: BODY CENTER
*******************************************************************************
Start time      : A.D. 2018-Jul-27 20:21:00.0003 TDB
Stop  time      : A.D. 2018-Jul-28 20:21:00.0003 TDB
Step-size       : 0 steps
*******************************************************************************
Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii    : (undefined)                                                  
Output units    : AU-D                                                         
Output type     : GEOMETRIC cartesian states
Output format   : 3 (position, velocity, LT, range, range-rate)
Reference frame : ICRF/J2000.0                                                 
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch                 
*******************************************************************************
JDTDB
   X     Y     Z
   VX    VY    VZ
   LT    RG    RR
*******************************************************************************
$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB 
 X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05
 VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05
 LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05
$$EOE
*******************************************************************************
Coordinate system description:

  Ecliptic and Mean Equinox of Reference Epoch

    Reference epoch: J2000.0
    XY-plane: plane of the Earth's orbit at the reference epoch
              Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76)
    X-axis  : out along ascending node of instantaneous plane of the Earth's
              orbit and the Earth's mean equator at the reference epoch
    Z-axis  : perpendicular to the xy-plane in the directional (+ or -) sense
              of Earth's north pole at the reference epoch.

  Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]:

    JDTDB    Julian Day Number, Barycentric Dynamical Time
      X      X-component of position vector (au)                               
      Y      Y-component of position vector (au)                               
      Z      Z-component of position vector (au)                               
      VX     X-component of velocity vector (au/day)                           
      VY     Y-component of velocity vector (au/day)                           
      VZ     Z-component of velocity vector (au/day)                           
      LT     One-way down-leg Newtonian light-time (day)                       
      RG     Range; distance from coordinate center (au)                       
      RR     Range-rate; radial velocity wrt coord. center (au/day)            

Geometric states/elements have no aberrations applied.

 Computations by ...
     Solar System Dynamics Group, Horizons On-Line Ephemeris System
     4800 Oak Grove Drive, Jet Propulsion Laboratory
     Pasadena, CA  91109   USA
     Information: http://ssd.jpl.nasa.gov/
     Connect    : telnet://ssd.jpl.nasa.gov:6775  (via browser)
                  http://ssd.jpl.nasa.gov/?horizons
                  telnet ssd.jpl.nasa.gov 6775    (via command-line)
     Author     : Jon.D.Giorgini@jpl.nasa.gov
*******************************************************************************


$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB 
 X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05
 VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05
 LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05
$$EOE

Чудесно! Теперь необходимо слегка обработать полученные данные напильником.

6. 38 попугаев и одно попугайское крылышко


Для начала определимся с масштабом, ведь наши уравнения движения (5) записаны в безразмерной форме. Данные, предоставленные NASA сами подсказывают нам, что за масштаб координат стоит взять одну астрономическую единицу. Соответственно в качестве эталонного тела, к которому мы будем нормировать массы других тел мы возьмем Солнце, а в качестве масштаба времени — период обращения Земли вокруг Солнца.

Все это конечно очень хорошо, но мы не задали начальные условия для Солнца. «Зачем?» — спросил бы меня какой-нибудь лингвист. А я бы ответил, что Солнце отнюдь не неподвижно, а тоже вращается по своей орбите вокруг центра масс Солнечной системы. В этом можно убедится, взглянув на данные NASA для Солнца

$$SOE
2458327.347916670 = A.D. 2018-Jul-27 20:21:00.0003 TDB 
 X = 6.520050993518213E+04 Y = 1.049687363172734E+06 Z =-1.304404963058507E+04
 VX=-1.265326939350981E-02 VY= 5.853475278436883E-03 VZ= 3.136673455633667E-04
 LT= 3.508397935601254E+00 RG= 1.051791240756026E+06 RR= 5.053500842402456E-03
$$EOE

Взглянув на параметр RG мы увидим, что Солнце вращается вокруг барицентра Солнечной системы, и на 27.07.2018 центр звезды находится от него на расстоянии в миллион километров. Радиус Солнца, для справки — 696 тысяч километров. То есть барицентр Солнечной системы лежит в полумиллионе километров от поверхности светила. Почему? Да потому что все остальные тела, взаимодействующие с Солнцем так же сообщают ему ускорение, главным образом, конечно тяжеленький Юпитер. Соответственно у Солнца тоже есть своя орбита.

Мы конечно можем выбрать эти данные в качестве начальных условий, но нет — мы же решаем модельную задачу трех тел, и Юпитер и прочие персонажи в неё не входят. Так что в ущерб реализму, зная положение и скорости Земли и Луны мы пересчитаем начальные условия для Солнца, так, чтобы центр масс системы Солнце — Земля — Луна находился в начале координат. Для центра масс нашей механической системы справедливо уравнение

$(m_1 + m_2 + m_3) \, \vec r_C = m_1 \, \vec r_1 + m_2 \, \vec r_2 + m_3 \, \vec r_3$



Поместим центр масс в начало координат, то есть зададимся $\vec r_C = 0$, тогда

$m_1 \, \vec r_1 + m_2 \, \vec r_2 + m_3 \, \vec r_3 = 0$


откуда

$\begin{align} & m_3 \, \vec r_3 = -m_1 \, \vec r_1 - m_2 \, \vec r_2 \\ & \vec r_3 = - \frac{m_1}{m_3} \vec r_1 - \frac{m_2}{m_3} \, \vec r_2 \end{align} $


Перейдем к безразмерным координатам и параметрам, выбрав $\mu = \mu_3$

$\vec \xi_3 = -\varkappa_1 \vec \xi_1 -\varkappa_2 \vec \xi_2\quad\quad(6)$


Дифференцируя (6) по времени и переходя к безразмерному времени получаем и соотношение для скоростей

$\vec u_3 = -\varkappa_1 \, \vec u_1 -\varkappa_2 \, \vec u_2$


где $\vec u_i = \cfrac{d\vec \xi_i}{d\tau}, \forall i=\overline{1,3}$

Теперь напишем программу, которая сформирует начальные условия в выбранных нами «попугаях». На чем будем писать? Конечно же на Питоне! Ведь, как известно, это самый лучший язык для математического моделирования.

Однако, если уйти от сарказма, то мы действительно попробуем для этой цели питон, а почему нет? Я обязательно приведу ссылку на весь код в моем профиле Github.

Расчет начальных условий для системы Луна - Земля - Солнце
#
# Исходные данные задачи
#

# Гравитационная постоянная
G = 6.67e-11

# Массы тел (Луна, Земля, Солнце)
m = [7.349e22, 5.792e24, 1.989e30]

# Расчитываем гравитационные параметры тел
mu = []

print("Гравитационные параметры тел")

for i, mass in enumerate(m):
    mu.append(G * mass)
    print("mu[" + str(i) + "] = " + str(mu[i]))

# Нормируем гравитационные параметры к Солнцу
kappa = []

print("Нормированные гравитационные параметры")

for i, gp in enumerate(mu):
    kappa.append(gp / mu[2])
    print("xi[" + str(i) + "] = " + str(kappa[i]))

print("\n")

# Астрономическая единица
a = 1.495978707e11

import math

# Масштаб безразмерного времени, c
T = 2 * math.pi * a * math.sqrt(a / mu[2])

print("Масштаб времени T = " + str(T) + "\n")

# Координаты NASA для Луны
xL = 5.771034756256845E-01
yL = -8.321193799697072E-01
zL = -4.855790760378579E-05

import numpy as np

xi_10 = np.array([xL, yL, zL])
print("Начальное положение Луны, а.е.: " + str(xi_10))

# Координаты NASA для Земли
xE = 5.755663665315949E-01
yE = -8.298818915224488E-01
zE = -5.366994499016168E-05

xi_20 = np.array([xE, yE, zE])
print("Начальное положение Земли, а.е.: " + str(xi_20))

# Расчитываем начальное положение Солнца, полагая что начало координат - в центре масс всей системы
xi_30 = - kappa[0] * xi_10 - kappa[1] * xi_20
print("Начальное положение Солнца, а.е.: " + str(xi_30))

# Вводим константы для вычисления безразмерных скоростей
Td = 86400.0
u = math.sqrt(mu[2] / a) / 2 / math.pi

print("\n")

# Начальная скорость Луны
vxL = 1.434571674368357E-02
vyL = 9.997686898668805E-03
vzL = -5.149408819470315E-05

vL0 = np.array([vxL, vyL, vzL])
uL0 = np.array([0.0, 0.0, 0.0])

for i, v in enumerate(vL0):
    vL0[i] = v * a / Td
    uL0[i] = vL0[i] / u

print("Начальная скорость Луны, м/с: " + str(vL0))
print(" -//- безразмерная: " + str(uL0))

# Начальная скорость Земли
vxE = 1.388633512282171E-02
vyE = 9.678934168415631E-03
vzE = 3.429889230737491E-07

vE0 = np.array([vxE, vyE, vzE])
uE0 = np.array([0.0, 0.0, 0.0])

for i, v in enumerate(vE0):
    vE0[i] = v * a / Td
    uE0[i] = vE0[i] / u

print("Начальная скорость Земли, м/с: " + str(vE0))
print(" -//- безразмерная: " + str(uE0))

# Начальная скорость Солнца
vS0 = - kappa[0] * vL0 - kappa[1] * vE0
uS0 = - kappa[0] * uL0 - kappa[1] * uE0

print("Начальная скорость Солнца, м/с: " + str(vS0))
print(" -//- безразмерная: " + str(uS0))


Выхлоп программы

Гравитационные параметры тел
mu[0] = 4901783000000.0
mu[1] = 386326400000000.0
mu[2] = 1.326663e+20
Нормированные гравитационные параметры
xi[0] = 3.6948215183509304e-08
xi[1] = 2.912016088486677e-06
xi[2] = 1.0


Масштаб времени T = 31563683.35432583

Начальное положение Луны, а.е.: [ 5.77103476e-01 -8.32119380e-01 -4.85579076e-05]
Начальное положение Земли, а.е.: [ 5.75566367e-01 -8.29881892e-01 -5.36699450e-05]
Начальное положение Солнца, а.е.: [-1.69738146e-06  2.44737475e-06  1.58081871e-10]


Начальная скорость Луны, м/с: [24838.98933473 17310.56333294   -89.15979106]
 -//- безразмерная: [ 5.24078311  3.65235907 -0.01881184]
Начальная скорость Земли, м/с: [2.40435899e+04 1.67586567e+04 5.93870516e-01]
 -//- безразмерная: [5.07296163e+00 3.53591219e+00 1.25300854e-04]
Начальная скорость Солнца, м/с: [-7.09330769e-02 -4.94410725e-02  1.56493465e-06]
 -//- безразмерная: [-1.49661835e-05 -1.04315813e-05  3.30185861e-10]

7. Интегрирование уравнений движения и анализ результатов


Собственно само интегрирование сводится к более-менее стандартной для SciPy процедуре подготовки системы уравнений: преобразованию системы ОДУ к форме Коши и вызову соответствующих функций-решателей. Для преобразования системы к форме Коши вспоминаем, что

$ \vec u_i = \frac{d\vec \xi_i}{d\tau}, \forall i=\overline{1,3}\quad\quad(7) $


Тогда введя вектор состояния системы

$\vec y = \left[\vec\xi_1, \vec\xi_2, \vec\xi_1, \vec u_1, \vec u_2, \vec u_3 \right]^T$


сводим (7) и (5) к одному векторному уравнению

$\frac{d\vec y}{d\tau} = \vec f(\tau, \vec y)\quad\quad(8)$


Для интегрирования (8) с имеющимися начальными условиями напишем немного, совсем немного кода

Интегрирования уравнений движения в задаче трех тел
#
#   Вычисление векторов обобщенных ускорений
#
def calcAccels(xi):
    k = 4 * math.pi ** 2

    xi12 = xi[1] - xi[0]
    xi13 = xi[2] - xi[0]
    xi23 = xi[2] - xi[1]

    s12 = math.sqrt(np.dot(xi12, xi12))
    s13 = math.sqrt(np.dot(xi13, xi13))
    s23 = math.sqrt(np.dot(xi23, xi23))

    a1 = (k * kappa[1] / s12 ** 3) * xi12 + (k * kappa[2] / s13 ** 3) * xi13
    a2 = -(k * kappa[0] / s12 ** 3) * xi12 + (k * kappa[2] / s23 ** 3) * xi23
    a3 = -(k * kappa[0] / s13 ** 3) * xi13 - (k * kappa[1] / s23 ** 3) * xi23

    return [a1, a2, a3]


#
#   Система уравнений в нормальной форме Коши
#
def f(t, y):
    n = 9

    dydt = np.zeros((2 * n))

    for i in range(0, n):
        dydt[i] = y[i + n]

    xi1 = np.array(y[0:3])
    xi2 = np.array(y[3:6])
    xi3 = np.array(y[6:9])

    accels = calcAccels([xi1, xi2, xi3])

    i = n
    for accel in accels:
        for a in accel:
            dydt[i] = a
            i = i + 1

    return dydt

# Начальные условия задачи Коши
y0 = [xi_10[0], xi_10[1], xi_10[2],
      xi_20[0], xi_20[1], xi_20[2],
      xi_30[0], xi_30[1], xi_30[2],
      uL0[0], uL0[1], uL0[2],
      uE0[0], uE0[1], uE0[2],
      uS0[0], uS0[1], uS0[2]]

#
# Интегрируем уравнения движения
#

# Начальное время
t_begin = 0
# Конечное время
t_end = 30.7 * Td / T;
# Интересующее нас число точек траектории
N_plots = 1000
# Шаг времени между точкими
step = (t_end - t_begin) / N_plots

import scipy.integrate as spi

solver = spi.ode(f)

solver.set_integrator('vode', nsteps=50000, method='bdf', max_step=1e-6, rtol=1e-12)
solver.set_initial_value(y0, t_begin)

ts = []
ys = []
i = 0

while solver.successful() and solver.t <= t_end:
    solver.integrate(solver.t + step)
    ts.append(solver.t)
    ys.append(solver.y)
    print(ts[i], ys[i])
    i = i + 1


Посмотрим что у нас получилось. Получилась пространственная траектория Луны на первые 29 суток от выбранной нами начальной точки


а так же её проекция в плоскость эклиптики.


«Эй, дядя, что ты нам впариваешь?! Это же окружность!».

Во-первых, таки не окружность — заметно смещение проекции траектории от начала координат вправо и вниз. Во-вторых — ничего не замечаете? Не, ну правда?


Обещаю подготовить обоснование того (на основе анализа погрешностей счета и данных NASA), что полученное смещение траектории не есть следствие ошибок интегрирования. Пока предлагаю читателю поверить мне на слово — это смещение есть следствие солнечного возмущения лунной траектории. Крутанем-ка еще один оборот



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

Можно сделать вывод, что солнечная гравитация влияет на орбиту Луны достаточно существенно — старушка не ходит по небу дважды одним и тем же путём. Картинка за полгода движения позволяет (по крайней мере качественно) убедится в этом (картинка кликабельна)

image

Интересно? Ещё бы. Астрономия вообще наука занятная.

Постскриптум


В вузе, где я учился и работал без малого семь лет — Новочеркасском политехе — ежегодно проводилась зональная олимпиада студентов по теоретической механике вузов Северного Кавказа. Трижды мы принимали и Всероссийскую олимпиаду. На открытии, наш главный «олимпиец», профессор Кондратенко А.И., всегда говорил: «Академик Крылов называл механику поэзией точных наук».

Я люблю механику. Всё то хорошее, чего я добился в своей жизни и карьере произошло благодаря этой науке и моим замечательным учителям. Я уважаю механику.

Поэтому, я никогда не позволю издеваться над этой наукой и нагло эксплуатировать её в своих целях никому, будь он хоть трижды доктор наук и четырежды лингвист, и разработал хоть миллион учебных программ. Я искренне считаю, что написание статей на популярном публичном ресурсе должно предусматривать их тщательную вычитку, нормальное оформление (формулы LaTeX — это не блажь разработчиков ресурса!) и отсутствие ошибок, приводящих к результатам нарушающим законы природы. Последнее вообще «маст хэв».

Я часто говорю своим студентам: «компьютер освобождает ваши руки, но это не значит, что при этом нужно отключать и мозг».

Ценить и уважать механику я призываю и вас, мои уважаемые читатели. Охотно отвечу на любые вопросы, а исходный текст примера решения задачи трех тел на языке Python, как и обещал, выкладываю в своем профиле Github.

Спасибо за внимание!
Поделиться публикацией
Комментарии 73
    –6
    если бы луна притягивалась к солнцу с силой большей чем притягивается к земле она бы уже притянулась))… или к земле бы притянулась… очевидно же они отталкиваются
      –15
      Такие публикации вызывают уныние. Претензий много в виде переписных из разных источников прописных истин. Ничего своего, кроме критики чужого, и огромное желание что-то доказать.Жалко Ваших студентов за внешним лоском пустота.
        +7
        Вы, Юрий Карлович, невежа. Причем, как я вижу — невежа воинствующий. Посмотрите сами, на то, какое дерьмо льете на ресурс да почитайте школьный учебник физики, для начала. А уж потом критикуйте на здоровье. А то я что-то не пойму, кто вы — механик-недучка, или лингвист-недомеханик.
          +1
          Никогда не думал, что классическая небесная механика может вызвать такое оживление в комментариях — и тут на тебе! Представляю заголовки завтрашних газет: Нью-Йорк Таймз, «Ученые из Днепра и Ростова ставят под сомнение классическую механику (друг друга)» и приписку более мелким шрифтом: «Мировое научное сообщество замерло в ожидании очередного комментария на хабре коммита в ipython notebook».
            0

            Дмитрий, только два лакмусовых вопроса: "Были американцы на Луне или нет? И если были то как вернулись обратно?" Катющик, кстати топит за отталкивание планет...

              +1
              В статье я уже на этот вопрос ответил
              Мы живем в мире, где человек ходил по Луне

              Я предлагаю сходить за этими данными к тем, кто собственно ходил по Луне, к NASA, а именно в Лабораторию реактивного движения, Пасадена, штат Калифорния.
            +3
            Претензий много в виде переписных из разных источников прописных истин

            По тексту, я где-то претендую на научную новизну? Да, я излагаю прописные истины, взятые из разных источников. Вы правы, черт возьми, и я не вижу в этом ничего криминального.
            огромное желание что-то доказать

            Да, доказать, что Луна вращается вокруг Земли не на палке с постоянной угловой скоростью независящей от расстояния до Земли, как это написано у Вас, а под действием сил гравитационного притяжения. Только и всего
            –16
            Да, кроме хамства от Вас услышать что-то другое было бы новостью. А вот школьный учебник физики- это предел Вашего образования. Жаль.
              +1
              А велико ли влияние на Землю и Луну газовых гигантов (в первую очередь, Сатурна и Юпитера)?
                +7
                Если подходить к вопросу с точки зрения критерия Лапласа, который приведен в статье — влияние незначительно, так как Луна далеко за пределами сфер влияния этих планет (особенно Сатурна).
                Однако Юпитер — та ещё сволочь — любит засасывать кометы и прочее, неосторожно пролетающее мимо космическое хозяйство. Но одно дело комета, движущаяся по вытянутому эллипсу и, в районе Юпитера имеющая существенно меньшую гелиоцентрическую скорость, чем в перигелии. Другое дело Луна, располагающаяся значительно дальше, внутри гравитационной сферы влияния Земли и имеющая гелиоцентрическую скорость сравнимую со скоростью своей хозяйки. Так что мой ответ — нет, для простых практических расчетов влиянием планет-гигантов можно пренебречь. Более развернутый ответ требует дополнительного исследования вопроса
                  +1
                  Спасибо! (плюсануть не могу)
                +3
                кроме хамства от Вас

                Расшаркиваться перед Вами я не собираюсь, зарубите себе это на носу. Уважения Вы не заслуживаете хотя бы по той причине, что у Вас хватает совести держать неисправленной статью с кошмарной для Вашего уровня образования ошибкой. Ваша сладкая жизнь на этом ресурсе закончилась — будьте уверены, стоит Вам оступится, я пну Вас дополнительно, для придачи большего ускорения. Где там очередной шедевр, про «необъявленную войну в атмосфере Земли»? Что так быстро у черновики-то убрали, а, профессор? Ждем же ж с нетерпением, попкорн уже закупили.

                Вы смешной дядька, Юрий Карлович)

                  +5
                  Статья в кэше у Гугла есть, если что. :)
                    +2
                    И на сохабре тоже.
                      +1
                      Из первых комментариев понятно, почему сей шедевр ушел в черновики.
                      Собственно спасибо за ссылки, и подтверждение того, что мне это не померещилось
                  +3
                  А как далеко Вы моделировали? При ограничении тремя телами и точностью используемых методов интегрирования СДУ — нет ли такого отрезка траектории, где Луна цепляет Землю или хотя бы атмосферу?
                    +1
                    А как далеко Вы моделировали?

                    В общем-то не далеко. Но вот результат на период сароса (18 лет). С Землей явно не сталкивается, но ширина «бублика» довольно подозрительна, видимо из-за накапливаемой погрешности интегрирования
                      +1
                      Прошу прощения за ляпсус — бублик образовался из-за того, что на все 18 лет я для графика взял только 1000 точек траектории, если присмотритесь — предыдущий график ломаный.
                      Собственно вот новый на тот же период но уже с 10000 точек и максимальным шагом интегрирования 1e-5

                        0
                        Тот же период, макс. шаг 1e-6, 1e5 точек на графике
                          +3
                          Просто оставлю здесь. Для более-менее точных расчетов в декартовой системе координат вам нужно порядка 1000 точек на виток.

                          За 18 лет 18*365.25/27.321661 = 2.40633247E+002 Почти 241 оборот. Т.е. всего точек надо ~2.4E+05
                            0
                            Спасибо, учту
                              +2

                              А разве методом Булирша-Штёра не получится интегрировать с шагом порядка суток? Траектория же гладкая, методом высокого порядка должна хорошо браться. Хотя вычислений силы там порядка тысячи на виток и потребуется, конечно.

                                +2
                                методом Булирша-Штёра

                                Я попробую. Тем более что есть желание прогнать эту задачку с хорошей точность на плюсах
                                  0
                                  Даже не знал о таком методе. Надо записать в таски на наш комплекс.
                                +1
                                Дмитрий, можно ли свести полученную систему каким-нибудь преобразованием к динамической системе с квадратичной правой частью по фазовым координатам?
                                  0
                                  Я думаю, свести правую часть этой системы к квадратичной форме по каким бы то ни было координатам не выйдет
                                    0
                                    {тер-мех забыл уже} смотреть Гамильтониан
                                      0
                                      Функцию Гамильтона для такой системы по сути полная механическая её энергия. Ну так если её составить, то квадратичной формы относительно канонических координат там не будет.

                                      Квадратичной формой относительно обобщенных скоростей является кинетическая энергия, причем всегда. А вот потенциальная энергия от выбранных обобщенных координат может зависеть совершенно по разному
                                +2
                                Ну что же, 18 лет можно за Луну не волноваться. Да я думаю, и 180 лет она нас не зацепит. Интересно, пожалуй, на тысячу или миллион лет вперёд заглянуть. Как-то, вообще, можно понять: это устойчивые колебания или некий переходный процесс?
                                  0
                                  Как-то, вообще, можно понять: это устойчивые колебания или некий переходный процесс?

                                  Применить второй метод Ляпунова для исследования устойчивости её движения
                                    0
                                    Уже все посчитали)
                                    +2

                                    Не очень удобный график, на самом деле. Не видно — это орбита сжимается или прецессирует. Его бы дополнить графиками для оскулирующих большой полуоси орбиты, эксцетриситета и наклонения.

                                      +1
                                      Хорошая идея для другой публикации
                                  +3
                                  Вы вбили мне клин в голову, пришлось перепроверять формулы по ресурсам… В общем спасибо за статью и крылатое выражение (про освобождённые руки...) с вашего позволения скомуниздю :) моим ученикам тоже надо.
                                    +1
                                    ИМХО для интегрирования астрономических задач на больших интервалах времени можно использовать специализированные методы высокого порядка точности. И, возможно, изменить уравнения (полу-аналитическое решение ...).
                                      +3
                                      Получается Луна притягивается к Солнцу силой, более чем в два раза превышающей силу её притяжения к Земле.
                                      Подобное возмущение уже нельзя не учитывать и оно определенно повлияет на конечную траекторию движения Луны.
                                      Замечание интересное, но ведь дело не в том, что сильней притягивает Луну: Солнце или Земля, а на сколько отличается ускорение свободного падения (на Солнце) в разных точках траектории Луны и Земли.
                                      Т.е. если поле силы тяжести Солнца было бы на пути Луны и Земли однородным (Солнце было бы бесконечно далеко), то влияния Солнца на орбиту Луны относительно Земли не было бы.

                                      Спасибо, хорошая статья!
                                        0
                                        Т.е. если поле силы тяжести Солнца было бы на пути Луны и Земли однородным (Солнце было бы бесконечно далеко), то влияния Солнца на орбиту Луны относительно Земли не было бы.

                                        Со всей очевидностью вы путаете потенциал грав поля и силу тяжести.
                                        Во-первых наличие любой силы меняет движение тела. Просто замените в программе автора возмущение от Солнца постоянной в декартовых координатах силой.

                                        Во-вторых, если поле однородно (потенциал постоянен), то сила равна нулю и действительно никакого эффекта нет.
                                          0
                                          Про потенциал я ничего не писал. Однородное поле — поле, в котором сила, действующая на тело, одинакова в любой точке этого поля. Потенциал при этом будет разный в разных точках (при движении вдоль линии действия силы он будет линейно уменьшаться/увеличиваться, эквипотенциальная поверхность будет плоскостью, перпендикулярной этой линии).
                                          В таком поле Земля и Луна будут получать одинаковое ускорение от силы притяжения к Солнцу и траектория Луны относительно Земли (т.е. в системе координат, связанных с Землёй (она получается не инерционная)) будет такой же, как если бы этого однородного поля не было.
                                          Так же как «невесомость» в космических аппаратах возникает потому, что в границах КА поле тяготения Земли можно считать однородным. Если габариты КА большие, то «невесомость» не будет столь полной (например, вытянутый КА будет «вставать» вертикально)
                                            0
                                            Забавная штука, когда в реальную физику движения вносится нечто абстрактно-оторванное: я прав и не прав одновременно.

                                            Я не прав, т.к. действительно вы не говорили о грав поле, но прав т.к. если сила (F) одна, а массы-то существенно разные, а значит вот это:
                                            В таком поле Земля и Луна будут получать одинаковое ускорение от силы притяжения к Солнцу

                                            неверно.
                                              0
                                              прав т.к. если сила (F) одна, а массы-то существенно разные

                                              Как же сила-то будет одна, если она равна напряжённости поля, умноженной на гравитационную массу, а массы-то существенно разные.

                                                0
                                                Вы попали в туже ловушку, что и я. Nokta_strigo предлагает умозрительно приравнять силу действия Солнца на Землю и Луну.

                                                Да, это не физично, согласен с вами. Но в качестве разминки для мозгов можно и заняться.
                                                  +1
                                                  Вы попали в туже ловушку, что и я. Nokta_strigo предлагает умозрительно приравнять силу действия Солнца на Землю и Луну.
                                                  Неееет!!! ;-) Я говорю про однородное гравитационное поле, где сила тяжести не зависит от положения тела, но это не значит, что она не зависит от массы тела. Т.е. постоянна напряжённость поля (или сила, действующая на пробную частицу, или ускорение свободного падения), а сила — пропорциональна массе тела.
                                                  Так что силы, действующие на Землю и Луну будут пропорциональны массам тел (разные), а ускорения — одинаковы.
                                                  Собственно, понятие однородного поля силы вполне устоявшееся (в теории поля и в классической механике), я ничего нового не придумал.
                                                    0
                                                    Прошу простить, что-то на меня нашло.
                                                0
                                                Ускорение материальной точки в гравитационном поле (именно в гравитационном!) не зависит от массы этой точки. По этому есть понятие «ускорение свободного падения», которое одинаково для тел разной массы.
                                                Так получается из-за того, что масса как мера инерции и масса как характеристика гравитационного взаимодействия точки во всех известных физике случаях совпадают с высокой точностью. Если расписать выражение для ускорения тела под действие силы тяжести, то масса тела будет и в числителе, и знаменателе и, соответственно, сократится.
                                          +1

                                          Я давно думал о моделировании n>2 тел, но не хватает опыта в астрономии.
                                          Что, если не вычислять притяжение каждого тела к каждому, а искривлять пространство? Алгоритм таков: создается доп. измерение. Дальше в нем подвешываются тела. Затем соответственно массам прогибается это измерение по формуле m*G под каждым телом (в разрезе этот прогиб выглядит как х^2). После расчета прогиба у нас получается некая кривая поверхность (типа простыни). В финале нужно просто катать тела по этой поверхности.
                                          П.С. для упрощения расчетов, прогибать измерение нужно только под телами, в пустом пространстве это можно игнорировать.

                                            0
                                            Ну так на каждом шаге придётся вычислять не просто n*(n-1)/2 сил, а пересчитывать искривление пространства под воздействием каждого тела и, потом, воздействие новой кривизны пространства на каждое тело. Т.е. «просто катать тела» можно только если эти тела будут вызывать пренебрежимо малое искривление «простыни», а большие тела будут зафиксированы.
                                            Задача с одним массивным телом и любым количеством «лёгких» имеет аналитическое решение. А если массивных тел больше одного, то они будут двигать друг друга и придётся на каждом шаге пересчитывать искривление.
                                            +1
                                            Спасибо. Очень интересно.
                                              0
                                              На сайте "Наброски для новой физики" в статье БИРЮЛЬКИ И ФИТЮЛЬКИ ВСЕМИРНОГО ТЯГОТЕНИЯ ) всё описано гораздо веселее:
                                              И – вот они, реалии: несмотря на то, что у Луны имеется собственное тяготение (действующее, впрочем, лишь в небольшой окололунной области), Луна движется в сфере тяготения Земли как болванка, которая Землю не притягивает. Ибо если бы Луна Землю притягивала, у Земли была бы полноценная динамическая реакция, т.е. обращение около общего центра масс. Но, как мы видели, вместо этого обращения имеют место одномерные колебания.
                                                0

                                                В каком смысле веселее?
                                                Насколько я помню, барицентр довольно далеко от центра Земли, возможно наибольшее расстояние во всей сол. системы

                                                  +1

                                                  В системе Плутон-Харон больше. Но там отношение масс 1/4 а не 1/81 как в системе

                                                    0
                                                    Веселее — в том смысле, что автор «Бирюлек...» утверждает, что Земля не вращается вокруг барицентра системы Земля-Луна, а совершает лишь колебания вдоль своей орбиты вокруг Солнца. Т.е. колебания к Солнцу/от Солнца регистрировать не удаётся, хотя техника, вроде, позволяет. Ну, или как-то так.
                                                      0
                                                      И как это проверить? Дальномером от Солнца? Сложно представить. А если и выйдет, то точность будет порядка самого расстояния к барицентру.
                                                      Получается, автор «Бирюлек...» утверждает, что Земля двигается по рельсам?
                                                        0
                                                        И как это проверить? Дальномером от Солнца? Сложно представить. А если и выйдет, то точность будет порядка самого расстояния к барицентру.

                                                        Не дальномером, и не по рельсам, но смысл не в этом. Вы сами «Бирюльки...» прочитайте, чтобы не обсуждать их голословно. Весело — не то слово. У автора там и гравитация не веществом вызывается, если что )
                                                    +1
                                                    БИРЮЛЬКИ И ФИТЮЛЬКИ ВСЕМИРНОГО ТЯГОТЕНИЯ

                                                    Выражаясь словами героев Стругацких — «всё это, конечно очень блаародно..», но как насчет 60 с лишним лет космических полетов? Аппараты летят туда, куда мы их направляем, руководствуясь неверной теорией? Сомневаюсь.

                                                    В упомянутой статье много рассуждизмов, но мало фактического материала в виде конкретных данных и количественных соотношений. Насчет астероида Эрос у меня два замечания:

                                                    1. Согласно приведенному выше критерию Лапласа, максимальный радиус сферы гравитационного влияния этого астероида — 432 километра. Вторая космическая скорость у поверхности — 10,3 м/с или около 36 километров в час. Солнце просто нагло сдирает всё, что пытается выйти на устойчивую орбиту вокруг него.
                                                    2. Вот эта каменюка имеется в виду?


                                                    Ничего, что это скорее картофелина, выросшая в сухой земле, чем тело формы близкой к сферической? Ничего что гравитационное поле вокруг неё точно не центральное?

                                                    Ну так пусть тот умник, что написал «бирюльки...» попробует в таких условиях поуправлять космическим аппаратом, я на это с удовольствием посмотрю
                                                      0
                                                      Ничего, что это скорее картофелина, выросшая в сухой земле, чем тело формы близкой к сферической? Ничего что гравитационное поле вокруг неё точно не центральное?
                                                      Так автор и говорит о том, что гравитационного поля там как раз и вовсе нет, и что вещество гравитации подчиняется, но никак не является для неё причиной.
                                                        0
                                                        вещество гравитации подчиняется, но никак не является для неё причиной

                                                        Автор так пишет и попадает пальцем в небо — природу гравитации до сих пор не могут понять, так и что теперь?
                                                        Да, закон всемирного тяготения — теория, основанная на концепции дальнодействия, очень «нефизична» по своей сути, но это то хорошее допущение, которое позволило человечеству добраться до Плутона включительно и уже 60 лет нормально подтверждается решением задач практической космонавтики.
                                                          0
                                                          Да, закон всемирного тяготения — теория, основанная на концепции дальнодействия, очень «нефизична» по своей сути, но это то хорошее допущение, которое позволило человечеству добраться до Плутона включительно и уже 60 лет нормально подтверждается решением задач практической космонавтики.

                                                          Так автор, как я понял, и не спорит с этим утверждением, просто приводит свои доказательства достаточно своеобразной теории о связи материи и тяготения.
                                                            0
                                                            Доказательства — это если бы он взял какой-нибудь вышедший из строя спутник, предсказал бы с помощью своей теории его траекторию на месяц или год вперёд, сопоставил бы с показаниями Stellarium и через этот самый месяц или год призвал всех посмотреть в телескоп, чтобы убедиться, что Stellarium давал неверные координаты этого спутника, а вот по теории Гришаева всё сходится.
                                                            Но так доказать не получится, потому что Stellarium даст правильные координаты, а теория Гришаева — нет.
                                                            Поэтому надо выкручиваться, придумывая заговор космических агентств о скрытии отсутствия следов человека притяжения на Луне.
                                                              0
                                                              Доказательства — это если бы он взял...

                                                              Всё верно, но я немного перефразирую Вас, коллега: доказательства — это фактический материал основанный на прямом (или косвенном, но основанном на проверенной гипотезе) измерении количественных показателей.

                                                              Нет цифр, нет доказательств. Иное — просто измышления
                                                                0
                                                                Доказательства — это если бы он взял какой-нибудь вышедший из строя спутник...
                                                                Спутник чего? Земли? Оно ничего про них не говорит, как бы. В пределах действия земной гравитации всё ровно.
                                                                Поэтому надо выкручиваться, придумывая заговор космических агентств о скрытии отсутствия притяжения на Луне.
                                                                У Гришаева не написано, что на Луне нет притяжения.

                                                                Могу ошибаться, но для подтверждения его теории нет смысла лететь в космос, достаточно поглубже погрузиться в недра Земли. Если он прав, то с погружением вес тела будет постоянно расти. Такое вот обывательское мнение…
                                                      0
                                                      И насколько это точно? Обычно уравнение Кеплера решается.
                                                        +2
                                                        Обычно уравнение Кеплера решается

                                                        Уравнение Кеплера обычно решатся для задачи двух тел. Уравнение Кеплера проистекает из допущения о движении точки в центральном поле тяготения. Уравнение Кеплер — готовое решение другой задачи. В случае с Луной поле тяготения не центрально, а определяется силами притяжения к Земле и Солнцу.
                                                        И насколько это точно?

                                                        В моем примере не слишком точно. На деле, используя подходящие методы решения ОДУ и соответствующую дискретизацию задачи, можно получить решение с любой практически требуемой точностью
                                                          0
                                                          Уравнение Кеплера обычно решатся для задачи двух тел. Уравнение Кеплера проистекает из допущения о движении точки в центральном поле тяготения. Уравнение Кеплер — готовое решение другой задачи. В случае с Луной поле тяготения не центрально, а определяется силами притяжения к Земле и Солнцу.
                                                          Это исправляется добавлением пертурбаций. Я видел формулы, которые таким образом дают точность до 1-2 арк минуты.
                                                            +1
                                                            Имеете в виду метод оскулирующих элементов?
                                                              0
                                                              Ага
                                                                0
                                                                Если речь про Луну, то, скорее всего, речь идет про элементы Делоне.
                                                          0
                                                          Для галактик во избежание расчета парных сил уже очень давно используются методы частиц в ячейках.
                                                            0
                                                            Просто дополню, что вы незаслуженно принижаете важность той самой как-бы сферы радиусом 920 тыс. км. Это граница нашей гравитационной ямы, на которой эффективная (с учетом притяжения Солнца) вторая космическая равна нулю.
                                                            Если вы нашли булыжник внутри этой сферы, то, возможно, он является спутником Земли и будет сопровождать ее примерно вечно.
                                                            Если булыжник вне этой сферы, то у него только один шанс остаться с нами навсегда — долбануться об Землю или Луну. Во всех остальных случаях он просвистит мимо, какая бы изначальная скорость у него ни была.
                                                            И в обратную сторону: если хоть как-то уползти от родной планеты на эти самые 925 000 км, хоть из пушки, хоть на лифте, хоть с помощью солнечного паруса — дальше вас унесет далеко-далеко и вы станете спутником Солнца (опять же за исключением варианта долбануться о поверхность).
                                                            Ну и о последнем вы сказали — любые маневры по изменению орбиты вокруг Земли имеют смысл только внутри этой сферы. По той причине, что снаружи никаких околоземных орбит нет.
                                                              0
                                                              силой, направленной вдоль прямой, соединяющей центры рассматриваемых небесных тел
                                                              прямая должна соединять не центры тел, а центры масс.
                                                                +1
                                                                Вы правы отчасти. Согласно формулировке закона, указанное выражение для модуля силы справедливо для точечных масс, вообще. Если взять тела, масса которых распределена в соответствии с центральной симметрией — шаров, то модуль равнодействующей силы притяжения так же будет соответствовать формулировке закона, а центры масс таких тел лежат в их геометрическом центре симметрии.

                                                                Для протяженных не сферических тел F ~ 1/r^2 уже не работает
                                                                0
                                                                Спасибо большое за интересную статью
                                                                  0
                                                                  К серии статей по «точному» численному решению задач хаотической динамики
                                                                  github.com/Libbum/cns-rs

                                                                Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                                                                Самое читаемое