Problem mit Program zum Lösen von Gleichungen bis 4. Grades
-
Guten Tag.
Ich wende mich nun an dieses Forum, da ich wirklich nciht mehr weiter weiß.
Mein Problem ist folgendes:Ich habe einen Raytracer, der bisher Würfel, dreiecke und Kugeln rendern kann.
Das macht er auch wirklich gut und schnell.
Nun wollte ich mich mal an eine etwas komplexere Form wagen: der Torus.Der Anfang ist auch ganz leicht nur, um den Schnittpunkt mit einer Geraden und einem Torus zu finden, muss man eine Gleichung der Form
ax^4 + bx^3 + cx^2 + dx + e = 0 lösen.
Die Koeffizienden a, b, c, d und e sind dabei bekann.Nun habe ich ein kleines in C geschriebenes Program gefunden, aus den Graphics Gems:
http://tog.acm.org/GraphicsGems/gems/Roots3And4.cIm wesentlichen geht es um die Funktion, die eine Gleichung 3. Grades Berechnet:
Die Zeile 44 ist die Interessante: wenn sqrtD < q, dann kann man auch nicht mehr die 3. Wurzel ziehen, obwohl eine reelle Lösung bei der Gleichung existiert.Gleichung bei der dieser Fehler z.B. auftritt: x3+4x2-3x+2
eigentliche reelle Nullstelle:
x1 = -4.724576727631849
siehe http://www.thkoehler.de/midnightblue/m_kdb.htmint rootFinder::solveDegree3(double c[4], double s[3]) { // Momentane Form: ax^3 + bx^2 + cx + d = 0 int i, num; double sub; double A, B, C; double sqA, p, q; double cbP, D; // Normalform: x^3 + Ax^2 + Bx + C = 0; A = c[ 2 ] / c[ 3 ]; B = c[ 1 ] / c[ 3 ]; C = c[ 0 ] / c[ 3 ]; // Substituieren mit x = y - A / 3 um quadratischen Term zu eliminieren // x^3 + px + q sqA = A * A; p = 1.0/3 * (-1.0/3 * sqA + B); q = 1.0/2 * (2.0/27 * A * sqA - 1.0/3 * A * B + C); // verwende Cardano`s Formeln cbP = p * p * p; D = q * q + cbP; if( isZero(D) ) { if( isZero(q) ) {// one tripel solution s[ 0 ] = 0; num = 1; } else { // one single and one double solution double u = cbrt(-q); s[ 0 ] = 2 * u; s[ 1 ] = -u; num = 2; } } else if( D < 0 ) { // casus irreducibilis: three real solutions double phi = 1.0/3 * acos( -q / sqrt( -cbP ) ); double t = 2 * sqrt( -p ); s[ 0 ] = t * cos(phi); s[ 1 ] = -t * cos(phi + M_PI / 3); s[ 2 ] = -t * cos(phi - M_PI / 3); num = 3; } else { // one real solution double sqrtD = sqrt(D); double u = cbrt( sqrtD - q ); double v = - cbrt( sqrtD + q ); s[ 0 ] = u + v; num = 1; // check auf NAN if( u != u ) { std::cout << "uncorrectable NaN detected: " << s[ 0 ] << " "; } } // Rücksubstituieren sub = 1.0 / 3 * A; for( i = 0; i < num; ++i ) s[ i ] -= sub; return num; }
Ich hoffe, jemand kann mir weiter helfen.
Wenn Ihr mehr code braucht, sagt bescheid.freundliche Grüße
devildeath
-
Hi,
devildeath0 schrieb:
Die Zeile 44 ist die Interessante: wenn sqrtD < q, dann kann man auch nicht mehr die 3. Wurzel ziehen, obwohl eine reelle Lösung bei der Gleichung existiert.
Doch, die dritte Wurzel kann gezogen werden. Aus deinem Sourcecode-Link:
#define cbrt(x) ((x) > 0.0 ? pow((double)(x), 1.0/3.0) : \ ((x) < 0.0 ? -pow((double)-(x), 1.0/3.0) : 0.0))
"pow" kann normalerweise keine Wurzel aus negativen Zahlen ziehen, deshalb der "
-pow(abs(x))
"-Workaround
-
DANKE!!!!!!!!!! genau das wars.
Ich hab in der Macro_Definition das Minus vor dem x vergessen...und da kamen Sachen raus kann ich nur sagen:
http://www.php-dir.de/torus.png
http://www.php-dir.de/torus1.png
http://www.php-dir.de/torus2.png
http://www.php-dir.de/torus3.pngDOCH DANN KAM DIE RETTUNG:
http://www.php-dir.de/torus4.pngBeleuchtung stimmt nur noch nciht ganz, doch das bekomm ich hin.
DANKEEEEE
und hier mal noch ein Bildchen, was er bisher bis vor dem Torus leistete...er wird noch erweitert:
-
Wenn Dich dieses Thema interessiert, schau auch mal bei:
http://www.povray.org/
Die gesamten Sourcecodes sind verfügbar, und man kann doch das eine oder andere herauslesen.
(Schön zu sehen, dass an diesem alten Dinosaurier immer noch dran gewerkelt wird)