Задача 7. Кубическое уравнение

Найти все корни (в том числе и комплексные) уравнения x³+ax²+bx+c= =0, где a, b, c - заданные вещественные коэффициенты. Использовать решение Кардано.

{ Решение задачи 7 на языке PASCAL }

TYPE Real = Extended;

Complex = RECORD Re,Im : Real; END;

Troot = ARRAY[1..3] OF Complex;

PROCEDURE ComplexRoot2(x:Real; VAR x1_2:Complex);

BEGIN WITH x1_2 DO

IF x>=0 THEN BEGIN Re:=Sqrt(x); Im:=0; END

ELSE BEGIN Re:=0; Im:=Sqrt(Abs(x)); END;

END;

PROCEDURE ComplexAdd(x1,x2:Complex; VAR y:Complex);

BEGIN WITH y DO BEGIN Re:=x1.Re+x2.Re; Im:=x1.Im+x2.Im; END;

END;

PROCEDURE ComplexRoot3(x:Complex; VAR x1_3:Complex);

VAR z,f : Real;

BEGIN

WITH x DO BEGIN

z:=Sqrt(Sqr(Re)+Sqr(Im));

IF Re=0 THEN BEGIN IF Im>=0 THEN f:=0.5*Pi ELSE f:=1.5*Pi END

ELSE BEGIN f:=Arctan(Im/Re);

IF Re<0 THEN f:=Pi+f ELSE IF Im<0 THEN f:=2*Pi+f;

END;

END;

IF z>0 THEN z:=Exp(Ln(z)/3);

WITH x1_3 DO BEGIN Re:=z*Cos(f/3); x1_3.Im:=z*Sin(f/3); END;

END;

PROCEDURE ComplexMult(x1,x2:Complex; VAR y:Complex);

BEGIN

WITH y DO BEGIN

Re:=x1.Re*x2.Re-x1.Im*x2.Im; Im:=x1.Re*x2.Im+x1.Im*x2.Re; END;

END;

PROCEDURE PrintComplex(x:Complex);

BEGIN Write(x.Re:10:5);

IF x.Im>=0 THEN Write('+',x.Im:8:5) ELSE Write(x.Im:9:5); Write('i');

END;

FUNCTION Min3(x1,x2,x3:Real):Byte;

BEGIN IF x1<x2 THEN IF x1<x3 THEN Min3:=1 ELSE Min3:=3

ELSE IF x2<x3 THEN Min3:=2 ELSE Min3:=3;

END;

VAR a,b,c,p,q,QQ : Real;

AA,BB,QQ1_2,Tmp,Tmp1,Tmp2,d,One_1_3_p,One_1_3_m : Complex;

y,x : TRoot;

i,j : Byte;

BEGIN

WITH One_1_3_p DO BEGIN Re:=-0.5; Im:= Sqrt(3)/2; END;

WITH One_1_3_m DO BEGIN Re:=-0.5; Im:=-Sqrt(3)/2; END;

Write('Введите коэффициенты a,b,c '); ReadLn(a,b,c);

p:=-Sqr(a)/3+b; q:=2*Sqr(a/3)*(a/3)-a*b/3+c;

QQ:=Sqr(p/3)*(p/3)+Sqr(q/2); ComplexRoot2(QQ,QQ1_2);

WITH AA DO BEGIN Re:=QQ1_2.Re-q/2; Im:=QQ1_2.Im; END;

WITH BB DO BEGIN Re:=-QQ1_2.Re-q/2; Im:=-QQ1_2.Im; END;

ComplexRoot3(AA,AA); ComplexRoot3(BB,BB); ComplexMult(AA,BB,Tmp);

CASE Min3(Abs(p/3+Tmp.Re),Abs(p/3-Tmp.Re/2-Sqrt(3)*Tmp.Im/2),

Abs(p/3-Tmp.Re/2+Sqrt(3)*Tmp.Im/2)) OF

2:ComplexMult(AA,One_1_3_p,AA);

3:ComplexMult(AA,One_1_3_m,AA);

END;

ComplexAdd(AA,BB,y[1]);

With Tmp1 DO BEGIN

Re:=(AA.Re+BB.Re)/2; Im:=(AA.Im+BB.Im)/2; END;

With Tmp2 DO BEGIN Re:=AA.Re-BB.Re; Im:=AA.Im-BB.Im; END;

With Tmp DO BEGIN Re:=0; Im:=Sqrt(3)/2; END;

ComplexMult(Tmp,Tmp2,Tmp2);

