Wie ein Lineares Gleichungssytem lösen?



  • Ich muss für ein Schulprojekt ein Lineares Gleichungssytem aufstellen, und lösen.
    Nun das aufstellen ist kein Problem und soweit fertig programmiert, nun hab ich jedoch eine n*(n+1) matrix und weiß nicht wie ich den Löser Programmiere.
    Das Problem ist, wir haben das gaussche Lösungsverfahren im LK einfach übergangen, weil wir nen Taschenrechner haben der das kann... 😞 Nun jetzt habe ich ein Lineares Gleichungssystem für Kubische Splines. Bei dem die erste Zeile dann so aussieht

    0 0 0 1 0 0 0 0 ... 0

    nun das Problem ist das das Gleichungssystem zwar überwiegend aus 0 werten besteht (liegt daran das kubische Splines damit berechnet werden sollen), ich mit meinen Mathekenntnissen aber keine Ahnung habe wie ich das anstellen soll.

    Ich bin dankbar für jede Hilfe
    Omti



  • Je nach dem, ob dein LGS dafür geeignet ist, würde sich die Cramersche Regel anbieten:
    http://de.wikipedia.org/wiki/Cramersche_Regel



  • Cramer ist im Regelfall viel zu ineffizient. Ich würde stattdessen eher was anderes benutzen. Der Gauss* ist eigentlich immer ein Standardverfahren, wenn man nichts besseres weiß. Gerne macht man sich auch die Struktur und sonstige Eigenschaften der Matrix zu Nutze (z.B. Bandmatrix, Dreiecksmatrix, Symmetrie, Definitheit, Dünnbesetztheit ...). Kannst Du was näheres zu der Struktur sagen? Von welchen Dimensionen n reden wir überhaupt maximal?

    * http://de.wikipedia.org/wiki/Gaußsches_Eliminationsverfahren



  • Walli schrieb:

    Cramer ist im Regelfall viel zu ineffizient.

    So wie ich das sehe, geht es hier nicht vorrangig um Effizienz. Der Vorteil an Cramer ist, dass es relativ einfach umzusetzen ist.



  • Wenn du es programmieren musst und nicht einmal von Hand berechnen und dann das Ergebnis hardcodieren kannst, dann würde ich dir zu einer LR-Zerlegung raten, das ist nicht so schwer zu berechnen, wenn du ohne Zeilentausch auskommst, d.h. wenn du auf der Diagonalen von 0 verschiedene Einträge hast.
    Andererseits ist mit Zeilentausch auch nicht viel schwerer, würd dir aber raten einfach Scilab/Octave für die Programmierung zu verwenden, dann sparst du dir die ganze Handarbeit und die Programmierung ist ein Kinderspiel 🙂



  • So sieht das Gleichungssystem in etwa aus.
    http://img368.imageshack.us/img368/9313/matrixkubsplinewo9.th.png

    Allerdings sollte es wesentlich größer werden können. Ich hatte eigentlich vor das es mindestens eine 60x61 matrix sein sollte.

    Wenn ich den Gaus richtig verstanden habe, dann subtrahiert man immer die erste Zeile von der zweiten, danach die zweite von der dritten, usw bis die matrix in dreiecks Form ist. Muss ich dann die erste Zeile mit einer vertauschen wo halt kein 0 an erster Stelle steht?

    Danke,
    Omti



  • Omti schrieb:

    So sieht das Gleichungssystem in etwa aus.
    http://img368.imageshack.us/img368/9313/matrixkubsplinewo9.th.png

    Allerdings sollte es wesentlich größer werden können. Ich hatte eigentlich vor das es mindestens eine 60x61 matrix sein sollte.

    Wenn ich den Gaus richtig verstanden habe, dann subtrahiert man immer die erste Zeile von der zweiten, danach die zweite von der dritten, usw bis die matrix in dreiecks Form ist. Muss ich dann die erste Zeile mit einer vertauschen wo halt kein 0 an erster Stelle steht?

    Danke,
    Omti

    muss es selbst implementiert sein, oder kann es schon was fertiges in form einer bibliothek sein? 🙂

    edit:

    wenns letzeres ist und C sein soll, kann ich die gnu scientific library empfehlen 🙂

    würde dann zum beispiel so aussehen:

    #include <gsl/gsl_linalg.h>
    
    double _A[60*60];
    double _Y[60];
    double a = 0.0;
    int s = 0;
    int i = 0;
    gsl_matrix_view A;
    gsl_vector_view Y;
    gsl_vector *_X = NULL;
    gsl_permutation *p = NULL;
    
    // _A und _Y initialisieren
    
    A = gsl_matrix_view_array(_A, 60, 60);
    Y = gsl_vector_view_array(_Y, 60);
    *_X = gsl_vector_alloc(60);
    *p = gsl_permutation_alloc(60);
    
    gsl_linalg_LU_decomp(&A, p, &s);
    gsl_linalg_LU_solve(&A.matrix, p, &Y.vector, _X);
    
    // mit gsl_vector_get lösungen erhalten
    for (i = 0; i < 60; i++)
    {
    	a = gsl_vector_get(_X, i);
    	// sonstwas damit machen
    }
    
    gsl_permutation_free(p);
    gsl_vector_free(_X);
    


  • Im (wirklich guten) Wikipedia-Artikel findest du auch einen passenden passenden Pseudocode.

    Wenn es auch was fertiges sein darf, kannst du boost.uBLAS benutzen.

    Omti schrieb:

    Das Problem ist, wir haben das gaussche Lösungsverfahren im LK einfach übergangen, weil wir nen Taschenrechner haben der das kann... 😞

    uff. Aus dem Grund bin ich dagegen solche Taschenrechner in Schulen einzusetzen. An deiner Stelle würde ich mir den Wikipedia-Artikel wirklich gut durchlesen.



  • Walli schrieb:

    Gerne macht man sich auch die Struktur und sonstige Eigenschaften der Matrix zu Nutze (z.B. Bandmatrix, Dreiecksmatrix, Symmetrie, Definitheit, Dünnbesetztheit ...). Kannst Du was näheres zu der Struktur sagen? Von welchen Dimensionen n reden wir überhaupt maximal?

    Bei kubischen Splines ergibt sich IIRC immer eine Tridiagonalmatrix, wenn man es geschickt anstellt, ist sie sogar symmetrisch. Das steht jedenfalls bei Sedgewick irgendwo.



  • Jan schrieb:

    Walli schrieb:

    Cramer ist im Regelfall viel zu ineffizient.

    So wie ich das sehe, geht es hier nicht vorrangig um Effizienz. Der Vorteil an Cramer ist, dass es relativ einfach umzusetzen ist.

    Naja, ich weiß nicht. Bei allen Fällen, wo ich mir den Einsatz von bikubischen Splines vorstellen kann, geht es schon irgendwie um Effizienz, auch wenn man die Koeffizienten nur einmal berechnen braucht falls die Stützpunkte konstant bleiben. Und der Gauß ist auch nicht viel schwerer zu implementieren.

    Bashar schrieb:

    Walli schrieb:

    Gerne macht man sich auch die Struktur und sonstige Eigenschaften der Matrix zu Nutze (z.B. Bandmatrix, Dreiecksmatrix, Symmetrie, Definitheit, Dünnbesetztheit ...). Kannst Du was näheres zu der Struktur sagen? Von welchen Dimensionen n reden wir überhaupt maximal?

    Bei kubischen Splines ergibt sich IIRC immer eine Tridiagonalmatrix, wenn man es geschickt anstellt, ist sie sogar symmetrisch. Das steht jedenfalls bei Sedgewick irgendwo.

    Ich weiß grad auch nicht genau, ob es unbedingt tridiagonal ist, aber zumindest Bandstruktur müsste sie haben. Die kann man ja relativ einfach in O(n) lösen.



  • Wenn es sich um einen Endomorphismus handelt kann man die Matrix ja eh immer auf die Jordansche Normalform bringen, also als Bidiagonal(abbildungs)matrix darstellen.

    Aber all diese Theorie hilft dir nichts, denn ohne Gauß-Verfahren kommst du nicht sehr weit. Aber du hast ja eh nur eine 60x60 Matrix da reicht auch einfach ein unoptimierter Gauß-Algorithmus, und im Moment bereitet der dir ja schon genug Probleme.

    Hast du mal deinen Lehrer darauf angesprochen ich bin mir sicher dass er dir offene Fragen dazu beantworten kann 🙂



  • Ich hab jetzt mein Progrämmchen so weitergeschrieben dass es die Zeilen vertauscht. Das funktioniert soweit auch ganz gut, es gibt nur noch ein Problem. Die oberste Zeile, in der Nur ein y wert und eine koeffiziente (so heißt da doch :D) stehenrutscht nach ganz unten durch.

    http://img208.imageshack.us/my.php?image=kubumsorthd8.png

    Der Code sieht so aus:

    for(int j=0;j<amt;j++)  //amt ist die Größe der Matrix
    	{
        if(kub [j][j]==0) //kub ist die Matrix
        	{
            for(int k=j;k<=amt;k++)
            	{
                if(kub [k][j]!=0)
                	{
                    for(int l=0;l<=amt;l++)
                    	{
                        tausch=0;
                        tausch=kub[j][l];
                        kub[j][l]=kub[k][l];
                        kub[k][l]=tausch;
                        }
                    }
                }
    		}
    

    Jetzt weiß ich zwar dank wikipedia in etwa wie es weitergeht... (ich habe das nach dem 10ten mal durchlesen endlich halbwegs verstanden ^^)
    Jetzt weiß ich halt leider nur nicht wie ich die Matrix richtig sortiert kriege.

    Naja... Wenn ihr mir dabei helfen könntet, sollte ich es auf die Reihe kriegen.
    Danke,
    Omti



  • So als Tipp Abstraktion erleichtert dir das Leben ich würde an deiner Stelle eine Funktion schreiben, welche zwei Zeilen vertauscht, das erhöht die Übersichtlichkeit und du verlierst nicht den Blick für das Wesentliche.



  • nunja ich habe jetzt meinen Löser fertiggestellt 😃

    for(int j=0;j<amt;j++) //amt ist die größe der matrix
    	{
         for(int k=j+1;k<amt;k++)
         	{
            double x=kub[k][j]/kub[j][j]; //kub ist die Matrix
            for(int l=0;l<=amt;l++)
            	{
                 tausch=kub[k][l]-(kub[j][l]*x);
                 if ((fabs(tausch) < 0.00001) && (l==j))
                    {
    
                     tausch = 0.0;
                     }
                 kub[k][l]=tausch;
                }
            }
    
        }
    
    for(int j=amt-1;j>=0;j--)
    	{
         for(int k=j-1;k>=0;k--)
         	{
            double y=kub[k][j]/kub[j][j];
            for(int l=0;l<=amt;l++)
            	{
                 tausch=kub[k][l]-(kub[j][l]*y);
                 if ((fabs(tausch) < 0.00001) && (l==j))
                    {
    
                     tausch = 0.0;
                     }
                 kub[k][l]=tausch;
                }
            }
    
        }
     for(int j=0;j<amt;j++)  //amt ist die Dimension der Matrix
    	{
        tausch=kub[j][j];
        kub[j][j]= kub[j][j]/kub[j][j];
        kub[j][amt]=kub[j][amt]/tausch;
        }
    

    Das ganze ist jetzt für ein bereits pivotisiertes Gleichungssystem...



  • Bei natuerlichen Randbedingungen fuer Interpolation mit kubischen Splines bekommst Du ein diagonaldominates, tridiagonales Gleichungssystem, welches keine Pivotsuche erfordert. Bei Hermite-Randbedingungen uebrigens auch.



  • Die Transformation in die Jordan-Normalform ist nicht stetig. Daher ist dieses Verfahren für numerische Anwendungen sowieso wertlos.


Log in to reply