Проектируем космическую ракету с нуля. Часть 5 — Первый закон Кеплера14.09.2019 16:48

Содержание


Часть 1 — Задача двух тел
Часть 2 — Полу-решение задачи двух тел
Часть 3 — Ужепочти-решение задачи двух тел
Часть 4 — Второй закон Кеплера

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

$ \ddot{\rho} = -\dfrac{\mu}{\rho^{2}} + \dfrac{h^{2}}{\rho^{3}} $


Это дифференциальное уравнение второго порядка, где в качестве неизвестной функции — длина радиуса вектора, зависящего от времени. Здесь $ h^{2}\geq0, \mu = G(m_{1} + m_{2})>0. $» /><img src= как мы помним, может равняться нулю в случае прямолинейного движения вдоль радиус-вектора. Этот случай слишком прост, его даже рассматривать не будем, а кто хочет может приравнять в уравнении к нулю и дальше его решить.
Здесь попытаемся решать для $ h\neq0. $ Первым что приходит в голову (в мою пришло, когда я впервые увидел это уравнение) то, что здесь нет $ \dot{\rho} $, ну, конечно, и $ t. $ А в таких случаях (специальных) можно проводить дальнейшую замену, которая понижает порядок уравнения до первого.

В общем случае уравнение второго порядка может быть записано в виде:

$ F(t,\rho,\dot{\rho},\ddot{\rho}) = 0 $


У нас же — уравнение попроще, когда:

$ F(\rho,\ddot{\rho}) = 0 $


А в таких случаях можно сделать замену, которая снизит порядок уравнения:

$ \dot{\rho} = z(\rho), $


где $ z $ — новая неизвестная функция, но которая зависит не напрямую от времени, а от $ \rho. $ Тогда:

$ \ddot{\rho} = \dfrac{dz}{d\rho}\dfrac{d\rho}{dt} = z^{'}\dot{\rho} = z^{'}z $


Здесь продифференцировали как сложную функцию, а потом штрихом обозначили производную $ z $ по $ \rho. $ Теперь всё готово, и можно подставить:

$ z^{'}z = -\dfrac{\mu}{\rho^{2}} + \dfrac{h^{2}}{\rho^{3}} $


Уравнение первого порядка, но относительно $ \rho, $ вместо времени. Причём с разделяющимися переменными:

$ \dfrac{dz}{d\rho}z = -\dfrac{\mu}{\rho^{2}} + \dfrac{h^{2}}{\rho^{3}} $

$ zdz = \left(-\dfrac{\mu}{\rho^{2}} + \dfrac{h^{2}}{\rho^{3}}\right)d\rho $


Проинтегрировать не составляет труда:

$ \dfrac{z^{2}}{2} = \dfrac{\mu}{\rho} - \dfrac{h^{2}}{2\rho^{2}} + \dfrac{C}{2} $


Ну добавили пол константы, какая разница? Зато потом жить проще:

$ z^{2} = \dfrac{2\mu}{\rho} - \dfrac{h^{2}}{\rho^{2}} + C = \dfrac{2\mu\rho - h^{2} + C\rho^{2}}{\rho^{2}} $


Пришло время вспомнить что такое $ z $:

$ \dot{\rho}^{2} = \dfrac{2\mu\rho - h^{2} + C\rho^{2}}{\rho^{2}} $


И извлечь корень:

$ \dot{\rho} = \pm\sqrt{\dfrac{2\mu\rho - h^{2} + C\rho^{2}}{\rho^{2}}} = \pm\dfrac{\sqrt{C\rho^{2} + 2\mu\rho - h^{2}}}{\rho} $


Уравнение опять с разделяющимися переменными, и даже интеграл вроде как берётся в элементарных функциях:

$ \int\dfrac{\rho d\rho}{\sqrt{C\rho^{2} + 2\mu\rho - h^{2}}} = \pm \int dt $


И всё бы хорошо, но проблема в том, что если мы и решим, то получим обратную зависимость, то есть времени от радиуса:

$ t = f(\rho) $


А хотелось бы наоборот:

$ \rho = f^{-1}(t) $


Да еще и этот $ \pm $ — думать какую ветвь выбирать. Но это не самое страшное, нужно будет рассматривать разные случаи соотношений постоянных величин под корнем: $ C>0, C<0, C = 0... $
Можно, конечно, и в wolframalpha вбить и прикинуть, что будет:
image

Это просто страх и ужас, а найти обратную элементарную функцию можно и не мечтать. А ведь нам еще нужно и угол искать:

$ \dot{\phi} = \dfrac{h}{\rho^{2}} $

$ \phi = h\int\dfrac{dt}{\rho^{2}} = h\int\dfrac{dt}{\left(f^{-1}(t)\right)^{2}} $


Слишком муторное дело. И даже умение считать интегралы от обратных функций нас не спасет скорее всего.

Умение считать интегралы от обратных функций

$ \int\dfrac{dt}{\left(f^{-1}(t)\right)^{2}} = \begin{vmatrix} y = f^{-1}(t) \\ f(y) = t \\ dt = f^{'}(y)dy \end{vmatrix} = \int\dfrac{f^{'}(y)dy}{y^{2}} $


