Programm zur Berechnung von Ax=b

[joh]

Hallo,
ich bin so langsam echt am verzweifeln. Und zwar soll ich ein Programm implementieren wie hier beschrieben:


Ich sollte vielleicht erwähnen dass ich kein blutiger Anfänger bin, aber auch nicht besonders versiert. Ich hab bisher folgenden Quellcode, aber ich bekomm einfach nicht das richtige heraus...

#include <iostream>
#include <cmath>

using namespace std;


const int n=5;


// Diagonalmatrix
void diagMatrix( double A[n][n] )
  {
    for( int i=0; i<n; ++i )
      {
        for( int k=0; k<n; ++k )
          {
            if( i==k )
              {
                A[i][k] = i+1.0;
              }
            else
              {
                A[i][k] = 0.0;
              }
          }
      }
  }


// Hilbertmatrix
void hilbertMatrix( double A[n][n] )
  {
    for( int i=0; i<n; ++i )
      {
        for( int k=0; k<n; ++k )
          {
            A[i][k] = 1/(i+k+1.0);
          }
      }
  }


// Multiplikation A * x = erg
void mult( const double A[n][n], const double x[n], double erg[n] )
  {
    for( int i=0; i<n; ++i )
      {
        erg[i] = 0.0;
        for( int k=0; k<n; ++k )
          {
            erg[i] += A[i][k] * x[i];
          }
      }
  }


// Funktion f
void f( const double x[n], const double tau, const double erg[n], const double b[n], double fX[n] )
  {
    for( int i=0; i<n; ++i )
      {
        fX[i] = x[i]-(tau*(erg[i]-b[i]));
      }
  }


// Residuum
void residuum( const double erg[n], const double b[n], double r[n] )
  {
    for( int i=0; i<n; ++i )
      {
        r[i] = erg[i] - b[i];
      }
  }


// Norm berechnen
double norm2( const double v[n] )
  {
    double summe = 0.0;
    for( int i=0; i<n; ++i )
      {
        summe += v[i]*v[i];
      }
    return sqrt(summe);
  }


// Hauptfunktion
int main()
  {
    double A[n][n], x[n], b[n], erg[n], fX[n], r[n];
    double tau;
    double eps = 0.00001;
    double e = 1.0;
    int m;
    int z=0;

    // Eingabe
    cout << "Bitte waehlen Sie die Matrix.\n1) Diagonalmatrix\n2) Hilbertmatrix" << endl; cin >> m; cout << endl;
    cout << "Bitte geben sie ein TAU > 0 ein: "; cin >> tau; cout << endl;

    // Vektor b erstellen
    for( int i=0; i<n; ++i )
      {
        b[i] = i;
      }

    // Vektor x initialisieren
    for( int i=0; i<n; ++i )
      {
        x[i] = 0.0;
      }

    // Matrix einfügen
    if( m==1 )
      {
        diagMatrix( A );
      }
    else if( m==2 )
      {
        hilbertMatrix( A );
      }
    else
      {
        cout << "Eingabe muss 1 oder 2 sein." << endl;
      }

    // Folge
    while( e>=eps )
      {
        mult( A, x, erg );
        f( x, tau, erg, b, fX );

        for( int i=0; i<n; ++i )
          {
            x[i] = fX[i];
          }

        mult( A, x, erg );
        residuum( erg, b, r );
        e = norm2( r )/norm2( b );
      }
    return 0;
  }

Könnt ihr mir vielleicht helfen?

mfg joh

[peschmae]

Was funktioniert denn gemäss dir nicht daran?

Ich hab das mal schnell durchgeguckt - sieht ok aus und funktioniert auch bei mir (zumindest der Fehler fällt gegen 0 nach ein paar Iterationen). Natürlich solange du tau richtig wählst...

MfG Peschmä

[joh]

naja... also fuer tau gleich 2 zum Beispiel bekomme ich wirklich sinnfreie werte. Also bezueglich der Diagonalmatrix sind die werte ja leicht zu berechnen.
Als Loesungsvektor sollte ja x=( 0, 1/2, 2/3, 3/4, 4/5 ) rauskommen.
Es kommt aber

Bitte waehlen Sie die Matrix.
1) Diagonalmatrix
2) Hilbertmatrix
1

Bitte geben sie ein TAU > 0 ein: 2

0
-1.93311e+154
-1.9507e+226
-4.86213e+273
-inf

heraus.
Gut es koennte sein, dass ich mein tau falsch waehle, aber fuer 2 sollte es doch funktionieren, oder nicht?

mfg joh

[joh]

also wenn man tau kleiner waehlt, zb 0.1 dann konvergiert das ganze...
haette ich vielleicht auch mal versuchen koennen ;-)
Aber danke fuer deine Bemuehungen.
Du koenntest mir vielleicht ein feedback geben, wie ich das ganze implementiert hab.
Das waere nett.

mfg joh

[peschmae]

also wenn man tau kleiner waehlt, zb 0.1 dann konvergiert das ganze...

Genau. Das ist der Witz des Algorithmus. Und der ganzen Numerik. Konvergiert halt nicht immer, man muss also wissen wann und sicherstellen dass...

Du koenntest mir vielleicht ein feedback geben, wie ich das ganze implementiert hab.
Das waere nett.

Ich weiss ja nicht genau was die Randbedingungen für deine Implementierung sind; für auf nichts aufbauend implementiert ists ganz ok.

Ich persönlich bin nicht so ein grosser Fan von der Verwendung von Referenzen für die Rückgabe von Resultaten wie z.B. in

void mult( const double A[n][n], const double x[n], double erg[n] )

Klar sieht man da am fehlenden Const und am Variablennamen dass da wohl das Ergebnis reinkommt, soo leserlich ist das trotzdem nicht.

