Witam, mam do napisania program ktory ma obliczac wielomian 4 stopnia za pomocą wzorów Cardano do postaci 3 stopnia a potem za pomoca wzorów Ferrari do 2 stopnia!

Może ktoś ma jakieś progaramy, albo przykłady jak to "ugryść"?

Nigdy nie zajmowałem się takim problemem jednak bym zaczął od:

http://pl.wikipedia.org/w...wartego_stopnia
http://pl.wikipedia.org/w...zeciego_stopnia

W artykule dotyczącym 3 stopnia jest zamieszczony algorytm implementacji w języku C jeśli chodzi o czwarty poziom to jest żmudne i czasochłonne implementowanie wzorów matematycznych.

Powodzenia

Ja umieściłem implementacje wzorów Cardano w języku C na stronie
wikipedii Obecnie jest ona usunięta .
Ostatnio napisałem program obliczający pierwiastki
równania czwartego stopnia z C->C (współczynniki i pierwiastki zespolone)

Oto odnośnik do mojego programu
Zaimplementowane tam są dwie metody
rozwiązywanie równań czwartego stopnia

http://download.4programm...oad.php?id=1763

1. Metoda analogiczna do metody Niccolo Fontany
dla równań trzeciego stopnia

Stosujemy dwa podstawienia
Pierwsze wyeliminuje składnik a[3]*x^3
Drugie da nam wzory Viete'a równania
rozwiązującego

2. Metoda Ferrari'ego opisana poniżej


Wybrałem język C++ ponieważ ma on jako tako zaimplementowany typ zespolony
(chyba klasę)

Aby rozwiązać równanie trzeciego stopnia trzeba umieć rozwiązać równanie kwadratowe
Aby rozwiązać równanie kwadratowe należy
1.Przenieść wyraz wolny na drugą stronę
2.Sprowadzić lewą stronę do postaci kwadratu zupełnego
a) stosując podstawienie y=x+a[1]/(2a[2])
b) dodając stronami a[1]^2/(4a[2])
3.Dzielimy równanie przez a[2] oraz wyciągamy obustronnie
pierwiastek kwadratowy przy czym po prawej stronie równania
bierzemy oba pierwiastki

Równanie trzeciego stopnia wymaga dwóch podstawień i
rozwiązania jednego równania kwadratowego
Pierwsze podstawienie wyeliminuje składnik a[2]*x**2
Po tym podstawieniu oraz po podzieleniu równana przez a[3]
otrzymamy równanie postaci
y^3+3py+2q=0
Drugie podstawienie da nam układ dwóch równań nieliniowych
który to układ sprowadza się do rozwiązania równania kwadratowego
Układ ten to wzory Viete'a dla równania kwadratowego
zwanego równaniem rozwiązującym
Równanie rozwiązujące jest postaci
t^2+2qt-p^3=0

Następnie trzeba wyciągnąć pierwiastki arytmetyczne trzeciego stopnia
z pierwiastków równania kwadratowego ,dobrać odpowiednie pierwiastki
trzeciego stopnia z jedynki i dodać
Pierwiastki arytmetyczne trzeciego stopnia należy wyciągnąć tak
aby ich iloczyn był równy -p
Wzory Cardano to
y[1]=t[1]^(1/3)+t[2]^(1/3)
y[2]=e[1]*t[1]^(1/3)+e[2]*t[2]^(1/3)
y[3]=e[2]*t[1]^(1/3)+e[1]*t[2]^(1/3)

Pierwiastki trzeciego stopnia z jedynki to e[k]=exp(2*k*I*Pi/3) k należy do Z3

Co do pierwiastków równania czwartego stopnia to

