Рівняння Кеплера

Анімація, що ілюструє справжню аномалію, ексцентричну аномалію, середню аномалію і розв'язок рівняння Кеплера (у правому верхньому куті), ексцентриситет — 0,6.

Рівня́ння Ке́плера описує рух тіла по еліптичній орбіті в задачі двох тіл і має вигляд:

$E-e\,\sin E=M$

де $E$ — ексцентрична аномалія, $e$ — ексцентриситет орбіти, а $M$ — середня аномалія.

Вперше це рівняння отримав астроном Йоганн Кеплер у 1619 році. Відіграє значну роль у небесній механіці.

Варіанти рівняння Кеплера

Рівняння Кеплера в класичній формі описує рух лише по еліптичних орбітах, тобто при $0\leq e<1$. Рух по гіперболічних орбітах $(e>1)$ підкоряється гіперболічному рівняння Кеплера, схожому за формою з класичним. Рух по прямій лінії $(e=-1)$ описує радіальне рівняння Кеплера. Нарешті, для опису руху по параболічній орбіті $(e=1)$ використовують рівняння Баркера. При $e<0$ орбіт не існує.

Завдання, що приводить до рівняння Кеплера

Розглянемо рух тіла по орбіті в полі іншого тіла. Знайдемо залежність положення тіла на орбіті від часу. З II закону Кеплера випливає, що

$r^{2}{\frac {d\vartheta }{dt}}={\rm {const}}={\sqrt {\mu a\left(1-e^{2}\right)}}$.

Тут $r$ — відстань від тіла до гравітуючого центра, $\vartheta$ — істинна аномалія — кут між напрямками на перицентр орбіти й на тіло, $\mu =GM_{0}$ — добуток гравітаційної сталої на масу гравітуючого тіла, $a$ — велика піввісь орбіти. Звідси можна отримати залежність часу руху по орбіті від істинної аномалії:

$t-t_{0}={\frac {1}{\sqrt {\mu a\left(1-e^{2}\right)}}}\int \limits _{0}^{\vartheta }r^{2}d\vartheta$.

Тут $t_{0}$ — час проходження через перицентр.

Подальше розв'язування задачі залежить від типу орбіти, по якій рухається тіло.

Еліптична орбіта

Рівняння еліпса в полярних координатах має вигляд

$r={\frac {a(1-e^{2})}{1+e\cos {\vartheta }}}$

Тоді рівняння часу набуває вигляду

$t-t_{0}={\frac {\left(a\left(1-e^{2}\right)\right)^{3/2}}{\sqrt {\mu }}}\int \limits _{0}^{\vartheta }{\frac {d\vartheta }{(1+e\cos {\vartheta })^{2}}}$

Для того, щоб взяти інтеграл вводять таку підстановку:

$\operatorname {tg} {\frac {\vartheta }{2}}={\sqrt {\frac {1+e}{1-e}}}\cdot \operatorname {tg} {\frac {E}{2}}$

Величина E називається ексцентричною аномалією. Завдяки такій підстановці інтеграл легко береться. Виходить таке рівняння:

$t-t_{0}={\sqrt {\frac {a^{3}}{\mu }}}\left(E-e\sin E\right)$

Величина ${\sqrt {\frac {\mu }{a^{3}}}}$ є середньою кутовою швидкістю руху тіла по орбіті. В небесній механіці для цієї величини використовується термін середній рух. Добуток середнього руху на час називається середньою аномалією M. Ця величина являє собою кут, на який повернувся б радіус-вектор тіла, якби воно рухалося по коловій орбіті з радіусом, рівним великій півосі орбіти тіла.

Таким чином отримуємо рівняння Кеплера для еліптичного руху:

$E-e\sin E=M$

Гіперболічна орбіта

