Численные методы нахождения корней многочленов с действительными и комплексными коэффициентами

Скляр Александр Яковлевич кандидат технических наук доцент; кафедра прикладная математика; Российский технологический университет (МИРЭА) 119602, Россия, г. Москва, Вернадского, 78 Sklyar Alexander Yakovlevich PhD in Technical Science Associate Professor; Department of Applied Mathematics; Russian Technological University (MIREA) 119602, Russia, Moscow, Vernadsky , 78 askliar@mail.ru Другие публикации этого автора

Введение

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

Методы решения уравнений до второй степени известны еще со времен древней Греции. В XVI веке получены аналитические выражения для многочленов 3 степени (формула Кордано) и 4 степени, полученные Сципионом дель Ферро, Тарталья, Феррари ^[3].

В как 20-х гг. Абель, а затем Галуа в 30-х гг. XIX в. доказали, что такие формулы для уравнений n-й степени в общем случае при любом n ≥ 5 заведомо не могут быть найдены.

Для численного приближенного решения уравнений высших степеней в настоящее время используются различные методы, такие как метод Лобачевского, метод Хичкока, схема Горнера для деления многочлена на двучлен и квадратный трехчлен ^{[1, 4]} и другие.

Другой тип алгоритмов, основанный на получении итерационных формул, за счет выделения из многочленов простых и квадратичных множителей, с последующим сопоставлением записей многочленов с остатком, когда корни являются приближенными, и без остатка, когда значения корней являются точными рассмотрен в статьях Чье Ен Ун и А.Б. Шеина ^[5,6,7].

Ряд алгоритмов решения полиномиальных задач и обширная библиография приведены у Г. П. Кутищева ^[8]. Стоит отметить также подходы для их решения в ^{[9, 10, 11]}.

Заметим, что само существование большого числа разнообразных методов в целом свидетельствует, что не существует ни одного «вполне удовлетворительного».

1. Удаление кратных корней

Многочлен степени n имеет в точности n корней и представим в виде

В общем случае не все корни x_k различны.

Нахождения корней многочлена наталкивается на определенные трудности в случаях, когда он имеет кратные корни.

Пусть имеется кратный корень x* с кратностью m, тогда P_n(x) можно представить в виде

Следовательно, P_n(x) и P'_n(x) будут иметь общий делитель (x-x*)^m-1.

Тогда, используя алгоритм Евклида можно найти наибольший делитель Q(x) многочленов P_n(x) и P'_n(x).

И исходный многочлен P_n(x) быть представлен, как

Многочлен H(x) при этом будет иметь те же корни, что и P_n(x), но не будет иметь кратных корней. Таким образом исходная задача сводится к нахождению корней многочлена, все корни которого различны.

2. Определение диапазона корней

Рассмотрим уравнение (a₀ будем полагать равным 1)

Согласно известной теореме ^[1,2] все корни x_k (k = 1, 2, …, n) многочлена в комплексной плоскости лежат в кольце

r <| x_k |< R, где

В то же время сразу отметим, что это утверждение можно несколько усилить.

Проведем замену переменной x=y/t. В этом случае исходное уравнение относительно y примет вид

Учитывая характер подстановки, можно записать

Меняя значение t, можно добиться уменьшения оценки верхней границы.

Пусть максимальное значение принимает модуль коэффициента при y_m. Тогда оценка по этому коэффициенту будет не менее, чем

t* находится из

Таким образом можно утверждать, что оценка лежит между значениями, соответствующими текущему значению t и t*.

В то же время при изменении t меняются и другие коэффициенты. Таким образом получаем ограничения на пересечение линий их изменения.

Если t_r не лежит внутри интервала t и t*, то его можно игнорировать. Из остальных выбираем значение ближайшее к t.

Если таковых нет, то оптимальное t_r=t*.

Если знаки производных

совпадают, то дальнейшее улучшение невозможно. Оптимальное t_r=t*.

Если знаки производных совпадают, то принимаем m равным k и повторяем процедуру.

Рассмотрим производную по k коэффициенту в точке t

Рассмотрим в качестве исходного многочлен четвертой степени

t⁴-10t³+35t²-50^t+24=(t-1)(t-2)(t-3)(t-4)

Исходная верхняя граница – 51, m=1, t=1,

