Слышал я, что Луна удаляется от Земли. (Espenak, F.; Meeus, J. Secular acceleration of the Moon. NASA (February 7, 2007))
Лагранжев формализм и немного терпения позволят нам составить уравнения для ограниченной круговой задачи трех тел.

– радиус орбиты Земли,

– гравитационный параметр Земли,

- гравитационный параметр Солнца,

- частота обращения Земли вокруг Солнца,

- геоцентрическая долгота Луны,

- фаза Луны,

- расстояние от Луны до Солнца,

- расстояние Земля-Луна.


Теперь, если задать скорость Луны в перигее скажем 1023 м/с, кривые

и

будут весьма похожи на правду, которая за окном и на сайте JPL NASA.

Уменьшая скорость Луны, мы увидим, что эксцентриситет ее геоцентрической орбиты уменьшается, чуть не до нуля, а потом снова начинает расти. Что и понятно. Допустим, в задаче двух тел пробное тело имело эллиптическую орбиту с энергией больше «круговой». Потеря энергии на торможение, например, сделает орбиту более округлой, затем круговой, а потом перигей станет апогеем, а новый перигей опустится ниже поверхности Земли и конец.
Возьмем невозмущенную задачу Кеплера и «возмутим» ее диссипативной силой вида

(это векторное равенство).

- скорость пробного тела. Переходя к обобщенным силам, получим нечто вроде:


где


Поправьте меня, если что не так. И есть ли надежда упростить и решить эти уравнения в аналитике?
Ладно. Решим численно. (И в Декарте, прошу прощения, пока не уверен в правильности формул для обобщенных сил). Картина движения следующая: эксцентриситет сначала падает, потом растет. При этом большая полуось почти не меняется, потому что апогей поднимается, а перигей опускается. Тем не менее, среднее расстояние сначала падает, потом, после прохождения «круговой стадии» растет.

Так вот. Луна пока не прошла круговую стадию. Среднее расстояние падает, получается…
Диссипативные силы приливного трения имеют неизвестную мне зависимость (кстати какую?). Может быть, картина движения будет другая?