With y[2] DO BEGIN

Re:=-Tmp1.Re+Tmp2.Re; Im:=-Tmp1.Im+Tmp2.Im; END;

With y[3] DO BEGIN

Re:=-Tmp1.Re-Tmp2.Re; Im:=-Tmp1.Im-Tmp2.Im; END;

WriteLn(' Корни уравнения и контроль :');

FOR i:=1 TO 3 DO BEGIN

x[i].Re:=y[i].Re-a/3; x[i].Im:=y[i].Im;

ComplexMult(x[i],x[i],Tmp1); ComplexMult(x[i],Tmp1,Tmp2);

d.Re:=Tmp2.Re+a*Tmp1.Re+b*x[i].Re+c;

d.Im:=Tmp2.Im+a*Tmp1.Im+b*x[i].Im;

PrintComplex(x[i]); Write(' ----> '); PrintComplex(d);WriteLn(' = 0');

END;

WriteLn('Нажмите Enter'); ReadLn;

END.

/*Решение задачи 7 на языке С*/

#include<math.h>

#include<stdio.h>

#include<conio.h>

struct s_complex{double Re,Im;};

typedef struct s_complex complex;

complex ComplexRoot2(double x) { complex x_1_2;

if(x>=0){x_1_2.Re=sqrt(x); x_1_2.Im=0;}

else{x_1_2.Re=0; x_1_2.Im=sqrt(fabs(x));}return x_1_2;}

complex ComplexAdd(complex x1,complex x2) {complex y;

y.Re=x1.Re+x2.Re; y.Im=x1.Im+x2.Im; return y;}

complex ComplexRoot3(complex x){complex x1_3; double z,f;

z=sqrt(pow(x.Re,2)+pow(x.Im,2)); f=atan2(x.Im,x.Re); if(z) z=pow(z,1./3);

x1_3.Re=z*cos(f/3); x1_3.Im=z*sin(f/3); return x1_3;}

complex ComplexMult(complex x1,complex x2){complex y;

y.Re=x1.Re*x2.Re-x1.Im*x2.Im; y.Im=x1.Re*x2.Im+x1.Im*x2.Re;

return y;}

int Min3(double x1,double x2,double x3)

{if(x1<x2)if(x1<x3)return 1;else return 3;

else if(x2<x3)return 2;else return 3;}

void PrintComplex(complex x){ printf("%10.5f",x.Re);

if(x.Im>0)printf("+%7.5fi",x.Im);else printf("%8.5fi",x.Im);}

void main(void) { double a,b,c,p,q,Q; int i,j; complex A,B,Q1_2,Tmp,Tmp1, Tmp2,y[3],x[3],d,One_1_3_p={-.5,.5},One_1_3_m={-.5,-.5};

One_1_3_p.Im*=sqrt(3); One_1_3_m.Im*=sqrt(3);

printf("\nВведите коэффициенты a,b,c "); scanf("%lg%lg%lg",&a,&b,&c);

p=-a*a/3+b; q=2*pow(a/3,3)-a*b/3+c; Q=pow(p/3,3)+pow(q/2,2);

Q1_2=ComplexRoot2(Q); A.Re=Q1_2.Re-q/2; A.Im=Q1_2.Im;

B.Re=-Q1_2.Re-q/2; B.Im=-Q1_2.Im; A=ComplexRoot3(A);

B=ComplexRoot3(B); Tmp=ComplexMult(A,B);

switch(Min3(fabs(p/3+Tmp.Re), fabs(p/3-Tmp.Re/2-sqrt(3)*Tmp.Im/2),

fabs(p/3-Tmp.Re/2+sqrt(3)*Tmp.Im/2))){

case 2:A=ComplexMult(A,One_1_3_p); break;

case 3:A=ComplexMult(A,One_1_3_m);}

y[0]=ComplexAdd(A,B);

Tmp1.Re=(A.Re+B.Re)/2; Tmp1.Im=(A.Im+B.Im)/2;

Tmp2.Re=A.Re-B.Re; Tmp2.Im=A.Im-B.Im;

Tmp.Re=0; Tmp.Im=sqrt(3)/2;

Tmp2=ComplexMult(Tmp,Tmp2);

y[1].Re=-Tmp1.Re+Tmp2.Re; y[1].Im=-Tmp1.Im+Tmp2.Im;

y[2].Re=-Tmp1.Re-Tmp2.Re; y[2].Im=-Tmp1.Im-Tmp2.Im;

printf("\nКорни уравнения и контроль :\n");

for(i=0;i<3;i++){x[i].Re=y[i].Re-a/3; x[i].Im=y[i].Im;

Tmp1=ComplexMult(x[i],x[i]); Tmp2=ComplexMult(x[i],Tmp1);

d.Re=Tmp2.Re+a*Tmp1.Re+b*x[i].Re+c;

d.Im=Tmp2.Im+a*Tmp1.Im+b*x[i].Im;

PrintComplex(x[i]);printf(" ----> ");PrintComplex(d);printf(" == 0\n");}

printf("\nНажмите любую клавишу\n"); getch();}

