Studium, Ausbildung und Beruf
 StudiumHome   FAQFAQ   RegelnRegeln   SuchenSuchen    RegistrierenRegistrieren   LoginLogin

Umformen einer monströsen Gleichung
Gehe zu Seite Zurück  1, 2, 3, 4  Weiter
Neues Thema eröffnen   Neue Antwort erstellen
Foren-Übersicht -> Mathe-Forum -> Umformen einer monströsen Gleichung
 
Autor Nachricht
algebrafreak
Senior Member
Benutzer-Profile anzeigen
Senior Member


Anmeldungsdatum: 28.10.2004
Beiträge: 4143
Wohnort: Passau

BeitragVerfasst am: 15 März 2005 - 23:58:06    Titel:

Zitat:
Eine Bemerkung am Rande: Flash kann keine kubischen Kurven rendern - die Kurve wird durch einen Spline aus mehreren quadratischen Kurven emuliert (mathematisch und algebraisch sicherlich völlig unhaltbar )


Nö, nö. Wie will man das sonst machen? Irgendwo liegt sowieso eine Methode, die alle Punkte nacheinander durchgeht und zeichnet. Diese kann man auch als Splines ansehen. Solange es noch keine kontinuierliche Ausgabegeräte gibt, wird das wohl auch so bleiben. Der Einwand von vorhin bezog sich auf die Komplexität und die Stabilität des Verfahrens.

Zitat:
Ich würde diesen Port ins ActionScript gerne in dieser Form der 'Flash Community' zur Verfügung stellen (natürlich unentgeltlich) - ist das ok für Dich?


Ja, klar. Der Einwand mit GPL sollte Leute davor abschrecken, den Code zu kopieren und in ihre Programme ohne nachzudenken zu pasten. Es gibt knifflige Argumente, warum der Code nachgebessert werden sollte, solange man nicht mit einer algebraischen Rational-Zahlen-Bibliothek arbeitet. Ich denke, die Fragen weiter unten beziehen sich unter anderem auf diese Eigenheiten. Ich verstecke meinen Namen nicht Smile , ich heiße Aless Lasaruk.

Zitat:
Frage: lassen sich diese Werte 'irgendwie' - (Einbezug der relativen Position von P2 und P3) so umrechnen, dass t mit der Länge der Kurve korreliert? (Im Ergebnis: t ist 0, 0.1, 0.2 etc.)


Ich glaube, ich bin der falsche Mann um über den Begriff "Länge" zu diskutieren (siehe meine andere Postings Smile ) Da ich aber Informatiker bin, muß ich auch mal das ganze aus der praktischen Sicht sehen. Ja, klar. Das geht natürlich. Dazu muß man die Kurve integrieren. Als Nebeneffekt kommt da auch die Länge raus. Ich bin mir nicht sicher, ob nicht dasselbe, wie t, dabei rauskommt, aber ich mache mal das morgen.

Zitat:
Wenn das möglich wäre liesse sich für manche (alle?) 'triviale Fälle' die jetzt -1 zurückgeben, ein effektiver Wert zurückgeben - z.B. obiges Beispiel:


Der Wert -1 hate ursprünglich eine andere Bedeutung. Der sollte semantisch ausdrücken: "der Punkt liegt nicht auf der kurve". Ich habe mich allerdings in ein paar Annahmen geirrt. Dabei ist das Programm "schlauer" als ich und gibt bis auf ein paar kleine Veränderungen doch die Lösungen raus, auch für degenerierte Fälle.

Ich habe, da ich die Kurven zum ersten mal im Leben sehe, natürlich, wie auch sonst immer, fälschlicherweise angenommen, daß die Funktion B(u) injektiv (injektiv bedeutet zu jedem Punkt gehört ein Wert von t) ist. Was weiß ich, wo ich das her habe. Nun, wie man anhand des Flash-Programms sieht, ist dem nicht so. D.h. für bestimmte Kurven (Schleifen) und Punkte gibt es zwei Antworten.

Zitat:
Wenn ich noch einen Wunsch frei hätte (rein hypothetisch nätürlich...) die gleiche Funktion für quadratische Bezier Kurven (und es wäre natürlich versprochen, dass ich dann nicht noch mit den quintischen und so weiter nachkomme...).


Das ist kein Problem. Das mit 5 Kontrollpunkte usw. wäre eins, denn da wird es richtig mathematisch.