In der Praxis würde ich für so Zeugs auf irgend eine Bibliothek mit netten Matrix & Vektorklassen zurückgreifen, für die dann auch die Multiplikations- und Additionsoperatoren überladen sind.
Das macht den Code dann einiges einfacher zu schreiben und lesen. Ein Beispiel wäre die Boost uBLAS Bibliothek.

Oder dann gleich zu Matlab/Octave greifen - der Riesenaufwand das Zeugs in C++ zu schreiben lohnt sich nicht wirklich; schon gar nicht für Uniübungen.

MfG Peschmä

[joh]

Ja, da hast du recht. Man sollte immer wissen was eigentlich rauskommen soll und den algorithmus erstmal mathematisch verstehen.
Bzgl. der Implementierung kann ich nur sagen, das ich mal ein Semester Informatik studiert hab und ich da JAVA "gelernt" hab. Und davon nicht wirklich viel. Kann halt Schleifen etc. Deswegen bin ich im Moment darauf angewiesen, was der Tutor in der Übung uns für Algorithmen gibt mit denen wir das machen sollen. Wie zb die Rückgabe der Werte. Also kann ich nicht groß was ändern, da ich es nicht kann.
Und die Numerik Übung baut auf C++ auf. Also wir müssen darin programmieren.
Aber danke für die Anregungen. Wenn ich mal Zeit aufbringen kann, werd ich mich mal etwas mehr damit beschäftigen, damit ich mal ne Spur von ner Ahnung bekomm über Bibliotheken etc.

mfg joh

[joh]

Jetzt hab ich leider noch ein Problem festgestellt.
Also für die Diagonalmatrix läuft alles einwandfrei. Für die Hilbertmatrix aber leider nicht... Da konvergiert das Ganze irgendwie nicht, wobei ich das seltsam finde da es ja der selbe algorithmus ist, oder habe ich da irgendwo nen fehler übersehen?
Mein fertiger Code sieht jetzt so aus:

#include <iostream>
#include <cmath>

using namespace std;


// Dimension
const int n=5;


// Diagonalmatrix
void diagMatrix( double A[n][n] )
  {
    for( int i=0; i<n; ++i )
      {
        for( int k=0; k<n; ++k )
          {
            if( i==k )
              {
                A[i][k] = i+1.0;
              }
            else
              {
                A[i][k] = 0.0;
              }
          }
      }
  }


// Hilbertmatrix
void hilbertMatrix( double A[n][n] )
  {
    for( int i=0; i<n; ++i )
      {
        for( int k=0; k<n; ++k )
          {
            A[i][k] = 1/(i+k+1.0);
          }
      }
  }


// Multiplikation A * x = erg
void mult( const double A[n][n], const double x[n], double erg[n] )
  {
    for( int i=0; i<n; ++i )
      {
        erg[i] = 0.0;
        for( int k=0; k<n; ++k )
          {
            erg[i] += A[i][k] * x[i];
          }
      }
  }


// Funktion f
void f( const double x[n], const double tau, const double erg[n], const double b[n], double fX[n] )
  {
    for( int i=0; i<n; ++i )
      {
        fX[i] = x[i]-(tau*(erg[i]-b[i]));
      }
  }


// Residuum
void residuum( const double erg[n], const double b[n], double r[n] )
  {
    for( int i=0; i<n; ++i )
      {
        r[i] = erg[i] - b[i];
      }
  }


// Norm berechnen
double norm2( const double v[n] )
  {
    double summe = 0.0;
    for( int i=0; i<n; ++i )
      {
        summe += v[i]*v[i];
      }
    return sqrt(summe);
  }


// Hauptfunktion
int main()
  {
    double A[n][n], x[n], b[n], erg[n], fX[n], r[n];
    double tau;
    double eps = 0.00001;
    double e = 1.0;
    int m;
    int z=0;

    // Eingabe
    cout << "Bitte waehlen Sie die Matrix.\n1) Diagonalmatrix\n2) Hilbertmatrix" << endl; cin >> m; cout << endl;
    cout << "Bitte geben sie ein TAU > 0 ein: "; cin >> tau; cout << endl;

    // Vektor b und x erstellen bzw initialisieren
    for( int i=0; i<n; ++i )
      {
        b[i] = i;
        x[i] = 0.0;
      }

    // Matrix einfügen
    if( m==1 )
      {
        diagMatrix( A );
      }
    else if( m==2 )
      {
        hilbertMatrix( A );
      }

    // Folge
    while( e>=eps )
      {
        mult( A, x, erg );
        f( x, tau, erg, b, fX );

        for( int i=0; i<n; ++i )
          {
            x[i] = fX[i];
          }

        mult( A, x, erg );
        residuum( erg, b, r );
        e = norm2( r )/norm2( b );
        ++z;
      }

      // Ausgabe
      cout << "--------------------------------------------------" << endl;
      cout << "Die Loesung des linearen Gleichungssystems" << endl << endl;

      for( int i=0; i<n; ++i )
        {
          cout << "| ";
          for( int k=0; k<n; ++k )
            {
              cout << A[i][k] << " ";
            }
          cout << "|" << "   ";
          cout << "| " << "x[" << i << "]" << " |";
          cout << "  =  ";
          cout << "| " << b[i] << " |" << endl;
        }

      cout << endl << "nach " << z << " Iterationsschritten, ist:" << endl << endl;

      for( int i=0; i<n; ++i )
        {
          cout << "x[" << i << "] = " << x[i] << endl;
        }
      cout << endl;
      cout << "--------------------------------------------------" << endl;

    return 0;
  }

Danke schonmal.

mfg joh