[Из песочницы] Алгоритм расчёта вещественных корней полиномов
Основополагающая идея этого алгоритма очень проста и может быть изложена двумя предложениями. Вещественный корень полинома всегда находится на участке монотонного изменения полинома, т.е. между корнями производной полинома. Но производная полинома — это тоже полином, однако, меньшей степени и, найдя его вещественные корни, надо искать корни исходного полинома между корнями производной методом деления пополам.
А теперь по порядку.
Проблема нахождения корней алгебраических полиномов известна давно, по крайней мере со средневековья. В школе учат решать квадратные уравнения. В википедии можно найти формулы Кардано для решения кубических уравнений и описание метода Феррари для решения в радикалах уравнения четвёртой степени. Там же описан метод Лобачевского для решения алгебраических уравнений произвольной степени. Суть метода Лобачевского вкратце сводится к следующему.
Нетрудно, имея некоторый исходный полином, построить полином2, имеющий те же по модулю корни, что и исходный полином, но с противоположным знаком. Перемножая исходный полином и полином2, получаем полином, корни которого равны квадратам корней исходного полинома.
Это преобразование (квадрирование) полезно повторить несколько раз. В результате, если корни исходного полинома не были равны друг другу, их многократно квадрированные значения оказываются далеко разнесёнными по величине, а их приближенные значения очень просто выражаются через коэффициенты соответствующего квадрированного полинома.
В частности, если коэффициент при старшей степени аргумента полинома равен единице, то следующий по старшинству коэффициент равен (с обратным знаком) сумме корней уравнения, а поскольку значения этих корней сильно разнесены, то приближенно можно считать эту сумму равной наибольшему по модулю корню.
Для конкретности сообщим, что для полинома 4-й степени с корнями 1, 2, 3, 4 метод Лобачевского уже после четвёртого квадрирования даёт правильные до второго знака после запятой значения корней. При этом для представления коэффициентов полиномов достаточно формата long double.
Бесспорно, этот метод является ценным инструментом в руках исследователя, наделённого интеллектом. Однако, его программирование для современной вычислительной техники вызывает серьёзные затруднения при необходимости строгой гарантии достоверности результата при всевозможных особых случаях расположения корней.
Теперь я начну описывать иной метод. В общедоступной печати упоминание о нём начинается с работы [1]. Какие-либо независимые публикации о применении такого метода мне неизвестны. Этот алгоритм сводится к последовательному исследованию интервалов монотонного изменения исходного полинома. Если на границах этого интервала монотонности значения полинома имеют разные знаки, то запускается процедура деления отрезка пополам для расчёта точного значения очередного корня. Границами интервалов монотонности являются точки, в которых значение производной полинома обращается в нуль, т.е. это корни производного полинома. Производный полином имеет степень на единицу меньшую, чем исходный полином, и процесс расчёта коэффициентов производных полиномов следует продолжить до полинома первой степени, корень которого находится непосредственно, без привлечения процедуры деления отрезка пополам. В результате этого шага получим два интервала монотонного изменения для производного полинома второй степени.
Теперь можно найти два вещественных корня производного полинома второй степени (если они существуют) и далее по лестнице из производных полиномов подниматься до корней исходного полинома. Остаётся пояснить, как технически реализуются границы «плюс, минус бесконечность» интервалов монотонности исходного и производных полиномов.
Нормируем полином так, чтобы коэффициент при старшей степени аргумента стал равным единице. Пусть M — наибольшее по модулю значение среди его остальных коэффициентов. Тогда значение полинома больше единицы для всех значений аргумента, больших, чем M+1.
Для доказательства рассмотрим расчёт полинома p (x)=x^n+k[n-1]*x^(n-1)+…+k[1]*x+k[0] по схеме Горнера.
На первом шаге вычисляется p[1]=k[n-1]+x и очевидно, что p[1]>1.
На втором шаге вычисляется p[2]=k[n-2]+x*p[1] и вновь очевидно, что p[2]>1.
Аналогичное имеет место на последующих шагах.
На последнем шаге вычисляется p (x)=k[0]+x*p[n-1] и окончательно получим p (x)>1.
Таким образом, если нужно определить знак полинома при бесконечном значении аргумента, следует взять аргумент равным M+1.
Прилагаемый текст соответствующей программы вполне заменяет занудное изложение отдельных технических подробностей описанного тут алгоритма.
Прокомментирую, наконец, не вполне очевидную особенность реализации алгоритма деления отрезка пополам.
Пробная точка pt, расположенная посередине между текущими концами ng и vg отрезка, вычисляется оператором pt=0.5*(ng+vg); , а цикл делений пополам прерывается оператором if (pt<=ng||pt>=vg)break; .
В силу конечной точности представления вещественных чисел в машине рано или поздно наступает состояние, при котором операция деления пополам вместо нового числа даёт значение одной из исходных границ. Именно в этом состоянии следует прекратить цикл делений пополам. Это состояние соответствует максимально достижимой точности результата.
Недавно мне удалось использовать этот алгоритм для решения задачи вычисления комплексного корня полинома, не имеющего вещественных корней. Но об этом я планирую рассказать на Хабре в следующей статье.
Ниже, как приложение, приведен полный текст файла polynomRealRoots.cpp, реализующего описанныйалгоритм.
//*************************************************************************
double polinom(int n,double x,double *k)
//полином вида x^n+k[n-1]*x^(n-1)+k[n-2]*x^(n-2)+...+k[1]*x+k[0]
//со старшим коэффициентом, равным единице
{
double s=1;
for(int i=n-1;i>=0;i--)
s=s*x+k[i];
return s;
}//polinom//расчёт корня полинома методом деления пополам отрезка, содержащего этот корень
double dihot(int degree,double edgeNegativ,double edgePositiv,double *kf)
{
for(;;)
{//цикл деления отрезка пополам
double x=0.5*(edgeNegativ+edgePositiv);
if(x==edgeNegativ||x==edgePositiv)return x;
if(polinom(degree,x,kf)<0)edgeNegativ=x;
else edgePositiv=x;
}//цикл деления отрезка пополам
}//dihot
//один шаг подъёма по лестнице из производных полиномов
void stepUp(
int level, //индекс достигаемой ступеньки лестницы
double **A, //вспомогательные массивы
double **B, //параметров производных полиномов
int *currentRootsCount //сформированные в вызывающей процедуре
)
{
//аргумент полинома, равносильный его бесконечно большому значению
double major=0;
for(int i=0;i
double s=fabs(A[level][i]);
if(s>major)major=s;
}//формирование major
major+=1.0;
currentRootsCount[level]=0; //рабочая инициализация
//основной цикл поиска корня уравнения
//A[level][0]+A[level][1]*x+...+A[level][level-1]*x^(level-1)+x^level=0
//на очередном интервале монотонности
for(int i=0;i<=currentRootsCount[level-1];i++)
{//очередной интервал монотонности
//signLeft signRight - знаки текущего A-полинома на левой и правой границе интервала монотонности
int signLeft,signRight;
//предварительная левая и правая границы интервала поиска
double edgeLeft,edgeRight;
//границы интервала монотонности, несущие информацию о знаке полинома на них
double edgeNegativ,edgePositiv;
//формирование левой границы поиска
if(i==0)edgeLeft=-major;
else edgeLeft=B[level-1][i-1];
//значение текущего A-полинома на левой границе
double rb=polinom(level,edgeLeft,A[level]);
if(rb==0)
{//маловероятный случай попадания в корень
B[level][currentRootsCount[level]]=edgeLeft;
currentRootsCount[level]++;
continue;
}//маловероятный случай попадания в корень
//запомнить знак текущего A-полинома на левой границе
if(rb>0)signLeft=1;else signLeft=-1;
//формирование правой границы поиска
if(i==currentRootsCount[level-1])edgeRight=major;
else edgeRight=B[level-1][i];
//значение текущего A-полинома на правой границе
rb=polinom(level,edgeRight,A[level]);
if(rb==0)
{//маловероятный случай попадания в корень
B[level][currentRootsCount[level]]=edgeRight;
currentRootsCount[level]++;
continue;
}//маловероятный случай попадания в корень
//запомнить знак текущего A-полинома на правой границе
if(rb>0)signRight=1;else signRight=-1;
//если знаки полинома на границах интервала монотонности совпадают,
//то корня нет
if(signLeft==signRight)continue;
//теперь можно определить плюс границу и минус границу поиска корня
if(signLeft<0){edgeNegativ=edgeLeft;edgePositiv=edgeRight;}
else {edgeNegativ=edgeRight;edgePositiv=edgeLeft;}
//всё готово для локализации корня методом деления пополам интервала поиска
B[level][currentRootsCount[level]]=dihot(level,edgeNegativ,edgePositiv,A[level]);
currentRootsCount[level]++;
}//очередной интервал монотонности
return;
}//stepUp
//процедура находит все вещественные корни полинома любой степени
void polynomRealRoots(
int n, //степень исходного полинома
double *kf, //массив коэффициентов исходного полинома
double *rootsArray, //выходной массив вычисленных корней
int &rootsCount //количество найденных корней
)
{
//используются вспомогательные массивы A и B, имеющие следующее содержание
//A это коэффициенты а B корни производных полиномов
//все A-полиномы нормируются так,
//чтобы коэффициент при старшей степени был равен единице
//A[n] - это массив нормированных коэффициентов исходного полинома
//B[n] - это массив корней исходного полинома
//A[n-1] - это массив нормированных коэффициентов производного полинома
//B[n-1] - это массив корней производного полинома
//аналогичным образом
//A[n-2] и B[n-2] - это коэффициенты и корни дважды производного полинома
//наконец A[1] - это массив коэффициентов последнего полинома
//в цепочке производных полиномов
//это линейный полином и B[1] - это массив его корней,
//представленный единственным значимым элементом
//выделение памяти для вспомогательных массивов
double **A=new double *[n+1];
double **B=new double *[n+1];
//количество вещественных корней для каждого из A-полиномов
int *currentRootsCount=new int[n+1];
for(int i=1;i<=n;i++)
{
A[i]=new double[i];
B[i]=new double[i];
}
//нормировка исходного полинома
for(int i=0;i
//расчёт производных A-полиномов
for(int i1=n,i=n-1;i>0;i1=i,i--)
{
for(int j1=i,j=i-1;j>=0;j1=j,j--)
{
A[i][j]=A[i1][j1]*j1/i1;
}
}
//формирование исходного корня последнего производного полинома
currentRootsCount[1]=1;
B[1][0]=-A[1][0];
//подъём по лестнице производных полиномов
for(int i=2;i<=n;i++)stepUp(i,A,B,currentRootsCount);
//формирование результата
rootsCount=currentRootsCount[n];
for(int i=0;i
//очистка памяти
for(int i=1;i<=n;i++)
{
delete[]A[i];
delete[]B[i];
}
delete[]A;
delete[]B;
delete[]currentRootsCount;
return;
}//polynomRealRoots
//*************************************************************************
Примите также текст заголовочного файла polynomRealRoots.h, позволяющего легко организовать ссылку на приведенный выше программный модуль.
//*************************************************************************
//процедура вычисления вещественных корней полинома
void polynomRealRoots(int n,double *kf,double *rootsArray,int &rootsCount);
//*************************************************************************
Литература
1. Костин И.К. Семейство алгоритмов расчета интервалов прохождения космического аппарата над круговым наземным объектом с учетом продольной ошибки определения параметров орбиты
Вопросы радиоэлектроники, сер. РЛТ, 2004 г., вып. 1
Эту ссылку можно найти в Яндексе поиском по закавыченной фразе «семейство алгоритмов расчета», но текст этой статьи в электронном виде, кажется, недоступен. Поэтому приведу здесь цитату из двух предложений этой статьи:
Вещественный корень полинома всегда находится на участке монотонного изменения полинома, т.е. между корнями производной полинома.
Но производная полинома — это тоже полином, однако, меньшей степени и, найдя его вещественные корни, надо искать корни исходного полинома между корнями производной методом деления пополам.