Das größte Problem, mit dem man hier zu kämpfen hat, ist das numerische Verhalten. In Wirklichkeit, und das wird Dir vermutlich nicht unbekannt sein, muß man Rundungsfehler berücksichtigen. Vor allem ist ein Test auf Gleichheit x == 0 problemmatisch, denn dieser ist normal syntaktisch. Dies schlägt sich deutlich im obigen Programm nieder

Code:

...
         for (var j=0, len2=yl.length; j<len2 ; j++){
            if (Math.abs(xl[i] - yl[j]) < 0.01){
                 ^ Hier sollte eigentlich xl[i] == yl[i] stehen.
               return xl[i]; 
            }
         }
...


Ersetzt man das so, wie es "sein sollte", kommt in etwa 30% der Fälle Mist raus, weil die Werte nicht exakt gleich sind. Der Wert 0.01 beschreibt dabei die "Ungenauigkeit", mit der man einen Punkt, der eigentlich nicht auf der Kurve liegt doch als einen Punkt auf der Kurve betrachtet.

Ich poste morgen mal den Code mit allen degenerierten Fällen, die Länge und den quadratischen Fall. Die Kommentare sollten nochmal überarbeitet werden.
Andreas Weber
Gast






BeitragVerfasst am: 16 März 2005 - 16:03:49    Titel:

Hallo Aless

Herzlichen Dank fürs 'Dranbleiben'!

Da ein kleines visuelles Modell das Verständnis doch manchmal erleichert hier noch die quadratische Kurve:

http://tinyurl.com/5uq5p

Im Gegensatz zur kubischen kann sich die quadratische nicht selber schneiden - mit einer Ausnahme:
wenn alle Punkte auf einer Linie liegen und wenn zusätzlich der Kontrollpunkt nicht auf dem Mittelpunkt der 'Basislinie' zwischen den beiden Ankern liegt - es überschneiden sich dann unendlich viele Punkte.

Auch typisch für die quadratischen Kurven ist was ich ganz untechnisch 'Ueberschiessen' nennen möchte. Man sieht das besonders gut wenn einer der Kontrollpunkte nahe der Basislinie liegt und nicht in der Mitte sondern nahe bei einem der Anker liegen: Die Kurve schiesst über einen der Anker hinaus, die von der Kurve gebildete Linie ist länger als die Basislinie.

Und hier noch einmal der Link zur Seite von Paul Bourke mit der mathematischen Beschreibung der Bézier Kurven:

http://astronomy.swin.edu.au/~pbourke/curves/bezier/

Andreas Weber
algebrafreak
Senior Member
Benutzer-Profile anzeigen
Senior Member


Anmeldungsdatum: 28.10.2004
Beiträge: 4143
Wohnort: Passau

BeitragVerfasst am: 16 März 2005 - 23:40:11    Titel:

Es dauert noch ein wenig. Ich bin so ein fauler Sack, eigentlich sollte ich an meiner Diplomarbeit arbeiten ...
Andreas Weber
Gast






BeitragVerfasst am: 17 März 2005 - 11:53:33    Titel:

Hallo Aless

Ich bin beim Arbeiten mit der Funktion auf ein Problem gestossen, dass sich beim 'Integrieren' der Kurve vielleicht ohnehin erledigt, auf das ich Dich aber doch aufmerksam machen möchte:
t kann für beide Punkte, durch die die Kurve führt (p2 und p3), nicht ermittelt werden.

Beispiel:

Wir haben 4 beliebige Punkte, die x oder y Koordinaten liegen nicht auf einer Linie, sind also kein 'trivialer Fall'
p1 = {x:0, y:300};
p2 = {x:50, y:130};
p3 = {x:70, y:55};
p4 = {x:250, y:200};

Das Problem taucht auf, wenn der Suchpunkt p0 einer der beiden Punkte ist, durch die die Kurve führt:
p0 = p2 oder p0 = p3;

Aufruf der Funktion getTOfPointOnCubicCurve(p0, p1, p2, p3, p4);

Resultat ist in beiden Fällen t == -1

Wie wär's wenn Du Deine Diplomarbeit über Bézier-Kurven schreiben würdest? Wink

Andreas Weber
algebrafreak
Senior Member
Benutzer-Profile anzeigen
Senior Member


Anmeldungsdatum: 28.10.2004
Beiträge: 4143
Wohnort: Passau

BeitragVerfasst am: 17 März 2005 - 13:31:45    Titel:

Das hat zunächst für Verwirrung bei mir gesorgt. Nach einigem hin und her ist die Antwort: Die Kurve geht im Allgemeinen durch keinen der Kontrollpunkte, außer p0 und pn. Die müssen also nicht auf der Kurve liegen. Dies ist auch im obigen Beispiel der Fall.

