Уравнение Кеплера

Анимация, иллюстрирующая истинную аномалию, эксцентрическую аномалию, среднюю аномалию и решение уравнения Кеплера (в правом верхнем углу), эксцентриситет — 0,6.

Уравне́ние Ке́плера описывает движение тела по эллиптической орбите в задаче двух тел и имеет вид:

[math]\displaystyle{ E- e \, \sin E = M }[/math]

где [math]\displaystyle{ E }[/math] — эксцентрическая аномалия, [math]\displaystyle{ e }[/math] — эксцентриситет орбиты, а [math]\displaystyle{ M }[/math] — средняя аномалия.

Впервые это уравнение было получено астрономом Иоганном Кеплером в 1619 году. Играет значительную роль в небесной механике.

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

Уравнение Кеплера в классической форме описывает движение только по эллиптическим орбитам, то есть при [math]\displaystyle{ 0 \le e \lt 1 }[/math]. Движение по гиперболическим орбитам [math]\displaystyle{ (e \gt 1) }[/math] подчиняется гиперболическому уравнению Кеплера, сходному по форме с классическим. Движение по прямой линии [math]\displaystyle{ (e = -1) }[/math] описывается радиальным уравнением Кеплера. Наконец, для описания движения по параболической орбите [math]\displaystyle{ (e = 1) }[/math] используют уравнение Баркера. При [math]\displaystyle{ e \lt 0 }[/math] орбит не существует.

Задача, приводящая к уравнению Кеплера

Рассмотрим движение тела по орбите в поле другого тела. Найдем зависимость положения тела на орбите от времени. Из II закона Кеплера следует, что

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

Здесь [math]\displaystyle{ r }[/math] — расстояние от тела до гравитирующего центра, [math]\displaystyle{ \vartheta }[/math] — истинная аномалия — угол между направлениями на перицентр орбиты и на тело, [math]\displaystyle{ \mu = G M_0 }[/math] — произведение постоянной тяготения на массу гравитирующего тела, [math]\displaystyle{ a }[/math] — большая полуось орбиты. Отсюда можно получить зависимость времени движения по орбите от истинной аномалии:

[math]\displaystyle{ t - t_0 = \frac{1}{\sqrt{\mu a \left(1 - e^2\right)}} \int\limits_0^\vartheta r^2 d\vartheta }[/math].

Здесь [math]\displaystyle{ t_0 }[/math] — время прохождения через перицентр.

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

Эллиптическая орбита

Уравнение эллипса в полярных координатах имеет вид

[math]\displaystyle{ r = \frac{a(1-e^2)}{1+e\cos{\vartheta}} }[/math]

Тогда уравнение для времени приобретает вид

[math]\displaystyle{ 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} }[/math]

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

[math]\displaystyle{ \operatorname{tg}\frac {\vartheta} {2} = \sqrt{\frac {1+e} {1-e}}\cdot\operatorname{tg}\frac {E} {2} }[/math]

Величина E называется эксцентрической аномалией. Благодаря такой подстановке интеграл легко берется. Получается следующее уравнение:

[math]\displaystyle{ t - t_0 = \sqrt{\frac{a^3}{\mu}}\left(E - e\sin E\right) }[/math]

Величина [math]\displaystyle{ \sqrt{\frac{\mu}{a^3}} }[/math] является средней угловой скоростью движения тела по орбите. В небесной механике для этой величины используется термин среднее движение. Произведение среднего движения на время называется средней аномалией M. Эта величина представляет собой угол, на которой повернулся бы радиус-вектор тела, если бы оно двигалось по круговой орбите с радиусом, равным большой полуоси орбиты тела.

Таким образом получаем уравнение Кеплера для эллиптического движения:

[math]\displaystyle{ E-e\sin E = M }[/math]

Гиперболическая орбита

Уравнение гиперболы в полярных координатах имеет тот же вид, что и уравнение эллипса. Значит, интеграл получается такой же по виду. Однако, использовать эксцентрическую аномалию в данном случае нельзя. Воспользуемся параметрическим представлением гиперболы: [math]\displaystyle{ x = -a\,\mathrm{ch}\,H }[/math], [math]\displaystyle{ y = a \sqrt{e^2 -1}\,\mathrm{sh}\,H }[/math]. Тогда уравнение для гиперболы принимает вид

[math]\displaystyle{ r = a \left(e\,\mathrm{ch}\,H - 1\right) }[/math],

а связь между [math]\displaystyle{ \vartheta }[/math] и [math]\displaystyle{ H }[/math]

[math]\displaystyle{ \mathrm{tg}\,\frac{\vartheta}{2} = \sqrt{\frac{e +1}{e -1}}\,\mathrm{th}\,\frac{H}{2} }[/math].

Благодаря такой подстановке интеграл приобретает ту же форму, что и в случае с эллиптической орбитой. После произведения преобразований получаем гиперболическое уравнение Кеплера:

[math]\displaystyle{ M = e\,\mathrm{sh}\,H - H }[/math]

Величина [math]\displaystyle{ H }[/math] называется гиперболической эксцентрической аномалией. Поскольку [math]\displaystyle{ \mathrm{sh}\,H = -i\sin{iH} }[/math], то последнее уравнение можно преобразовать следующим образом:

[math]\displaystyle{ M = -ei\sin{iH} -H = i\left(iH - e\sin{iH}\right) = i\left(E-e\sin E\right) }[/math].

