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

в 10:24, , рубрики: diy или сделай сам, задача двух тел, космонавтика, математика, Научно-популярное, первый закон кеплера, ракета, ракетостроение, физика, эллипс

Содержание

Часть 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. $ $ h, $ как мы помним, может равняться нулю в случае прямолинейного движения вдоль радиус-вектора. Этот случай слишком прост, его даже рассматривать не будем, а кто хочет может приравнять в уравнении к нулю и дальше его решить.

Здесь попытаемся решать для $ hneq0. $ Первым что приходит в голову (в мою пришло, когда я впервые увидел это уравнение) то, что здесь нет $ 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}{drho}dfrac{drho}{dt}=z^{'}dot{rho}=z^{'}z $

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

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

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

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

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

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

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

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

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

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

$ dot{rho}^{2}=dfrac{2murho - h^{2} + Crho^{2}}{rho^{2}} $

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

$ dot{rho}=pmsqrt{dfrac{2murho - h^{2} + Crho^{2}}{rho^{2}}}=pmdfrac{sqrt{Crho^{2} + 2murho - h^{2}}}{rho} $

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

$ intdfrac{rho drho}{sqrt{Crho^{2} + 2murho - 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=hintdfrac{dt}{rho^{2}}=hintdfrac{dt}{left(f^{-1}(t)right)^{2}} $

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

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

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

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

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

$ rho^{2}>0 $ — эта штука всегда больше нуля. Правильней с точки зрения математики сказать: больше либо равно нулю, но с физической точки зрения — только больше. Ведь тела у нас всё таки — не материальные точки. Значит и сблизиться на нулевое расстояние никогда не смогут. То бишь центры масс их никогда не совпадут, иначе им пришлось бы пройти друг сквозь друга. Так что кто боится деления на ноль — не боитесь.

Так о чём это я, ах да. $ h $ тоже всегда больше либо меньше нуля, либо ноль. Ведь это постоянная величина (кстати, когда ноль — тогда угол постоянен и движение вдоль радиуса вектора, еще раз убедились). А это значит, что производная угла постоянна по знаку на протяжении всего движения:

$ dot{phi} > 0 $

или

$ dot{phi} < 0 $

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

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

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

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

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

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

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

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

$ rho=dfrac{l}{1 + ecos(phi)}, $

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

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

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

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

image
Эллипс

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

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

image
Парабола

Гипербола ($ e > 1 $):

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

image
Гипербола

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

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

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

$ dfrac{dphi}{dt}dfrac{dt}{drho}=pmdfrac{h}{rho^{2}}dfrac{rho}{sqrt{Crho^{2} + 2murho - h^{2}}}, $

$ dfrac{dphi}{drho}=dfrac{h}{rhosqrt{Crho^{2} + 2murho - h^{2}}} $

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

$ int dphi=hintdfrac{drho}{rhosqrt{Crho^{2} + 2murho - h^{2}}} $

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

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

$ rho=dfrac{l}{1 + ecos(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 + bcos(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( -2dfrac{1}{u^{3}}dot{u}^{2} + dfrac{1}{u^{2}}ddot{u} right)=- dfrac{ddot{u}}{u^{2}} + 2dfrac{dot{u}^{2}}{u^{3}} $

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

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

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

$ u=u(phi), $

$ dot{u}=dfrac{du}{dphi}dfrac{dphi}{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^{'})=hleft( dfrac{du^{2}}{dt}u^{'} + u^{2}dfrac{du^{'}}{dt} right)=$

$=h left( 2udfrac{du}{dt}u^{'} + u^{2}dfrac{du^{'}}{dphi}dfrac{dphi}{dt} right)=h left( 2udot{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}} + 2dfrac{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}} + 2dfrac{(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 + bcos(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)=Acos(phi - omega), $

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

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

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

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

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

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

$ u^{'}=-Asin(phi - omega) $

И система для нахождения произвольных констант:
begin{equation*}
begin{cases}
dfrac{mu}{h^{2}} + Acos(0 — omega) = u_{0},
\
-Asin(0 — omega) = u^{'}_{0}.
end{cases}
end{equation*}
begin{equation*}
begin{cases}
Acos(omega) = u_{0} — dfrac{mu}{h^{2}} = a,
\
Asin(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{Asin(omega)}{Acos(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}{dphi}=dfrac{d}{dphi}left( dfrac{1}{rho} right)=-dfrac{1}{rho^{2}}dfrac{drho}{dphi}=-dfrac{1}{rho^{2}}dfrac{drho}{dt} / dfrac{dphi}{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 = rhocos(phi),
\
y = rhosin(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) — rhosin(phi)dot{phi}
\
dot{y} = dot{rho}sin(phi) + rhocos(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): 0x9234016395e0e6ef7cf6c0aa0f6f48f91ab39239
Ripple (XRP): rLW9gnQo7BQhU6igk5keqYnH3TVrCxGRzm (адрес), 270547561 (тег)
Bitcoin Cash (BCH): bitcoincash:qzxfz2hdcl0hv23a3hlcefsy07mglssjtgwrckhyg8

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

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

Автор: Серый

Источник


* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js