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.c

    Im 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.htm

    int
    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.png

    DOCH DANN KAM DIE RETTUNG:
    http://www.php-dir.de/torus4.png

    Beleuchtung 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:

    http://www.php-dir.de/ronja.model.sdf.fullHD.png



  • 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)


Log in to reply