Invertieren einer Matrix



  • Hallo

    Ich soll für meine Informatikvorlesung ein Programm schreiben, das u.a. eine untere Dreiecksmatrix invertiert, ohne eine zusätzliche Matrix zu verwenden. Ich habe mir den passenden Algorithmus nach Gauss-Jordan-Verfahren erarbeitet und versuche diesen momentan zu programmieren. Allerdings stimmt da irgendwas nicht, und ich finde den Fehler einfach nicht.

    der Quellcode lautet:

    for(int j=0; j<MAX; j++)
    {
    	for (int k=j; k<MAX; k++)
    	{
    		double x=a[k][k];
    		a[k][k]=1;
    
    		for (int i=0; i<=k; i++)
    		{
    	        if (j==0) a[k][i]=a[k][i]/x;
    			double y=a[k][j-1];
    			if (j>0 && j==i+1) a[k][i]=0-a[j-1][i]*y;
    			else if (j>0 && j!=i+1) a[k][i]=a[k][i]-a[j-1][i]*y;
    		}
    	}
    
    }
    

    Die schleifen gehen die Matrix mehrfach durch und sollen so alle Einträge berechnen. Ohne die k-Schleife (in den Gleichungen k durch j ersetzen) bekomme ich die Hauptdiagonale und die erste Nebendiagonale korrekt raus. Die k-Schleife soll jetzt für die anderen Diagonalen korrekte Ergebnisse bringen, indem sie die Berechnung einer Zeile immer mit allen Zeilen darunter wiederholt. Jetzt bekomme ich die erste Spalte korrekt raus. Allerdings stimmen jetzt Haupt- und Nebendiagonale nicht mehr.

    In der k-Schleife wird die Hauptdiagonale gespeichert und auf 1 gesetzt, in der i-Schleife wird jeder Wert der Zeile durch den gespeicherten Wert dividiert. Seid ich die k-Schleife hinzugefügt habe, wird aber nur die Erste Spalte dividiert, die Hauptdiagonale bleibt auf 1 gesetzt und der Rest wird ohne die Division berechnet. Erst wenn das wieder Stimmt kann ich prüfen, ob der Rest auch korrekte Ergebnisse liefert.

    Vielleicht bin ich gerade einfach nur blind und übersehe einen offensichtlichen Fehler. Kann mir jemand sagen, wo das Problem liegt? Ich denke, ich habe alles Nötige beschrieben, falls noch Fragen sind werde ich diese aber gerne beantworten.

    Grüße

    Thorsten



  • Ich habe 5 Minuten auf Deinen Code geschaut und gebe es jetzt auf. Ich habe keine Ahnung, was Du da machst.
    War Gauß-Jordan nicht einfach, eine Zeile von der anderen so oft abziehen, dass ein bestimmtes Element verschwindet, um somit die untere Dreieckmatrix zu leeren?
    Ich sehe bei Dir bloß komische Indizes und Tests und Inkremente...
    Da müsste doch folgender aufbau sein:

    for( unsigned int zeile = 0; zeile < zeilen; ++zeile ) {
        // Zeile normalisieren, sodass matrix[zeile][zeile] == 1
    
        // Für alle folgenden Zeilen
        for( unsigned int zeile2 = zeile + 1; zeile2 < zeilen; ++zeile2 ) {
           // matrix[zeile] so oft von matrix[zeile2] abziehen, dass matrix[zeile2][zeile1] == 0 (also matrix[zeile2][zeile1] mal)
        }
    }
    

    Na gut, man könnte jetzt auch die obere Dreieckmatrix mit leeren, aber es ging doch bei euch nur um untere Dreiecksmatrizen?



  • Hi,

    Du hast dann aber nur ein halbes Gauss-Jordan-Verfahren. Du hast recht, man zieht eine Zeile von den anderen ab, um Einträge zu eliminieren. Um eine Matrix zu invertieren muss man das selbe aber gleichzeitig mit einer Einheitsmatrix machen. Man schreibt üblicherweise beide Nebeneinander und macht aus der Matrix nach Gauss-Jordan eine Einheitsmatrix, und wenn man alle Rechenschritte auch an der Einheitsmatrix daneben durchführt wird diese zur invertierten Matrix.

    Die i-Schleife geht die Zeile nacheinander durch, die k-Schleife zieht die Zeile von den anderen ab. die j-Schleife geht die Zeilen durch, die von den andeeren abgezogen werden sollen. Der Algorithmus sieht etwas anders aus als das übliche invertieren, weil die Aufgabe ist, das Invertieren ohne die zusätzliche Einheitsmatrix zu machen.



  • Soooo, ich habe mal auf dem Zettel rumgefrickelt und dann implementiert und komme auf folgendes:

    void invert_lowertri( float* mat, unsigned int size )
    {
    	for( unsigned int zeile = 0; zeile < size; ++zeile ) {
    		float diagelem = mat[zeile*size+zeile];
    		mat[zeile*size+zeile] = 1.0f;
    
    		// normalisieren
    		for( unsigned int spalte=0; spalte<=zeile; ++spalte ) {
    			mat[zeile*size+spalte] *= 1.0f/diagelem;
    		}
    
    		// abziehen
    		for( unsigned int zeile2 = zeile+1; zeile2 < size; ++zeile2 ) {
    			float factor = mat[zeile2*size+zeile];
    			mat[zeile2*size+zeile] = 0;
    
    			for( unsigned int spalte = 0; spalte <= zeile; ++spalte ) {
    				mat[zeile2*size+spalte] = mat[zeile2*size+spalte] - mat[zeile*size+spalte] * factor;
    			}
    		}
    	}
    }
    

    Die zwei Randfälle für 1/x und 0-x könnte man natürlich noch aus den Schleifen rausziehen (aber so finde ich es klarer, an welcher Stelle ein Element die Interpretation von unterer Quellmatrix zu oberer Zielmatrix wechselt). In meinen Tests und nach meiner Logik stimmte das Ergebnis.



  • Thorsten90 schrieb:

    Hi,

    Du hast dann aber nur ein halbes Gauss-Jordan-Verfahren. Du hast recht, man zieht eine Zeile von den anderen ab, um Einträge zu eliminieren. Um eine Matrix zu invertieren muss man das selbe aber gleichzeitig mit einer Einheitsmatrix machen. Man schreibt üblicherweise beide Nebeneinander und macht aus der Matrix nach Gauss-Jordan eine Einheitsmatrix, und wenn man alle Rechenschritte auch an der Einheitsmatrix daneben durchführt wird diese zur invertierten Matrix.

    Die i-Schleife geht die Zeile nacheinander durch, die k-Schleife zieht die Zeile von den anderen ab. die j-Schleife geht die Zeilen durch, die von den andeeren abgezogen werden sollen. Der Algorithmus sieht etwas anders aus als das übliche invertieren, weil die Aufgabe ist, das Invertieren ohne die zusätzliche Einheitsmatrix zu machen.

    Das Invertieren einer Matrix mit dem Gauß-Jordan-Algorithmus funktioniert mit Hilfe einer Einheitsmatrix.
    Eine zusätzliche brauchst du nicht. 🤡


Anmelden zum Antworten