Zitat:
Wie wär's wenn Du Deine Diplomarbeit über Bézier-Kurven schreiben würdest?


Neee. Das wäre doch etwas, was man gebrauchen könnte. Ich bin ein Informatiker, der sich nur mit Sachen befaßt, die kein Mensch je verwenden wird Smile Ich schreibe über uniforme Quantorenelimination in Presburger-Arithmetik. Das Verfahren ist inhärent 3-fach exponentiell und die Ausgabe schon bei einfachen Beispielen so kompliziert, daß ich deren Korrektheit mit Papier und Stift kaum überprüfen kann. Ich sitze (schlafe besser gesagt) seit 3 Wochen an einem Beweis über 3 Seiten, in dem ich sändig Fehler entdecke... Naja.
algebrafreak
Senior Member
Benutzer-Profile anzeigen
Senior Member


Anmeldungsdatum: 28.10.2004
Beiträge: 4143
Wohnort: Passau

BeitragVerfasst am: 17 März 2005 - 15:09:45    Titel:

Eine gute Nachricht und eine schlechte. Die gute ist, daß ich denke, die degenerierten Fälle behandelt zu haben und den Code auch an quadratische Fälle angepasst zu haben. Hier der Code

Code:

// Der Code ist ein Prototyp und sollte auf keinen Fall so wie er ist,
// schon allein wegen der Rechte und GPL (erstellt in GPL Emacs) in
// komerzieller Software verwendet werden. Als weiteren Grund ist
// zu nennen, daß hier keinerlei Gedanken bzgl. der Numerischen
// Stabilität verschwendet wurden, was auch als Grund dienen kann :)

#include <vector>
#include <iostream>
#include <cmath>

class v2D {

  // Zweidimensionale vektoren

public:
 
  // Es ist dämlich auf koordinaten über accessor methoden zuzugreifen
  double x,y;

  v2D::v2D(double nx,double ny) {
    x = nx;
    y = ny;
  };

};

// Kubikwurzel
double crt(double x) {

    // Symmetrie verwenden
    if (x < 0) return -pow(-x,1.0/3.0);
    return pow(x,1.0/3.0);

}

// Gibt eine liste der reellen Nullstellen eines kubischen Polynoms zurück,
// falls nicht (alle Koeffizienten 0) gilt. Sonst liefert der Algorithmus
// eine Liste mit einem Element NAN mit der semantik R
std::vector<double> solve(double a,double b,double c,double d) {

    //std::cout << a << "," << b << "," << c << "," << d << std::endl;

  std::vector<double> r;

  // Triviale Fälle
 
  // Das Nullpolynom
  if (a == 0 && b == 0 && c == 0 && d == 0) {
      r.push_back(NAN);
      return r;
  }
  // Lineares Polynom c x + d
  if (a == 0 && b == 0 && c != 0) {
    r.push_back(-d/c);
    return r;
  }
  // Quadratisches Polynom
  if (a == 0 && b == 0) {
    double diskr = c*c-4.0*b*d;
    // Keine Nullstellen
    if (diskr < 0) return r;
    // Eine Nullstelle
    if (diskr == 0) {
      r.push_back(-c/(2.0*b));
      return r;
    };
    // Zwei Nullstellen
    r.push_back((-c + sqrt(diskr))/(2.0*b));
    r.push_back((-c - sqrt(diskr))/(2.0*b));
    return r;   
  }

  // Cardansche Formeln
  double p = (3.0*a*c - b*b)/(3.0*a*a);
  double q = (2.0*b*b*b-9.0*a*b*c+27.0*a*a*d)/(27.0*a*a*a);

  // Gelöst wird nun y^3 + py + q = 0

  double diskr = q*q/4.0+p*p*p/27.0;
 
  if (diskr > 0) {
    // Eine reelle Lösung
      double x = crt(-q/2.0+sqrt(diskr))+
     crt(-q/2.0-sqrt(diskr))-b/(3.0*a);
      r.push_back(x);
    return r;
  }

  if (diskr == 0) {
    // Zwei reelle Lösungen
    double x1 = crt(q/2.0)-b/(3.0*a);
    double x2 = -crt(4.0*q)-b/(3.0*a);
    r.push_back(x1);
    r.push_back(x2);
    return r;
  }

  // 3 reelle Lösungen

  double pi = 3.1415927;
 
  double x1 = 2.0*sqrt(-p/3.0)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p))))
    - b/(3.0*a);
  double x2 = -2.0*sqrt(-p/3.0)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p)))+
              pi/3.0) - b/(3*a);
  double x3 = -2.0*sqrt(-p/3.0)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p)))-
              pi/3.0) - b/(3.0*a);
  r.push_back(x1);
  r.push_back(x2);
  r.push_back(x3);
  return r;
}

