Gausscher Algorithmus mit Pivotisierung und Inverse einer Matrix:
Hilfsprogramme MatIn und MatOut im File MatrixIO.c
Testumgebung im File MatrixTest.c
Header File Matrix.h zu Matrix.c
Testdaten
// Matrix.c // Modul mit einigen Matrix-Methoden #includeBeispiele#include <stdio.h> #include "Matrix.h" // Kopiert die quadratische n x n Matrix a und speichert sie in b // b = a void Copy(double a[][NMAX], double b[][NMAX], int n) { int i, j; for (i = 0; i < n; i++) for (j = 0; j < n; j++) b[i][j] = a[i][j]; } // Gauss'sches Eliminationsverfahren mit Zeilenpivotisierung // Argumente: // double a[N][N+1] erweiterte Koeffizientenmatrix Read/Write // int n Anzahl der Gleichungen Read // double x[N] Loesungen Write // Resultat: // int Fehlercode 0 fuer Fehler, 1 fuer Erfolg int Gaussalg (double a[][NMAX], int n, double x[]) { int i, j; // Zeile, Spalte int s; // Elimininationsschritt int pzeile; // Pivotzeile int fehler = 0; // Fehlerflag double f; // Multiplikationsfaktor const double Epsilon = 0.01; // Genauigkeit double Maximum; // Zeilenpivotisierung extern FILE *fout; s = 0; do { // die einzelnen Eliminationsschritte fprintf(fout, "Schritt %2i von %2i\n", s+1, n-1); Maximum = fabs(a[s][s]); // groesstes Element pzeile = s ; // suchen for (i = s+1; i < n; i++) if (fabs(a[i][s]) > Maximum) { Maximum = fabs(a[i][s]) ; pzeile = i; } fehler = (Maximum < Epsilon); if (fehler) break; // nicht loesbar if (pzeile != s) // falls erforderlich, Zeilen tauschen { double h; for (j = s ; j <= n; j++) { h = a[s][j]; a[s][j] = a[pzeile][j]; a[pzeile][j]= h; } } // Elimination --> Nullen in Spalte s ab Zeile s+1 for (i = s + 1; i < n; i++ ) { f = -(a[i][s]/a[s][s]); // Multiplikationsfaktor a[i][s] = 0.0; for (j = s+1; j <= n ; j++) // die einzelnen Spalten a[i][j] += f*a[s][j]; // Addition der Zeilen i, s } #if DEBUG MatOut (stdout, a, n, n+1); #endif s++; } while ( s < n-1 ) ; if (fehler) { fprintf (fout, "gauss: Gleichungssystem nicht loesbar\n"); return 0; } else { // Berechnen der Loesungen aus der entstandenen Dreiecksmatrix // letzte Zeile x[n-1] = a[n-1][n] / a[n-1][n-1]; // restliche Zeilen for (i = n-2 ; i >= 0; i-- ) { for (j = n-1 ; j > i ; j-- ) { a[i][n] -= x[j]*a[i][j]; // rechte Seite berechnen } x[i] = a[i][n] / a[i][i]; // Loesung } return 1; } } // Inverse einer Matrix // Berechung nach dem Gauss-Jordan Verfahren // Lit. Helmut Selder, Einführung in die numerische Mathematik für Ingenieure, HANSER // // Das Verfahren löst die n Gleichungssysteme (für je eine Spalte der Einheitsvektoren) // in einem gemeinsamen Eliminationsverfahren. // Statt nur einem Vektor für die rechte Seite, erweitert man die Matrix um alle n Einheitsvektoren // a11 a12 a13 | 1 0 0 // a21 a22 a23 | 0 1 0 // a31 a32 a33 | 0 0 1 // und eliminiert in der Matrix a alle Elemente ausserhalb der Diagonalen // zusätzlich sorgt man durch eine Division fÜr lauter 1 in der Diagonalen // Das schaut dann so aus: // x1 x2 x2 // 1 0 0 | b11 b12 b13 // 0 1 0 | b21 b22 b23 // 0 0 1 | b31 b32 b33 // Dadurch fällt die Berechnung der Lösungen aus der Dreiecksmatrix (xn, xn-1, ... x1) // weg, weil die Lösungen sofort abgelesen werden können. // x1 = b11 bzw. b12 oder b13 // x2 = b21 bzw. b22 oder b23 // x3 = b31 bzw. b23 oder b33 // Das bedeutet, die Einheitsmatrix ist bei der Elimination in die Matrix B übergegangen, // und das ist die Inverse Matrix. int Inverse (double a[][NMAX], double ainv[][NMAX], int n) { int i, j; // Zeile, Spalte int s; // Elimininationsschritt int pzeile; // Pivotzeile int fehler = 0; // Fehlerflag double f; // Multiplikationsfaktor const double Epsilon = 0.01; // Genauigkeit double Maximum; // Zeilenpivotisierung extern FILE *fout; int pivot = 1; // ergänze die Matrix a um eine Einheitsmatrix (rechts anhängen) for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { a[i][n+j] = 0.0; if (i == j) a[i][n+j] = 1.0; } } #if DEBUG MatOut (stdout, a, n, 2*n); #endif // die einzelnen Eliminationsschritte s = 0; do { // Pivotisierung vermeidet unnötigen Abbruch bei einer Null in der Diagnonalen und // erhöht die Rechengenauigkeit Maximum = fabs(a[s][s]); if (pivot) { pzeile = s ; for (i = s+1; i < n; i++) if (fabs(a[i][s]) > Maximum) { Maximum = fabs(a[i][s]) ; pzeile = i; } } fehler = (Maximum < Epsilon); if (fehler) break; // nicht lösbar if (pivot) { if (pzeile != s) // falls erforderlich, Zeilen tauschen { double h; for (j = s ; j < 2*n; j++) { h = a[s][j]; a[s][j] = a[pzeile][j]; a[pzeile][j]= h; } } } // Eliminationszeile durch Pivot-Koeffizienten f = a[s][s] dividieren f = a[s][s]; for (j = s; j < 2*n; j++) a[s][j] = a[s][j] / f; // Elimination --> Nullen in Spalte s oberhalb und unterhalb der Diagonalen // durch Addition der mit f multiplizierten Zeile s zur jeweiligen Zeile i for (i = 0; i < n; i++ ) { if (i != s) { f = -a[i][s]; // Multiplikationsfaktor for (j = s; j < 2*n ; j++) // die einzelnen Spalten a[i][j] += f*a[s][j]; // Addition der Zeilen i, s } } #if DEBUG fprintf(stdout, "Nach %1i-tem Eliminationschritt:\n", s+1); MatOut (stdout, a, n, 2*n); #endif s++; } while ( s < n ) ; if (fehler) { fprintf (fout, "Inverse: Matrix ist singulär\n"); return 0; } // Die angehängte Einheitsmatrix Matrix hat sich jetzt in die inverse Matrix umgewandelt // Umkopieren auf die Zielmatrix for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { ainv[i][j] = a[i][n+j]; } } return 1; } // Produkt // Multipliziert die Matrix a mit der Matrix b , Resultat in Matrix c // n ... Anzahl der Zeilen und Spalten void Produkt (double a[][NMAX], double b[][NMAX], double c[][NMAX], int n) { int i, j, k; for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { c[i][j] = 0.0; for (k = 0; k < n; k++) c[i][j] += a[i][k]*b[k][j]; } } } // end Produkt