// Решение задачи 7 на языке C++

#include<stdio.h>

#include<iostream.h>

#include<iomanip.h>

#include<conio.h>

#include<math.h>

#include<complex.h>

void main(void){

double a,b,c,p,d1,d2,d3; int i,j;

complex q,Q,tmp,A,B,y[3],x[3],

One_1_3_p=complex(-.5,.5*sqrt(3)),One_1_3_m=complex(-.5,-.5*sqrt(3));

cout<<setfill(' ' )<<setprecision(5)<<

setiosflags(ios::showpoint | ios::fixed | ios::right);

cout<<"Введите коэффициенты a,b,c "; cin>>a>>b>>c;

p=-a*a/3+b; q=2*pow(a/3,3)-a*b/3+c; Q=pow(p/3,3)+pow(q/2,2);

A=pow(-q/2+sqrt(Q),complex(1./3,0));

B=pow(-q/2-sqrt(Q),complex(1./3,0));

tmp=A*B; d1=fabs(p/3+real(tmp));

d2=fabs(p/3-real(tmp)/2-sqrt(3)*imag(tmp)/2);

d3=fabs(p/3-real(tmp)/2+sqrt(3)*imag(tmp)/2);

if(d2<d3){if(d2<d1)A*=One_1_3_p;}else if(d3<d1)A*=One_1_3_m;

y[0]=A+B; y[1]=-(A+B)/2+complex(0,sqrt(3)/2)*(A-B);

y[2]=-(A+B)/2-complex(0,sqrt(3)/2)*(A-B);

for(i=0;i<3;i++)x[i]=y[i]-a/3;

cout <<"\nКорни уравнения и контроль:\n";

for(i=0;i<3;i++){ tmp=pow(x[i],3)+a*x[i]*x[i]+b*x[i]+c;

cout<<setw(10)<<real(x[i])<<setw(10)<<imag(x[i])<<"i ---->"<<setw(10)

<<real(tmp)<<setw(10)<<imag(tmp)<<"i = 0"<<endl;}

cout <<"\nНажмите любую клавишу\n"; getch();}

C Решение задачи 7 на языке FORTRAN

DOUBLE PRECISION a,b,c,p,q,QQ,d1,d2,d3

COMPLEX*16 tmp,AA,BB,y(3),x(3),One_1_3_p,One_1_3_m

PRINT*,'Введите коэффициенты a,b,c '

READ*,a,b,c

One_1_3_p=(-.5,0)+(0,.5)*DSQRT(3)

One_1_3_m=CONJG(One_1_3_p)

p=-a**2/3+b

q=2*(a/3)**3-a*b/3+c

QQ=(p/3)**3+(q/2)**2

AA=(-q/2+CDSQRT(QQ))**(1./3)

BB=(-q/2-CDSQRT(QQ))**(1./3)

tmp=AA*BB

d1=DABS(p/3+DREAL(tmp))

d2=DABS(p/3-DREAL(tmp)/2-DSQRT(3)*DIMAG(tmp)/2)

d3=DABS(p/3-DREAL(tmp)/2+DSQRT(3)*DIMAG(tmp)/2)

IF(d2.EQ.MIN(d1,d2,d3)) AA=AA*One_1_3_p

IF(d3.EQ.MIN(d1,d2,d3)) AA=AA*One_1_3_m

y(1)=AA+BB

y(2)=-(AA+BB)/2+(0,.5)*DSQRT(3)*(AA-BB)

y(3)=-(AA+BB)/2-(0,.5)*DSQRT(3)*(AA-BB)

x=y-a/3

y=x**3+a*x**2+b*x+c

PRINT*,'Корни уравнения и контроль'

PRINT 999,(x(i),y(i),i=1,3)

999 FORMAT(1X,SP,2G12.6,'i ---->',2G12.6,'i = 0')

PAUSE 'Нажмите ENTER'

END