// B(u) = p0 * (1-u)^3 + 3 * p1 * u * (1-u)^2 + 3 * p2 u^2 * (1-u) + p3 * u^3
v2D bezier(double u,v2D& p0,v2D& p1,v2D& p2, v2D& p3) {

  double v = 1-u;

  double x = p0.x*v*v*v+p1.x*3*u*v*v+p2.x*3*u*u*v+p3.x*u*u*u;
  double y = p0.y*v*v*v+p1.y*3*u*v*v+p2.y*3*u*u*v+p3.y*u*u*u;

  return(v2D(x,y));
}

// B(u) = p0 * (1-u)^2 + 2 * p1 * u * (1-u) + p2 u^2
v2D bezier(double u,v2D& p0,v2D& p1,v2D& p2) {

  double v = 1-u;

  double x = p0.x*v*v+p1.x*2*u*v+p2.x*u*u;
  double y = p0.y*v*v+p1.y*2*u*v+p2.y*u*u;

  return(v2D(x,y));
}

std::vector<double> reizeb(v2D& p,v2D& p0,v2D& p1,v2D& p2, v2D& p3) {
 
  // Berechnen von koeffizienten
  double ax = p3.x-3.0*p2.x+3.0*p1.x-p0.x;
  double bx = 3.0*p2.x-6.0*p1.x+3.0*p0.x;
  double cx = 3.0*p1.x-3.0*p0.x;
  double dx = p0.x - p.x;

  double ay = p3.y-3.0*p2.y+3.0*p1.y-p0.y;
  double by = 3.0*p2.y-6.0*p1.y+3.0*p0.y;
  double cy = 3.0*p1.y-3.0*p0.y;
  double dy = p0.y - p.y;

  std::vector<double> xl = solve(ax,bx,cx,dx);
  std::vector<double> yl = solve(ay,by,cy,dy);
 
  std::vector<double> r;

  /*std::cout << "X[i] = ";
  for (int i = 0; i < xl.size(); i++) std::cout << xl[i] << ",";
  std::cout << std::endl;
  std::cout << "Y[i] = ";
  for (int i = 0; i < yl.size(); i++) std::cout << yl[i] << ",";
  std::cout << std::endl;*/

  // Degenerierte Fälle
  if (xl.size() == 1 && xl[0] == NAN)
      // Erfüllungsmenge der ersten Gleichung sind die ganzen reellen Zahlen
      return yl;
  // Degenerierte Fälle
  if (yl.size() == 1 && yl[0] == NAN)
      // Erfüllungsmenge der zweiten Gleichung sind die ganzen reellen Zahlen
      return xl;

  for (int i = 0; i < xl.size(); i++)
      for (int j = 0; j < yl.size(); j++) {
     double diff = xl[i] - yl[j];
     if (fabs(diff) < 0.01) r.push_back(xl[i]);
      }

  return r;
}

std::vector<double> reizeb(v2D& p,v2D& p0,v2D& p1,v2D& p2) {
 
    // Berechnen von koeffizienten
    double ax = 0;
    double bx = p0.x-2*p1.x+p2.x;
    double cx = 2*p1.x-2*p0.x;
    double dx = p0.x;

    double ay = 0;
    double by = p0.y-2*p1.y+p2.y;
    double cy = 2*p1.y-2*p0.y;
    double dy = p0.y;

    std::vector<double> xl = solve(ax,bx,cx,dx);
    std::vector<double> yl = solve(ay,by,cy,dy);

    std::vector<double> r;

    /*for (int i = 0; i < xl.size(); i++) std::cout << xl[i] << ",";
    std::cout << std::endl;
    for (int i = 0; i < yl.size(); i++) std::cout << yl[i] << ",";
    std::cout << std::endl;*/

    // Degenerierte Fälle
    if (xl.size() == 1 && xl[0] == NAN)
   // Erfüllungsmenge der ersten Gleichung sind die ganzen reellen Zahlen
   return yl;
    // Degenerierte Fälle
    if (yl.size() == 1 && yl[0] == NAN)
   // Erfüllungsmenge der zweiten Gleichung sind die ganzen reellen Zahlen
   return xl;

    for (int i = 0; i < xl.size(); i++)
   for (int j = 0; j < yl.size(); j++)
       if (fabs(xl[i] - yl[j]) < 0.01)
      r.push_back(xl[i]);
   
    return r;
}

