Алгоритм обчислення комплексного кореня полінома довільного ступеня

Це завершення статті habrahabr.ru/post/303342

Дякую коментаторам, які зробили більш ясним моє надто вже конспективное виклад методу Лобачевского.

Насправді, мені слід було явно написати, що квадрированный поліном треба розглядати як поліном від аргументу x^2, де x — аргумент вихідного полінома.

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

Тепер на цьому фундаменті буде побудований цілком елементарний алгоритм обчислення комплексного кореня многочлена, не має дійсних коренів.

Спочатку нормуємо поліном так, щоб його вільний член став рівним одиниці.

Тоді при значенні аргументу x=0 значення полінома буде речовим і рівним одиниці.
Це широко відома математикам відправна точка міркувань.

Уявімо аргумент x полінома у вигляді
x=r*exp(i*fi)


Тут i це уявна одиниця.

При зміні fi від нуля до 2*pi (=6.28318...) значення полінома p опише деяку замкнену криву.
Назвемо цю криву эпициклоидой.

При дуже малому r эпициклоида буде схожа на маленьке кільце навколо точки p=1.
Не будемо відволікатися на ситуації, коли лінійний член многочлена дорівнює нулю.

При нескінченно великому значенні r эпициклоида опише кілька петель навколо точки p=1.
Кількість петель одно ступеня полінома, форма близька до кола великого радіусу.

Було б корисно навчитися обчислювати всі значення fi, при яких уявна частина полінома звертається в нуль.
Це дозволить кожному значенню r конструктивно поставити у відповідність мінімальне значення координати x, при якій эпициклоида перетинає вісь абсцис.

Це значення позитивно при малих r і негативно при великих r.

Методом ділення відрізка навпіл можна знайти таке r, при якому це значення буде в точності дорівнює нулю.
Відповідні значення r і fi визначать шуканий комплексний корінь полінома.

Перейдемо до реалізації наміченого плану.

Уявна частина многочлена з цілими коефіцієнтами kf0, kf1, kf2,..., kfn представиться у вигляді:

Im=kf1*r*sin(fi)+kf2*r^2*sin(2*fi)+...+kfn*r^n*sin(n*fi)


Хотілося б представити її у вигляді якого-небудь полінома, оскільки речові корені поліномів ми вміємо рахувати.
І справді, цю уявну частину можна представити як добуток синуса fi, на деякий поліном від косинуса fi.
При невеликому ступені поліному це можна безпосередньо перевірити.

Використовуємо метод математичної індукції для доведення цього твердження в загальному вигляді.

Нехай


cos(n*x)=polynomC(cos(x))
sin(n*x)=sin(x)*polynomS(cos(x))


Тоді


cos(n*x+x)=cos(n*x)*cos(x)-sin(n*x)*sin(x)=
=polynomC(cos(x))*cos(x)-sin^2(x)*polynomS(cos(x))


Оскільки
sin^2(x)=1-cos^2(x)
, ясно, що останній вираз є деяким поліномом від cos(x).

Аналогічно

sin(n*x+x)=sin(n*x)*cos(x)+cos(n*x)*sin(x)=
=sin(x)*(polynomS(cos(x)*cos(x)+polynomC(cos(x))))


І тут очевидно, що співмножник при sin(x) в останньому виразі є поліномом від cos(x).

Зрозуміло, що ретельне виписування наведених співвідношень дозволяє отримати рекурентні формули для коефіцієнтів потрібних нам поліномів.

Подальші подробиці реалізації наміченого плану зручніше простежити за наведеними нижче анотований текстів відповідної програми.

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

Текст файлу polynomComplexRoot.hВміст


//*************************************************************************
class polinomComplexRoot
{
public:
//тіло класу визначається ступенем досліджуваного полінома і набором його коефіцієнтів
polinomComplexRoot(int _degree,double* _kf);
~polinomComplexRoot();

//основна процедура зовнішнього інтерфейсу класу
//повертає код помилки
//при успішному завершенні код помилки дорівнює нулю,
//а значення параметрів rootRe і rootIm стануть дорівнюють відповідно
//дійсної та уявної частини комплексного кореня досліджуваного полінома
int run(double& rootRe,double& rootIm);

private:
int degree,errorCode;
double *kf;
double *ki;
double *cosRootsArray;
double *rpw;
int **ck;
int **sk;

//основна робоча процедура
//при заданому r знаходить мінімальне значення речової частини
//досліджуваного полінома від аргументу r*exp(i*fi)
//при тому значенні fi, яке забезпечує рівність нулю
//його уявної частини. Це значення fi 
//буде присвоєно другого аргументу процедури по закінченні виклику
double minAxeXcrossVal(double r,double& fi);
};//class polinomComplexRoot
//*************************************************************************