Кстати, можно заметить некоторое свойство для угла $ \phi $ вот из этого равенства:

$ \dot{\phi} = \dfrac{h}{\rho^{2}} $


$ \rho^{2}>0 $» /> — эта штука всегда больше нуля. Правильней с точки зрения математики сказать: больше либо равно нулю, но с физической точки зрения — только больше. Ведь тела у нас всё таки — не материальные точки. Значит и сблизиться на нулевое расстояние никогда не смогут. То бишь центры масс их никогда не совпадут, иначе им пришлось бы пройти друг сквозь друга. Так что кто боится деления на ноль — не боитесь. </p> <p>Так о чём это я, ах да. <img src= тоже всегда больше либо меньше нуля, либо ноль. Ведь это постоянная величина (кстати, когда ноль — тогда угол постоянен и движение вдоль радиуса вектора, еще раз убедились). А это значит, что производная угла постоянна по знаку на протяжении всего движения:

$ \dot{\phi} > 0 $» /></p> <p><br />или </p> <p><img src=


в зависимости от знака $ h $.

Иными словами сам угол $ \phi $ — строго монотонная функция. Кто не помнит как это, вот наглядная картинка:
image
Монотонная и немонотонная функция

Этот факт подтверждает то, что мы видели при численном моделировании. Тело всегда движется в одном направлении. Никогда не останавливается (производная равна нулю) и не начинает крутиться в обратном направлении (производная изменила знак).

Соответственно траектории могут и не могут быть такими:
image
Синие — возможные траектории, красные — невозможные (в полярных координатах)

Ну и анимации лишними не бывают:
image
Производная $ \dot{\phi} $ меняет знак. Таких решений у нас точно не будет.

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

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

В полярных координатах коники можно задать так (с центром в одном из фокусов и нулевым направлением вдоль главной оси):

$ \rho = \dfrac{l}{1 + e\cos(\phi)}, $


где $ e $ обозначает эксцентриситет, а $ l $ фокальный параметр.

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

Эллипс ($ e < 1 $):

$ \rho = \dfrac{1}{1 + 0.5\cos(\phi)} $


image
Эллипс

Парабола ($ e = 1 $):

$ \rho = \dfrac{1}{1 + \cos(\phi)} $


image
Парабола

Гипербола ($ e > 1 $» />): <br /></p> <p><img src=


image
Гипербола

Очень всё похоже на траектории движения спутника, полученные численным моделированием…

Как видим, здесь $ \rho = \rho(\phi) $ — радиус как функция от угла. И выражения довольно таки простые. Так может имеет смысл найти сначала зависимость радиуса от угла, а потом угла как функции времени? Можно попробовать и проверить нашу гипотезу, что решения являются кониками.

Вспомним что у нас есть:
\begin{equation*}
\begin{cases}
\dot{\phi} = \dfrac{h}{\rho^{2}},
\\
\dot{\rho} = \pm\dfrac{\sqrt{C\rho^{2} + 2\mu\rho — h^{2}}}{\rho}.
\end{cases}
\end{equation*}
Поделив одно на другое, можно избавиться от $ dt $:

$ \dfrac{d\phi}{dt}\dfrac{dt}{d\rho} = \pm\dfrac{h}{\rho^{2}}\dfrac{\rho}{\sqrt{C\rho^{2} + 2\mu\rho - h^{2}}}, $

$ \dfrac{d\phi}{d\rho} = \dfrac{h}{\rho\sqrt{C\rho^{2} + 2\mu\rho - h^{2}}} $


Но опять неудобно получается, наоборот:

$ \int d\phi = h\int\dfrac{d\rho}{\rho\sqrt{C\rho^{2} + 2\mu\rho - h^{2}}} $


Этот интеграл конечно можно взять и даже, скорее всего, потом найти обратную функцию труда не составит. Но есть способ попроще найти зависимость $ \rho = \rho(\phi). $

Мы предполагаем, что решение будет иметь такой вид:

$ \rho = \dfrac{l}{1 + e\cos(\phi)}, $


Вот этого уравнения:

$ \ddot{\rho} = -\dfrac{\mu}{\rho^{2}} + \dfrac{h^{2}}{\rho^{3}} $


Но производная по времени, всё же. А нам бы по углу. Ничего не мешает просто сделать сложную функцию:

$ \rho = \rho(\phi(t)) $


и дальше продифференцировать пару раз, при этом используя уже полученные нами факты ($ \rho^{2}\dot{\phi} = h $).

Но перед этим само уравнение и предполагаемое решение само просится сделать такую замену:

$ \rho = \dfrac{1}{u} $


Тогда сама функция $ u $ будет такой:

$ u = a + b\cos(\phi) $


А это решения вот таких элементарных диффуров:

$ \ddot{u} + k u = f $


То есть сделав последовательную замену:

$ \rho = \dfrac{1}{u}, $

$ u = u(\phi) $


Можно надеяться, что уравнение будет сведено к очень простенькому.

Приступим:

$ \rho = \dfrac{1}{u}, $

$ \dot{\rho} = \dfrac{d}{dt}u^{-1} = -\dfrac{1}{u^{2}}\dot{u} = -\dfrac{\dot{u}}{u^{2}}, $

$ \ddot{\rho} = -\dfrac{d}{dt}\left(\dfrac{1}{u^{2}}\dot{u}\right) = -\left( \dfrac{du^{-2}}{dt}\dot{u} + \dfrac{1}{u^{2}}\ddot{u} \right) = $

$ = -\left( -2\dfrac{1}{u^{3}}\dot{u}^{2} + \dfrac{1}{u^{2}}\ddot{u} \right) = - \dfrac{\ddot{u}}{u^{2}} + 2\dfrac{\dot{u}^{2}}{u^{3}} $


Подставляя в исходное будет так:

$ - \dfrac{\ddot{u}}{u^{2}} + 2\dfrac{\dot{u}^{2}}{u^{3}} = -\mu u^{2} + h^{2} u^{3} $


И продолжим, пусть теперь:

$ u = u(\phi), $

$ \dot{u} = \dfrac{du}{d\phi}\dfrac{d\phi}{dt} = u^{'}\dot{\phi} = \begin{vmatrix} \dot{\phi} = \dfrac{h}{\rho^{2}} \\ \dot{\phi} = hu^{2} \end{vmatrix} = hu^{2}u^{'}, $