int main() {
   
    v2D p0(0.0,300.0);
    v2D p1(50.0,130.0); 
    v2D p2(70.0,55.0);
    v2D p3(250.0,200.0);
   
    v2D p(50.0,130.0);
   
    std::vector<double> r = reizeb(p,p0,p1,p2,p3);
   

    std::cout << "Result: " << std::endl;
    for (int i = 0; i < r.size(); i++)
   std::cout << r[i] << std::endl;
 
    double t = 0;
   
    while (t <= 1) {
   
   v2D v = bezier(t,p0,p1,p2,p3);
   std::vector<double> r = reizeb(v,p0,p1,p2,p3);
   
   //std::cout << "Testwert : " << t << " Antwort : " << r << std::endl;
   
   if (r.size() != 0 && fabs(r[0] - t) > 0.01)
       std::cout << "Fuck " << std::endl;
   
   t = t + 0.00001;
     
   }
   
}


Schlechte Nachricht: Was die Länge der Kurve anbetrifft, ist es schwieriger als ich dachte. Maple wird mit der Integration der Formel nicht fertig. Es gibt mehrere Lösungen, die allerdings alle recht problematisch sind:

(i) Einbinden bzw. Umformulieren einer numerischen Integrationsbibliothek. Wäre am einfachsten, wenn Du C++ benutzen würdest. Bei ActionScript weiß ich nicht ob es so eine gibt. Die Bibliotheken sind meistens hochoptimiert und relativ lang. Müsste man suchen.

(ii) Explizite Approximation der Formel durch ein Polynom und anschließende Integration. Es besteht die Möglichkeit eine Integrationsbibliothek zu "simulieren", indem man die zu integrierende Funktion durch ein Polynom approximiert und dann integriert. Problem ist: man bekommt möglicherweise sehr näherungsweise Werte.

Eure Lösung zu behalten erscheint mir zunächst einmal am aller einfachsten. Ich habe hier

http://www.uni-protokolle.de/foren/viewt/17931,0.html

die Frage nach der Längenbestimmung für alle gepostet. Vielleicht ist jemand, der schlauer als ich ist, in der Lage das schön zu lösen.

Am sonsten können wir über die Möglichkeiten diskutieren.
Andreas Weber
Gast






BeitragVerfasst am: 17 März 2005 - 16:02:54    Titel:

Hallo Aless

Herzlichen Dank für die Ueberarbeitung der Funktion!

Ich werd sie nach ActionScript umschreiben und dann ausgiebig testen - ich melde mich allenfalls erst morgen wieder.

Betreffend Längenproblem:
Danke für das saubere Formieren des Problems im anderen Thread!

Zitat:
(i) Einbinden bzw. Umformulieren einer numerischen Integrationsbibliothek. Wäre am einfachsten, wenn Du C++ benutzen würdest. Bei ActionScript weiß ich nicht ob es so eine gibt. Die Bibliotheken sind meistens hochoptimiert und relativ lang. Müsste man suchen.


Das gibt es in ActionScript ziemlich sicher noch nicht - also allenfalls Umschreiben? Was mich etwas skeptisch macht, ist, dass Du schreibst die Bibliotheken seinen relativ lang - wenn das heisst, dass die Funktion zur Bestimmung der Länge sehr viele Operationen ausführen muss, wird das in ActionScript - das _viel_ langsamer ist als z.B. C++ - zu langsam. Für den Moment werde ich auf alle Fälle den von Dir eröffneten Thread aufmerksam mitverfolgen - wer weiss?

Mein Problem vom letzten Post, betreffend unauffindbarkeit von P2, P3:
Zitat:
Das hat zunächst für Verwirrung bei mir gesorgt.


Sorry - meinerseits war die Verwirrtheit sehr gross und sehr hartnäckig - aus meinem Kommentar zur Verwendung der Funktion:
Zitat:
// - p2 and p3 are points _on_ the curve, _not_ the control points


Richtig ist - wie fast immer - das Gegenteil. Sorry!

Bevor Du den Schwung, in den Du mit diesen Kurven gekommen bist, wieder verlierst - und bevor sie allenfalls in Vergessenheit gerät: die arme quadratische harrt noch ihrer mathematischen Erlösung...

Very Happy

Herzlichen Dank!

Andreas Weber
algebrafreak
Senior Member
Benutzer-Profile anzeigen
Senior Member


Anmeldungsdatum: 28.10.2004
Beiträge: 4143
Wohnort: Passau