Dzielimy równanie przez a[4]
Przenosimy trójmian kwadratowy na drugą stronę równania
Sprowadzamy lewą stronę równania do postaci kwadratu zupełnego
a) podstawiając np y=x+a[3]/(4*a[4])
b) dodając stronami (b[3]*x/2)**2, b[3]=a[3]/a[4]
Wprowadzamy nową zmienną (np u) tak aby lewa strona równania nadal była
kwadratem zupełnym
Prawa strona równania będzie kwadratem zupełnym gdy wyróżnik będzie
równy zero .Należy więc tak wyznaczyć zmienną (np u) aby wyróżnik był
równy zero. W tym celu trzeba rozwiązać równanie trzeciego stopnia względem
wprowadzonej zmiennej (np u).

Tyle wstępu teraz kod

#include<iostream.h>
#include<iomanip.h>
#include<complex.h>
#include<conio.h>

complex horner(int,complex[],complex);

int linear(complex[],complex[]);
int quadratic(complex[],complex[]);
int cubic(complex[],complex[]);
int quartic(complex[],complex[]);

const double eps=1e-12;

int main(){
int i,n;
complex a[5];
complex x[5];
char ch;
clrscr();
do{
for(i=4;i>=0;i--){
cout<<"a["<<i<<"]=";
cin>>a;
}
n=quartic(a,x);
for(i=1;i<=n;i++)
cout<<"x["<<i<<"]="<<setprecision(12)<<x<<endl;
ch=getch();
}
while(ch!=27);
return 0;
}

complex horner(int n,complex a[],complex x){
complex s=complex(0,0);
int i;
for(i=n;i>=0;i--)
s=s*x+a;
return s;
}

int linear(complex a[],complex x[]){
int n=0;
if (a[1]!=complex(0,0)){
x[1]=complex(-a[0]/a[1]);

n=1;
}
return n;
}

int quadratic(complex a[],complex x[]){
int n;
complex delta;
if (a[2]==complex(0,0)) return linear(a,x);
else{
delta=pow(a[1],2)-4.0*a[2]*a[0];
x[1]=complex((-a[1]-sqrt(delta))/(2.0*a[2]));
x[2]=complex((-a[1]+sqrt(delta))/(2.0*a[2]));

n=2;
return n;
}
}

int cubic(complex a[],complex x[]){

int j,n,k=0;

complex p,q;
complex b[3],t[3];

if(a[3]==complex(0.0,0.0)) return quadratic(a,x);
else{

for(j=1;j<=3;j++)
x[j]=complex(-a[2]/(3.0*a[3]));

p=complex((3.0*a[3]*a[1]-pow(a[2],2))/(9.0*pow(a[3],2)));
q=complex((2.0*pow(a[2],3)-9.0*a[3]*a[2]*a[1]+27.0*pow(a[3],2)*a[0])/(54.0*pow(a[3],3)));

b[2]=complex(1.0,0.0);
b[1]=2*q;
b[0]=-pow(p,3);

quadratic(b,t);

t[1]=pow(t[1],1./3.);
t[2]=pow(t[2],1./3.);

for(j=0;j<3&&k==0;j++){
if (abs(horner(3,a,exp(complex(0,2*j*M_PI/3))*t[1]+t[2]-complex(a[2]/(3*a[3]))))<abs(eps)) k=1;
}
t[1]*=exp(complex(0,2*(j-1)*M_PI/3));

x[1]+=(t[1]+t[2]);
x[2]+=(exp(complex(0.0,2.0*M_PI/3.0))*t[1]+exp(complex(0.0,-2.0*M_PI/3.0))*t[2]);
x[3]+=(exp(complex(0.0,-2.0*M_PI/3.0))*t[1]+exp(complex(0.0,2.0*M_PI/3.0))*t[2]);

n=3;
return n;
}

}

