Уравнение Кеплера без осечек
Программа-планетарий показывает кометы на небе. Чтобы узнать, где комета на орбите в нужный момент, ей нужно решить уравнение Кеплера
$$E - e \sin E = M$$относительно эксцентрической аномалии E. Здесь e — эксцентриситет орбиты (0 ≤ e < 1), M — средняя аномалия, то есть время, прошедшее от перигелия, в долях оборота, умноженное на 2π. Решение на [0, 2π) единственное. Планетарий решает уравнение для каждой кометы много раз в секунду и не имеет права ошибиться ни разу: комета, прыгнувшая не туда, видна сразу.
Вход
В первой строке число запросов q (1 ≤ q ≤ 40 000). В каждой из следующих q строк — два вещественных числа e и M (0 ≤ e ≤ 0,9999, 0 ≤ M < 2π).
Выход
Для каждого запроса — E из [0, 2π). Ответ принимается с абсолютной или относительной погрешностью 10⁻⁹.
Примеры
Пример 1
4 0.0000000000 1.0000000000 0.0167000000 2.0000000000 0.9679000000 3.3602000000 0.9920000000 0.4084070450
1.000000000000 2.015078752618 3.252791912317 1.382957944904
Ограничения
- Время · Python
- 5000 мс
- Время · Go
- 1000 мс
- Память
- 256 МБ
- Проверка
- 6 тестов
Разбор
Открывайте по ступеням: каждая следующая подсказывает больше. Можно сразу к решению — но задача, решённая с одной подсказкой, запоминается лучше.
1Подсказка 1
Стандартный рецепт — метод Ньютона из E₀ = M. Попробуйте его на четвёртой строке примера, e = 0,992, M = 0,13π, и напечатайте первые десять приближений.
2Подсказка 2
Функция f(E) = E − e sin E − M возрастает (её производная 1 − e cos E ≥ 1 − e > 0), f(0) ≤ 0 и f(2π) > 0. Какой метод с такими условиями не может промахнуться?
3Идея
Корень всегда лежит на [0, 2π], и знаки на концах разные, поэтому бисекция находит его гарантированно — за 43 шага до 10⁻¹². Ньютон быстрее, но при e около 0,99 и малых M производная у старта почти ноль, и касательная выбрасывает приближение на сотни радиан. Надёжное решение соединяет оба: держим отрезок [lo, hi] с разными знаками на концах, делаем шаг Ньютона, а если он вылетел за отрезок — шаг бисекции.
4Разбор
На каждом шаге вычисляем f(E) и по её знаку сдвигаем одну из границ: отрицательная — корень правее, lo = E; положительная — hi = E. Отрезок только сужается и всегда содержит корень. Следующее приближение — шаг Ньютона, если он попал строго внутрь (lo, hi); иначе середина отрезка. Останавливаемся, когда шаг меньше 10⁻¹⁵ или отрезок стал короче.
Старт берём π при e > 0,8, как советуют справочники по небесной механике: в наших проверках из π метод Ньютона не подвёл ни разу. Но полагаться на это не нужно — страховка бисекцией делает решение надёжным при любом старте.
Третий тест — настоящие кометы из базы JPL на 1 января каждого года с 2000 по 2060. Голый Ньютон из M ошибается на нём одиннадцать раз — у комет Хейла — Боппа и NEOWISE, чьи эксцентриситеты 0,995 и 0,9992.
4+Решение
5Сложность
Подстрахованный Ньютон почти всегда делает 4–12 шагов, а там, где касательная выбрасывает за отрезок, его выручают шаги бисекции, каждый из которых вдвое сужает отрезок. Итого O(q) с небольшой константой.
6Типичные ошибки
- Ньютон из E₀ = M без страховки. Сходится почти всегда, а при e ≥ 0,99 и M до 0,5 иногда не сходится вовсе. Четвёртый и пятый тесты собраны из таких случаев.
- Фиксированные пять шагов Ньютона. Для круговых орбит хватает, у почти параболических — нет.
- Итерация Кеплера E ← M + e sin E. Сходится всегда, но при e = 0,9999 ошибка уменьшается за шаг лишь в 1,0001 раза: на 10⁻⁹ уйдут сотни тысяч шагов на запрос.
- Ответ вне [0, 2π). E + 2π не решение: левая часть уравнения тоже вырастет на 2π.