Точки пересечения

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

3. Поиск действительных корней

Рассмотрим отдельно поиск положительных и отрицательных корней.

Пусть исходный многочлен имеет вид (старший коэффициент a₀ не снижая общности можно полагать равным 1).

Введем b_k=(a_k+|a_k|)/2 и c_k=(a_k-|a_k|)/2. Для a_k>0 b_k=a_k, c_k=0 для a_k≤0 b_k=0. c_k=-a_k. В этих обозначениях

Многочлены P₊(x) и P_-(x) при x>0 представляют собой непрерывные неотрицательные монотонно возрастающие функции.

Используя введенные многочлены P₊(x) и P_-(x) и верхнюю границу значений корней R будем искать корни в диапазоне [A, B], где A=0, B=R. Нижнюю границу можно уточнить до величины r*, получив диапазон [r*, R*] вместо диапазона [r, R], но это не имеет принципиального значения.

Рассмотрим подробнее алгоритм нахождения корней.

Вычисляем значения P₊(x) и P_-(x) на концах интервала и вызываем функцию Root вычисления корней с параметрами A, B, P₊(A), P_-(A) и P₊(B), P_-(B). Функция Root возвращает либо найденный корень (значение большее 0), либо -1, если на данном интервале многочлен P(x) не имеет корней.

Алгоритм функции Root имеет следующий вид.

1. Вычисляем значения исходного многочлена на концах интервала P(A)=P₊(A)-P_-(A) и P(B)= P₊(B)-P-(B)

2. Если P(A)=0 возвращаем найденный корень A, если P(B)=0 возвращаем найденный корень B, если P(A)P(B)вызываем стандартную процедуру STROOT нахождения корней на заданном интервале [A, B] (например, методом дихотомии или хорд) и возвращаем, найденное ей значение.

3. Вычисляем Q₁=P₊(A)-P_-(B) – нижняя граница значений P(x) на интервале [A, B] и Q₂=P₊(B)-P_-(A) – верхняя граница значений P(x) на интервале [A, B].

4. Если Q₁Q₂≥0, то на данном интервале корней нет и возвращаем значение -1.

5. Вычисляем C=(A+B)/2.

Вызываем функцию Root вычисления корней с параметрами A, С, P₊(A), P_-(A) и P₊(С), P_-(С).

6. Если вычисленное значение x>0, то возвращаем найденный корень x.

7. Вызываем функцию Root вычисления корней с параметрами С, B, P₊(C), P_-(C) и P₊(B), P_-(B). Возвращаем найденный корень x (если корней нет, то -1).

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

По нахождении корня u переходим к поиску следующего корня, заменив исходный многочлен P_n(x) многочленом меньшей степени P_n-1(x), P_n(x)=(x-u)P_n-1(x). Многочлен P_n-1(x) имеет все оставшиеся корни P_n(x). Повторяя приведенный выше алгоритм, находим все положительные корни.

Для нахождения отрицательных корней заменим исходный многочлен P_n(x) многочленом (-1)ⁿP_n(-x), корни которого имеют те же значения, но с обратным знаком. Найдя все положительные корни второго многочлена, мы тем самым найдем все действительные корни исходного многочлена.

Поскольку других действительных корней нет, то оставшиеся корни представляют пары комплексно-сопряженных чисел и оставшийся многочлен имеет четную степень.

4. Поиск комплексных корней

Пусть многочлен P_2n(x) с действительными коэффициентами не имеет действительных корней, тогда он представим в виде

Рассмотрим представление многочлена степени 2n в виде

Для стандартизации расчетов для всех m введем q_-2=q_-1=q_n-1=q_n≡0, q_n-2=1.

Нахождение корней будем проводить на основе оптимизации аппроксимации исходного многочлена Pn(x) многочленом Rn(x).

Введем функцию

При заданных значениях p₀, p₁ коэффициенты q_i находятся из требования минимизации F(p,q).

Или

При заданных значениях q_i коэффициенты p₀, p₁ находятся из требования минимизации F(p,q).

Или

Для нахождения многочлена x²+p₁x+p₀ можно воспользоваться следующим итеративным алгоритмом.

1. Задаем начальное значение коэффициентов p₀, p₁.