BeitragVerfasst am: 17 März 2005 - 16:19:57    Titel:

Der quadratische Fall ist ja schon im Code oben implementiert. Dazu werden die bezier und die reizeb Methoden überladen. Mathematisch gesehen ist die Umformung noch einfacher, als beim ersten Fall

B(u) = x0 (1-u)^2 + 2 x1 (1-u) u + x2 u^2

Ausmultiplizieren ergibt

u^2 (x0-2*x1+x2)+
u (2*x1-2*x0)
+ p0.x;

Die quadratische Gleichung ist als Speziallfall einer kubischen in der Methode solve enthalten.

Für alle Fälle: Du musst aufpassen, wenn die Methoden zwei Werte zurückliefern. Ich bin mir nicht sicher, welche Folgerungen speziell das Ganze hat, aber da ist auf jeden Fall Vorsicht geboten.

Der Wert 0.01 sollte irgendwo als Parameter definiert sein (vermutlich abgeleitet von der Qualität der Kurve), denn für verschiedene Auflösungen und Rasterungen der Kurve kann es notwendig sein gröbere Werte bzw. feinere Werte zu nehmen.

Eine numerische Integrations-Bibliothek wird vermutlich mit irgend einer Form von Interpolation arbeiten (z.B. Newton-Cotes-Formeln als stückweise Polynominterpolation) und benötigt meistens viele Schritte, vor allem aber solche, die die Fließkommaberechnungen involvieren. Es ist daher abzuwägen, ob die Addition von 100 Werten nicht besser ist, als ein Aufruf einer solchen Bibliothek.

Wie gesagt, man könnte eine bestimmte Approximation bzw. Interpolationsart vorberechnen und dann nur dort Einsetzen. Ich müsste mir aber dafür deutlich mehr Zeit nehmen, wie für dieses Beispiel, da die Theorie dazu sauber hergeleitet werden müsste.
Andreas Weber
Gast






BeitragVerfasst am: 17 März 2005 - 23:45:56    Titel:

Hallo Aless

Bei der quadratischen Kurve steh ich noch irgendwo auf der Leitung...