Рівняння гіперболи в полярних координатах має такий самий вигляд, як і рівняння еліпса. Отже, інтеграл виходить такий самий на вигляд. Однак, використовувати ексцентричну аномалію в цьому випадку не можна. Скористаємося параметричним поданням гіперболи: $x=-a\,\mathrm {ch} \,H$, $y=a{\sqrt {e^{2}-1}}\,\mathrm {sh} \,H$. Тоді рівняння гіперболи набуває вигляду

$r=a\left(e\,\mathrm {ch} \,H-1\right)$,

а зв'язок між $\vartheta$ і $H$

$\mathrm {tg} \,{\frac {\vartheta }{2}}={\sqrt {\frac {e+1}{e-1}}}\,\mathrm {th} \,{\frac {H}{2}}$.

Завдяки такій підстановці інтеграл набуває такої ж форми, що й у випадку з еліптичною орбітою. Після проведення перетворень отримуємо гіперболічне рівняння Кеплера:

$M=e\,\mathrm {sh} \,H-H$

Величину $H$ називають гіперболічною ексцентричною аномалією. Оскільки $\mathrm {sh} \,H=-i\sin {iH}$, то останнє рівняння можна перетворити таким чином:

$M=-ei\sin {iH}-H=i\left(iH-e\sin {iH}\right)=i\left(E-e\sin E\right)$.

Звідси видно, що $E=iH$.

Параболічна орбіта

Рівняння параболи в полярних координатах має вигляд

$r={\frac {2\,r_{\pi }}{1+\cos \vartheta }}$

де $r_{\pi }$ — відстань до перицентра. Другий закон Кеплера для випадку руху по параболічній траєкторії

$r^{2}\,{\frac {d\vartheta }{dt}}={\rm {const}}={\sqrt {2\,\mu \,r_{\pi }}}$

Звідки одержуємо інтеграл, що визначає час руху

$t-t_{0}=2\,r_{\pi }\,{\sqrt {\frac {2\,r_{\pi }}{\mu }}}\int \limits _{0}^{\vartheta }{\frac {d\vartheta }{(1+\cos \vartheta )^{2}}}$

Вводимо універсальну тригонометричну заміну

$z={\rm {tg}}\,{\frac {\vartheta }{2}},\quad \vartheta =2\,{\rm {arctg}}\,z,\quad d\vartheta ={\frac {2\,dz}{1+z^{2}}},\quad \cos \vartheta ={\frac {1-z^{2}}{1+z^{2}}}$

і перетворюємо інтеграл

$t-t_{0}=4\,r_{\pi }\,{\sqrt {\frac {2\,r_{\pi }}{\mu }}}\int \limits _{0}^{{\rm {tg\,}}{\frac {\vartheta }{2}}}{\frac {\cfrac {dz}{1+z^{2}}}{\left(1+{\cfrac {1-z^{2}}{1+z^{2}}}\right)^{2}}}=r_{\pi }\,{\sqrt {\frac {2\,r_{\pi }}{\mu }}}\int \limits _{0}^{{\rm {tg\,}}{\frac {\vartheta }{2}}}(1+z^{2})\,dz=r_{\pi }\,{\sqrt {\frac {2\,r_{\pi }}{\mu }}}\left.\left(z+{\frac {z^{3}}{3}}\right)\right|_{0}^{{\rm {tg\,}}{\frac {\vartheta }{2}}}$

остаточно одержуємо

$t-t_{0}=r_{\pi }\,{\sqrt {\frac {2\,r_{\pi }}{\mu }}}\left({\rm {tg\,}}{\frac {\vartheta }{2}}+{\frac {1}{3}}{\rm {tg^{3}\,}}{\frac {\vartheta }{2}}\right)$

Останнє співвідношення відоме в небесній механіці як рівняння Баркера.

Радіальна орбіта

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

$v={\frac {dr}{dt}}$

Зв'язок між положенням тіла на орбіті і часом знайдемо з енергетичних міркувань

${\frac {m\,v^{2}}{2}}-{\frac {m\,\mu }{r}}={\rm {const}}$