Отсюда видно, что [math]\displaystyle{ E = i H }[/math].

Параболическая орбита

Уравнение параболы в полярных координатах имеет вид

[math]\displaystyle{ r = \frac{2\,r_{\pi}}{1+\cos \vartheta} }[/math]

где [math]\displaystyle{ r_{\pi} }[/math] — расстояние до перицентра. Второй закон Кеплера для случая движения по параболической траектории

[math]\displaystyle{ r^2 \, \frac{d\vartheta}{dt} = {\rm const} = \sqrt{2\, \mu \, r_{\pi}} }[/math]

Откуда получаем интеграл, определяющий время движения

[math]\displaystyle{ t - t_0 = 2 \, r_{\pi} \, \sqrt{\frac{2\,r_{\pi}}{\mu}} \int\limits_{0}^{\vartheta} \frac{d\vartheta}{(1+\cos\vartheta)^2} }[/math]

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

[math]\displaystyle{ 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} }[/math]

и преобразуем интеграл

[math]\displaystyle{ 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}} }[/math]

получаем окончательно

[math]\displaystyle{ 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) }[/math]

Последнее соотношение известно в небесной механике как уравнение Баркера.

Радиальная орбита

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

[math]\displaystyle{ v = \frac{dr}{dt} }[/math]

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

[math]\displaystyle{ \frac{m \, v^2}{2} - \frac{m \, \mu}{r} = {\rm const} }[/math]

[math]\displaystyle{ v^2 = \frac{2\mu}{r} + h }[/math]

— интеграл энергии. Отсюда имеем дифференциальное уравнение

[math]\displaystyle{ \frac{dr}{dt} = \pm \sqrt{\frac{2\mu}{r} + h} }[/math]

Разделяя переменные в этом уравнении, приходим к интегралу

[math]\displaystyle{ \mp \left(t_1 - t_0\right) = \int\limits_{r_0}^{r_1} \frac{dr}{\sqrt{\cfrac{2\mu}{r} + h}} }[/math]

способ вычисления которого определяется знаком константы [math]\displaystyle{ h }[/math]. Выделяют три случая

• [math]\displaystyle{ h \lt 0 }[/math] прямолинейно-эллиптическая орбита

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

[math]\displaystyle{ \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 }[/math]

вычисляем интеграл

[math]\displaystyle{ \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} }[/math]

Полагая [math]\displaystyle{ E = 2\,u }[/math], запишем результат

[math]\displaystyle{ \mp \left(t_1 - t_0\right) = \frac{\mu}{-h\sqrt{-h}} \left(E_1 - E_0 - \sin E_1 + \sin E_0 \right) }[/math]

приняв в качестве (недостижимого в реальности) условного перицентра [math]\displaystyle{ r_0 = 0 }[/math], и направление начальной скорости от притягивающего центра, получим так называемое радиальное уравнение Кеплера, связывающее расстояние от притягивающего центра со временем движения

[math]\displaystyle{ t_1 - t_0 = \frac{\mu}{-h\sqrt{-h}} \left(E - \sin E\right) }[/math]

где [math]\displaystyle{ E = 2 \arcsin \sqrt{\frac{-h\,r}{2\,\mu}} }[/math].

• [math]\displaystyle{ h = 0 }[/math] прямолинейно-параболическая орбита

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

[math]\displaystyle{ \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) }[/math]

Принимая начальные условия первого случая, получаем явный закон движения

[math]\displaystyle{ r(t) = \left[3 \sqrt{\frac{\mu}{2}} \left(t_1 - t_0 \right) \right]^{\frac{2}{3}} }[/math]

• [math]\displaystyle{ h \gt 0 }[/math] прямолинейно-гиперболическая орбита

Соответствует уходу от притягивающего центра на бесконечность. На бесконечности тело будет иметь скорость, [math]\displaystyle{ v_{\infty} = \sqrt{h} }[/math]. Вводим замену

[math]\displaystyle{ \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 }[/math]

и вычисляем интеграл

[math]\displaystyle{ \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} }[/math]

Полагая [math]\displaystyle{ H = 2\,u }[/math], получаем

[math]\displaystyle{ \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) }[/math]

Полагая начальные условия аналогичными первому случаю, имеем гиперболическое радиальное уравнение Кеплера

[math]\displaystyle{ t_1 - t_0 = \frac{\mu}{h\sqrt{h}} \left({\rm sh} H - H\right) }[/math]

где [math]\displaystyle{ H = 2 \, {\rm arcsh} \sqrt{\frac{h\,r}{2\,\mu}} }[/math]

Решение уравнения Кеплера

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

[math]\displaystyle{ E = M + 2\cdot\sum_{n=1}^{n} \frac{1}{n}J_n\left(n\, e \, \right)\cdot\sin{nM} }[/math],

где

[math]\displaystyle{ J_m\left(x\right) = \frac{1}{\pi}\int\limits_0^\pi \cos\left(mE - x\sin{E}\right) dE }[/math]

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

Этот ряд сходится, когда величина ε не превышает значения предела Лапласа.

Приближённые методы

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

[math]\displaystyle{ E_{n+1}=e \, \sin E_n+M }[/math]

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

[math]\displaystyle{ H_{n+1}=\operatorname{Ar sh}\frac{H_n+M}{e} }[/math]

Литература

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