А теперь по порядку.
Проблема нахождения корней алгебраических полиномов известна давно, по крайней мере со средневековья. В школе учат решать квадратные уравнения. В википедии можно найти формулы Кардано для решения кубических уравнений и описание метода Феррари для решения в радикалах уравнения четвёртой степени. Там же описан метод Лобачевского для решения алгебраических уравнений произвольной степени. Суть метода Лобачевского вкратце сводится к следующему.
Нетрудно, имея некоторый исходный полином, построить полином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<level;i++)
{//формирование major
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<n;i++)A[n][i]=kf[i]/kf[n];
//расчёт производных 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<rootsCount;i++)rootsArray[i]=B[n][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
Эту ссылку можно найти в Яндексе поиском по закавыченной фразе «семейство алгоритмов расчета», но текст этой статьи в электронном виде, кажется, недоступен. Поэтому приведу здесь цитату из двух предложений этой статьи:
Вещественный корень полинома всегда находится на участке монотонного изменения полинома, т.е. между корнями производной полинома.
Но производная полинома — это тоже полином, однако, меньшей степени и, найдя его вещественные корни, надо искать корни исходного полинома между корнями производной методом деления пополам.
Комментарии (25)
mamkaololosha
15.06.2016 14:49-17А теперь по порядку.
Где применяется? Как? Зачем? В университетах учат учиться? Покажут на заводе когда попадете по распределению? Альё, уже 25 лет никого никуда не распределяют. Я могу реализовать Ваш метод где-то на конкретной задаче? Где? Выложить на гитхаб и написать в резюме, что я его знаю и понимаю? А когда спросят «как вы его примените для оптимизации сети или базы данных» мне ссылку на эту статью кидать?diller61
15.06.2016 21:19+2хм, мне кажется это применяется в широчайшем спектре задач
распространение нагрева в твердом теле, предсказание траектории движения некой цели, любой график скорости или ускорения —
все это прекрасно описывается полиномами различных степеней(первое что пришло в голову из реальной практики), и полиномы же могут использоваться для расчетов
сплайны, кривые Безье(где применяются надеюсь писать не надо?) — все это тоже полиномы
mamkaololosha
16.06.2016 10:29-1> в широчайшем спектре задач
Вы философ-фундаменталист? У вас спектр задач на столько широк, что вы не можете решить ни одну из них? Можно применить этот метод в ракетостроении чтобы ракеты не падали? Как? Вот конкретная задача. Её могут попросить закодировать в гугл-доке на собеседовании в NASA.
Кривые Безье? Геимдев, 2 курс унивеситета, 10 класс физмата. Ничего оригинального, вот серьезно.diller61
16.06.2016 14:17я практик программист, приведенные мной примеры из реальных рабочих проектов которые я успешно выполнил
и да, в алгоритме наведения ракеты можно использовать и используются полиномы
diller61
16.06.2016 14:25>Кривые Безье? Геимдев, 2 курс унивеситета, 10 класс физмата.
еще автомобилестроение
кстати чтоб ракеты не падали первоначально необходимо определить причину паденияdiller61
16.06.2016 14:33+1может причина в том что кто-то из разработчиков забыл как рассчитывать вещественные корни полиномов)
diller61
16.06.2016 14:31+2когда начинаешь программировать что-то отличное от «оптимизации сети или базы данных» внезапно, иногда приходится вспоминать многое из того чему учили в «универе» и в «10 классе физмата». Физика, геометрия, матан, сопромат, техмех, теория вероятности, социология…
может потребоваться вспомнить все и раскопать новое, что не входило в учебную программу
veontomo
15.06.2016 14:50+6Хотел бы прояснить этот момент:
> Нетрудно, имея некоторый исходный полином, построить полином2, имеющий те же по модулю корни,
> что и исходный полином, но с противоположным знаком.
> Перемножая исходный полином и полином2,
> получаем полином, корни которого равны квадратам корней исходного полинома.
Перемножая полиномы, получим полином множество корней которого равно объединению множеств корней исходных полиномов. И совсем необязательно элементы из объединения равны квадратам элементов исходных множеств.
Пример: полином p1 = x — 2 с одним корнем равным 2, строим полином2 p2 = x + 2, который имеет один корень -2. Произведение полиномов p1*p2 = x^2 — 4 имеет корни +2 и -2. Который из них равен квадрату корней исходых полиномов?
tsul
15.06.2016 21:10Опущена операция замены переменной: x^2 -> z. Т.е. квадрированным называется полином*полином2 при условии замены. Имеет ту же степень, что и исходный полином.
Нашёл здесь: stu.sernam.ru/book_dig_m.php?id=69
artemoz
15.06.2016 21:21Если я правильно понял, квадрированный полином предполагается решать относительно x^2. То бишь несколько раз квадрируем, потом берем корень соответствующей степени из второго коэффициента.
ezz666
15.06.2016 21:21Нужно ещё сделать замену x^2 на x. Очевидно, что получившийся полином уже будет иметь корни равные квадратам корней исходных полиномов.
Но есть другая проблема: если у полинома есть комплексные корни этот метод уже не сработае, замену сделать не удастся.veontomo
16.06.2016 14:15Предлагаемый метод претендует на нахождение только вещественных корней, а при при квадрировании комплексные корни появиться не могут.
Timka21213
15.06.2016 19:03http://www.sciencedirect.com/science/article/pii/S0747717197901905 < — как вариант
BalinTomsk
15.06.2016 21:06А можно вставить какой-нибудь реальный пример? Типа такого?
int main(int argc, char** argv) { double kf = 5; double rootsArray[] = {1,3,6,9}; int rootsCount = 2; polynomRealRoots(15, &kf, rootsArray, rootsCount); return 0; }
Заодно и красивый график для привлечения внимания?
pchelintsev_an
15.06.2016 22:36Есть неплохая книга по этой теме Г.П. Кутищева «Решение алгебраических уравнений произвольной степени: Теория, методы, алгоритмы» (URSS, 2010). Аналитически для уравнения третьей степени — мой топик.
vadimr
16.06.2016 06:12Не вполне понятна практическая применимость описанного. Давно известны простые численные методы нахождения вещественных корней произвольных дифференцируемых функций. В частности, можно взять в библиотеке старинную книжку «Библиотека алгоритмов. Алгоритмы 1б-50б» и посмотреть там алгоритм 25б. Помнится, я в детстве ещё искал с его помощью корни уравнений вида sin(x)=ln(x) и т.п.
lostmsu
16.06.2016 18:45+1> Это преобразование (квадрирование) полезно повторить несколько раз. В результате, если корни исходного полинома не были равны друг другу, их многократно квадрированные значения оказываются далеко разнесёнными по величине
Это как они окажутся далеко разнесёнными по величине, если они лежат на [-1; 1]?chersanya
17.06.2016 01:32Далеко — не обязательно в абсолютных значениях. 0.001 отличается от 1 так же, как 1 от 1000.
rafuck
17.06.2016 02:00А можно, наверное, глупый вопрос. Зачем на практике находить корни полиномов высокой степени?
Calvrack
17.06.2016 10:33ну масса есть задач. Из бытовой практики — 5ти точечный метод калибровки стерео-пары. Или — есть образмеренный чертеж в Solidworks с наложенными ограничениями, и надо его отобразить — степень там может быть более-менее большой. Вот любой кому нужна такая функциональность для его собственных нужд вынужден искать корни полиномов (или систем полиномов).
Вообще поиск собственных чисел и векторов можно сводить к полиномам, хотя на практике бывает что и полиномы раскладывают записывая матрицу с характеристическим многочленом совпадающих с данным полиномом.rafuck
17.06.2016 17:15Спасибо. Про такой вариант решения задачи на СЗ знал, но не пришло в голову почему-то…
P.S. Немного не в тему, навеяло полиномами. Вспомнил замечательную фразу, не помню, чью: «Полиномы Чебышева всюду плотны в вычислительной математике».
SteveNers
Теоремой Штурма + бин. поиском вроде побыстрее будет
Sirion
Написать собственный велосипед — бесценно.