$v^{2}={\frac {2\mu }{r}}+h$

— інтеграл енергії. Звідси маємо диференціальне рівняння

${\frac {dr}{dt}}=\pm {\sqrt {{\frac {2\mu }{r}}+h}}$

Розділяючи змінні в цьому рівнянні, приходимо до інтегралу

$\mp \left(t_{1}-t_{0}\right)=\int \limits _{r_{0}}^{r_{1}}{\frac {dr}{\sqrt {{\cfrac {2\mu }{r}}+h}}}$

спосіб обчислення якого визначається знаком константи $h$. Виділяють три випадки

• $h<0$ — прямолінійно-еліптична орбіта

Відповідає випадку, коли повна механічна енергія тіла від'ємна, і, віддалившись на деяку найбільшу відстань, від притягального центра, воно почне рухатися у зворотному напрямі. Це аналогічно руху по еліптичній орбіті. Для обчислення інтеграла введемо заміну

${\frac {2\,\mu }{r}}={\frac {-h}{\sin ^{2}u}},\quad u_{0}={\rm {arcsin}}{\sqrt {\frac {-h\,r_{0}}{2\,\mu }}},\quad u_{1}={\rm {arcsin}}{\sqrt {\frac {-h\,r_{1}}{2\,\mu }}},\quad dr={\frac {4\,\mu }{-h}}\sin u\cos u\,du$

обчислюємо інтеграл

$\mp \left(t_{1}-t_{0}\right)={\frac {4\,\mu }{-h{\sqrt {-h}}}}\int \limits _{u_{0}}^{u_{1}}\sin ^{2}u\,du={\frac {2\,\mu }{-h{\sqrt {-h}}}}\int \limits _{u_{0}}^{u_{1}}\left(1-\cos {2u}\right)\,du={\frac {\mu }{-h{\sqrt {-h}}}}\left.\left(2\,u-\sin 2u\right)\right|_{u_{0}}^{u_{1}}$

Вважаючи $E=2\,u$, запишемо результат

$\mp \left(t_{1}-t_{0}\right)={\frac {\mu }{-h{\sqrt {-h}}}}\left(E_{1}-E_{0}-\sin E_{1}+\sin E_{0}\right)$

прийнявши за (недосяжний в реальності) умовний перицентр $r_{0}=0$, і напрямок початкової швидкості від притягального центра, отримаємо так зване радіальне рівняння Кеплера, що зв'язує відстань від притягального центра з часом руху

$t_{1}-t_{0}={\frac {\mu }{-h{\sqrt {-h}}}}\left(E-\sin E\right)$

де $E=2\arcsin {\sqrt {\frac {-h\,r}{2\,\mu }}}$.

• $h=0$ — прямолінійно-параболічна орбіта

Занедбане радіально тіло піде на нескінченність від притягального центра, маючи на нескінченності швидкість, рівну нулю. Відповідає випадку руху з параболічною швидкістю. Найпростіший випадок, бо не вимагає заміни в інтегралі

$\mp \left(t_{1}-t_{0}\right)=\int \limits _{r_{0}}^{r_{1}}{\frac {dr}{\sqrt {\cfrac {2\mu }{r}}}}={\frac {1}{3}}{\sqrt {\frac {2}{\mu }}}\left(r_{1}{\sqrt {r_{1}}}-r_{0}{\sqrt {r_{0}}}\right)$

Беручи початкові умови першого випадку, отримуємо явний закон руху

$r(t)=\left[3{\sqrt {\frac {\mu }{2}}}\left(t_{1}-t_{0}\right)\right]^{\frac {2}{3}}$

• $h>0$ — прямолінійно-гіперболічна орбіта

Відповідає віддаленню від притягального центра на нескінченність. На нескінченності тіло буде мати швидкість, $v_{\infty }={\sqrt {h}}$. Вводимо заміну