2. На основе (4.5) рассчитываем коэффициенты q_i. Задача сводится к решению системы линейных уравнений (СЛАУ) для 5-диагональной матрицы. Учитывая вид матрицы, задача сводится к ряду подготовительных операций трудности O(n²) и решению СЛАУ размерности 3∙3.

3. Полученный набор коэффициентов q_i используем для получения в соответствии с (4.6) скорректированных значений p₀, p₁.

4. Оцениваем величину погрешности вычисления в соответствии с (4.3). Если погрешность меньше заданного порога, то заканчиваем работу алгоритма. В противном случае переходим к шагу 2 алгоритма.

В результате получаем многочлена x₂+p₁x+p₀ дающий первую пару корней. Коэффициенты qi дают многочлен

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

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

Пусть на i-ом и i+1-ом шагах получили значения p_i,1, p_i+1,1 и p_i,0, p_i+1,0 с невязками F(p,q) – Fa, Fb соответственно.

Вычислим значения p_c1=p_i,1+2(p_i+1,1-p_i+1,1) и p_c0= p_i+1,0+2(p_i+1,0-p_i,0) и соответствующее ему значение Fc.

Если F_c<F_b проводим замену точек и значений в них a=b, b=c и повторяем расчет. Если F_c<F_b, то квадратичной интерполяцией находим точку c*, принимая ее за субоптимальное значение и возвращаемся в основной алгоритм.

В точке оптимума должны одновременно выполняться (4.5 и 4.6)

В матричном виде совокупность (4.5, 4.6) представима в виде

При |A|≠0 алгоритм сходится к решению задачи.

Кроме того, с точки зрения практической реализации значительно эффективнее выглядит гибридный алгоритм.

Несколько итераций базового алгоритма приводят нас в окрестность корня. Учитывая, что при отсутствии кратных корней |P'_n(x)| гарантированно не обращается в 0, при P_n(x)=0 схема с использованием значений производных (метод Ньютона), либо более тонких его модификаций ^[12]) обеспечивает достаточно быструю сходимость.

Пусть получено приближение s=x²+p₁x+p₀. Его корень будет

Значение P_n(z) в окрестности z₀ представимо в виде

Это позволяет рассчитать приближение корня z*.

При z₀, близких к z* алгоритм сходится достаточно быстро. В результате получаем для s=x²+p₁x+p₀

Полученное выражение используем для получения многочлена Q_n-2(x) из P_n(x)=(x²+p₁x+p₀)Q_n-2(x), после чего ищем следующие корни, многочлена, совпадающие с корнями Q_n-2(x).

5. Поиск корней многочленов с комплексными коэффициентами

Рассмотрим уравнение (a₀ будем полагать равным 1). Коэффициенты a_k – комплексные числа.

Пусть его корни – x₁, x₂, …, x_n, тогда корни уравнения

будут

Построим многочлен R_2n(x)=P_n(x)Q_n(x). Этот многочлен имеет корни

Величины b_k представимы как суммы и произведения комплексно-сопряженных чисел и, следовательно являются действительными числами. Таким образом многочлен R_2n(x) является многочленом с действительными коэффициентами.

Соответственно поиск корней многочлена P_n(x) с комплексными коэффициентами сводится к поиску корней многочлена R_2n(x) с действительными коэффициентами. Отметим, что все действительные его корни будут кратными, а из полученных комплексных сопряженных корней часть будет посторонней.

Заключение

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

Для нахождения корней многочлена степени n с действительными коэффициентами строится алгоритм, включающий следующие основные этапы:

· определение кратных корней,

· выделение диапазона корней;

· нахождение интервалов, гарантированно содержащих корни (за время O(ln(R/H)), где R – оценка общего диапазона корней, а H длина интервала, гарантированно содержащего корень);

· итеративное построение трехчленов, служащих оценкой значений пар комплексно - сопряженных корней (количество итераций не обязательно должно соответствовать требованиям точности решения, достаточно выполнения требования изоляции корней, после чего решение достигается традиционными методами, например методом Ньютона).

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

Программная реализация предлагаемых алгоритмов выполнена на языке C++ в среде операционной системы Windows. Учитывая вычислительный характер алгоритмов перенос их в иную операционную среду не вызывает каких-либо проблем.