Хабрахабр

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

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

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

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

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

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

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

где m1, m2 — массы, соответственно Луны и Земли; G = 6,67e-11 — гравитационная постоянная; 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}$, м — среднее расстояние от Земли до Солнца. Подставим в эту формулу реальные параметры

99 \cdot 10^{30}}{5. $\frac{F_{1,3}}{F_{1,2}} = \frac{1. 844\cdot10^{8}}{1. 98\cdot10^{24}} \, \left( \frac{3. 844\cdot10^{8}}\right)^2 = 2. 496\cdot10^{11} + 3. 19$

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

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

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

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

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

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

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

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

Лапласом. Оценим возмущение, вносимое Солнцем в траекторию Луны относительно Земли, а так же возмущение, вносимое Землей в траекторию Луны относительно Солнца, воспользовавшись критерием, предложенным П. Рассмотрим три тела: Солнце (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 — расстояние между центрами тел. В нашем случае

496\cdot10^{11} \, \left(\frac{5. $r = a \, \left(\frac{m_2}{m_3} \right)^{\frac{2}{5}} = 1. 99\cdot10^{30}} \right)^{\frac{2}{5}} = 925000,\, км$ 98\cdot10^{24}}{1.

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

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

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

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

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

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

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

Все эти величины, в силу законов механики, связаны соотношением
Итак, введем некое абстрактное небесное тело с гравитационным параметром $\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) дает нам девять скалярных уравнений второго порядка, что при переходе к нормальной форме Коши повышает порядок системы ещё в 2 раза. Как я уже писал ранее, прежде чем начинать численное интегрирование, следует озаботится расчетом начальных условий для решаемой задачи. Где мы возьмем данные о положении интересующих нас небесных тел? То есть нам необходимо рассчитать целых 18 параметров — начальные положения и компоненты начальной скорости всех точек системы. Мы живем в мире, где человек ходил по Луне — естественно человечество должно обладать информацией, как эта самая Луна движется и где она находится.

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

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

Полный вывод для эфемерид Луны на 27.07.2018 20:21 (начало координат в центре Земли)

*******************************************************************************
Revised: Jul 31, 2013 Moon / (Earth) 301

Mean Radius, km = 1737. GEOPHYSICAL DATA (updated 2018-Aug-13):
Vol. 03 Mass, x10^22 kg = 7. 53+-0. 0 Surface emissivity = 0. 349
Radius (gravity), km = 1738. 4 GM, km^3/s^2 = 4902. 92
Radius (IAU), km = 1737. 3437 GM 1-sigma, km^3/s^2 = +-0. 800066
Density, g/cm^3 = 3. 21 Surface accel., m/s^2 = 1. 0001
V(1,0) = +0. 3005690769 Farside crust. 62
Earth/Moon mass ratio = 81. = ~80 - 90 km
Mean crustal density = 2. thick. 07 g/cm^3 Nearside crust. 97+-. 1+-. thick.= 58+-8 km
Heat flow, Apollo 15 = 3. 024059
Heat flow, Apollo 17 = 2. 6 mW/m^2 k2 = 0. 5 mW/m^2 Rot. 2+-. 0000026617
Geometric Albedo = 0. Rate, rad/s = 0. 12

2" Orbit period = 27. Mean angular diameter = 31'05. 67 deg Eccentricity = 0. 321582 d
Obliquity to orbit = 6. 145 deg
Mean motion, rad/s = 2. 05490
Semi-major axis, a = 384400 km Inclination = 5. 38 d
Apsidal period = 3231. 6616995x10^-6 Nodal period = 6798. of inertia C/MR^2= 0. 50 d Mom. 310213 gamma (B-A/C), x10^-4 = 2. 393142
beta (C-A/B), x10^-4 = 6. 277317

2 5. 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
******************************************************************************* 2 5.

D. *******************************************************************************
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. 0003 TDB
Stop time : A. 2018-Jul-27 20:21:00. 2018-Jul-28 20:21:00. D. 00000000,0. 0003 TDB
Step-size : 0 steps
*******************************************************************************
Center geodetic : 0. 0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0. 00000000,0. 00000000,0. 00000000,0. 1 x 6378. 0000000 {E-lon(deg),Dxy(km),Dz(km)}
Center radii : 6378. 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. 1 x 6356. 347916670 = A. 0
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch
*******************************************************************************
JDTDB
X Y Z
VX VY VZ
LT RG RR
*******************************************************************************
$$SOE
2458327. 2018-Jul-27 20:21:00. D. 537109094089627E-03 Y =-2. 0003 TDB
X = 1. 112037386426180E-06
VX= 4. 237488447258137E-03 Z = 5. 187527302531735E-04 VZ=-5. 593816208618667E-04 VY= 3. 567825598846416E-05 RG= 2. 183707711777675E-05
LT= 1. 707898607099066E-06
$$EOE
*******************************************************************************
Coordinate system description: 714605874095336E-03 RR=-2.

Ecliptic and Mean Equinox of Reference Epoch

0
XY-plane: plane of the Earth's orbit at the reference epoch
Note: obliquity of 84381. Reference epoch: J2000. 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.

700 km, 1 day= 86400. Symbol meaning [1 au= 149597870. 0 s]:

center (au/day) 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.

Geometric states/elements have no aberrations applied.

D. 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. Giorgini@jpl.nasa.gov
*******************************************************************************

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

347916670 = A. $$SOE
2458327. 2018-Jul-27 20:21:00. D. 537109094089627E-03 Y =-2. 0003 TDB
X = 1. 112037386426180E-06
VX= 4. 237488447258137E-03 Z = 5. 187527302531735E-04 VZ=-5. 593816208618667E-04 VY= 3. 567825598846416E-05 RG= 2. 183707711777675E-05
LT= 1. 707898607099066E-06
$$EOE
714605874095336E-03 RR=-2.

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

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

Полный вывод эфемерид Земли на 27.07.2018 20:21 (начало координат в центре масс Солнечной системы)

*******************************************************************************
Revised: Jul 31, 2013 Earth 399

Mean Radius (km) = 6371. GEOPHYSICAL PROPERTIES (revised Aug 13, 2018):
Vol. 02 Mass x10^24 (kg)= 5. 01+-0. 0006
Equ. 97219+-0. 137 Mass layers:
Polar axis, km = 6356. radius, km = 6378. 1 x 10^18 kg
Flattening = 1/298. 752 Atmos = 5. 4 x 10^21 kg
Density, g/cm^3 = 5. 257223563 oceans = 1. 6 x 10^22 kg
J2 (IERS 2010) = 0. 51 crust = 2. 043 x 10^24 kg
g_p, m/s^2 (polar) = 9. 00108262545 mantle = 4. 835 x 10^24 kg
g_e, m/s^2 (equatorial) = 9. 8321863685 outer core = 1. 675 x 10^22 kg
g_o, m/s^2 = 9. 7803267715 inner core = 9. 435436 Inner core rad = 1215 km
GM 1-sigma, km^3/s^2 = 0. 82022 Fluid core rad = 3480 km
GM, km^3/s^2 = 398600. 186 km/s
Rot. 0014 Escape velocity = 11. 00007292115 Surface Area:
Mean sidereal day, hr = 23. Rate (rad/s) = 0. 48 x 10^8 km
Mean solar day 2000. 9344695944 land = 1. 002 sea = 3. 0, s = 86400. 0, s = 86400. 62 x 10^8 km
Mean solar day 1820. 3308 Love no., k2 = 0. 0
Moment of inertia = 0. pressure = 1. 299
Mean Temperature, K = 270 Atm. mag. 0 bar
Vis. 86 Volume, km^3 = 1. V(1,0) = -3. 367 Magnetic moment = 0. 08321 x 10^12
Geometric Albedo = 0. 6 (mean), 1414 (perihelion), 1322 (aphelion) 61 gauss Rp^3
Solar Constant (W/m^2) = 1367.

4392911 Sidereal orb period = 1. ORBIT CHARACTERISTICS:
Obliquity to orbit, deg = 23. 79 Sidereal orb period = 365. 0000174 y
Orbital speed, km/s = 29. 9856474 Hill's sphere radius = 234. 25636 d
Mean daily motion, deg/d = 0. 9
*******************************************************************************

D. *******************************************************************************
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. 0003 TDB
Stop time : A. 2018-Jul-27 20:21:00. 2018-Jul-28 20:21:00. D. 00000000,0. 0003 TDB
Step-size : 0 steps
*******************************************************************************
Center geodetic : 0. 0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0. 00000000,0. 00000000,0. 00000000,0. 0
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch
*******************************************************************************
JDTDB
X Y Z
VX VY VZ
LT RG RR
*******************************************************************************
$$SOE
2458327. 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. D. 347916670 = A. 0003 TDB
X = 5. 2018-Jul-27 20:21:00. 298818915224488E-01 Z =-5. 755663665315949E-01 Y =-8. 388633512282171E-02 VY= 9. 366994499016168E-05
VX= 1. 429889230737491E-07
LT= 5. 678934168415631E-03 VZ= 3. 009940888883960E+00 RR=-3. 832932117417083E-03 RG= 1. 947237246302148E-05
$$EOE
*******************************************************************************
Coordinate system description:

Ecliptic and Mean Equinox of Reference Epoch

0
XY-plane: plane of the Earth's orbit at the reference epoch
Note: obliquity of 84381. Reference epoch: J2000. 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.

700 km, 1 day= 86400. Symbol meaning [1 au= 149597870. 0 s]:

center (au/day) 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.

Geometric states/elements have no aberrations applied.

D. 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. Giorgini@jpl.nasa.gov
*******************************************************************************

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

347916670 = A. $$SOE
2458327. 2018-Jul-27 20:21:00. D. 755663665315949E-01 Y =-8. 0003 TDB
X = 5. 366994499016168E-05
VX= 1. 298818915224488E-01 Z =-5. 678934168415631E-03 VZ= 3. 388633512282171E-02 VY= 9. 832932117417083E-03 RG= 1. 429889230737491E-07
LT= 5. 947237246302148E-05
$$EOE
009940888883960E+00 RR=-3.

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

Полный вывод эфемерид Луны на 27.07.2018 20:21 (начало координат в центре масс Солнечной системы)

*******************************************************************************
Revised: Jul 31, 2013 Moon / (Earth) 301

Mean Radius, km = 1737. GEOPHYSICAL DATA (updated 2018-Aug-13):
Vol. 03 Mass, x10^22 kg = 7. 53+-0. 0 Surface emissivity = 0. 349
Radius (gravity), km = 1738. 4 GM, km^3/s^2 = 4902. 92
Radius (IAU), km = 1737. 3437 GM 1-sigma, km^3/s^2 = +-0. 800066
Density, g/cm^3 = 3. 21 Surface accel., m/s^2 = 1. 0001
V(1,0) = +0. 3005690769 Farside crust. 62
Earth/Moon mass ratio = 81. = ~80 - 90 km
Mean crustal density = 2. thick. 07 g/cm^3 Nearside crust. 97+-. 1+-. thick.= 58+-8 km
Heat flow, Apollo 15 = 3. 024059
Heat flow, Apollo 17 = 2. 6 mW/m^2 k2 = 0. 5 mW/m^2 Rot. 2+-. 0000026617
Geometric Albedo = 0. Rate, rad/s = 0. 12

2" Orbit period = 27. Mean angular diameter = 31'05. 67 deg Eccentricity = 0. 321582 d
Obliquity to orbit = 6. 145 deg
Mean motion, rad/s = 2. 05490
Semi-major axis, a = 384400 km Inclination = 5. 38 d
Apsidal period = 3231. 6616995x10^-6 Nodal period = 6798. of inertia C/MR^2= 0. 50 d Mom. 310213 gamma (B-A/C), x10^-4 = 2. 393142
beta (C-A/B), x10^-4 = 6. 277317

2 5. 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
******************************************************************************* 2 5.

D. *******************************************************************************
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. 0003 TDB
Stop time : A. 2018-Jul-27 20:21:00. 2018-Jul-28 20:21:00. D. 00000000,0. 0003 TDB
Step-size : 0 steps
*******************************************************************************
Center geodetic : 0. 0000000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 0. 00000000,0. 00000000,0. 00000000,0. 0
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch
*******************************************************************************
JDTDB
X Y Z
VX VY VZ
LT RG RR
*******************************************************************************
$$SOE
2458327. 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. D. 347916670 = A. 0003 TDB
X = 5. 2018-Jul-27 20:21:00. 321193799697072E-01 Z =-4. 771034756256845E-01 Y =-8. 434571674368357E-02 VY= 9. 855790760378579E-05
VX= 1. 149408819470315E-05
LT= 5. 997686898668805E-03 VZ=-5. 012655462859054E+00 RR=-3. 848610189172283E-03 RG= 1. 979984423450087E-05
$$EOE
*******************************************************************************
Coordinate system description:

Ecliptic and Mean Equinox of Reference Epoch

0
XY-plane: plane of the Earth's orbit at the reference epoch
Note: obliquity of 84381. Reference epoch: J2000. 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.

700 km, 1 day= 86400. Symbol meaning [1 au= 149597870. 0 s]:

center (au/day) 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.

Geometric states/elements have no aberrations applied.

D. 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. Giorgini@jpl.nasa.gov
*******************************************************************************

347916670 = A. $$SOE
2458327. 2018-Jul-27 20:21:00. D. 771034756256845E-01 Y =-8. 0003 TDB
X = 5. 855790760378579E-05
VX= 1. 321193799697072E-01 Z =-4. 997686898668805E-03 VZ=-5. 434571674368357E-02 VY= 9. 848610189172283E-03 RG= 1. 149408819470315E-05
LT= 5. 979984423450087E-05
$$EOE
012655462859054E+00 RR=-3.

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

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

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

347916670 = A. $$SOE
2458327. 2018-Jul-27 20:21:00. D. 520050993518213E+04 Y = 1. 0003 TDB
X = 6. 304404963058507E+04
VX=-1. 049687363172734E+06 Z =-1. 853475278436883E-03 VZ= 3. 265326939350981E-02 VY= 5. 508397935601254E+00 RG= 1. 136673455633667E-04
LT= 3. 053500842402456E-03
$$EOE

Взглянув на параметр RG мы увидим, что Солнце вращается вокруг барицентра Солнечной системы, и на 27. 051791240756026E+06 RR= 5. 2018 центр звезды находится от него на расстоянии в миллион километров. 07. То есть барицентр Солнечной системы лежит в полумиллионе километров от поверхности светила. Радиус Солнца, для справки — 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))

