Алгоритм Дженкинс-Трауб для полиномиальных нулей является быстрым глобально сходящимся итерационным Полиномом корня ознакомительного метода опубликован в 1970 годом Майкл А. Дженкинс и Джозеф Ф. Трауб . Они предложили два варианта: один для общих многочленов с комплексными коэффициентами, широко известный как алгоритм «CPOLY», и более сложный вариант для частного случая многочленов с действительными коэффициентами, широко известный как алгоритм «RPOLY». Последнее является «практически стандартом для полиномиальных кортежей черного ящика». [1]
В этой статье описан сложный вариант. Учитывая многочлен P ,
с комплексными коэффициентами вычисляет приближения к n нулямиз Р ( г ), по одному за раз в примерно порядке возрастания величины. После вычисления каждого корня его линейный множитель удаляется из полинома. Использование этого дефлятирования гарантирует, что каждый корень вычисляется только один раз и что все корни найдены.
Реальный вариант следует тому же шаблону, но вычисляет два корня одновременно, либо два действительных корня, либо пару сопряженных комплексных корней. Избегая сложной арифметики, реальный вариант может быть быстрее (в 4 раза), чем сложный вариант. Алгоритм Дженкинса – Трауба стимулировал значительные исследования теории и программного обеспечения для методов этого типа.
Обзор
Алгоритм Дженкинса – Трауба вычисляет все корни многочлена с комплексными коэффициентами. Алгоритм начинается с проверки полинома на наличие очень больших или очень маленьких корней. При необходимости коэффициенты масштабируются путем изменения масштаба переменной. В алгоритме правильные корни находятся один за другим и, как правило, в увеличивающемся размере. После того, как каждый корень найден, полином дефлируется путем деления соответствующего линейного множителя. Действительно, факторизация полинома в линейный множитель и оставшийся спущенный полином уже является результатом процедуры поиска корня. Процедура поиска корня состоит из трех этапов, соответствующих различным вариантам итерации обратной степени . См. Дженкинс и Трауб . [2] Описание также можно найти у Ральстона и Рабиновица [3] с. 383. Алгоритм по духу аналогичен двухэтапному алгоритму, исследованному Траубом. [4]
Процедура поиска корней
Начиная с текущего полинома P ( X ) степени n , вычисляется наименьший корень P (x) . С этой целью строится последовательность так называемых H- полиномов. Все эти многочлены имеют степень n - 1 и должны сходиться к множителю P ( X ), содержащему все оставшиеся корни. Последовательность многочленов H встречается в двух вариантах: ненормализованный вариант, который позволяет легко получить теоретическое понимание, и нормализованный вариант полиномы, которые удерживают коэффициенты в числовом значении.
Построение полиномов H зависит от последовательности комплексных чисел называется смены. Сами эти сдвиги зависят, по крайней мере на третьем этапе, от предыдущих многочленов H. В H полиномы определяются как решение неявного рекурсии
- а также
Прямым решением этого неявного уравнения является
где полиномиальное деление точное.
Алгоритмически можно было бы использовать, например, схему Хорнера или правило Руффини для вычисления многочленов ви одновременно получить частные. Используя полученные частные p ( X ) и h ( X ) в качестве промежуточных результатов, следующий полином H получается как
Поскольку коэффициент наивысшей степени получается из P (X) , старший коэффициент является . Если это разделить, нормализованный многочлен H равен
Этап первый: бессменный процесс
Для набор . Обычно M = 5 выбирается для полиномов умеренных степеней вплоть до n = 50. Этот этап не является необходимым только из теоретических соображений, но полезен на практике. Он подчеркивает в полиномах H кофактор (линейного множителя) наименьшего корня.
Этап второй: процесс фиксированной смены
Сдвиг для этого этапа определяется как некоторая точка, близкая к наименьшему корню многочлена. Он квазислучайно расположен на окружности с внутренним радиусом корня, который, в свою очередь, оценивается как положительное решение уравнения
Поскольку левая часть является выпуклой функцией и монотонно возрастает от нуля до бесконечности, это уравнение легко решить, например, методом Ньютона .
Теперь выберите на окружности этого радиуса. Последовательность многочленов, , генерируется с фиксированным значением сдвига . Во время этой итерации текущее приближение для корня
прослеживается. Второй этап успешно завершен, если условия
- а также
одновременно встречаются. Если после некоторого количества итераций успеха не было, пробуется другая случайная точка на круге. Обычно используется 9 итераций для полиномов умеренной степени со стратегией удвоения в случае множественных отказов.
Этап третий: процесс переменной смены
В теперь генерируются с использованием переменных сдвигов которые генерируются
является последней оценкой корня второго этапа и
- где - нормированный многочлен H , т. е. делится на его старший коэффициент.
Если размер шага на третьем этапе не падает достаточно быстро до нуля, второй этап запускается заново с другой случайной точки. Если это не удается после небольшого количества перезапусков, количество шагов на втором этапе удваивается.
Конвергенция
Можно показать , что, при условии , L выбран достаточно большой, ево λ всегда сходится к корню Р .
Алгоритм сходится для любого распределения корней, но может не найти все корни многочлена. Более того, сходимость немного быстрее, чем квадратичная сходимость итерации Ньютона – Рафсона, однако она использует как минимум вдвое больше операций на шаг.
В чем сила алгоритма?
Сравните с итерацией Ньютона – Рафсона
Итерация использует заданные P и. В отличие от третьего этапа Дженкинса – Трауба
в точности итерация Ньютона – Рафсона, выполняемая над некоторыми рациональными функциями . Точнее, Ньютон – Рафсон выполняется над последовательностью рациональных функций
Для достаточно большой,
как можно ближе к полиному первой степени
где один из нулей . Несмотря на то, что стадия 3 представляет собой итерацию Ньютона – Рафсона, дифференцирование не выполняется.
Анализ полиномов H
Позволять - корни P ( X ). Так называемые факторы Лагранжа P (X) являются кофакторами этих корней,
Если все корни различны, то множители Лагранжа образуют базис пространства многочленов степени не выше n - 1. Анализ процедуры рекурсии показывает, что многочлены H имеют координатное представление
Каждый фактор Лагранжа имеет старший коэффициент 1, так что старший коэффициент полиномов H является суммой коэффициентов. Таким образом, нормированные полиномы H имеют вид
Приказы конвергенции
Если условие выполняется почти для всех итераций, нормализованные полиномы H будут сходиться по крайней мере геометрически к .
При условии, что
получаем асимптотические оценки для
- этап 1:
- для этапа 2, если s достаточно близко к:
- а также
- и для этапа 3:
- а также
- приводя к более высокому порядку сходимости, чем квадратичный , где это золотое сечение .
Интерпретация как итерация обратной мощности
Все этапы комплексного алгоритма Дженкинса – Трауба можно представить как задачу линейной алгебры определения собственных значений специальной матрицы. Эта матрица является координатным представлением линейного отображения в n- мерном пространстве многочленов степени n - 1 или меньше. Основная идея этой карты - интерпретировать факторизацию
с корнем а также оставшийся множитель степени n - 1 в качестве уравнения на собственный вектор для умножения с переменной X с последующим вычислением остатка с делителем P ( X ),
Это отображает полиномы степени не выше n - 1 в полиномы степени не выше n - 1. Собственные значения этого отображения являются корнями P ( X ), поскольку уравнение для собственных векторов имеет вид
откуда следует, что , это, является линейным множителем в P ( X ). В мономиальном базисе линейное отображениепредставлена сопутствующей матрицей полинома P , как
результирующая матрица коэффициентов
К этой матрице применяется итерация обратной мощности в трех вариантах: без сдвига, постоянного сдвига и обобщенного рэлеевского сдвига на трех этапах алгоритма. Более эффективно выполнять операции линейной алгебры в полиномиальной арифметике, а не с помощью матричных операций, однако свойства обратной степенной итерации остаются прежними.
Реальные коэффициенты
Алгоритм Дженкинса – Трауба, описанный ранее, работает для многочленов с комплексными коэффициентами. Те же авторы создали трехэтапный алгоритм для многочленов с действительными коэффициентами. См. Дженкинс и Трауб Трехэтапный алгоритм для вещественных многочленов с использованием квадратичной итерации . [5] Алгоритм находит линейный или квадратичный фактор, полностью работающий в реальной арифметике. Если комплексный и реальный алгоритмы применяются к одному и тому же действительному многочлену, реальный алгоритм будет примерно в четыре раза быстрее. Реальный алгоритм всегда сходится, и скорость сходимости выше второго порядка.
Связь со сдвинутым QR-алгоритмом
Есть удивительная связь со сдвинутым QR-алгоритмом для вычисления собственных значений матрицы. См. Деккер и Трауб «Сдвинутый QR-алгоритм для эрмитовых матриц» . [6] Сдвиги снова можно рассматривать как итерацию Ньютона-Рафсона на последовательности рациональных функций, сходящихся к полиному первой степени.
Программное обеспечение и тестирование
Программное обеспечение для алгоритма Дженкинса – Трауба было опубликовано как Алгоритм Дженкинса и Трауба 419: нули комплексного многочлена . [7] Программное обеспечение для реального алгоритма было опубликовано как Алгоритм Дженкинса 493: нули действительного многочлена . [8]
Методы были протестированы многими людьми. Как и предполагалось, они пользуются более быстрой, чем квадратичная сходимость для всех распределений нулей.
Однако существуют многочлены, которые могут привести к потере точности, как показано в следующем примере. Все нули многочлена лежат на двух полукругах разного радиуса. Уилкинсон рекомендует, чтобы для стабильного дефляции сначала вычислялись меньшие нули. Сдвиги второй ступени выбираются таким образом, чтобы сначала находились нули на меньшем полукруге. Известно, что после дефляции многочлен с нулями на полукруге плохо обусловлен, если степень велика; см. Wilkinson, [9] p. 64. Исходный полином был 60-й степени и страдал от сильной дефляционной нестабильности.
Рекомендации
- ^ Press, WH, Teukolsky, SA, Vetterling, WT и Flannery, BP (2007), Численные рецепты: Искусство научных вычислений, 3-е изд., Cambridge University Press, стр. 470.
- ^ Дженкинс, М.А. и Трауб, Дж. Ф. (1970), Трехэтапная итерация сдвига переменных для полиномиальных нулей и ее связь с обобщенной итерацией Рэлея , Numer. Математика. 14, 252–263.
- ↑ Ральстон, А. и Рабиновиц, П. (1978), Первый курс численного анализа, 2-е изд., Макгроу-Хилл, Нью-Йорк.
- ^ Трауб, Дж. Ф. (1966), Класс глобально сходящихся итерационных функций для решения полиномиальных уравнений , Math. Comp., 20 (93), 113–138.
- ^ Дженкинс, М.А. и Трауб, Дж. Ф. (1970), Трехэтапный алгоритм для вещественных многочленов с использованием квадратичной итерации , SIAM J. Numer. Анализ., 7 (4), 545–566.
- ^ Деккер, Т. Дж. И Трауб, Дж. Ф. (1971), Сдвинутый QR-алгоритм для эрмитовых матриц , Lin. Алгебра и приложения, 4 (2), 137–154.
- ^ Дженкинс, М.А., Трауб, Дж. Ф. (1972), Алгоритм 419: нули комплексного многочлена , Comm. АКМ, 15, 97–99.
- ^ Дженкинс, Массачусетс (1975), Алгоритм 493: нули действительного многочлена , ACM TOMS, 1, 178–189.
- ^ Уилкинсон, JH (1963), ошибки округления в алгебраических процессах, Prentice Hall, Englewood Cliffs, NJ
Внешние ссылки
- Бесплатное загружаемое приложение для Windows, использующее метод Дженкинса – Трауба для многочленов с действительными и комплексными коэффициентами.
- RPoly ++ Оптимизированная для SSE реализация C ++ алгоритма RPOLY.