Wielomian 3 stopnia

kolumb · Post autor: **kolumb** » 30 cze 2006, o 11:30

Witam

szukam jakiegos uniwersalnego sposobu rozwiazywania wielomianow stopnia 3, tak zeby dalo sie to wykorzytsac przy pisaniu programu w C++ do obliczania tych pierwiastkow.

pozdr

juzef · Post autor: **juzef** » 30 cze 2006, o 11:59

Wzory Cardano nie wystarczają?

bisz · Post autor: **bisz** » 30 cze 2006, o 14:04

Kod: Zaznacz cały

#include <stdio.h>
#include <conio.h>
#include <iostream.h>
#include <math.h>

main()
{
       
    double long a,b,c,d,aa,bb,cc,lol,dd,ee,t1,t2,k1,k2,k3,k4,k5,k6,k7,k8,eu1,eu2,eu3,n,t3,w1,w2,w3,w4,w5,w6,w7,w8,t4,re,im,x1,x2,x3,y1,y2,y3,p,qq,pp,rr,q,z2,z1,z3,k,uu,u,vv,v,arg,delta;
     double long pi=3.1415926535897932384626433832795;
      cout<<"
Wielomian ax^4+bx^3+cx^2+dx+e=0";
      cout<<"
Podaj wspolczynnik a:";
      cin>>aa;
      cout<<"
Podaj wspolczynnik b:";
       cin>>bb;
      cout<<"
Podaj wspolczynnik c:";
       cin>>cc;
      cout<<"
Podaj wspolczynnik d:";
       cin>>dd;
       cout<<"
Podaj wspolczynnik e:";
       cin>>ee;
              
   
       pp=(-3.0*bb*bb)/(8.0*aa*aa)+(cc/aa);
       qq=(bb*bb*bb)/(8.0*aa*aa*aa)-(bb*cc)/(2.0*aa*aa)+(dd/aa);
       rr=-3.0*bb*bb*bb*bb/(256.0*aa*aa*aa*aa)+bb*bb*cc/(16.0*aa*aa*aa)-bb*dd/(4.0*aa*aa)+ee/aa;

    
               
      a=1;
      b=-2*pp;
      c=pp*pp-4*rr;
      d=qq*qq;
      
      p=c/a-b*b/(3*a*a);
      q=2*b*b*b/(27*a*a*a)+d/a-b*c/(3*a*a);
     delta=(p/3)*(p/3)*(p/3)+(q/2)*(q/2);
     if (sqrt(delta*delta)<0.00000001)
     {delta=0;}
      
        
     if(delta==0)
      { 
      if (q<0)
      {q=-q;
      
      y1=-pow((q/2),(1.0/3.0));
      y2=2*pow((q/2),(1.0/3.0));
      
       x1=y1-b/(3*a);
      x2=y2-b/(3*a);
      cout<<"
x= "<<x1;
      cout<<"
x= "<<x2;
      }
       else
       {          
      y1=pow((q/2),(1.0/3.0));
      y2=-2*pow((q/2),(1.0/3.0));
      
       x1=y1-b/(3*a);
      x2=y2-b/(3*a);
       }
       z1=x1;
       z2=x1;
       z3=x2;
      
      cout<<"
t= "<<t2;
      cout<<"
t= "<<t3;
      cout<<"
t= "<<t4;
      cout<<"
x1= "<<x1;
      cout<<"
x2= "<<x2;
       
             }
      if(delta<0)
     {
     arg=acos(3*q/(2*p*sqrt(-p/3)));
     y1=2*sqrt(-p/3)*cos((arg)/3);
     y2=2*sqrt(-p/3)*cos((arg+2*pi)/3);
     y3=2*sqrt(-p/3)*cos((arg+4*pi)/3);
     x1=y1-b/(3*a);
     if (sqrt(x1*x1)<0.0000001)
      {x1=0;}     
     x2=y2-b/(3*a);
     x3=y3-b/(3*a);
     
       z1=x1;
       
       z2=x2;
       z3=x3;
      a=0;
      b=0;
       
      if(z1>0)
       {eu1=-sqrt(z1)*0.5;
       ;
       b=b+eu1;
       }
       else
       {eu1=-sqrt(-z1)*0.5;
       a=a+eu1;
       ;}
        if(z2>0)
       {eu2=-sqrt(z2)*0.5;
       b=b+eu2;
       
       }
       else
       {eu2=-sqrt(-z2)*0.5;
       a=a+eu2;
       }
        if(z3>0)
       {eu3=-sqrt(z3)*0.5;
       b=b-eu3;
       }
       else
       {eu3=-sqrt(-z3)*0.5;
       a=a-eu3;
     }
     
     
      arg=atan(b/a);
      re=(sqrt(a*a+b*b)*sqrt(a*a+b*b)*sqrt(a*a+b*b)*sqrt(a*a+b*b))*(cos(arg*4))+pp*(sqrt(a*a+b*b)*sqrt(a*a+b*b))*(cos(arg*2))+qq*sqrt(a*a+b*b)*cos(arg)+rr;
      im=(sqrt(a*a+b*b)*sqrt(a*a+b*b)*sqrt(a*a+b*b)*sqrt(a*a+b*b))*(sin(arg*4))+pp*(sqrt(a*a+b*b)*sqrt(a*a+b*b))*(sin(arg*2))+qq*sqrt(a*a+b*b)*sin(arg);
      
      cout<<"
re="<<re;
      cout<<"
re="<<im;
      
      cout<<"
"<<a<<"+"<<b;
       
      k1=+eu1+eu2+eu3;
      k2=+eu1+eu2-eu3;
      k3=+eu1-eu2+eu3;
      k4=+eu1-eu2-eu3;
      k5=-eu1+eu2+eu3;
      k6=-eu1+eu2-eu3;
      k7=-eu1-eu2+eu3;
      k8=-eu1-eu2-eu3;
      
      w1=k1*k1*k1*k1+pp*k1*k1+qq*k1+rr;
      w2=k2*k2*k2*k2+pp*k2*k2+qq*k2+rr;
      w3=k3*k3*k3*k3+pp*k3*k3+qq*k3+rr;
      w4=k4*k4*k4*k4+pp*k4*k4+qq*k4+rr;
      w5=k5*k5*k5*k5+pp*k5*k5+qq*k5+rr;
      w6=k6*k6*k6*k6+pp*k6*k6+qq*k6+rr;
      w7=k7*k7*k7*k7+pp*k7*k7+qq*k7+rr;
      w8=k8*k8*k8*k8+pp*k8*k8+qq*k8+rr;
      
       if(sqrt(w1*w1)<0.0001)
      {cout<<"
x1= "<<k1-bb/(4.0*aa);
      } 
        if(sqrt(w2*w2)<0.0001)
      {cout<<"
x2= "<<k2-bb/(4.0*aa);
      } 
        if(sqrt(w3*w3)<0.0001)
      {cout<<"
x3= "<<k3-bb/(4.0*aa);
      }  
     if(sqrt(w4*w4)<0.0001)
      {cout<<"
x4= "<<k4-bb/(4.0*aa);
      } 
      if(sqrt(w5*w5)<0.0001)
      {cout<<"
x5= "<<k5-bb/(4.0*aa);
      } 
        if(sqrt(w6*w6)<0.0001)
      {cout<<"
x6= "<<k6-bb/(4.0*aa);
      } 
       if(sqrt(w7*w7)<0.0001)
      {cout<<"
x7= "<<k7-bb/(4.0*aa);
      } 
        if(sqrt(w8*w8)<0.0001)
      {cout<<"
x8= "<<k8-bb/(4.0*aa);
      } 

     
       }   
            
    if(delta>0)
    { 
    lol=(-108*q+12*sqrt(12*p*p*p+81*q*q));
      x1=  1.0/6.0* pow((-108*q+12*sqrt(12*p*p*p+81*q*q)),(1.0/3.0))-2*p/pow((-108*q+12*sqrt(12*p*p*p+81*q*q)),1.0/3.0)-b/(3*a);
      if(lol<0)
      {
      x1=-( pow((-(-q/2.0+sqrt(delta))),(1.0/3.0))+ pow((-(-q/2.0-sqrt(delta))),(1.0/3.0))+b/(3*a));
     }
      cout<<"
x= "<<x1;
      cout<<"
x= "<<-(a*x1+b)/(2*a)<<"+"<<sqrt(-((a*x1+b)*(a*x1+b)-4*a*(a*x1*x1+b*x1+c)))/(2*a)<<"i";
      cout<<"
x= "<<-(a*x1+b)/(2*a)<<"-"<<sqrt(-((a*x1+b)*(a*x1+b)-4*a*(a*x1*x1+b*x1+c)))/(2*a)<<"i";
      
      z1=x1;
       z2=x2;
       z3=x3;
       
       t4=0.5*(-sqrt(-z1)+sqrt(-z2)+sqrt(-z3))-bb/(4.0*aa);
       t1=0.5*(sqrt(-z1)-sqrt(-z2)+sqrt(-z3))-bb/(4.0*aa);
       t2=0.5*(sqrt(-z1)+sqrt(-z2)-sqrt(-z3))-bb/(4.0*aa);
       t3=0.5*(-sqrt(-z1)-sqrt(-z2)-sqrt(-z3))-bb/(4.0*aa);
       
       
      cout<<"
t= "<<t1;
      cout<<"
t= "<<t2;
      cout<<"
t= "<<t3;
      cout<<"
t= "<<t4;
      
      
       }     
getch();
}