Текст файлу polynomComplexRoot.cppВміст


//*************************************************************************
#include <math.h>
#include <polynomComplexRoot.h>
#include <polynomRealRoots.h>

const double cPI=3.14159265358979323846;

polinomComplexRoot::polinomComplexRoot(int _degree,double* _kf)
{
degree=_degree;errorCode=0;

//попередня ініціалізація
kf=0;ck=0;sk=0;rpw=0;ki=0;cosRootsArray=0;

//відбраковування неприйнятних значень коефіцієнтів вихідного полінома
if(_kf[0]==0){errorCode=1;return;}
if(_kf[degree]==0){errorCode=2;return;}
if(degree<2){errorCode=3;return;}

kf=new double[1+degree];
ki=new double[1+degree];
cosRootsArray=new double[1+degree];
rpw=new double[degree+1];
ck=new int*[1+degree];
sk=new int*[1+degree];
for(int i=0;i<=degree;i++)
{
ck[i]=new int[1+degree];
sk[i]=new int[1+degree];
}
//відбраковування вихідного полінома, що має речові коріння
polynomRealRoots(degree,_kf,kf,errorCode);
if(errorCode>0){errorCode=4;return;}

for(int i=0;i<=degree;i++)kf[i]=_kf[i]/_kf[0];


for(int i=0;i<=degree;i++)
for(int j=0;j<=degree;j++)
{ck[i][j]=0;sk[i][j]=0;}

//коефіцієнти косинусного поліномів для подання
//синусів і косинусів кратних за формулами
//cos(n*x)=ck[n,0]+ck[n,1]*cos(x)+ck[n,2]*cos^2(x)+...+ck[n,n]*cos^n(x)
//sin(n*x)=sin(x)*(sk[n,0]+sk[n,1]*cos(x)+sk[n,2]*cos^2(x)+...+sk[n,n]*cos^n(x))
ck[1][1]=1;
sk[1][0]=1;

//розрахунок ck і sk по рекуррентным формулами
//тут np1=n+1
for(int n=1,np1=2;n<degree;n=np1,np1++)
for(int k=0;k<=np1;k++)
{//реалізація рекурентних формул
//ck[n+1,k]=ck[n,k-1]+sk[n,k-2]-sk[n,k];
//sk[n+1,k]=sk[n,k-1]+ck[n,k];
if(k>=1)ck[np1][k]+=ck[n][k-1];
if(k>=2)ck[np1][k]+=sk[n][k-2];
if(k>=1)sk[np1][k]+=sk[n][k-1];
if(k < =n){ck[np1][k]-=sk[n][k];sk[np1][k]+=ck[n][k];}
}//реалізація рекурентних формул
return;
}//constructor

polinomComplexRoot::~polinomComplexRoot()
{
if(kf)delete[] kf;
if(ki)delete[] ki;
if(cosRootsArray)delete[] cosRootsArray;
if(rpw)delete[] rpw;
if(ck)
{
for(int i=0;i<=degree;i++)delete[] ck[i];
delete[] ck;
}

if(sk)
{
for(int i=0;i<=degree;i++)delete[] sk[i];
delete[] sk;
}

}//destructor

int polinomComplexRoot::run(double& rootRe,double& rootIm)
{
rootRe=0;rootIm=0;
if(errorCode>0)return errorCode;

//верхня rup і нижня rdn кордону пошуку для r 
double fidn=0,rdn=0,fiup=0,rup=1;

//подвоюємо rup поки воно недостатньо велика
for(;minAxeXcrossVal(rup,fiup)>0;)rup+=rup;

for(;;)
{//цикл ділення інтервалу навпіл (rdn, rup)
double fit,rt=0.5*(rdn+rup);
if(rt==rdn||rt==rup)break;
if(minAxeXcrossVal(rt,fit)>0){rdn=rt;fidn=fit;}else {rup=rt;fiup=fit;}
}//цикл ділення інтервалу навпіл (rdn, rup)

//формальний вибір для видачі результату одного з двох
//практично однакових рішень (rdn, fidn) або (rup, fiup)
//буде обчислений модуль досліджуваного полінома в кожній з цих двох точок
//і як результату рішення буде видана та з цих точок,
//якої відповідає обчислювально менше значення цього модуля
double minmod; //мінімальне значення модуля полінома
bool mindn=true; //мінімум досягнуто у точці з суфіксом dn

for(int c=0;c<2;c++)
{//контрольний розрахунок в точках up і dn
double rc,fic;
if(c==0){rc=rdn;fic=fidn;}
else {rc=rup;fic=fiup;}

//rec - речова частина досліджуваного полінома в пробній точці
//imc - уявна частина досліджуваного полінома в пробній точці
double rec=kf[degree]*cos(degree*fic),
imc=kf[degree]*sin(degree*fic);
for(int i=degree-1;i>0;i--)
{//gorner
rec=rc*rec+kf[i]*cos(i*fic);
imc=rc*imc+kf[i]*sin(i*fic);
}//gorner
rec=rc*rec+kf[0];
imc=rc*imc;

//mc - квадрат модуля досліджуваного полінома в пробній точці
double mc=rec*rec+imc*imc;
if(c==0)minmod=mc;
else if(mc<minmod)
{//точка up краще точки dn
minmod=mc;
mindn=false;
}//точка up краще точки dn
}//контрольний розрахунок в точках up і dn

//формування результату
if(mindn)
{//вибір точки dn в якості результату
rootRe=rdn*cos(fidn);
rootIm=rdn*sin(fidn);
}//вибір точки dn в якості результату
else
{//вибір точки up в якості результату
rootRe=rup*cos(fiup);
rootIm=rup*sin(fiup);
}//вибір точки up в якості результату

return errorCode;
}//run