$ \ddot{u} = \dfrac{d}{dt}(hu^{2}u^{'}) = h\left( \dfrac{du^{2}}{dt}u^{'} + u^{2}\dfrac{du^{'}}{dt} \right) = $

$ = h \left( 2u\dfrac{du}{dt}u^{'} + u^{2}\dfrac{du^{'}}{d\phi}\dfrac{d\phi}{dt} \right) = h \left( 2u\dot{u}u^{'} + u^{2}u^{''}\dot{\phi} \right) = $

$ = h \left( 2uhu^{2}u^{'}u^{'} + u^{2}u^{''}hu^{2} \right) = h^{2} \left( 2u^{3}(u^{'})^{2} + u^{4}u^{''} \right) $


Да, немножко поднапрячься нужно, чтобы всё продифференцировать и не ошибиться, но всё-таки дифференцировать это не интегрировать. Да и вообще, всё можно машине поручить, она правильно возьмет производную. Машины, конечно, могут и интегрировать, но не всегда. А вот дифференцировать всё что угодно. Интегрирование всё-таки — это творчество и всегда ним останется IMHO (Ломать — не строить).

В итоге, что у нас есть:

$ - \dfrac{\ddot{u}}{u^{2}} + 2\dfrac{\dot{u}^{2}}{u^{3}} = -\mu u^{2} + h^{2} u^{3}, $

$ - \dfrac{h^{2} \left( 2u^{3}(u^{'})^{2} + u^{4}u^{''} \right)}{u^{2}} + 2\dfrac{(hu^{2}u^{'})^{2}}{u^{3}} = -\mu u^{2} + h^{2} u^{3}, $

$ -2h^{2}u(u^{'})^{2} - h^{2}u^{2}u^{''} + 2h^{2}u(u^{'})^{2} = -\mu u^{2} + h^{2} u^{3}, $

$ - h^{2}u^{2}u^{''} = -\mu u^{2} + h^{2} u^{3}, $

$ - h^{2}u^{''} = -\mu + h^{2} u, $

$ u^{''} = \dfrac{\mu}{h^{2}} - u, $

$ u^{''} + u = \dfrac{\mu}{h^{2}}. $


Ну что, будем решать это уравнение?

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

Стандартный метод решения, но для красивости, всё запишем, ничего не пропустим:

$ u^{''} + u = 0, $

$ \lambda^{2} + 1 = 0, $

$ \lambda = \pm i, $


А значит решение будет иметь вид:

$ u_{o} = C_{1}sin(\phi) + C_{2}cos(\phi), $


Общее решение будет как сумма однородного и частного. Частное можно искать в виде константы:

$ u_{\text{ч}} = K, $


Тогда общее решение:

$ u = C_{1}sin(\phi) + C_{2}cos(\phi) + K $


Дифференцируя два раза и подставляя в исходное уравнение, находим константу частного решения:

$ u^{''} = -C_{1}sin(\phi) - C_{2}cos(\phi) $

$ u^{''} + u = -C_{1}sin(\phi) - C_{2}cos(\phi) + C_{1}sin(\phi) + C_{2}cos(\phi) + K = \dfrac{\mu}{h^{2}} $

$ K = \dfrac{\mu}{h^{2}}$


Ну вот, так совпало, что константа это и есть правая часть нашего уравнения. В итоге:

$ u = C_{1}sin(\phi) + C_{2}cos(\phi) + \dfrac{\mu}{h^{2}} $


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

$ u = a + b\cos(\phi) $


Так же? Ну ничего страшного, кто хорошо знаком с тригонометрией, тот легко всё запишет как пожелает (известный трюк):

$ C_{1}sin(\phi) + C_{2}cos(\phi) = \sqrt{C_{1}^{2} + C_{2}^{2}} \left( \dfrac{C_{1}}{\sqrt{C_{1}^{2} + C_{2}^{2}}}sin(\phi) + \dfrac{C_{2}}{\sqrt{C_{1}^{2} + C_{2}^{2}}}cos(\phi) \right) = $

$ = \sqrt{C_{1}^{2} + C_{2}^{2}} \left( sin(\omega)sin(\phi) + cos(\omega)cos(\phi) \right) = A\cos(\phi - \omega), $


где $ A, \omega $ — две произвольные константы, выраженные через старые $ C_{1}, C_{2} $ несложным образом. И нам вообще это неважно, ибо это константы.

Таким образом имеем решение, в более-менее ожидаемом виде:

$ u = \dfrac{\mu}{h^{2}} + A\cos(\phi - \omega) $


А неизвестные постоянные находятся тоже ожидаемо, для этого нам нужно значение функции в начальной точке, а также производная:
\begin{equation*}
\begin{cases}
u (0) = u_{0},
\\
u^{'}(0) = u^{'}_{0}.
\end{cases}
\end{equation*}
Не стоит также забывать, что зависимость здесь от угла, а не от времени:

$ u(0) = u(\phi = 0) $


Нужна производная, пожалуйста:

$ u^{'} = -A\sin(\phi - \omega) $


И система для нахождения произвольных констант:
\begin{equation*}
\begin{cases}
\dfrac{\mu}{h^{2}} + A\cos (0 — \omega) = u_{0},
\\
-A\sin (0 — \omega) = u^{'}_{0}.
\end{cases}
\end{equation*}
\begin{equation*}
\begin{cases}
A\cos (\omega) = u_{0} — \dfrac{\mu}{h^{2}} = a,
\\
A\sin (\omega) = u^{'}_{0} = b.
\end{cases}
\end{equation*}
Для удобства обозначили за a и b правые части. Теперь легко найти $ A $ возведя в квадраты равенства и сложив их:

$ A^{2}\cos^{2}(\omega) + A^{2}\sin^{2}(\omega) = a^{2}+b^{2} $

$ A = \pm \sqrt{a^{2} + b^{2}} $


Вероятно нужно будет выбрать +, ведь $ A $ — это эксцентриситет. Ну да Бог с ним, потом если что и пере-обозначить можно будет в $ e $.
Начальный угол $ \omega $ можно найти поделив одно уравнение на другое:

$ \dfrac{A\sin(\omega)}{A\cos(\omega)} = \dfrac{b}{a}, $

$ \tan(\omega) = \dfrac{b}{a}, $

$ \omega = \arctan(\dfrac{b}{a}), $


Вроде всё нашли, но теперь бы распутать клубок с $ a, b $:
\begin{equation*}
\begin{cases}
a = u_{0} — \dfrac{\mu}{h^{2}},
\\
b = u^{'}_{0}.
\end{cases}
\end{equation*}
Разберем и вспомним что здесь кто:

$ \mu = G(m_{1} + m_{2})$


Чтобы остальное найти, нужно определиться с начальными условиями в полярной системе, они будут такие (здесь уже в начальный момент времени $ t = 0 $):
\begin{equation*}
\begin{cases}
\rho (0) = \rho_{0},
\\
\dot{\rho}(0) = \dot{\rho}_{0},
\\
\phi (0) = \phi_{0},
\\
\dot{\phi}(0) = \dot{\phi}_{0}. \\
\end{cases}
\end{equation*}
Тогда в любой момент времени:

$ h = \rho^{2}\dot{\phi} $


В частности нулевой момент позволит определить постоянную $ h $:

$ h = \rho^{2}_{0}\dot{\phi}_{0} $


Момент импульса в нулевой момент времени (игра слов).

И наконец:

$ u_{0} = \dfrac{1}{\rho_{0}} $


А с производной нужно немноооожко повозиться:

$ u^{'} = \dfrac{du}{d\phi} = \dfrac{d}{d\phi}\left( \dfrac{1}{\rho} \right) = -\dfrac{1}{\rho^{2}}\dfrac{d\rho}{d\phi} = -\dfrac{1}{\rho^{2}}\dfrac{d\rho}{dt} / \dfrac{d\phi}{dt} = -\dfrac{1}{\rho^{2}} \dfrac{\dot{\rho}}{\dot{\phi}} $


То что нам нужно:

$ u^{'}_{0} = -\dfrac{1}{\rho^{2}_{0}} \dfrac{\dot{\rho}_{0}}{\dot{\phi}_{0}} $


Вроде бы разобрались. Но мы то начинали не с полярной системы координат. Не с неё… Еще возни немного будет.

Как мы в полярную попали. Так:
\begin{equation*}
\begin{cases}
x = \rho\cos (\phi),
\\
y = \rho\sin (\phi).
\end{cases}
\end{equation*}
Обратно как? Так (сразу в нулевой момент):
\begin{equation*}
\begin{cases}
\rho_{0} = \sqrt{x^{2}_{0} + y^{2}_{0}},
\\
\phi_{0} = \arctan (\dfrac{y_{0}}{x_{0}}).
\end{cases}
\end{equation*}
Но если кто помнит, мы в прошлые разы направили ось $ x $ вдоль $ \vec{r}_{0} $, а потому должно быть так:

$ x_{0} = |\vec{r}_{0}|, y_{0} = 0 $


Соответственно:

$ \rho_{0} = x_{0}, \phi_{0} = 0 $


Ну и не зря же мы до этого начальный угол брали ноль, в смысле здесь $ u(0) = u(\phi = 0). $

Со скоростями тоже всё просто (копируем формулы с прошлой статьи):
\begin{equation*}
\begin{cases}
\dot{x} = \dot{\rho}\cos (\phi) — \rho\sin (\phi)\dot{\phi}
\\
\dot{y} = \dot{\rho}\sin (\phi) + \rho\cos (\phi)\dot{\phi}
\end{cases}
\end{equation*}
Правда нужно наоборот, тут линейная система уравнений, решаем:
\begin{equation*}
\begin{cases}
\dot{\rho} = \dot{x}\cos (\phi) + \dot{y}\sin (\phi)
\\
\dot{\phi} = \dfrac{\dot{y}\cos (\phi) — \dot{x}\sin (\phi)}{\rho}
\end{cases}
\end{equation*}
Этим можно пользоваться даже когда начальный угол не ноль, но в нашем случае:
\begin{equation*}
\begin{cases}
\dot{\rho}_{0} = \dot{x}_{0}
\\
\dot{\phi}_{0} = \dfrac{\dot{y}_{0}}{\rho_{0}}
\end{cases}
\end{equation*}
Ну, а дальше я не стану расписывать как из начальных данных $ \vec{r}_{10}, \vec{r}_{20}, \dot{\vec{r}}_{10}, \dot{\vec{r}}_{20} $ (мы ведь с этого начинали) вычислить $ x_{0}, y_{0}, \dot{x}_{0}, \dot{y}_{0} $. Скажу лишь что следующий шаг будет применение обратной матрицы преобразования, она где то в прошлых статьях затерялась. И таким образом мы обратно выйдем из плоскости в трехмерное пространство. Как ни уютно было в двумерном, но мы живем в трехмерном…

На сегодня всё. Продолжение следует…

хлеб наш насущный дай нам на сей день.
PayPal ($): what.is.truth.19@gmail.com
Bitcoin (BTC): 1AodAFYCbwrwTiZb5JVsQjv37G5toBcyQ
Ethereum Classic (ETC): 0×9234016395e0e6ef7cf6c0aa0f6f48f91ab39239
Ripple (XRP): rLW9gnQo7BQhU6igk5keqYnH3TVrCxGRzm (адрес), 270547561 (тег)
Bitcoin Cash (BCH): bitcoincash: qzxfz2hdcl0hv23a3hlcefsy07mglssjtgwrckhyg8

или webmoney (Ниже: Поддержать автора → Отправить деньги)

Нет столь святаго, как Господь; ибо нет другого, кроме Тебя; и нет твердыни, как Бог наш.
Не умножайте речей надменных; дерзкие слова да не исходят из уст ваших; ибо Господь есть Бог ведения, и дела у Него взвешены.
Лук сильных преломляется, а немощные препоясываются силою;
сытые работают из хлеба, а голодные отдыхают; даже бесплодная рождает семь раз, а многочадная изнемогает.
Господь умерщвляет и оживляет, низводит в преисподнюю и возводит;
Господь делает нищим и обогащает, унижает и возвышает.
Из праха подъемлет Он бедного, из брения возвышает нищего, посаждая с вельможами, и престол славы дает им в наследие; ибо у Господа основания земли, и Он утвердил на них вселенную.
Стопы святых Своих Он блюдет, а беззаконные во тьме исчезают; ибо не силою крепок человек.
Господь сотрет препирающихся с Ним; с небес возгремит на них. Господь будет судить концы земли, и даст крепость царю Своему и вознесет рог помазанника Своего.
1-я Царств 2

© Habrahabr.ru