Как найти комплексный корень многочлена

Содержимое

//*************************************************************************
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
//*************************************************************************

Содержимое

//*************************************************************************
#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
//*************************************************************************

//*************************************************************************
//демонстрация расчёта комплексного корня полинома
#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были странные ситуацииn");
printf("возможно из-за недостаточной точности аппаратного представления вещественных чиселn");
}
printf("кореньВещт=%fnкореньМним=%fn",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.1en",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("#код ошибки=%dn",ret);
}//печать краткой диагностики произошедшей ошибки
return 0;
}//main
//*************************************************************************

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

Пусть 1 = 1 + 1 и 2 = 2 + 2. Тогда основные арифметические операции с комплексными числами определяются так:

1)1 ± 2 = ( 1 ± 2) + ( 1 ± 2);

2)1 2 = ( 1 + 1)( 2 + 2) = ( 1 2 − 1 2) + ( 1 2 + 2 1);

3)1 = 1 + 1 = ( 1 + 1)( 2 − 2) =2 2 + 2 ( 2 + 2)( 2 − 2)

Имеет место равенство ( + ) = + . Два комплексных числа равны тогда и только тогда, когда равны их действительные и мнимые части:

та выражается через тригонометрические функции и наоборот:

Отсюда еще одна форма представления комплексного числа

–экспоненциальная:

= + = (cos + · sin ) = = ln + .

Из формулы Эйлера следует интересное отношение

= −1.

На множестве вещественных чисел мы лишены возможности видеть эти связи, как человек, рассматривающий с берега моря живописные острова, не подозревает, что они всего лишь вершины сложной горной системы подводного царства. На комплексной плоскости «в порядке вещей» многое из того, что категорически запрещено на вещественной оси. Здесь существуют решения уравнений ( ) = 5, 2 = −3

и т. д. Конечно, и теперь мы рассматриваем комплексное число как некую абстракцию, но не в большей мере, чем число «три». Ведь числа «три» также в природе не существует! Естественно возникает вопрос: не придется ли для решения других уравнений, например 4 + 1 = 0, и даль-

144 ГЛАВА 3. УРАВНЕНИЯ СТАРШЕГО ПОРЯДКА

Таким образом, при перемножении комплексных чисел их модули перемножаются, а аргументы складываются (рис. 28).

Рис. 28. Умножение комплексных чисел

Аналогично

Следовательно,

| 1 2|= | 1|| 2|, ( 1 2) = ( 1) + ( 2),

В частности,

1

cos + · sin = cos − · sin .

Из графика многочлена видно, что у него три корня: 1, 2 и 3.