Gauss-Algo

Gauss-Algo

Beitragvon CommanderTomalak » 15.05.2010 13:11

Hi Kinners,

habt ihr auch das Problem, dass euer Progi das zweite LGS nicht richtig lösen kann, obwohl alle verwendeten Variablen schon als long double deklariert sind?? Die 5e-17 wird bei mir intern immer falsch behandelt.

Source:

Code: Alles auswählen
void task13()
{
    // Open file
    char fileName[] = "a13-lgs2.dat", pivot;

    ifstream fin (fileName);
    if (fin == 0)
    {
        cout << "File not found..." << endl ;
        return;
    }

    cout << endl << "Use pivot mode? (y, n)";
    cin >> pivot;

    // Get size
    int size;
    fin >> size;


    /**********************************************************/
    // Declaration
    long double matrix[size][size], b[size], matrixC[size][size], bC[size], a, fE;
    long double res[size], diff, cmp, maxVal = 0;
    int i, j, k, maxRow;

    /**********************************************************/
    // Read matrix
    for(i = 0; i < size; i++)
        for(j = 0; j < size; j++)
            fin >> matrix[i][j];

    // Read vector b
    for(i = 0; i < size; i++)
        fin >> b[i];

    // Close file
    fin.close();

    // Show system
    cout << endl << "The system is:" << endl << endl;

    for(i = 0; i < size; i++)
    {
        for(j = 0; j < size; j++)
        {
            if(j == 0) cout << "(";

            cout << setw(5) << matrix[i][j];

            if(j == size - 1)
            {
                cout << setw(5) << ")";
                if(i == (size-1)/2) cout << setw(5) << " * x = ";
            }
        }

        cout << setw(i == (size-1)/2 ? 3 : 10) << "(  " << b[i] << endl;
    }

    /**********************************************************/
    // Transform system to triangle form
    for(i = 1; i < size; i++)
    {
        // Sort matrix
        if(pivot == 'y')
        {
            // Create copies
            memcpy(matrixC, matrix, sizeof matrix);
            memcpy(bC, b, sizeof b);

            maxVal = 0;
            maxRow = 0;

            // Find row with maximum value
            for(k = i-1; k < size; k++)
            {
                cmp = (matrix[k][i-1] < 0 ? -matrix[k][i-1] : matrix[k][i-1]);

                if(cmp > maxVal)
                {
                    maxVal = cmp;
                    maxRow = k;
                }
            }

            // Switch rows
            memcpy(matrix[i-1], matrix[maxRow], sizeof matrix[maxRow]);
            memcpy(matrix[maxRow], matrixC[i-1], sizeof matrixC[i-1]);
            b[i-1] = b[maxRow];
            b[maxRow] = bC[i-1];
        }

        a = matrix[i-1][i-1];

        for(k = i; k < size; k++)
        {
            fE = matrix[k][i-1];

            for(j = i-1; j < size; j++)
                matrix[k][j] -= matrix[i-1][j] / a*fE;

            b[k] -= b[i-1]/a*fE;
        }

        cout << endl;
       
        // Debug output
        for(k = 0; k < size; k++)
        {
            for(j = 0; j < size; j++) cout << matrix[k][j] << " ";
            cout << b[k] << endl;
        }
    }

    /**********************************************************/
    // Calculate result vector
    for(i = size-1; i >= 0; i--)
    {
        diff = 0;

        for(j = i+1; j < size; j++) diff += matrix[i][j]*res[j];

        res[i] = (b[i] - diff)/matrix[i][i];
    }

    /**********************************************************/
    // Output result vector
    cout << endl;

    for(i = 0; i < size; i++)
    {
        if(i == (size-1)/2) cout << "x = ";
        cout << setw(i == (size-1)/2 ? 4 : 8) << "( " << res[i] << setw(5) << " )" << endl;
    }

    return;
}


Ausgabe mit Pivot:
Code: Alles auswählen
Use pivot mode? (y, n)y

The system is:

(5e-17    1    ) * x = (  0.5
(    1    1    )       (  1

1 1 1
0 1 0.5

x =   ( 0.5    )
      ( 0.5    )



Ausgabe ohne Pivot:
Code: Alles auswählen
Use pivot mode? (y, n)n

The system is:

(5e-17    1    ) * x = (  0.5
(    1    1    )       (  1

5e-17 1 0.5
0 -2e+16 -1e+16

x =   ( 0.499817    )
      ( 0.5    )


mit derm ersten LGS kommt das Programm problemlos klar^^
"Das Volk hat das Vertrauen der Regierung verscherzt. Wäre es da nicht doch einfacher, die Regierung löste das Volk auf und wählte ein anderes?"
- Bertolt Brecht
CommanderTomalak
 
Beiträge: 204
Registriert: 19.01.2009 00:42
Wohnort: Karlsruhe

Re: Gauss-Algo

Beitragvon Kazaar » 15.05.2010 21:32

Bei mir kommt ohne Pivotisierung (0; 0,5) raus, mit Pivotisierung wie bei dir (0,5; 0,5). Ich runde aber auch sehr unbehutsam und benutze nur double.

Das Ergebnis für das LGS 1 ist bei mir: 2.7 1 8 2 8 1 8 2.8

Falls jemand noch ein bisschen Code lesen will, hier meiner:
Code: Alles auswählen
/*
* source.cpp
*
*  Created on: 15.05.2010
*      Author: Kazaar
*       Title: Aufgabe 13 - Gauß-Algorithmus mit Spalten-Pivotisierung, nur für eindeutig lösbare LGS
*/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>

using namespace std;

int main() {

   // Aufgabe des Programms: Löse das LGS A x = b mit A aus R^(n X n) und x, b aus R^n
   // A = (a)_ij

   cout << "Löse das LGS A x = b mit A aus R^(n X n) und x, b aus R^n.\n\n";

   // Variablen und Konstanten deklarieren
   char dateiname[30] = "lgs1.txt"; // Dateiname der Daten
   const double eps_double = 0.000000001; // genauer Wert ist noch zu suchen <<<
   int n; // Dimension des betrachteten Vektorraums
   bool pivot; // Pivotisierung?


   // -------------- Daten aus Datei einlesen -------------------------------------------

   ifstream einlesen;
   einlesen.open(dateiname);
   einlesen >> n; // n einlesen

   double a[n][n], b[n], x[n]; // Arrays für Matrix A und Vektoren b, x deklarieren
   cout << "Dimension des Vektorraums: " << n << endl;

   // Matrix einlesen
   for (int i = 0; i < n; i++) { // Zeilen der Matrix einlesen
      for (int j = 0; j < n; j++) { // in einer Zeile die Elemente spaltenweise einlesen
         einlesen >> a[i][j];
      }
   }

   // Vektor b einlesen
   for (int i = 0; i < n; i++) {
      einlesen >> b[i];
   }
   einlesen.close(); // Datei schließen, ab hier wird nichts mehr daraus eingelesen

   // Matrix ausgeben
   cout << "Mit Vektor b erweiterte Matrix A: " << endl;
   for (int i = 0; i < n; i++) { // Zeilen der Matrix ausgeben
      for (int j = 0; j < n; j++) { // in einer Zeile die Elemente spaltenweise ausgeben
         cout << a[i][j] << "\t";
      }
      cout << "|  " << b[i] << "\n";
   }


   // Abfrage nach Pivotisierung
   cout << "Soll mit Spalten-Pivotisierung gearbeitet werden (0/1)? "; // 0 entspricht false, 1 entspricht true
   cin >> pivot;
   //pivot = false;



   // ----------- Lösung berechnen mit Gauß-Algoritmus, evtl. Spalten-Pivotisierung -----------------
   cout << "Berechne Lösung...\n";
   // Gauß
   for (int k = 0; k < n; k++) { // Die aktuell bearbeitete Untermatrix beginnt linksoben bei Element a_kk.
         // Das heißt, die erste zu betrachtende Zeile und Spalte haben den Index k.

      // Pivotisierung
      if (pivot) {

         // Maximumsnorm bestimmen der Elemente in der Spalte mit/unter a_kk
         int maxind = k; // Anfangswert des Elements der Maximumsnorm
         for (int i = k+1; i < n; i++) {
            if ( abs(a[i][k]) > abs(a[maxind][k]) ) { // besser so was wie abs(abs(a[i][k] - a[maxind][k])) > eps_double
               maxind = i;
            }
         }

         if (k != maxind) { // Nur tauschen, wenn nicht oberste Zeile sowieso Pivotzeile ist
               cout << "Pivotisierung: Tausche Zeile " << k+1 << " (oberste in aktueller Untermatrix) mit Zeile " << maxind+1 << ".\n";
               double zeile_tmp[n+1]; // Hilfs-Array zum Tauschen
               for (int i = 0; i < n ; i++) { // Die Pivotzeile zwischenspeichern
                  zeile_tmp[i] = a[maxind][i]; // Die Pivotzeile zwischenspeichern
                  a[maxind][i] = a[k][i]; // Die Zeile k (oberste in aktueller Untermatrix) in Zeile mit Pivotelement verschieben
                  a[k][i] = zeile_tmp[i]; // Die Zeile k (oberste in aktueller Untermatrix) mit Pivotzeile überschreiben
               }
               // Im Vektor b auch die Zeilen tauschen
               zeile_tmp[n] = b[maxind]; // Die Pivotzeile zwischenspeichern
               b[maxind] = b[k]; // Die Zeile k (oberste in aktueller Untermatrix) in Zeile mit Pivotelement verschieben
               b[k] = zeile_tmp[n]; // Die Zeile k (oberste in aktueller Untermatrix) mit Pivotzeile überschreiben
         }
      } // Pivotisierung


      // Eliminiere Einträge in Spalte k unter a_kk
      // Ziehe a_i1/a_kk von Zeile i ab. i = k+1, ..., n-1 (letzte Spalte)
      for (int i = k+1; i < n; i++) { // Zeilen durchgehen beginnend bei Zeile k+1
         double a_ik = a[i][k]; // Dieser Wert (zu nullendes Element) ändert sich gleich, ist aber wichtig. Alternativ: Von rechts durch die Zeile gehen.
         for (int j = 0; j < n; j++) { // Alle Spalten durchgehen
            a[i][j] = a[i][j] - a[k][j] * (a_ik / a[k][k]);
            if (abs(a[i][j]) < eps_double) a[i][j] = 0; // Zu Null machen, was praktisch Null ist.
         }
         b[i] = b[i] - b[k] * (a_ik / a[k][k]); // Vektor b auch ändern
         if (abs(b[i]) < eps_double) b[i] = 0; // Zu Null machen, was praktisch Null ist.

      }

   } // Gauß


   // Lösungsvektor bestimmen durch Rückwärts-Einsetzen
   x[n-1] = b[n-1] / a[n-1][n-1];
   for (int i = n-1; i >= 0; i--) {
      double sum = 0;
      for (int j = i+1; j < n; j++) sum += a[i][j] * x[j];
      x[i] = 1 / a[i][i] * (b[i] - sum);
   }


   // Matrix ausgeben
   cout << "Mit Vektor b erweiterte Matrix A - geändert: " << endl;
   for (int i = 0; i < n; i++) { // Zeilen der Matrix ausgeben
      for (int j = 0; j < n; j++) { // in einer Zeile die Elemente spaltenweise ausgeben
         cout << setprecision(4) << a[i][j] << "\t\t"; // setprecision legt die Zahl der auszugebenden signifikanten Stellen fest
      }
      cout << "|  " << b[i] << "\n";
   }

   // Lösungsvektor ausgeben
   cout << "\nDer Lösungsvektor ist: \n";
   for (int i = 0; i < n; i++) cout << x[i] << "\t";
   cout << endl;



   cout << endl << "-- Beendet --" << endl;
   return 0;
}
Kazaar
 
Beiträge: 326
Registriert: 29.10.2008 21:35

Re: Gauss-Algo

Beitragvon M.A. » 16.05.2010 23:12

Kazaar hat geschrieben:Das Ergebnis für das LGS 1 ist bei mir: 2.7 1 8 2 8 1 8 2.8

Soll es auch. Der Hinweis auf dem Blatt mit der irrationalen Zahl bezieht sich auf die eulersche. e = 2,718281828

Bei mir wird das zweite LGS auch nicht gelöst. Das haben die wohl nur so gewählt, um uns zu ärgern.
Bild
Benutzeravatar
M.A.
 
Beiträge: 366
Registriert: 25.10.2008 16:37

Re: Gauss-Algo

Beitragvon mabl » 17.05.2010 21:02

Kazaar hat geschrieben:Falls jemand noch ein bisschen Code lesen will, hier meiner:
Code: Alles auswählen
   const double eps_double = 0.000000001; // genauer Wert ist noch zu suchen <<<


In float.h existieren DBL_MIN und FLT_MIN Definitionen, evt. macht es auch Sinn einfach nur DBL_EPSILON zu nutzen.
Benutzeravatar
mabl
Site Admin
 
Beiträge: 741
Registriert: 25.10.2008 11:28
Wohnort: Ettlingen, Karlsruhe

Re: Gauss-Algo

Beitragvon mabl » 18.05.2010 11:32

Das Ganze mal mit standard lib vector Klassen. Hat den Vorteil, dass die Matrizen beliebig groß sind, und das Vertauschen der Zeilen sehr einfach zu erledigen ist. :mrgreen:

Code: Alles auswählen
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <float.h>
#include <math.h>

using namespace std;

typedef double element;
typedef vector<element> Vec;
typedef vector<Vec> Matrix;

void print_matrix(Matrix m, bool solution_vec = false) {
   for (unsigned int i = 0; i < m.size(); i++) {
      for (unsigned int j = 0; j < m[i].size(); j++) {
         if (solution_vec && j == m[i].size() - 1)
            cout << '|';

         cout << setw(7) << setprecision(3) << m[i][j];
      }
      cout << endl;
   }
}

int main() {
   unsigned int dimension;
   Matrix m;

   // Read in file
   ifstream file_in;
   file_in.open("/home/mabl/Desktop/c++/a13-lgs1.dat");

   file_in >> dimension;
   cout << "Reading in Matrix with " << dimension << " dimensions..." << endl;

   // Read in matrix
   Vec row(dimension);
   for (unsigned int i = 0; i < dimension; i++) {
      for (unsigned int j = 0; j < dimension; j++)
         file_in >> row[j];
      m.push_back(row);
   }

   // Read in solution vector and append it to matrix
   for (unsigned int i = 0; i < dimension; i++) {
      element tmp;
      file_in >> tmp;
      m[i].push_back(tmp);
   }

   // Print out matrix
   cout << "Matrix and vector read:" << endl;
   print_matrix(m, true);
   cout << endl;

   // Okay, now we attempt to solve it
   cout << "Solving.." << endl;
   for (unsigned int n = 0; n < m.size(); n++) {

      // pivot table function

      // there is a maximum iterator in the std class, but it does not help here as we
      // can only iterate over rows, not columns... so we use brute force and do it by hand.
      unsigned int max_row = n;
      unsigned int max_value =0;
      for  (unsigned int i = n; i < m.size(); i++){
         if (fabs(m[i][n]) > max_value){
            max_value = fabs(m[i][n]);
            max_row = i;
         }
      }

      // now that we have the row index of the maximal value => swap
      m[n].swap(m[max_row]);

      // And do the magic (TM)
      for (unsigned int i = n + 1; i < m.size(); i++) {
         element factor = m[i][n] / m[n][n];
         for (unsigned int j = n; j < m[i].size(); j++) {
            if (j == n)
               m[i][j] = 0;
            else
               m[i][j] -= m[n][j] * factor;

            if (fabs(m[i][j]) < DBL_MIN)
               m[i][j] = 0;
         }
      }
   }

   // Print out matrix
   cout << "Row-Echelon Form of  matrix:" << endl;
   print_matrix(m, true);
   cout << endl;

   // Ok now solve reverse
   for (unsigned int n = dimension; n > 0; n--) {
      element value = m[n-1][dimension];

      for (unsigned int i=n; i<dimension; i++){
         value -= row[i] * m[n-1][i];
      }


      row[n-1] = value / m[n-1][n-1];
   }

   // print out solution
   cout << "Solution vector:" << endl;
   for (unsigned int i=0; i<dimension; i++)
      cout << row[i] << endl;


   return 0;
}


Benutzeravatar
mabl
Site Admin
 
Beiträge: 741
Registriert: 25.10.2008 11:28
Wohnort: Ettlingen, Karlsruhe

Re: Gauss-Algo

Beitragvon CommanderTomalak » 18.05.2010 12:28

Die Matrizen können bei mir auch so groß werden wie sie wollen, ich leg die Größe ja dynamisch fest. Auch wenn der Gieseke behauptet hat, das ginge nicht^^
"Das Volk hat das Vertrauen der Regierung verscherzt. Wäre es da nicht doch einfacher, die Regierung löste das Volk auf und wählte ein anderes?"
- Bertolt Brecht
CommanderTomalak
 
Beiträge: 204
Registriert: 19.01.2009 00:42
Wohnort: Karlsruhe

Re: Gauss-Algo

Beitragvon mabl » 18.05.2010 13:36

Hmmmm. Ich habe gerade mal ein dynamisches

Code: Alles auswählen
cin >> size;
int test[size];


ausprobiert, und komischer weise geht es. Dabei sollte sowas nach meinem Verständnis definitiv nicht funktionieren?! Nicht mal nen compiler warning kommt...

Eigentlich müsste man expizit Speicher auf dem Heap allozieren:

Code: Alles auswählen
cin >> size;
int* test = new int[size];


Ich habe das böse Gefühl dass du im obrigen Beispiel in Speicherbereiche schreibst, in die man nicht schreiben darf.
Benutzeravatar
mabl
Site Admin
 
Beiträge: 741
Registriert: 25.10.2008 11:28
Wohnort: Ettlingen, Karlsruhe

Re: Gauss-Algo

Beitragvon Kazaar » 20.05.2010 07:20

mabl hat geschrieben:
Kazaar hat geschrieben:Falls jemand noch ein bisschen Code lesen will, hier meiner:
Code: Alles auswählen
   const double eps_double = 0.000000001; // genauer Wert ist noch zu suchen <<<


In float.h existieren DBL_MIN und FLT_MIN Definitionen, evt. macht es auch Sinn einfach nur DBL_EPSILON zu nutzen.

Danke.
Kazaar
 
Beiträge: 326
Registriert: 29.10.2008 21:35


Zurück zu Programmieren

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 1 Gast

cron