Это завершение моей первой статьи: "Алгоритм расчёта вещественных корней полиномов".

Спасибо комментаторам, сделавшим более ясным мое слишком уж конспективное изложение метода Лобачевского. В самом деле, мне следовало явно написать, что квадрированный полином надо рассматривать как полином от аргумента 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были странные ситуации\n");
printf("возможно из-за недостаточной точности аппаратного представления вещественных чисел\n");
}
printf("кореньВещт=%f\nкореньМним=%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.1e\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
//*************************************************************************

Поделиться с друзьями
-->

Комментарии (3)


  1. Playa
    17.06.2016 03:50

    Добавьте, пожалуйста, файл «polynomRealRoots.h»


  1. icoz
    19.06.2016 09:21

    И хотя бы одну картинку с иллюстрацией того, что происходит…


  1. luarviq
    19.06.2016 12:11

    Вызывает сомнения вот этот участок кода:
    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