Code:
std::vector<double> reizeb(v2D& p,v2D& p0,v2D& p1,v2D& p2) {
 
    // Berechnen von koeffizienten
    double ax = 0;
    double bx = p0.x-2*p1.x+p2.x;
    double cx = 2*p1.x-2*p0.x;
    double dx = p0.x;

    double ay = 0;
    double by = p0.y-2*p1.y+p2.y;
    double cy = 2*p1.y-2*p0.y;
    double dy = p0.y;


Kannst Du rasch bestätigen, dass das so stimmt (ich vermisse eigentlich den Einbezug des Suchpunkts 'p')?
Falls ja, muss ich bei der Uebersetzung ins ActionScript noch irgendwo eine Unebenheit haben...
Danke für ein kurzes Feedback!

Andreas Weber
algebrafreak
Senior Member
Benutzer-Profile anzeigen
Senior Member


Anmeldungsdatum: 28.10.2004
Beiträge: 4143
Wohnort: Passau

BeitragVerfasst am: 18 März 2005 - 00:04:53    Titel:

Ich muss was vor'm Pasten durcheinandergebracht haben. Sorry.
Ich habe einen interessanten Ansatz zur Längenapproximation gefunden und schreibe gerade dran. Das dauert aber ein wenig...

Da ist der "korrekte" Code. Beachte in der solve Methode eine Änderung.
El stupido grande Smile Ich bin irgendwie zu garnix zu gebrauchen.

Code:

// Der Code ist ein Prototyp und sollte auf keinen Fall so wie er ist,
// schon allein wegen der Rechte und GPL (erstellt in GPL Emacs) in
// komerzieller Software verwendet werden. Als weiteren Grund ist
// zu nennen, daß hier keinerlei Gedanken bzgl. der Numerischen
// Stabilität verschwendet wurden, was auch als Grund dienen kann :)

#include <vector>
#include <iostream>
#include <cmath>

class v2D {

  // Zweidimensionale vektoren

public:
 
  // Es ist dämlich auf koordinaten über accessor methoden zuzugreifen
  double x,y;

  v2D::v2D(double nx,double ny) {
    x = nx;
    y = ny;
  };

};

// Kubikwurzel
double crt(double x) {

    // Symmetrie verwenden
    if (x < 0) return -pow(-x,1.0/3.0);
    return pow(x,1.0/3.0);

}

// Gibt eine liste der reellen Nullstellen eines kubischen Polynoms zurück,
// falls nicht (alle Koeffizienten 0) gilt. Sonst liefert der Algorithmus
// eine Liste mit einem Element NAN mit der semantik R
std::vector<double> solve(double a,double b,double c,double d) {

    //std::cout << a << "," << b << "," << c << "," << d << std::endl;

  std::vector<double> r;

  // Triviale Fälle
 
  // Das Nullpolynom
  if (a == 0 && b == 0 && c == 0 && d == 0) {
      r.push_back(NAN);
      return r;
  }
  // Lineares Polynom c x + d
  if (a == 0 && b == 0 && c != 0) {
    r.push_back(-d/c);
    return r;
  }
  // Quadratisches Polynom
  if (a == 0) {
   // ^ Hier -------------------------------------------------- watch out
    double diskr = c*c-4.0*b*d;
    // Keine Nullstellen
    if (diskr < 0) return r;
    // Eine Nullstelle
    if (diskr == 0) {
      r.push_back(-c/(2.0*b));
      return r;
    };
    // Zwei Nullstellen
    r.push_back((-c + sqrt(diskr))/(2.0*b));
    r.push_back((-c - sqrt(diskr))/(2.0*b));
    return r;   
  }

  // Cardansche Formeln
  double p = (3.0*a*c - b*b)/(3.0*a*a);
  double q = (2.0*b*b*b-9.0*a*b*c+27.0*a*a*d)/(27.0*a*a*a);

  // Gelöst wird nun y^3 + py + q = 0

  double diskr = q*q/4.0+p*p*p/27.0;
 
  if (diskr > 0) {
    // Eine reelle Lösung
      double x = crt(-q/2.0+sqrt(diskr))+
     crt(-q/2.0-sqrt(diskr))-b/(3.0*a);
      r.push_back(x);
    return r;
  }

  if (diskr == 0) {
    // Zwei reelle Lösungen
    double x1 = crt(q/2.0)-b/(3.0*a);
    double x2 = -crt(4.0*q)-b/(3.0*a);
    r.push_back(x1);
    r.push_back(x2);
    return r;
  }

  // 3 reelle Lösungen

  double pi = 3.1415927;
 
  double x1 = 2.0*sqrt(-p/3.0)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p))))
    - b/(3.0*a);
  double x2 = -2.0*sqrt(-p/3.0)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p)))+
              pi/3.0) - b/(3*a);
  double x3 = -2.0*sqrt(-p/3.0)*cos(1.0/3.0*acos(-q/2.0*sqrt(-27.0/(p*p*p)))-
              pi/3.0) - b/(3.0*a);
  r.push_back(x1);
  r.push_back(x2);
  r.push_back(x3);
  return r;
}

// B(u) = p0 * (1-u)^3 + 3 * p1 * u * (1-u)^2 + 3 * p2 u^2 * (1-u) + p3 * u^3
v2D bezier(double u,v2D& p0,v2D& p1,v2D& p2, v2D& p3) {

  double v = 1-u;

  double x = p0.x*v*v*v+p1.x*3*u*v*v+p2.x*3*u*u*v+p3.x*u*u*u;
  double y = p0.y*v*v*v+p1.y*3*u*v*v+p2.y*3*u*u*v+p3.y*u*u*u;

  return(v2D(x,y));
}

// B(u) = p0 * (1-u)^2 + 2 * p1 * u * (1-u) + p2 u^2
v2D bezier(double u,v2D& p0,v2D& p1,v2D& p2) {

  double v = 1-u;

  double x = p0.x*v*v+p1.x*2*u*v+p2.x*u*u;
  double y = p0.y*v*v+p1.y*2*u*v+p2.y*u*u;

  return(v2D(x,y));
}

