Simples Lineares Gleichungssystem lösen



  • Huhu,

    ich möchte ein Gleichungssystem mit zwei unbekannten und drei Gleichungen lösen. Genauer geht es um die Berechnung des Schnittpunktes zweier Geraden im Raum. Zur Verfügung habe ich die beiden Stützvektoren und die beiden Richtungsvektoren (als structs).

    Leider habe ich keinen wirklich Ansatz für dieses Problem. Pakete wie LAPACK oder ublas aus der Boost verwirren mich eher, als dass sie was nützen. Auf dem Papier klappt das lösen natürlich nach dem normalen Lösungsverfahren, ich habe aber ganz ehrlich keinen Plan, wie ich sowas in C++ implementieren könnte.

    Hat jemand einen geeigneten Ansatz oder kann mir einen netten Denkanstoß geben?

    Grüße,

    DeBlob



  • Löse es auf Papier allgemein, nicht für konkrete Zahlenwerte. Dann hast du zwei Formeln, die du 1:1 in C++ abtippen kannst.



  • Wenn es allgemein sein soll: Stelle die Matrix des Gleichungssystems auf und implementiere das Gaußsche Eliminationsverfahren. Ziemlich trivial.



  • Meine Empfehlung:

    (1) Reduziere den Geraden-Schnitt auf den Ebenen-Schnitt. Jede Gerade lässt sich als Schnitt zweier Ebenen darstellen. Das führt je Gerade zu zwei Gleichungen mit drei Unbekannten. Gegeben sind Stützvektor s mit Richtung r einer Geraden. Dann kannst du mit Hilfe des Kreuzproduktes zwei Vektoren p und q bestimmen (die Normalvektoren der beiden Ebenen), die jeweils die Länge 1 haben und mit r zusammen alle senkrecht aufeinander stehen. Das geht in etwa so:

    int dimension_mit_kleinstem_betrag(vec3d const& x);
    
    vec3d naeherungsweise_senkrecht_zu(vec3d const& x)
    {
      vec3d result (0,0,0);
      result[dimension_mit_kleinstem_betrag(x)] = 1;
      return result;
    }
    
    pair<vec3d,vec3d> orthogonales_system(vec3d const& r)
    {
      vec3d p = naeherungsweise_senkrecht_zu(r);
      vec3d q = kreuzprodukt(r,p);
      q /= euklidische_norm(q);
      p = kreuzprodukt(q,r);
      p /= euklidische_norm(p);
      return make_pair(p,q);
    }
    

    Dann sind die zwei Ebenen-Normalgleichungen gegeben durch
    pTx=pTsp^T x = p^T s
    qTx=qTsq^T x = q^T s
    wobei x=(x1;x2;x3) der gesuchte 3D-Punkt ist. Der Winkel zwischen diesen beiden Ebenen ist noch 90° und das ist gut so[TM].

    (2) Für jede der Ebenen hast du jetzt so einen Normalvektor n und einen vorzeichenbehafteten Abstand zum Ursprung a (die rechte Seite). Dann stelle folgendes Gleichungssystem auf:
    Cx=yC x = y
    mit
    C=i=1nn_iTn_iC = \sum_{i=1}^n n\_i^T n\_i
    und
    y=i=1nn_iTa_iy = \sum_{i=1}^n n\_i^T a\_i
    Das ist das sogenannte "Normalgleichungssystem" und das hat sich der liebe Herr Gauß überlegt. In diesem Fall ist C eine 3x3-Matrix ganz egal, wieviele Ebenen und/oder Geraden du hattest.

    Ab drei Ebenen, deren Normalvektoren zusammen den ganzen Vektorraum R^3 aufspannen, gibt es eine eindeutige Lösung. Und die Lösung minimiert dann die Summe der quadratischen Abstände zu allen Ebenen. Und so, wie wir aus jeweils einer Geraden zwei Ebenen gemacht haben mit einem 90° Winkel zwischen beiden Ebenen-Normalvektoren führt das dann auch zu einer Minimierung der quadratischen Abstände zu allen Geraden. Toll, oder?

    (3) Die Matrix C des Normalgleichungssystem C*x=y ist symmetrisch. Und sie ist mindestens positiv semidefinit. Das heißt, du kannst das System relativ leicht über eine Cholesky-Zerlegung lösen. Bibliotheken wie Eigen haben so etwas eingebaut. Damit ist das Lösen echt ein Einzeiler. Eigen nennt das u.a. LLT oder LDLT Zerlegung oder so etwas. Das ging irgendwie so:

    vec3d loesung = C.ldlt().solve(y);
    

    oder so ähnlich. Also wirklich sehr einfach.

    Viel Spaß bei der Umsetzung!



  • Also wirklich sehr einfach

    lol, du hast es hier vermutlich mit einem Schueler zu tun und kommst mit Normalgleichungssystem, machst aus Geraden Ebenen, kommst mit Cholesky-Zerlegung und willst Eigen mit ins Boot holen beim Schnittpunkt zweier geraden ...

    Hier fuer 2D, nD musste selbst daraus ableiten: http://paulbourke.net/geometry/pointlineplane/



  • Gut, an einen Schüler habe ich jetzt nicht gedacht. Aber einem Schüler traue ich schon zu, nachzuvollziehen, dass und wie sich eine Gerade im 3D-Raum über zwei Gleichungen mit drei Unbekannten darstellen lässt und was jede dieser zwei Gleichungen in Isolation für eine Lösungsmenge hat (eine Ebene, in der die Gerade liegt).

    Den Vorschlag habe ich deswegen gemacht, weil man damit auch gleich all die anderen (linearen) Schnitt-Probleme in 3D erschlagen kann, also alle möglichen Kombinationen von Geraden und Ebenen und das im Sinne der kleinsten Fehlerquadrate. 😃



  • DeBlob schrieb:

    ich möchte ein Gleichungssystem mit zwei unbekannten und drei Gleichungen lösen.

    Das ist überbestimmt. Wenn die Geraden windschief sind gibt es da keine Lösung.

    Wenn man überbestimmte Gleichungssysteme im Sinne der kleinsten Fehlerquadrate lösen will, reduziert man sie auf "normale" Gleichungssysteme per Givens-Rotationen, Householder-Transformationen oder durch Aufstellen des Normalgleichungssystem. Numerisch, also in der Genauigkeit, gibt es da schon Unterschiede. Wenn eine LA-Bibliothek nicht direkt so etwas anbietet, dann muss man sich da mit einer QR-Zerlegung behelfen. In Eigen z.B. so:

    dieMatrix.householderQr().solve(rechteSeite)
    

    Das Normalgleichungssystem in dem Fall (bei zwei Unbekannten) hätte nur eine 2x2-Matrix. Und das ist dann auch recht einfach invertierbar.



  • Ich würde hier einfach mal den Normalabstand zwischen den beiden Geraden berechnen - inklusive zumindest eines Fusspunktes.

    Dann kann man gucken ob der Normalabstand klein genug ist, um behaupten zu können dass die beiden Geraden sich schneiden (weil ja Gleitkommaberechnungen nie wirklich genau sind).

    Und wenn man beschlossen hat dass sie sich schneiden, dann nimmt man einfach den Fusspunkt als Schnittpunkt. Bzw. falls man beide Fusspunkte hat die Mitte zwischen den beiden Fusspunkten.

    => Problem gelöst.

    Ohne Householder oder Givens oder 1000 anderer Sachen die sowieso keiner Checkt der danach fragt wie man den Schnittpunkt zweier Geraden ermittelt.



  • hustbaer schrieb:

    Ich würde hier einfach mal den Normalabstand zwischen den beiden Geraden berechnen - inklusive zumindest eines Fusspunktes.

    Den Abstand bekommt man übrigens recht einfach:

    abs((stuetz1-stuetz2).dot(richtung1.cross(richtung2).normalized()))
    

    (Eigen-Syntax)

    hustbaer schrieb:

    => Problem gelöst.

    Ohne Householder oder Givens oder 1000 anderer Sachen die sowieso keiner Checkt der danach fragt wie man den Schnittpunkt zweier Geraden ermittelt.

    So'n Fußpunkt oder "Schnittpunkt" (Mittel der beiden Fußpunkte) zu berechnen ist schon schwieriger. Ich gehe mal davon aus, dass der OP folgendes aufstellte:

    \begin{align} & s\_1 + \alpha \, r\_1 \approx s\_2 + \beta \, r\_2 \\ \; \; \Leftrightarrow & \alpha \, r\_1 - \beta \, r\_2 \approx s\_2 - s\_1\\ \end{align}

    alpha und beta unbekannt, 3 Gleichungen (für x,y,z). Und das ist jetzt "einfach" oder wie? So einfach ist das auch nicht. Paul Bourke hat die Rechnungen weggelassen und zeigt nur seine Formeln, die er als Endergebnis bekommen hat. Wenn

    αr_1βr_2(s_2s_1)2\; \; \left|\left| \alpha \, r\_1 - \beta \, r\_2 - \left( s\_2 - s\_1 \right) \right|\right|^2

    über alpha und beta minimiert werden soll, dann schreit das eben nach QR-Zerlegung oder Normalgleichungssystem. Das Normalgleichungssystem erhält man, indem man die Ableitung obigen Ausdrucks bzgl alpha und beta bestimmt und dann Null gleichsetzt, weil im gesuchten Minimum ja die Ableitung 0 ist. Ob der OP den Anspruch hat, jedes Detail verstanden haben zu wollen, weiß ich nicht. Ich wäre zufrieden mit "hier, so kannst du das überbestimmte Gleichungssystem lösen, was du aufgestellt hast:

    Eigen::Matrix<double,3,2> meineMatrix; // 3x2 Matrix
    meineMatrix.col(0) =  r1;
    meineMatrix.col(1) = -r2;
    Eigen::Vector2d alpha_beta = dieMatrix.householderQr().solve(s2-s1);
    

    "

    Wie A.householderQr().solve(b) jetzt implementiert ist, muss man nicht wissen, so lange man weiß, was dabei herauskommt. Solche Bibliotheken gibt es ja nicht zum Spaß.

    Diese Variablen/Index-Schlacht im Quellcode von Paul Bourke gefällt mir ehrlich gesagt nicht. Das ist IMHO schlecht lesbar, zumindest schlechter als

    // returns a and b so that ||s1+a*r1-(s2+b*r2)|| is minimized
    Vector2d intersect2lines(Vector3d s1, Vector3d r1, Vector3d s2, Vector3d r2)
    {
      // minimize  ||s1+a*r1-(s2+b*r2)||  over a;b
      //      <=>  r1*a - r2*b \approx s2 - s1
      Matrix<double,3,2> mat;
      mat.col(0) =  r1;
      mat.col(1) = -r2;
      return mat.householderQr().solve(s2-s1);
    }
    

    🙂



  • Das ist IMHO schlecht lesbar, zumindest schlechter als

    Es ist einfach nur direct geloest. Natuerlich ist die Indexschlacht schlecht lesbar. Aber es ist ebenfalls nicht dumm, diese mit einem Aufruf von Bibliotheksfunktionen zu vergleichen. Du stellst die Funktion 'solve_my_problem' als Bibliothek dem tatsaechlichem Loesen des Problems gegenueber. Beide operieren auf ganz unterschiedlichen Abstraktionsebenen.



  • Naja, "solve_my_problem" trifft es nicht ganz. Hier ist ein Problem auf ein anderes reduziert worden, dessen Lösung man nicht neu erfinden muss sondern einfach verwenden kann.

    Ich hoffe mal, die Mehrheit schreibt sich nicht ständig Matrix-, Vektor-Klassen und Gaußelimination selbst ...



  • krümelkacker schrieb:

    Ich hoffe mal, die Mehrheit schreibt sich nicht ständig Matrix-, Vektor-Klassen und Gaußelimination selbst ...

    Macht doch Spaß und wenn man nicht allzuviel mit linearer Algebra am Hut hat echt kein Ding. 😉



  • krümelkacker schrieb:

    DeBlob schrieb:

    ich möchte ein Gleichungssystem mit zwei unbekannten und drei Gleichungen lösen.

    Das ist überbestimmt. Wenn die Geraden windschief sind gibt es da keine Lösung.

    Wenn die Geraden parallel sind, gibt es auch keine Lösung.:)



  • VS500 schrieb:

    krümelkacker schrieb:

    DeBlob schrieb:

    ich möchte ein Gleichungssystem mit zwei unbekannten und drei Gleichungen lösen.

    Das ist überbestimmt. Wenn die Geraden windschief sind gibt es da keine Lösung.

    Wenn die Geraden parallel sind, gibt es auch keine Lösung.:)

    Oder unendlich viele. 😃


Log in to reply