//основна робоча процедура
//при заданому r знаходить мінімальне значення речової частини
//досліджуваного полінома від аргументу r*exp(i*fi)
//при тому значенні fi, яке забезпечує рівність нулю її уявної частини
//значення fi буде присвоєно другого аргументу процедури по закінченні виклику
double polinomComplexRoot::minAxeXcrossVal(double r,double& fi)
{
//попередньо формуємо допоміжний масив rpw[k]=r^k
double rb=1;
for(int i=0;i<=degree;i++){rpw[i]=rb;rb*=r;}

//значення досліджуваного полінома від дійсного аргументу r
rb=kf[0];for(int i=1;i < =degree;i++)rb+=rpw[i]*kf[i];
double rez=rb;fi=0;

//значення досліджуваного полінома від дійсного аргументу -r
rb=kf[0];
for(int i=1,od=1;i < =degree;i++,od=1-od)if(od)rb=rpw[i]*kf[i];else rb+=rpw[i]*kf[i];

//ми шукаємо мінімальне значення абсциси эпициклоиды
if(rb<rez){rez=rb;fi=cPI;}

//уявна частина досліджуваного полінома від комплексного аргументу,
//представленого параметрами r і fi, виражається так
//im=r*kf[1]*sin(fi)+r^2*kf[2]*sin(2*fi)+...+r^n*kf[n]*sin(n*fi)
//вона буде представлена з участю деякого полінома від cos(fi)
//коефіцієнти цього полінома позначені ідентифікатором ki
//im=sin(fi)*(ki[0]+ki[1]*cos(fi)+ki[2]*cos^2(fi)+...+ki[n]*cos^n(fi))
//коефіцієнти ki виражаються через коефіцієнти досліджуваного полінома kf
//та коефіцієнти ck і sk, обчислені в конструкторі класу за формулами
//ki[0]=r*kf[1]*sk[1][0]+r^2*kf[2]*sk[2][0]+...+r^n*kf[n]*sk[n][0]
//ki[1]=r*kf[1]*sk[1][1]+r^2*kf[2]*sk[2][1]+...+r^n*kf[n]*sk[n][1]
//...
//ki[n]=r*kf[1]*sk[1][n]+r^2*kf[2]*sk[2][n]+...+r^n*kf[n]*sk[n][n]
for(int i=0;i<=degree;i++)
{//обчислення коефіцієнтів ki
rb=0;
for(int j=i+1;j < =degree;j++)rb+=(rpw[j]*kf[j]*sk[j][i]);
ki[i]=rb;
}//обчислення коефіцієнтів ki

//cosDegree це ступінь косинусного полінома, представленого коефіцієнтами ki
//страхуємося від можливості рівності нулю старших коефіцієнтів
int cosDegree=0,cosRootsCount=0;
for(int i=degree-1;i>0;i--)if(fabs(ki[i])>0){cosDegree=i;break;}

//інтерпретуємо ситуацію виродження ki-полінома як внутрішню помилку
if(cosDegree<1){errorCode=5;return rez;}

//знаходимо всі речові коріння ki-полінома
polynomRealRoots(cosDegree,ki,cosRootsArray,cosRootsCount);

//обстеження знайдених коренів ki-полінома
for(int i=0;i<cosRootsCount;i++)if(fabs(cosRootsArray[i])<1)
{//розрахунок fi і корекція rez
double x=acos(cosRootsArray[i]),
im=0,re=kf[0];

//розрахунок дійсної (re) і уявної (im) частини досліджуваного полінома
//при черговому значенні знайденого кореня
for(int j=1;j < =degree;j++)
{
re+=(kf[j]*rpw[j]*cos(j*x));
im+=(kf[j]*rpw[j]*sin(j*x));
}

//істотно ненульове значення im інтерпретуємо як внутрішню помилку
if(fabs(im)>1e-6)errorCode+=6;

//ми шукаємо мінімальне значення абсциси эпициклоиды
if(re<rez){rez=re;fi=x;}
}//розрахунок fi і корекція rez

//інтерпретуємо неймовірний випадок fi=0 або fi=cPI
//як внутрішню помилку
if(fi==0||fi==cPI)errorCode+=7;
return rez;
}//minAxeXcrossVal
//*************************************************************************



