Алгоритм розрахунку дійсних коренів поліномів

Основна ідея цього алгоритму дуже проста і може бути викладена двома пропозиціями. Дійсний корінь полінома завжди знаходиться на ділянці монотонного зміни полінома, тобто між корінням похідної полінома. Але похідної полінома — це теж поліном, однак, меншою мірою і, знайшовши його речові коріння, треба шукати коріння вихідного полінома між корінням похідної методом ділення навпіл.

А тепер по порядку.

Проблема знаходження коренів алгебраїчних поліномів відома давно, принаймні з середньовіччя. В школі вчать вирішувати квадратні рівняння. У вікіпедії можна знайти формули Кардано для вирішення кубічних рівнянь і опис методу Феррарі для розв'язування радикалів рівняння четвертого ступеня. Там описано метод Лобачевського для вирішення алгебраїчних рівнянь довільного ступеня. Суть методу Лобачевского коротко зводиться до наступного.

Неважко, маючи певний вихідний поліном, побудувати полином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, що реалізує описанныйалгоритм.

файл 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 j 1=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

Посилання можна знайти в Яндексі пошуком по закавыченной фразі «сімейство алгоритмів розрахунку», але текст цієї статті в електронному вигляді, здається, недоступний. Тому наведу тут цитату з двох речень цієї статті:
Дійсний корінь полінома завжди знаходиться на ділянці монотонного зміни полінома, тобто між корінням похідної полінома.

Але похідної полінома — це теж поліном, однак, меншою мірою і, знайшовши його речові коріння, треба шукати коріння вихідного полінома між корінням похідної методом ділення навпіл.
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.