to byl moj 1 prog w c++ wiec troche kijowo napisany mam tez do 4 stopnia ktory wykorzystuje klase liczb zespolonych i opiera sie tez na wielomianie 3 stopnia no i jest troche czytelniejszy

Kod: Zaznacz cały

#include <iostream>
#include <complex.h>
#include <math.h>
#define eps 0.0000001

using namespace std; 
typedef complex<double> cmpx;

ostream& operator<<(ostream& s, cmpx &z)
{
if(fabs(z.imag())<eps)
{return s << z.real() ;}
if(fabs(z.real())<eps)
{return s << z.imag() <<" i" ;}

if(z.imag()<0)
{return s << z.real() <<" - " <<-z.imag() <<" i";}
if(z.imag()==0)
{return s << z.real() ; }
if(z.imag()>0 )
{return s << z.real() <<" + " << z.imag() <<" i"; }

}




void druk(cmpx z);

cmpx i(0,1.0);
cmpx x1,x2,x3,z1,z2,z3,z4,il,w;
cmpx z5,z6,z7,z8;
double a, b, c, d,aa,bb,cc,dd,ee,pp,qq,rr,p, q;
cmpx z[8];
cmpx wielomian(cmpx x);

main()
{
cout<<"Rownanie wyglada :
 ax^4+bx^3+cx^2+dx+e=0
";
cout<<"Podaj a: ";
cin>>aa;
cout<<"Podaj b: ";
cin>>bb;
cout<<"Podaj c: ";
cin>>cc;
cout<<"Podaj d: ";
cin>>dd;
cout<<"Podaj e: ";
cin>>ee;

pp=-3.0*bb*bb/(8.0*aa*aa)+1.0*cc/aa;
qq=bb*bb*bb/(8.0*aa*aa*aa)-bb*cc/(2.0*aa*aa)+1.0*dd/aa;
rr=-3.0*bb*bb*bb*bb/(256.0*aa*aa*aa*aa)+cc*bb*bb/(16.0*aa*aa*aa)-1.0*bb*dd/(4.0*aa*aa)+1.0*ee/aa;


a=1;
b=2.0*pp;
c=pp*pp-4.0*rr;
d=-1.0*(qq*qq);

p=1.0*c/a-b*b/(3.0*a*a);
q=2*b*b*b/(27.0*a*a*a)+1.0*d/a-b*c/(3.0*a*a);


x1=pow(-q/2.0+sqrt(q*q/4.0+p*p*p/27.0+0.0*i) ,1/3.0)+pow(-q/2.0-sqrt(q*q/4.0+p*p*p/27.0+0.0*i) ,1/3.0)-b/(3.0*a);
x2=0.5*(-1.0+sqrt(3)*1.0*i)*pow(-q/2.0+sqrt(q*q/4.0+p*p*p/27.0+0.0*i) ,1/3.0)+0.5*(-1.0-sqrt(3)*1.0*i)*pow(-q/2.0-sqrt(q*q/4.0+p*p*p/27.0+0.0*i) ,1/3.0)-b/(3.0*a);
x3=0.5*(-1.0-sqrt(3)*1.0*i)*pow(-q/2.0+sqrt(q*q/4.0+p*p*p/27.0+0.0*i) ,1/3.0)+0.5*(-1.0+sqrt(3)*1.0*i)*pow(-q/2.0-sqrt(q*q/4.0+p*p*p/27.0+0.0*i) ,1/3.0)-b/(3.0*a);

if(fabs(x1.real())<eps) x1=0;
if(fabs(x2.real())<eps) x2=0;
if(fabs(x3.real())<eps) x3=0;

il=sqrt(x1)*sqrt(x2)*sqrt(x3);

cout<<"iloczn to "<<il;

z1=0.5*(+sqrt(x1)+sqrt(x2)+sqrt(x3))-bb/(4.0*aa);
z2=0.5*(+sqrt(x1)-sqrt(x2)-sqrt(x3))-bb/(4.0*aa);
z3=0.5*(-sqrt(x1)+sqrt(x2)-sqrt(x3))-bb/(4.0*aa);
z4=0.5*(-sqrt(x1)-sqrt(x2)-sqrt(x3))-bb/(4.0*aa);
z5=0.5*(-sqrt(x1)-sqrt(x2)+sqrt(x3))-bb/(4.0*aa);
z6=0.5*(-sqrt(x1)+sqrt(x2)+sqrt(x3))-bb/(4.0*aa);
z7=0.5*(+sqrt(x1)-sqrt(x2)+sqrt(x3))-bb/(4.0*aa);
z8=0.5*(+sqrt(x1)+sqrt(x2)-sqrt(x3))-bb/(4.0*aa);


if(fabs((wielomian(z1)).real())<eps && fabs((wielomian(z1)).imag())<eps) cout<<z1<<endl;
if(fabs((wielomian(z2)).real())<eps && fabs((wielomian(z2)).imag())<eps) cout<<z2<<endl;
if(fabs((wielomian(z3)).real())<eps && fabs((wielomian(z3)).imag())<eps) cout<<z3<<endl;
if(fabs((wielomian(z4)).real())<eps && fabs((wielomian(z4)).imag())<eps) cout<<z4<<endl;
if(fabs((wielomian(z5)).real())<eps && fabs((wielomian(z5)).imag())<eps) cout<<z5<<endl;
if(fabs((wielomian(z6)).real())<eps && fabs((wielomian(z6)).imag())<eps) cout<<z6<<endl;
if(fabs((wielomian(z7)).real())<eps && fabs((wielomian(z7)).imag())<eps) cout<<z7<<endl;
if(fabs((wielomian(z8)).real())<eps && fabs((wielomian(z8)).imag())<eps) cout<<z8<<endl;


cout<<"

x1="<<z1<<endl;
cout<<"x2="<<z2<<endl;
cout<<"x3="<<z3<<endl;
cout<<"x4="<<z4<<endl;
cout<<"x5="<<z5<<endl;
cout<<"x6="<<z6<<endl;
cout<<"x7="<<z7<<endl;
cout<<"x8="<<z8<<endl;

cout<<wielomian(z1).real()<<endl;
cout<<wielomian(z2)<<endl;
cout<<wielomian(z3)<<endl;
cout<<wielomian(z4)<<endl;
cout<<wielomian(z5)<<endl;
cout<<wielomian(z6)<<endl;
cout<<wielomian(z7)<<endl;
cout<<wielomian(z8)<<endl;

 return 0;
} 