std::vector<double> reizeb(v2D& p,v2D& p0,v2D& p1,v2D& p2, v2D& p3) {
 
  // Berechnen von koeffizienten
  double ax = p3.x-3.0*p2.x+3.0*p1.x-p0.x;
  double bx = 3.0*p2.x-6.0*p1.x+3.0*p0.x;
  double cx = 3.0*p1.x-3.0*p0.x;
  double dx = p0.x - p.x;

  double ay = p3.y-3.0*p2.y+3.0*p1.y-p0.y;
  double by = 3.0*p2.y-6.0*p1.y+3.0*p0.y;
  double cy = 3.0*p1.y-3.0*p0.y;
  double dy = p0.y - p.y;

  std::vector<double> xl = solve(ax,bx,cx,dx);
  std::vector<double> yl = solve(ay,by,cy,dy);
 
  std::vector<double> r;

  /*std::cout << "X[i] = ";
  for (int i = 0; i < xl.size(); i++) std::cout << xl[i] << ",";
  std::cout << std::endl;
  std::cout << "Y[i] = ";
  for (int i = 0; i < yl.size(); i++) std::cout << yl[i] << ",";
  std::cout << std::endl;*/

  // Degenerierte Fälle
  if (xl.size() == 1 && xl[0] == NAN)
      // Erfüllungsmenge der ersten Gleichung sind die ganzen reellen Zahlen
      return yl;
  // Degenerierte Fälle
  if (yl.size() == 1 && yl[0] == NAN)
      // Erfüllungsmenge der zweiten Gleichung sind die ganzen reellen Zahlen
      return xl;

  for (int i = 0; i < xl.size(); i++)
      for (int j = 0; j < yl.size(); j++) {
     double diff = xl[i] - yl[j];
     if (fabs(diff) < 0.01) r.push_back(xl[i]);
      }

  return r;
}

std::vector<double> reizeb(v2D& p,v2D& p0,v2D& p1,v2D& p2) {
 
    // Berechnen von koeffizienten
    double ax = 0;
    double bx = p0.x-2*p1.x+p2.x;
    double cx = 2*p1.x-2*p0.x;
    double dx = p0.x-p.x;

    double ay = 0;
    double by = p0.y-2*p1.y+p2.y;
    double cy = 2*p1.y-2*p0.y;
    double dy = p0.y-p.y;

    std::vector<double> xl = solve(ax,bx,cx,dx);
    std::vector<double> yl = solve(ay,by,cy,dy);

    std::vector<double> r;

    /*for (int i = 0; i < xl.size(); i++) std::cout << xl[i] << ",";
    std::cout << std::endl;
    for (int i = 0; i < yl.size(); i++) std::cout << yl[i] << ",";
    std::cout << std::endl;*/

    // Degenerierte Fälle
    if (xl.size() == 1 && xl[0] == NAN)
   // Erfüllungsmenge der ersten Gleichung sind die ganzen reellen Zahlen
   return yl;
    // Degenerierte Fälle
    if (yl.size() == 1 && yl[0] == NAN)
   // Erfüllungsmenge der zweiten Gleichung sind die ganzen reellen Zahlen
   return xl;

    for (int i = 0; i < xl.size(); i++)
   for (int j = 0; j < yl.size(); j++)
       if (fabs(xl[i] - yl[j]) < 0.01)
      r.push_back(xl[i]);
   
    return r;
}

int main() {
   
    v2D p0(0.0,300.0);
    v2D p1(50.0,130.0); 
    v2D p2(70.0,55.0);
    v2D p3(250.0,200.0);
   
    v2D p(50.0,130.0);
   
    std::vector<double> r = reizeb(p,p0,p1,p2,p3);
   

    std::cout << "Result: " << std::endl;
    for (int i = 0; i < r.size(); i++)
   std::cout << r[i] << std::endl;
 
    double t = 0;

    int counter = 0;
    int k = 0;
   
    while (t <= 1) {
   
   v2D v = bezier(t,p0,p1,p2);
   std::vector<double> r = reizeb(v,p0,p1,p2);
   
   //std::cout << "Testwert : " << t << " Antwort : " << r << std::endl;
   
   if (r.size() == 0 || (r.size() != 0 && fabs(r[0] - t) > 0.01))
       counter++;
   
   t = t + 0.00001;
   k++;
     
   }

    std::cout << counter << " Fehler von " << k << " Versuchen" << std::endl;
   
}



P.S. Kennst Du den folgenden Informatiker-Witz:

Betrachte folgendes Axiomensystem:

- Jedes Programm enthält mindestens einen Fehler
- Jedes Programm kann eine Zeile kürzer geschrieben Werden

Mit Induktion folgt daraus:

Jedes Programm kann bis auf eine Zeile reduziert werden, die nicht funktioniert

Ha.Ha. So ist es wohl bei meinem Code oben.

Ich hoffe alles wendet sich zum besten Smile
Beiträge der letzten Zeit anzeigen:   
Foren-Übersicht -> Mathe-Forum -> Umformen einer monströsen Gleichung
Neues Thema eröffnen   Neue Antwort erstellen Alle Zeiten sind GMT + 1 Stunde
Gehe zu Seite Zurück  1, 2, 3, 4  Weiter
Seite 3 von 4

 
Gehe zu:  
Du kannst keine Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum nicht antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.

Chat :: Nachrichten:: Lexikon :: Bücher :: Impressum