${\frac {2\,\mu }{r}}={\frac {h}{{\rm {sh}}^{2}u}},\quad u_{0}={\rm {arcsh}}{\sqrt {\frac {h\,r_{0}}{2\,\mu }}},\quad u_{1}={\rm {arcsh}}{\sqrt {\frac {h\,r_{1}}{2\,\mu }}},\quad dr={\frac {4\,\mu }{h}}\,{\rm {sh}}u\,{\rm {ch}}u\,du$

і обчислюємо інтеграл

$\mp \left(t_{1}-t_{0}\right)={\frac {4\,\mu }{h{\sqrt {h}}}}\int \limits _{u_{0}}^{u_{1}}{\rm {sh}}^{2}u\,du={\frac {2\,\mu }{h{\sqrt {h}}}}\int \limits _{u_{0}}^{u_{1}}\left({\rm {ch}}{2u}-1\right)\,du={\frac {\mu }{h{\sqrt {h}}}}\left.\left({\rm {sh}}2u-2\,u\right)\right|_{u_{0}}^{u_{1}}$

Вважаючи $H=2\,u$, отримуємо

$\mp \left(t_{1}-t_{0}\right)={\frac {\mu }{h{\sqrt {h}}}}\left({\rm {sh}}\,H_{1}-{\rm {sh}}\,H_{0}-H_{1}+H_{0}\right)$

Вважаючи початкові умови аналогічними першому випадку, маємо гіперболічне радіальне рівняння Кеплера

$t_{1}-t_{0}={\frac {\mu }{h{\sqrt {h}}}}\left({\rm {sh}}H-H\right)$

де $H=2\,{\rm {arcsh}}{\sqrt {\frac {h\,r}{2\,\mu }}}$

Розв'язок рівняння Кеплера

Розв'язок рівняння Кеплера в еліптичному і гіперболічному випадках існує і єдиний за будь-яких дійсних M. Для колової орбіти (e = 0) рівняння Кеплера набуває тривіального вигляду М = E. В загальному вигляді рівняння Кеплера трансцендентне. Воно не розв'язується в алгебраїчних функціях. Однак, його розв'язок можна знайти різними способами за допомогою збіжних рядів. Загальний розв'язок рівняння Кеплера можна записати за допомогою рядів Фур'є:

$E=M+2\cdot \sum _{n=1}^{n}{\frac {1}{n}}J_{n}\left(n\,e\,\right)\cdot \sin {nM}$,

де

$J_{m}\left(x\right)={\frac {1}{\pi }}\int \limits _{0}^{\pi }\cos \left(mE-x\sin {E}\right)dE$

— функція Бесселя.

Цей ряд збігається, коли величина ε не перевищує значення границі Лапласа.

Наближені методи

Серед чисельних методів розв'язування рівняння Кеплера часто використовують метод нерухомої точки («метод простої ітерації») і метод Ньютона. Для еліптичного випадку в методі нерухомої точки за початкове значення E₀ можна взяти M, а послідовні наближення мають такий вигляд:

$E_{n+1}=e\,\sin E_{n}+M$

В гіперболічному випадку метод нерухомої точки подібним чином використовувати не можна, однак цей метод дає можливість вивести для такого випадку іншу формулу наближень (з гіперболічним арксинусом):

$H_{n+1}=\operatorname {Arsh} {\frac {H_{n}+M}{e}}$

Література

• Д. Е. Охоцимский, Ю. Г. Сихарулидзе. Основы механики космического полета. — Москва : «Наука», 1990.
• В. Е. Жаров. Сферическая астрономия. — Фрязино, 2006. — С. 480. — ISBN ISBN 5-85099-168-9.
• Г. М. Фихтенгольц. Курс дифференциального и интегрального исчисления. Том 3.
• Лукьянов Л.Г., Ширмин Г. И. Лекции по небесной механике. — Алматы, 2009. — С. 276.