cmpx wielomian(cmpx x)
{w=aa*x*x*x*x+bb*x*x*x+cc*x*x+dd*x+ee;
return w;
}
void druk(cmpx x)
{
if(fabs(x.imag())<eps) printf("%2.12lf
",x.real());
else if(x.imag()>0) printf("%2.12lf + %2.12lf i
",x.real(),x.imag());
if(x.imag()<0 && fabs(x.real())>eps) printf("%2.12lf - %2.12lf i
",x.real(),-x.imag());
if(fabs(x.real())<eps) printf("%2.12lf i
",x.imag());

}

z tym ze ten drugi w niektorych przypadkach liczy zle (pierwiastki w st 4) czego do konca nie rozgryzlem...

guzik15 · Post autor: **guzik15** » 30 cze 2006, o 22:05

o w morde..., jestem pod wrażeniem

bisz · Post autor: **bisz** » 1 lip 2006, o 09:06

ja jestem pod wrazeniem jak nieudolnym programistą wtedy bylem ;]

kolumb · Post autor: **kolumb** » 1 lip 2006, o 11:22

ja wpisałem wzory Cardana i wyrzuca mi taki wyniki dla 1 x: -1.#IND
wielomian jest prosty powinien miec peirwiastki calkowite

oto kod

#include
#include
using namespace std;

main()
{

double long a=1;
double long b=3;
double long c=(-4);
double long d=(-12);

double long p;
p=((c/a)-((b*b)/(3*a*a))) ;
double long q;
q=(((2*b*b*b)/(27*a*a*a))+(d/a)-((b*c)/(3*a*a))) ;

double long w, u, k, s, e, o;
s=(((q*q)/4)+((p*p*p)/27));
w=(((-q)/2)+sqrt(s));
u=(((-q)/2)-sqrt(s));

e=pow(w,1.0/3.0); //pierwiastek 3 stopnia
o=pow(u,1.0/3.0);

k=e+o-(b/(3*a));

cout

Rogal · Post autor: **Rogal** » 1 lip 2006, o 13:08

Bo delta może być całkiem śmiało ujemna, bo wtedy są trzy piewiastki rzeczywiste.
Musisz to uwzględnić, rozbijając najlepiej na przypadki trzy: kiedy D > 0, D = 0 i D < 0.
Ponadto wypadałoby od razu zmodyfikować delikatnie wzory na pierwiastki zespolone, gdy Delta > 0, żeby część rzeczywistą mieć ładnie osobno, a osobno zespoloną.

kolumb · Post autor: **kolumb** » 1 lip 2006, o 14:17

Mowisz ze jak delta bedzie ujemna to sa 3 pierwiastki rzeczywiste. Ale jak je wtedy policzyc jak z tego jednego wzoru ktory wpisalem na x1 wychodzi mi juz bledny wynik w programie ?

[ Dodano: 1 Lipiec 2006, 14:21 ]
probowalem tez uruchomic programy bisza ale niezbyt u mnie dzialaja

Rogal · Post autor: **Rogal** » 1 lip 2006, o 22:54

No tak, musisz to przewidzieć i wprowadzić sobie i, czyli jednostkę urojoną (i=sqrt(-1)).
Gdy masz D < 0, to wtedy \(\displaystyle{ \sqrt{D} = i\sqrt{-D}}\)

[ Dodano: Sob Lip 01, 2006 10:55 pm ]
A dupa, korzysta się z postaci trygonometrycznej. Zajrzyj do Kompendium 2+2 -> wzory na pierwiastki równań wyższych stopni, tam jest wszystko wyjaśnione.

bisz · Post autor: **bisz** » 3 lip 2006, o 15:49

Herezje tu tworzycie, więc napisałem coś co u mnie działa, w razie problemów z kompilacją odpowiadam na pytania.

Kod: Zaznacz cały

#include <iostream>
#include <complex.h>
#include <math.h>
#define eps 0.00000001

using namespace std; 
typedef complex<double> cmpx;

ostream& operator<<(ostream& s, cmpx &z) 
{ 
if(fabs(z.imag())<eps)
{return s << z.real() ;}
if(fabs(z.real())<eps)
{return s << z.imag() <<" i" ;}

if(z.imag()<0)
{return s << z.real() <<" - " <<-z.imag() <<" i";}
if(z.imag()==0)
{return s << z.real() ; }
if(z.imag()>0 )
{return s << z.real() <<" + " << z.imag() <<" i"; }

}

void druk(cmpx z);

cmpx i(0,1.0);
cmpx x1,x2,x3;
double p,q,a, b, c, d;

main()
{
cout<<"Rownanie wyglada :
 ax^3+bx^2+cx+d=0
";
cout<<"Podaj a: ";
cin>>a;
cout<<"Podaj b: ";
cin>>b;
cout<<"Podaj c: ";
cin>>c;
cout<<"Podaj d: ";
cin>>d;
p=1.0*c/a-b*b/(3.0*a*a); 
q=2*b*b*b/(27.0*a*a*a)+1.0*d/a-b*c/(3.0*a*a);


x1=pow(0.0*i-q/2.0+sqrt(q*q/4.0+p*p*p/27.0+0.0*i),1/3.0)+pow(0.0*i-q/2.0-sqrt(q*q/4.0+p*p*p/27.0+0.0*i),1/3.0)-b/(3.0*a);
x2=0.5*(-1.0+sqrt(3)*1.0*i)*pow(0.0*i-q/2.0+sqrt(q*q/4.0+p*p*p/27.0+0.0*i),1/3.0)+0.5*(-1.0-sqrt(3)*1.0*i)*pow(0.0*i-q/2.0-sqrt(q*q/4.0+p*p*p/27.0+0.0*i),1/3.0)-b/(3.0*a);
x3=0.5*(-1.0-sqrt(3)*1.0*i)*pow(0.0*i-q/2.0+sqrt(q*q/4.0+p*p*p/27.0+0.0*i),1/3.0)+0.5*(-1.0+sqrt(3)*1.0*i)*pow(0.0*i-q/2.0-sqrt(q*q/4.0+p*p*p/27.0+0.0*i),1/3.0)-b/(3.0*a);


druk(x1);
druk(x2);
druk(x3);

}
void druk(cmpx x)
{
if(fabs(x.imag())<eps) printf("%2.12lf
",x.real());
else if(x.imag()>0) printf("%2.12lf + %2.12lf i
",x.real(),x.imag());
if(x.imag()<0 && fabs(x.real())>eps) printf("%2.12lf - %2.12lf i
",x.real(),-x.imag());
if(fabs(x.real())<eps) printf("%2.12lf i
",x.imag());

}

Rogal · Post autor: **Rogal** » 3 lip 2006, o 16:40

W C jest pierwiastkowanie sześcienne liczb zespolonych? Na przykład ile by to było w przybliżeniu \(\displaystyle{ \sqrt[3]{2+i}}\)?
A tak w ogóle, sprawdzałeś to na dużej liczbie przykładów?

bisz · Post autor: **bisz** » 3 lip 2006, o 21:29

w C - chyba ze napiszesz sobie typ zespolony na biednych strukturach i wymyslisz do niego funkcje, w C++ jest gotowa klasa complex z wszystkimi funkcjami z math.h przystosowanymi do typu zespolonego , przyklad , wstawic za main w poprzednim :

Kod: Zaznacz cały

main() 
{ 

x1=pow(2.0+1.0*i,1/3.0);
druk(x1); 


getch();
}

zwraca naturalnie pierwszy pierwiastek z wszystkich 3, dlaczego, chyba nie musze mowic

na ekranie ukazuje sie 1.292074512673103 + 0.20129431282890373 i

jest dobrze ?:]

wzory eulera itp TEZ działają... jeszcze jakies pytania ?

Rogal · Post autor: **Rogal** » 3 lip 2006, o 21:40

Niezłe to C++ : )