int quartic(complex a[],complex x[]){

int j,n;
complex p,q,r;
complex ua[4],ux[4];
complex a1[3],x1[3];
complex a2[3],x2[3];

if(a[4]==complex(0.0,0.0)) return cubic(a,x);
else{

p=complex((8*a[4]*a[2]-3*pow(a[3],2))/(8*pow(a[4],2)));
q=complex((pow(a[3],3)-4*a[4]*a[3]*a[2]+8*pow(a[4],2)*a[1])/(8*pow(a[4],3)));
r=complex((16*a[4]*pow(a[3],2)*a[2]+256*pow(a[4],3)*a[0]-3*pow(a[3],4)-64*pow(a[4],2)*a[3]*a[1])/(256*pow(a[4],4)));

for(j=1;j<=4;j++){
x[j]=complex(-a[3]/(4*a[4]));
}

ua[3]=complex(1.0,0.0);
ua[2]=-p;
ua[1]=-4*r;
ua[0]=4*p*r-pow(q,2);

cubic(ua,ux);


if(abs(ux[1]-p)<abs(eps)){

a1[2]=complex(1.0,0.0);
a1[1]=p;
a1[0]=r;

quadratic(a1,x1);

x[1]+=sqrt(x1[1]);
x[2]-=sqrt(x1[1]);
x[3]+=sqrt(x1[2]);
x[4]-=sqrt(x1[2]);

}
else{

a1[2]=complex(1.0,0.0);
a1[1]=-sqrt(ux[1]-p);
a1[0]=complex(q*sqrt(ux[1]-p)/(2*(ux[1]-p)))+complex(ux[1]/2);

a2[2]=complex(1.0,0.0);
a2[1]=sqrt(ux[1]-p);
a2[0]=complex(-q*sqrt(ux[1]-p)/(2*(ux[1]-p)))+complex(ux[1]/2);

quadratic(a1,x1);
quadratic(a2,x2);

x[1]+=x1[1];
x[2]+=x1[2];
x[3]+=x2[1];
x[4]+=x2[2];

}

n=4;
return n;
}
}



Skompilowane w dosowym środowisku

A teraz wersja dla Dev Cpp

#include<iostream.h>
#include<iomanip.h>
#include<complex.h>
#include<conio.h>

const double eps=1e-12;

complex<double> horner(int,complex<double>[],complex<double>);

int linear(complex<double>[],complex<double>[]);
int quadratic(complex<double>[],complex<double>[]);
int cubic(complex<double>[],complex<double>[]);
int quartic(complex<double>[],complex<double>[]);


int main(){
char ch;
int i,n;
complex<double> a[5],x[5];
do{
for(i=4;i>=0;i--){
cout<<"a["<<i<<"]=";
cin>>a;
}
n=quartic(a,x);
for(i=1;i<=n;i++)
cout<<"x["<<i<<"]="<<setprecision(12)<<x<<endl;
ch=getch();
}
while(ch!=27);
return 0;
}

complex<double> horner(int n,complex<double> a[],complex<double> x){
int i;
complex<double> s=complex<double>(0.0,0.0);
for(i=n;i>=0;i--)
s=s*x+a;
return s;
}

int linear(complex<double> a[],complex<double> x[]){
int n=0;
if(a[1]!=complex<double>(0.0,0.0)){
x[1]=complex<double>(-a[0]/a[1]);
n=1;
}
return n;
}

int quadratic(complex<double> a[],complex<double> x[]){
int i,n;
complex<double> delta;
if(a[2]==complex<double>(0.0,0.0)) return linear(a,x);
else{
delta=pow(a[1],2)-4.0*a[2]*a[0];
for(i=1;i<=2;i++)
x=complex<double>(-a[1]/(2.0*a[2]));

x[1]-=complex<double>(sqrt(delta)/(2.0*a[2]));
x[2]+=complex<double>(sqrt(delta)/(2.0*a[2]));

n=2;
return n;
}
}

