Бисекция и метод Ньютона

Комета Галлея сейчас в 35 астрономических единицах от Солнца — дальше Нептуна — и медленно разворачивается: к Солнцу она вернётся в 2061 году. Её орбиту знают с огромной точностью, но чтобы сказать, где на этой орбите комета в нужный день, нужно решить уравнение Кеплера: M = E − e · sin E. Формулы для его решения нет, и больше четырёхсот лет его решают подбором — уточняют E, пока равенство не выполнится. Ниже пять настоящих комет и четыре способа подбора. Двигайте время, меняйте комету и способ.

В уравнении e — эксцентриситет орбиты: у почти круглой орбиты Земли он 0,017, у кометы Галлея 0,968, у кометы NEOWISE 0,9992 — её эллипс почти парабола. M — средняя аномалия, она растёт с постоянной скоростью от перигелия до перигелия, это просто время в других единицах. E — эксцентрическая аномалия, угол, по которому положение на эллипсе находится сразу: расстояние до Солнца равно a(1 − e cos E). На втором графике — f(E) = E − e · sin E − M: корень там, где кривая пересекает ось. Бисекция находит его всегда за 43 шага, Ньютон — за считанные: у кометы Галлея сегодня за четыре. Но переключитесь на комету Хейла — Боппа или NEOWISE и посмотрите на полосу под ползунком: красные отметки — дни, когда Ньютон со стандартного старта E₀ = M так и не сходится.

Идея

Бисекция: если непрерывная функция на концах отрезка имеет разные знаки, корень внутри. Делим отрезок пополам и оставляем ту половину, где знаки по-прежнему разные. Каждый шаг вдвое сужает место, где прячется корень, — гарантированно, какой бы ни была функция. Метод Ньютона: вблизи точки график почти совпадает с касательной. Где касательная пересекает ось, найти легко — туда и шагаем: x ← x − f(x)/f′(x). Рядом с корнем ошибка на каждом шаге возводится в квадрат, и число верных знаков удваивается. Но только рядом с корнем.

У двух методов противоположные характеры, и это главное, что о них нужно знать. Бисекции нужны отрезок с разными знаками на концах и умение вычислять знак f; производная не нужна, гладкость не нужна, плохого старта не бывает. Зато скорость одна и та же — бит за шаг — и для простой функции, и для сложной. Ньютону нужны производная и хорошая начальная точка, и за это он платит скоростью, которой у бисекции нет. Школьная сторона обоих — теорема о промежуточном значении и доказательство квадратичной сходимости — в «Царице наук»: глава о пределах и глава об экстремумах.

Тысяча знаков

Возьмём √2, корень x² − 2. Бисекция на отрезке [1, 2] и Ньютон из x₀ = 1, шаг за шагом:

step   bisection           digits   Newton              digits
   1   1.250000000000000       0   1.500000000000000       1
   2   1.375000000000000       1   1.416666666666667       2
   3   1.437500000000000       1   1.414215686274510       5
   4   1.406250000000000       2   1.414213562374690      11
   5   1.421875000000000       2   1.414213562373095      16
bisection needs 50 steps for all 16 digits, Newton 5

Пять шагов Ньютона — и все 16 знаков, которые помещаются в число с плавающей точкой. Бисекции для того же нужно 50 шагов. Откуда удвоение? Пусть r — корень и x = r + ε. Разложим f около корня: f(x) ≈ f′(r)ε + ½f″(r)ε². Подставив в шаг Ньютона, получим новую ошибку

Для √2 множитель равен 1/(2√2) ≈ 0,35, и ошибка 10⁻⁵ превращается в 3,5 · 10⁻¹¹: пять верных знаков становятся десятью-одиннадцатью. Строгое доказательство с оценкой остатка — в «Царице наук».

В числах с плавающей точкой удвоению некуда расти дальше шестнадцати знаков. А в длинной арифметике есть. Ниже Ньютон считает √2, ∛2 и 1/7 с тысячей знаков после запятой, а рядом бисекция делает столько же шагов.

У Ньютона для √2 верных знаков 1, 2, 5, 11, 24, 48, 97, 195, 391, 783 — каждый шаг примерно вдвое больше, и на одиннадцатом шаге их больше тысячи. Бисекция за те же одиннадцать шагов добывает три-четыре знака: шаг даёт ей ровно один двоичный знак, а десятичный стоит log₂ 10 ≈ 3,32 шага. Обратите внимание на 1/7: шаг Ньютона для уравнения 1/x − 7 = 0 — это x ← x(2 − 7x), в нём нет ни одного деления. Так делят числа там, где деление дорогое: в длинной арифметике и даже в обычных программах — компилятор GCC с флагом -mrecip заменяет деление чисел одинарной точности на быструю приближённую обратную величину плюс один шаг Ньютона.