0
mu[1] = 386326400000000. Выхлоп программы
Гравитационные параметры тел
mu[0] = 4901783000000. 326663e+20
Нормированные гравитационные параметры
xi[0] = 3. 0
mu[2] = 1. 912016088486677e-06
xi[2] = 1. 6948215183509304e-08
xi[1] = 2. 0

35432583 Масштаб времени T = 31563683.

77103476e-01 -8. Начальное положение Луны, а.е.: [ 5. 85579076e-05]
Начальное положение Земли, а.е.: [ 5. 32119380e-01 -4. 29881892e-01 -5. 75566367e-01 -8. 69738146e-06 2. 36699450e-05]
Начальное положение Солнца, а.е.: [-1. 58081871e-10] 44737475e-06 1.

98933473 17310. Начальная скорость Луны, м/с: [24838. 15979106]
-//- безразмерная: [ 5. 56333294 -89. 65235907 -0. 24078311 3. 40435899e+04 1. 01881184]
Начальная скорость Земли, м/с: [2. 93870516e-01]
-//- безразмерная: [5. 67586567e+04 5. 53591219e+00 1. 07296163e+00 3. 09330769e-02 -4. 25300854e-04]
Начальная скорость Солнца, м/с: [-7. 56493465e-06]
-//- безразмерная: [-1. 94410725e-02 1. 04315813e-05 3. 49661835e-05 -1. 30185861e-10]

Для преобразования системы к форме Коши вспоминаем, что
Собственно само интегрирование сводится к более-менее стандартной для 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. Ценить и уважать механику я призываю и вас, мои уважаемые читатели.

Спасибо за внимание!

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть