Решение линейных уравнений 3-й и 4-й
степеней в комплексных числах на C++
|
| !!! Внимание !!! Предложенный код не был
использован в реальных приложениях - также не проводились исследования
на
корректность источников формул на основании которых был построен код. |
Здесь далее приведены тексты файла кодов на С++ и его заголовочного файла (RadicalsX4.cpp и RadicalsX4.h соответственно).
Итак, скопировав нижеприведенный текст в какой нибудь "блокнот" и сохранив с указанными именами, получаем файлы которые можно подключать к своим проектам. Возможно придется убрать в строке " #include <complex.h>" расширение файла ".h" - по крайней мере в справке к среде разработки (я пользовался "Borland C++Builder 4") имя библиотечного файла было указано без расширения - но у меня оно потребовалось. Функции решают уравнения 2, 3, и 4-й степеней в
комплексных числах. Исходные коэффициенты уравнений a, b, c и d -
задаются в том виде - как они находятся в закомментированных
пояснениях. Корни уравнений содержатся в массивах x
которые должны быть обязательно предварительно объявлены.
Необязательная переменная ToI -указывает количество нулей после запятой
когда комплексная и действительная части округляются до нуля.
Получить действительную и комплексную составляющие
чисел можно при помощи функций real и imag например real(x[0]) и
imag(x[0]).
Также имеется функция нахождения всех трех
кубических корней из комплексного числа - просто применяется для
решения по формуле Кардано.
Комментарий окончен - далее код:
| // Есьман Юрий
Григорьевич 2006г //--------------------------------------------------------------------------- #pragma hdrstop #include <complex.h> #include "RadicalsX4.h" //--------------------------------------------------------------------------- #pragma package(smart_init) //--------------------------------------------------------------------------- inline __fastcall int RadicalsX3(double a,double b,double c,double d,radicalsX3 &x,int ToI) {//решение уравнения третьей степени a*x3+b*x2+c*x+d=0 radicalsX2 y; radicalsX3 A,B; int r;//счетчик действительных корней if (a==0)//уравнение не является уравнением третьей степени {r=RadicalsX2((complex<double>)b,(complex<double>)c,(complex<double>)d,y); x[0]=y[0];x[1]=y[1];if (r!=2)return r;//одно или нет решений else {if (imag(y[0])!=0) return 0;//нет действительных корней else return 2;}}//два действительных корня d=d/a;c=c/a;b=b/a;//приведение к виду x3+b*x2+c*x+d=0 double p=c-b*b/3,q=d+(b*b/13.5-c/3)*b;//приведение к виду x3+p*x+q=0 complex<double> DL=-27*q*q-4*p*p*p;//дискриминант уравнения DL=sqrt(-DL/(complex<double>)108); double Roundoff=exp(log(0.1)*ToI);//минимальная учитываемае величина комплексной составляющей числа try{ complex<double> ab=-(complex<double>)q/(complex<double>)2; sqrt3(ab+DL, &A); sqrt3(ab-DL, &B); ab=A[0]*B[0]; if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff) &&(imag(ab)>-Roundoff)) x[0]=A[0]+B[0]; else {ab=A[0]*B[1]; if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff) &&(imag(ab)>-Roundoff)) x[0]=A[0]+B[1]; else {ab=A[0]*B[2]; if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff) &&(imag(ab)>-Roundoff)) x[0]=A[0]+B[2]; else {ab=A[1]*B[1]; if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff) &&(imag(ab)>-Roundoff)) x[0]=A[1]+B[1]; else {ab=A[1]*B[2]; if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff) &&(imag(ab)>-Roundoff)) x[0]=A[1]+B[2]; else {ab=A[2]*B[2]; if ((real(ab)<(-p/3+Roundoff))&&(real(ab)>(-p/3-Roundoff))&&(imag(ab)<Roundoff) &&(imag(ab)>-Roundoff)) x[0]=A[2]+B[2]; }}}}} }catch ( ... ){return 0;}//уравнение не разрешилось x[0]=x[0]-b/3;//обратная подстановка r=RadicalsX2((complex<double>)1,(complex<double>)(b+x[0]),(complex<double>)((b+x[0])*x[0]+c),y); x[1]=y[0]; x[2]=y[1]; if (r>2) return 999999;//решений бесконечное множество if (r==0) return 0;//уравнение не разрешилось r=0; for (int i=0;i<3;i++) if (((imag(x[i]))<Roundoff)&&((imag(x[i]))>-Roundoff)) {x[i]=real(x[i]);//округляем мелкие комплексные части r++;}//считаем действительные корни return r; } //--------------------------------------------------------------------------- inline __fastcall int RadicalsX2(complex<double> a,complex<double> b,complex<double> c,radicalsX2 &x) {//решение уравнения второй степени a*x2+b*x+c=0 if (a==(complex<double>)0)//линейное уравнение b*x+c=0 {if (b==(complex<double>)0)//нет переменных {if (c==(complex<double>)0) return 999999;//бесконечное множиство решений else return 0;}//нет решений x[0]=x[1]=-c/b;return 1;}//едиственное решение x[1]=sqrt(b*b-((complex<double>)4*a*c)); x[0]=(-b+x[1])/((complex<double>)2*a); x[1]=(-b-x[1])/((complex<double>)2*a); return 2;//уравнение действительно квадратное } //--------------------------------------------------------------------------- inline __fastcall void sqrt3(complex<double> a,radicalsX3 &x) {//три корня третьей степени из комплексного числа а double X=real(a),y=imag(a),R=sqrt(y*y+X*X); if (R==0) {x[0]=x[1]=x[2]=0; return;} double Q=acos(X/R); complex<double> i(0.0, 1.0); R=exp(log(R)/3); if (y<0) Q=6.28318530717958647692528676655901-Q; x[0]=R*cos(Q/3)+R*sin(Q/3)*i; x[1]=R*cos((Q+3.1415926535897932384626433832795)/3)+R*sin((Q+3.1415926535897932384626433832795)/3)*i; x[2]=R*cos((Q+6.28318530717958647692528676655901)/3)+R*sin((Q+6.28318530717958647692528676655901)/3)*i; } //--------------------------------------------------------------------------- inline __fastcall int RadicalsX4(double a,double k,double t,double p,double q,radicalsX4 &x,int ToI) {//решение уравнения четвертой степени a*x4+k*x3+t*x2+p*x+q=0 int r;//счетчик действительных корней radicalsX3 A;//корни уравнения третьей степени radicalsX2 x01,x23;//корни урравнений второй степени double Roundoff=exp(log(0.1)*ToI);//минимальная учитываемае величина комплексной составляющей числа if (a==0)//уравнение не является уравнением четвертой степени, решаем урвнение третьей степени {r=RadicalsX3(k,t,p,q,A,ToI);x[0]=A[0];x[1]=A[1];x[2]=A[2];return r;} q=q/a;p=p/a;t=t/a;k=k/a;//приведение к виду x4+k*x3+t*x2+p*x+q=0 k=k/4; q=q+((t-k*k*3)*k-p)*k; p=p+k*(k*k*8-t*2); t=t-k*k*6;//приведение к виду x4+t*x2+p*x+q=0 r=RadicalsX3(1,t*2,t*t-q*4,-p*p,A,ToI);//решение приведенного уравнения третьей степени if ((real(A[0])<0)&&(real(A[0])>-Roundoff))A[0]=0;//приближение малых отрицательных чисел к нулю (избегаем комплексных решений) else A[0]=0; // первый член всегда положителен -комплексная часть результат погрешностей complex<double> A12=sqrt(A[0]),AT12=sqrt((A[0]+t)*(A[0]+t)-q*4); RadicalsX2(1,-A12,(A[0]+t+AT12)*0.5,x01); RadicalsX2(1,A12,(A[0]+t-AT12)*0.5,x23); x[0]=x01[0];x[1]=x01[1];x[2]=x23[0];x[3]=x23[1]; x[0]=x[0]-k;x[1]=x[1]-k;x[2]=x[2]-k;x[3]=x[3]-k; //используем только первый из трех корней A - две другие четверки корней примерно равны - с точностью вычислений r=0; for (int i=0;i<4;i++) {if (((imag(x[i]))<Roundoff)&&((imag(x[i]))>-Roundoff)) {x[i]=real(x[i]);//округляем мелкие комплексные части A12=x[r];x[r]=x[i];x[i]=A12;//перемещаем целые части вперед r++;}//считаем действительные корни if (((real(x[i]))<Roundoff)&&((real(x[i]))>-Roundoff)) x[i]=x[i]-real(x[i]);//округляем мелкие целые части } return r; } //--------------------------------------------------------------------------- |
| // Есьман Юрий
Григорьевич 2006г //--------------------------------------------------------------------------- #ifndef RadicalsX4H #define RadicalsX4H //--------------------------------------------------------------------------- typedef complex<double> radicalsX4[4];//корни уравнения 4-й степени typedef complex<double> radicalsX3[3];//корни уравнения 3-й степени typedef complex<double> radicalsX2[2];//корни уравнения 2-й степени //решение уравнения третьей степени a*x3+b*x2+c*x+d=0 inline __fastcall int RadicalsX3(double a,double b,double c,double d,radicalsX3 &x,int ToI=10); //решение уравнения второй степени a*x2+b*x+c=0 inline __fastcall int RadicalsX2(complex<double> a,complex<double> b,complex<double> c,radicalsX2 &x); //три корня третьей степени из комплексного числа а inline __fastcall void sqrt3(complex<double> a,radicalsX3 &x); //решение уравнения четвертой степени a*x4+k*x3+t*x2+p*x+q=0 inline __fastcall int RadicalsX4(double a,double k,double t,double p,double q,radicalsX4 &x,int ToI=10); //--------------------------------------------------------------------------- #endif |
|
|
|