Почему бисекция не промахивается

Число шагов не зависит от функции — только от ширины отрезка и нужной точности. Отрезок [0, 2π] для уравнения Кеплера сжимается до 10⁻¹² за ⌈log₂(2π · 10¹²)⌉ = 43 шага — для кометы Чурюмова — Герасименко и для NEOWISE одинаково, вы это видели в сцене. Это тот же двоичный поиск, только по вещественным числам вместо индексов массива, и тот же приём, что двоичный поиск по ответу.

Код

Оба метода в общем виде: бисекции нужны функция и отрезок, Ньютону — функция, производная и старт. Программа печатает таблицу для √2, которую вы видели выше, и два примера, где Ньютон теряет скорость или уходит прочь, — о них ниже.

Остановка — отдельный вопрос. Бисекцию останавливают по ширине отрезка: она известна точно. Ньютона — по величине шага: если поправка меньше допуска, следующие будут ещё меньше. Сравнивать x с предыдущим на равенство нельзя: в числах с плавающей точкой Ньютон может качаться между двумя соседними числами. Здесь так и выходит — на шестом шаге √2 отличается от правильного double на одну единицу последнего разряда.

А это уравнение Кеплера на настоящих данных: элементы орбит пяти комет из базы JPL и четыре способа решения на 16 октября 2026 года.

16 Oct 2026            e        M°       E°   r, au   bisect  N(E0=M)  N(E0=pi)  Kepler
1P              0.967936  192.925  186.575   35.17      43        4         4    667
C/2020 F3       0.999178    0.333   18.473   18.75      43       14        10    466

Последний столбец — способ самого Кеплера из «Сокращения коперниковой астрономии» 1621 года: подставлять E в правую часть E = M + e · sin E снова и снова. Он сходится всегда, но медленно — ошибка каждый раз умножается примерно на e · |cos E|. Комете Галлея сейчас, вблизи афелия, где cos E ≈ −1, нужно 667 шагов: множитель 0,96 съедает ошибку очень неторопливо. Ньютону хватает четырёх.

Где Ньютон ломается

Обещание «знаки удваиваются» дано с оговоркой «вблизи простого корня». Каждое слово в ней важно.

Далеко от корня. У уравнения Кеплера производная 1 − e · cos E около перигелия близка к нулю, если e близко к единице: касательная почти горизонтальна и отбрасывает очень далеко. Классический пример из статьи Чарлза и Тейтума 1997 года — e = 0,992, M = 0,13π:

e = 0.992, M = 0.13pi, Newton from M: 0.4084 4.806 -1.13 -0.0178 50.07 -1761 202.2
from pi: 7 steps; bisection: 43 steps; E = 1.382957944863

Корень — 1,383, а итерации прыгают по сотням и тысячам радиан. Дальше шестого шага программа путь не печатает: он зависит от последнего бита синуса, и Python с Go, считая одно и то же, расходятся. Старт E₀ = π решает задачу за семь шагов. Программа ниже проверяет это на всех пяти кометах: берёт 4000 равномерно расставленных значений M и пускает Ньютона с двух стартов, не больше 50 шагов.

Для первых трёх комет оба старта сходятся всегда. У Хейла — Боппа (e = 0,995) старт E₀ = M в нашем запуске на Python не сошёлся в 58 случаях из 4000, у NEOWISE — в 91; старт π — ни разу и не дольше чем за 12 шагов. Версия на Go насчитала 70 и 94. Код тот же, но синус и косинус в Go и в системной библиотеке C, которой пользуется Python, — разные программы, и они изредка расходятся в последнем бите, а в этой области траектория настолько чувствительна, что последний бит решает, сойдётся ли старт. Поэтому у этой программы нет эталонного вывода, и на другой машине числа могут немного отличаться. Отсюда практическое правило астрономических программ: при большом e начинать с π, а надёжнее — страховать Ньютона бисекцией.

Кратный корень. Если график касается оси, а не пересекает её, f′(r) = 0, и квадратичное сжатие ошибки пропадает. Для (x − 1)² из x₀ = 2 ошибка лишь делится пополам: 0,5, 0,25, 0,125, … — та же скорость, что у бисекции, только без её гарантий.

Нет гладкости. У ∛x корень простой, но производная в нём бесконечна. Касательная из точки x пересекает ось в −2x, и метод уходит прочь: −2, 4, −8, 16, …

Циклы. У x³ − 2x + 2 Ньютон из нуля попадает в единицу, из единицы — обратно в ноль, и так без конца. Сколько таких ловушек, лучше всего видно на комплексной плоскости.

Фракталы Ньютона

Метод Ньютона работает и для комплексных чисел: та же формула z ← z − p(z)/p′(z). У z³ − 1 три корня: 1 и (−1 ± i√3)/2. Покрасим каждую точку плоскости в цвет корня, к которому из неё приходит метод, и сделаем ярче там, где он приходит быстро.

Около каждого корня — спокойная область своего цвета: оттуда метод сходится быстро и туда, куда ждёшь. А границы между областями устроены бесконечно сложно. Увеличьте любую: там, где встречаются два цвета, всегда есть и третий, и так на любом масштабе. Это не особенность картинки, а теорема: у всех трёх областей одна и та же граница, множество Жюлиа отображения z ↦ z − p(z)/p′(z), и в любой окрестности любой её точки есть старты, ведущие к каждому из трёх корней. Переключитесь на z³ − 2z + 2: чёрные острова вокруг 0 и 1 — те самые старты, которые затягивает в цикл.

Та же картинка в виде текста — программа ставит в каждую клетку знак корня, к которому из неё приходит метод:

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

Сколько это стоит

Бисекции нужно ⌈log₂((b − a)/ε)⌉ вычислений знака f, и это число известно до начала работы. Ньютону для d верных знаков из удачного старта нужно около log₂ d шагов: 5 шагов на 16 знаков, 11 — на тысячу. Но шаг Ньютона дороже: он вычисляет и f, и f′. Если производную считать дорого, берут метод секущих — вместо касательной прямую через две последние точки; знаки у него растут не вдвое, а в 1,618 раза за шаг, зато вычисление одно.

В длинной арифметике есть ещё одна хитрость. На первых шагах верных знаков мало, и считать их с тысячей знаков бессмысленно. Поэтому точность наращивают вместе с шагами: 2 знака, 4, 8, …, 1000. Последний шаг стоит больше, чем все предыдущие вместе, и весь корень обходится в несколько умножений тысячезначных чисел.

На практике оба характера соединяют. Надёжные решатели держат отрезок с разными знаками, как бисекция, но шагают быстрыми методами, а если быстрый шаг выпрыгивает из отрезка или сжимает его слишком медленно — делают шаг бисекции. Так устроен метод Брента 1973 года; в SciPy это функция brentq, которую документация называет лучшим из своих решателей с отрезком.

Где это работает

• Корни и деление в длинной арифметике. В Go big.Int.Sqrt повторяет шаг Ньютона z ← ⌊(z + ⌊x/z⌋)/2⌋ из заведомо большого старта, а big.Float.Sqrt решает уравнение 1/t² − x = 0, чтобы обойтись без деления. В справочной реализации модуля decimal на самом Python (_pydecimal) целую часть корня тоже находят методом Ньютона.
• Быстрый обратный квадратный корень. В игре Quake III Arena 1999 года 1/√x получали из особой константы 0x5F3759DF и одного шага Ньютона; второй шаг в исходниках был, но закомментирован — точности хватало.
• Небесная механика. Положение планет, комет и спутников по элементам орбиты — это уравнение Кеплера на каждый момент времени, решённое Ньютоном или его родственниками.
• Финансы. Внутреннюю норму доходности — ставку, при которой платежи окупают вложение, — функция IRR в Excel ищет итерациями: до точности 0,00001 % и не больше 20 попыток, иначе ошибка #NUM!. Первая задача ниже — про то же самое.
• Подбор параметра. Любое «при какой скорости, ставке, температуре выполнится условие» с монотонной зависимостью — бисекция; именно поэтому решатели вроде brentq требуют отрезок с разными знаками.

История

Само уравнение вывел Иоганн Кеплер в 1609 году в «Новой астрономии», а в 1621 году предложил решать его подстановкой — той самой итерацией из последнего столбца таблицы. Кто придумал «метод Ньютона», спорят до сих пор. Ньютон в рукописи 1669 года, напечатанной только в 1711-м, уточнял корень многочлена поправками; в печать этот приём впервые попал в «Трактате об алгебре» Джона Уоллиса 1685 года. Джозеф Рафсон в 1690 году записал его как повторяемую формулу — отсюда второе имя, метод Ньютона — Рафсона. В нынешнем виде, с производной и для любой функции, его опубликовал Томас Симпсон в 1740 году, и историк Ник Коллерстром в 1992 году доказывал, что метод правильнее считать симпсоновским. Сам способ извлекать квадратный корень, x ← (x + 2/x)/2, гораздо старше: он есть в «Метрике» Герона, и о нём подробно рассказывает глава о корне из двух «Царицы наук».

Задачи

В первой задаче корень ищут для денег: какую на самом деле ставку платит покупатель в рассрочку. Во второй — уравнение Кеплера с тестами, где старт E₀ = M не сходится.