Текст файлу main.cpp


//*************************************************************************
//демонстрація розрахунку комплексного кореня многочлена
#include < stdio.h>
#include <math.h>
#include <polynomComplexRoot.h>

int main()
{
//завдання мірою і коефіцієнтів вихідного полінома
int degree=4;
//double kf[5]={24,24,12,4,1};
//double kf[5]={2,-2,3,-2,1};
double kf[5]={4,0,0,0,1};

//запуск процесу розрахунку комплексного кореня
double re,im;
int ret=polinomComplexRoot(degree,kf).run(re,im);

//друк результату розрахунку кореня
if(ret==0||ret>4)
{//успішне завершення процедури розрахунку кореня
//друк результату розрахунку кореня
if(ret==0)printf("успішне рішення\n");
else 
{
printf("ПОПЕРЕДЖЕННЯ\пбыли дивні ситуації\n");
printf("можливо із-за недостатньої точності апаратного подання дійсних чисел\n");
}
printf("кореньВещт=%f\пкореньМним=%f\n",re,im);
//перевірка достовірності знайденого кореня
//використовується наступна формула для полінома від суми аргументів:
//p0(x+y)=p0(x)+p1(x)*y+p2(x)*y^2/2!+...+pn*y^n/n!
//що представляє собою кінцеву суму ряду Тейлора по збільшенню y
//для вихідного полінома в точці p0 x
//тут p1, p2, ... , pn - похідні 1-го, 2-го, ... , n-го порядку від вихідного полінома
//у тексті програми використаний далі зворотний порядок індексації вихідної і похідних поліномів
//індекс многочлена дорівнює його ступеня
double** kfx=new double*[degree+1];
for(int i=0;i<=degree;i++)kfx[i]=new double[degree+1];
double* kfy=new double[degree+1];

for(int i=degree;i>=0;i--)
{//вихідні присвоєння
for(int j=0;j<=degree;j++)kfx[i][j]=0;
kfy[i]=0;
kfx[degree][i]=kf[i];
}//вихідні присвоєння

//розрахунок коефіцієнтів похідних поліномів
for(int i=degree;i>0;i--)
for(int j=i;j>0;j--)kfx[i-1][j-1]=kfx[i][j]*j;

double fact=1;
for(int i=0,dmi=degree;i<=degree;i++,dmi--)
{//розрахунок коефіцієнтів полінома по y
//обчислюється за схемою Горнера значення похідного полінома від аргументу re
kfy[i]=kfx[dmi][dmi];
for(int j=dmi-1;j>=0;j--)kfy[i]=re*kfy[i]+kfx[dmi][j];
kfy[i]/=fact;
fact*=(i+1);
}//розрахунок коефіцієнтів полінома по y

//масив ступенів уявної одиниці
const int ipw[4]={1,1,-1,-1};

double ypw=1;

//реальна та уявна частини значення вихідного полінома
//при значенні комплексного аргументу re+i*im
double fre=0,fim=0;
for(int i=0,od=0;i<=degree;i++,od=1-od)
{//розрахунок значення вихідного полінома
if(od==0)fre+=(ypw*kfy[i]*ipw[i%4]);
else fim+=(ypw*kfy[i]*ipw[i%4]);
ypw*=im;
}//обчислення значення вихідного полінома

//друк похибки розрахунку вихідного кореня многочлена
//це модуль полінома в знайденому корені,
//нормований на значення вільного члена
printf("похибка=%7.1 e\n",sqrt(fre*fre+fim*fim)/kf[0]);

//заключне звільнення займаної пам'яті
for(int i=0;i<=degree;i++)delete[] kfx[i];
delete[] kfx;
delete[] kfy;
}//успішне вирішення
if(ret>0)switch(ret)
{//друк короткою діагностики сталася помилки
case 1:
printf("#помилка 1:нульове значення вільного члена вихідного полінома\n");
break;
case 2:
printf("#помилка 2:нульове значення старшого коефіцієнта полінома\n");
break;
case 3:
printf("#помилка 3:ступінь вихідного полінома занадто мала\n");
break;
case 4:
printf("#помилка 4:вихідний поліном має речові коріння\n");
break;
default:printf("#код помилки=%d\n",ret);
}//друк короткою діагностики сталася помилки
return 0;
}//main
//*************************************************************************

Джерело: Хабрахабр

0 коментарів

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