int cubic(complex<double> a[],complex<double> x[]){
int i,j,n;
complex<double> p,q;
complex<double> b[3],t[3];
complex<double> e[3];

if (a[3]==complex<double>(0.0,0.0)) return quadratic(a,x);
else{
p=complex<double>((3.0*a[3]*a[1]-pow(a[2],2))/(9.0*pow(a[3],2)));
q=complex<double>((2.0*pow(a[2],3)-9.0*a[3]*a[2]*a[1]+27.0*pow(a[3],2)*a[0])/(54.0*pow(a[3],3)));

for(i=1;i<=3;i++)
x=complex<double>(-a[2]/(3.0*a[3]));

for(i=0;i<3;i++)
e=exp(complex<double>(0.0,2.0*i*M_PI/3.0));

b[2]=complex<double>(1.0,0.0);
b[1]=2.0*q;
b[0]=-pow(p,3);

quadratic(b,t);

t[1]=pow(t[1],1.0/3.0);
t[2]=pow(t[2],1.0/3.0);

j=0;
while(abs(e[j/3]*t[1]*e[j%3]*t[2]+p)>=eps)
j++;

t[1]*=e[j/3];
t[2]*=e[j%3];

x[1]+=(t[1]+t[2]);
x[2]+=(e[1]*t[1]+e[2]*t[2]);
x[3]+=(e[2]*t[1]+e[1]*t[2]);

n=3;
return n;
}
}

int quartic(complex<double> a[],complex<double> x[]){
int i,n;
complex<double> p,q,r;
complex<double> ua[4],ux[4];
complex<double> b1[3],x1[3];
complex<double> b2[3],x2[3];

if(a[4]==complex<double>(0.0,0.0)) return cubic(a,x);
else{

p=complex<double>((8.0*a[4]*a[2]-3.0*pow(a[3],2))/(8.0*pow(a[4],2)));
q=complex<double>((pow(a[3],3)-4.0*a[4]*a[3]*a[2]+8.0*pow(a[4],2)*a[1])/(8.0*pow(a[4],3)));
r=complex<double>((16.0*a[4]*pow(a[3],2)*a[2]+256.0*pow(a[4],3)*a[0]-3.0*pow(a[3],4)-64.0*pow(a[4],2)*a[3]*a[1])/(256.0*pow(a[4],4)));

for(i=1;i<=4;i++)
x=complex<double>(-a[3]/(4.0*a[4]));

ua[3]=complex<double>(1.0,0.0);
ua[2]=-p;
ua[1]=-4.0*r;
ua[0]=4.0*p*r-pow(q,2);

cubic(ua,ux);

if(abs(ux[1]-p)<eps){

b1[2]=complex<double>(1.0,0.0);
b1[1]=p;
b1[0]=r;

quadratic(b1,x1);

x[1]-=sqrt(x1[1]);
x[2]+=sqrt(x1[1]);
x[3]-=sqrt(x1[2]);
x[4]+=sqrt(x1[2]);

}
else{

b1[2]=complex<double>(1.0,0.0);
b1[1]=-sqrt(ux[1]-p);
b1[0]=complex<double>(ux[1]/2.0)+complex<double>(q*sqrt(ux[1]-p)/(2.0*(ux[1]-p)));

b2[2]=complex<double>(1.0,0.0);
b2[1]=sqrt(ux[1]-p);
b2[0]=complex<double>(ux[1]/2.0)-complex<double>(q*sqrt(ux[1]-p)/(2.0*(ux[1]-p)));

quadratic(b1,x1);
quadratic(b2,x2);

x[1]+=x1[1];
x[2]+=x1[2];
x[3]+=x2[1];
x[4]+=x2[2];

}

n=4;
return n;
}
}

[ Dodano: 9 Styczeń 2009, 09:32 ]
Powyższy kod źródłowy jest to kod aplikacji konsolowej
Napisałem także kod źródłowy tej aplikacji (pod okienka)
z użyciem C++ WinAPI oraz w Javie

[ Dodano: 31 Styczeń 2009, 22:29 ]
Dla niektórych danych może zwracac Quiet Not a Number wystarczy wtedy zmniejszyc dokładnośc

---

Quartics.zip Kod aplikacji w Javie

Pobierz Plik ściągnięto 103 raz(y) 92.64 KB

Quartic.zip Kod aplikacji w C++ z użyciem WinAPI

Pobierz Plik ściągnięto 118 raz(y) 3.5 KB