Laden...

Gauss-Seidel Algorithmus

Erstellt von merc vor 17 Jahren Letzter Beitrag vor 17 Jahren 11.584 Views
merc Themenstarter:in
15 Beiträge seit 2006
vor 17 Jahren
Gauss-Seidel Algorithmus

Hallo,

Ich habe mich im Rahmen einer Seminararbeit mit den beiden iterativen Verfahren linearer Gleichungssysteme und deren Implementierung beschäftigt. Der Jacobi-Algorithmus funktioniert soweit einwandfrei. Nachfolgend habe ich aber mal den Gauss-Seidel Algorithmus (Komponentenschreibweise) und meinen C#-Code aufgelistet. Habe auch ein Beispiel. Jedoch stimmt das Ergebnis noch nicht wirklich.

Hier erstmal die Komponentenschreibweise:

Dann der Algorithmus in Pseudocode: Link zu Wikipedia (engl.)

Und zuletzt mein Code:


/* Es werden 3 Arrays der Methode übergeben: für die Koeffizientenmatrix (matArr), für den
 * b-Vektor (vecArr) und für den Ergebnisvektor (resArr) sowie die Matrixgröße (a)
 * und die Iterationsschritte (iSteps) */

public void solvedByGaussSeidel(decimal[,] matArr, decimal[] vecArr, decimal[] resArr, int iSteps, int a)
        {
            decimal[] tmpArr = new decimal[a]; //Array zum Zwischenspeichern
            
            //c=convergence; i=row; j=column of matrix
            for (int c = 0; c < iSteps; c++)
            {
                for (int i = 0; i < a; i++)
                {
                    resArr[i] = 0;
                    for (int j = 0; j < i-1; j++)
                    {                        
                        resArr[i] = resArr[i] + (matArr[i, j] * resArr[j]);                        
                    }                    
                    for (int j = i + 1; j < a; j++)
                    {
                        resArr[i] = resArr[i] + (matArr[i, j] * tmpArr[j]);
                    }
                    resArr[i] = (vecArr[i] - resArr[i]) / matArr[i, i];
                    tmpArr = resArr;
                }                
            }
            showResult(resArr, gbStartResult, a); //für die Ausgabe
        }

Vermute mal, dass ich im Code irgendwo mit den Indizes falsch gespielt habe oder eine Zuweisung vergessen habe. Nur wo??? 🤔

So, und nun hier mal eine Beispiel für ein LGS um meine Vermutung zu untermauern:

Bekomme mit 10 Iterationsschritten aber nicht ganz folgende Werte raus. Von der Tendenz und dem Vorzeichen stimmt es zwar. Doch je mehr Iterationsschritte ich benutze, desto weniger konvergiert es. Es divergiert eher. Dabei soll das Beispiel (und zahlreiche andere, die probiert habe) stabil und konvergent verlaufen. 🤔

Außerdem ist mir aufgefallen, dass sich die Ergebnisse für den x-Vektor auch beim funktionierenden Jacobi-Verfahren in den ersten 3-4 Iterationsschritten ab der 3. Nachkommastelle im Vergleich zu den Ergebnissen, die mit dem Mathematik-Programm Mathematica ermittelt wurden, leicht unterscheiden. Welchen Datentyp sollte man denn für solche iterative Verfahren am besten in C# benutzen? Habe gelesen, dass decimal mit 28 Kommastellen gut zurecht kommt. Oder ist es vielmehr ein Rundungsfehler-Problem?

T
210 Beiträge seit 2006
vor 17 Jahren

1.: Ich bin mir nicht sicher, aber könnte es sein, daß Du das tmpArray initialisieren (mit Nullen füllen) mußt?
Das Array an sich ist zwar Initialisiert, aber die Felder nicht, oder?

2.: Mußt Du nicht statt

tmpArray = resArray;

tmpArray = Array.Copy(resArray, tmpArray, a);

benutzen?

Gruß,
T-Man

merc Themenstarter:in
15 Beiträge seit 2006
vor 17 Jahren

Hi T-Man,

Habe deine beiden Tipps mal gleich umgesetzt. Leider kein Unterschied zu vorher. 🙁

T
512 Beiträge seit 2006
vor 17 Jahren

Also der Algorithmus scheint mir eine Annäherung zu sein.
In deinem Code iterierst du über eine feste Anzahl an Schritten. Der Algorithmus auf Wikipedia bricht erst dann ab, wenn die Änderung im aktuellen Schritt unter eine bestimmte Grenze fällt.

Außerdem stört mich auch diese Zeile:
tmpArr = resArr;

In der Schleife über i berechnest du doch sowieso nur die i-te Komponente des Vektors. Das Array zu kopieren ist also völlig unnötig, es kann sich nur an der Stelle i verändert haben.

Ich würde also folgendes daran ändern:

  1. iSteps musst du uminterpretieren als eine maximale Anzahl an Schritten. Sprich im schlechtesten Falle macht er iSteps Iterationen, aber er kann auch schon vorher aufhören wenn eine Konvergenz erreicht ist. Dementsprechend hoch müsstest du iSteps anlegen um ein gutes Ergebniss zu bekommen.

  2. tmpArr = resArr; kopiert nicht den Inhalt, sondern nur die Referenzen. Änderungen an tmpArray ändern also auch resArray. Deswegen musst du das kopieren. Ich würde dazu resArr.CopyTo( tmpArr, 0 ); benutzen.

  3. Kopiere das Array erst hinter der Zählschleife über i, nicht darin. Und bevor du kopierst musst du die Differenz zum vorhergehenden Vektor tmpArray ausrechnen. Wikipedia schweigt sich darüber aus. Ich würde da aber spontan das Betragsquadrat des Differenzvektors nehmen. Also die Summe über (tmpArr_-resArr_)*(tmpArr_-resArr_)

e.f.q.

Aus Falschem folgt Beliebiges

B
1.529 Beiträge seit 2006
vor 17 Jahren

Und falls du dich wirklich mit den Indize verhauen hast, kannst du das einfach durch Stift, Zettel, mitrechnen und F11 herausfinden...

merc Themenstarter:in
15 Beiträge seit 2006
vor 17 Jahren

Erstmal vielen Dank für die Tipps und Hilfen! 🙂

Original von Traumzauberbaum
3. Kopiere das Array erst hinter der Zählschleife über i, nicht darin. Und bevor du kopierst musst du die Differenz zum vorhergehenden Vektor tmpArray ausrechnen. Wikipedia schweigt sich darüber aus. Ich würde da aber spontan das Betragsquadrat des Differenzvektors nehmen. Also die Summe über (tmpArr_-resArr_)*(tmpArr_-resArr_)
Vom mathematischen Prinzip hab ich's verstanden. Es wird so lange iteriert, bis die Konvergenz erreicht wird, also die Differenz zwischen der neu Iterierten und der zuvor ermittelten Iterierten beispielsweise kleiner als ein bestimmer Wert (0,00001) ist. So habe ich mir das jedenfalls vorgestellt.

Die Frage ist jetzt, wie erhalte ich die Vektordifferenz? Und wieso würdest du das Betragsquadrat benutzen? Zwischen den beiden sehe ich keinen Zusammenhang... 🤔

merc Themenstarter:in
15 Beiträge seit 2006
vor 17 Jahren
Abbruchkriterium

Hi,

Wann sollte man denn die Iterationsschleife abbrechen? Wenn ALLE Elemente des Differenzvektors (Vektor der neu Iterierten - Vektor der alten Iterierten) einen gewissen Grenzwert unterschritten haben, oder nur EINIGE von denen?

X(

49.485 Beiträge seit 2005
vor 17 Jahren

Hallo merc,

ohne den ganzen Thread gelesen zu haben, würde ich pauschal sagen: wenn die "Länge" des Differenzverktors einen gewissen Grenzwert unterschritten hat.

herbivore

T
512 Beiträge seit 2006
vor 17 Jahren

Original von merc
Hi,

Wann sollte man denn die Iterationsschleife abbrechen? Wenn ALLE Elemente des Differenzvektors (Vektor der neu Iterierten - Vektor der alten Iterierten) einen gewissen Grenzwert unterschritten haben, oder nur EINIGE von denen?

X(

Die euklidische Norm. Dort spielen alle Komponenten des Vektors eine Rolle. Und weil es unsinnig ist eine Wurzel zu ziehen, wenn man das andere einfach Quadrieren kann, betrachtet man da meist das Normquadrat.

Sqrt( a ) < e

ist äquivalent zu (zumindest in unserem Fall wo a ≥ 0 und e ≥ 0)
a < e*e

Zweiteres ist nur eben leichter und genauer zu berechnen.

Und wenn du das geschaft hast, geb ich dir noch nen Tipp wie du tmpArr komplett sparen kannst 😉 Aber bring erstmal eine Version zum laufen.

e.f.q.

Aus Falschem folgt Beliebiges

merc Themenstarter:in
15 Beiträge seit 2006
vor 17 Jahren

Hallo Traumzauberbaum, Danke für den Hinweis!! 🙂 Hab deinen Tipp mit der "2-Norm" gleich in meinen anderen Algorithmus des Jacobi-Verfahrens eingebaut und hat auch super geklappt! Das gesamte Jacobi-Verfahren funktioniert einwandfrei, wie ich an einem Beispiel feststellen konnte. Bleibt jetzt nur noch mein Sorgenkind "Gauß-Seidel" übrig (konvergiert immer noch nicht). 🙁

Hier mal der Code der gesamten Methode:

public void solvedByGaussSeidel(decimal[,] matArr, decimal[] vecArr, decimal[] resArr, int iSteps, int dim, int loadedDim, GroupBox gbStartResult)
        {            
            if (loaded)
                dim = loadedDim;

            decimal[] tmpArr = new decimal[dim];
            dif = 1;

            for (int i = 0; i < dim; i++)
                tmpArr[i] = 0;
            
        //> ALGORITHMUS START:

            //Iterationsschleife bis Konvergenz erreicht:
            for (int c = 0; c < iSteps; c++)
            {
                /* Normwerte := 0 damit Differenz bei jedem
                 * neuen Schritt neu berechnet werden kann:*/
                normNew = 0;
                normOld = 0;

                //Falls Grenzwert kleiner Differenzwert:
                if (limit < dif)
                {
                    //Zeilenschleife:
                    for (int i = 0; i < dim; i++)
                    {
                        resArr[i] = 0;

                        //Spaltenschleifen:
                        for (int j = 0; j < i - 1; j++)
                        {
                            resArr[i] = resArr[i] + (matArr[i, j] * resArr[j]);
                        }
                        for (int j = i + 1; j < dim; j++)
                        {
                            resArr[i] = resArr[i] + (matArr[i, j] * tmpArr[j]);
                        }
                        resArr[i] = (vecArr[i] - resArr[i]) / matArr[i, i];
                    }
                }

                //Berechnung der Euklidischen Norm:
                for (int k = 0; k < dim; k++)
                {
                    normNew += Math.Abs(resArr[k] * resArr[k]);
                    normOld += Math.Abs(tmpArr[k] * tmpArr[k]);
                }

                //Berechnung der Differenz:
                dif = normNew - normOld;

                MessageBox.Show(dif.ToString());

                resArr.CopyTo(tmpArr, 0);                
            }            
            showResult(resArr, gbStartResult, dim);
        }

Woran liegt's? Müsste nicht beim gleichen Matrizen-Beispiel wie beim Jacobi-Verfahren genau derselbe Lösungsvektor rauskommen? Gauß-Seidel konvergiert doch nur schneller, oder?

Und ja, wie kann ich das tmpArr weglassen? Aber mich würde erst mal interessieren, weshalb der Algorithmus immer noch keine brauchbaren Werte ausspuckt. 🤔

Gruß,
merc

T
512 Beiträge seit 2006
vor 17 Jahren

Imo müsste die eine Schleife:

for (int j = 0; j < i - 1; j++)

eher so aussehen:

for (int j = 0; j < i; j++)

Damit wird es für alle j die kleiner als i sind erledigt. So ist es ja auch gedacht. Für alle j kleiner i wird der neue Wert genommen, für alle j größer i wird der alte Wert genommen.

Ich hab jetzt meine Lösung nicht mehr da, aber die hat auf jeden Fall funktioniert, und sah ziemlich ähnlich aus.

Und du machst einen Fehler beim Berechnen der Differenz. Du nimmst die Differenz der Normen, du musst aber die Norm der Differenz nehmen. Mal zur veranschaulichung 2 Vektoren: a = (10,10), b = (-10,-10)
Die sind offensichtlich verflucht unterschiedlich. Nach deiner Formel käme aber als Differenz der Normen 0 raus, was heißt es wäre konvergiert.
Der Differenzvektor wäre aber (20,20) oder (-20,-20). Die Norm davon ist 800. Was ziemlich weit von 0 entfernt ist. Deines könnte also durch dumme Zufälle viel eher konvergieren. Konvergenz nennt man es ja dann, wenn die echte Lösung nicht weit von der aktuellen Lösung entfernt liegt. Und das ist die Norm des Differenzvektors, nicht die Differenz der Normen.

Du musst also im Prinzip Summieren über (resArr_-tmpArr_)²

Außerdem würde ich empfehlen, bei eintreten der Konvergenz ( norm < dif ) einfach ein break einzuwerfen. Das ändert jetzt nichts an der Korrektheit, aber die Schleife wird eben verlassen, weil es ja nix mehr zu tun gibt.

e.f.q.

Aus Falschem folgt Beliebiges

merc Themenstarter:in
15 Beiträge seit 2006
vor 17 Jahren

Hi, stimmt. Lag an der Schleifenvariablen. 🙂

Danke für das Bsp. zur Berechnung des Differenzvektors. Habe das nun wie folgt implementiert:

for (int k = 0; k < dim; k++)
{
   dif += Math.Abs((resArr[k] - tmpArr[k]) * (resArr[k] - tmpArr[k]));                 
}

Hab mir mal in einer MessageBox zu jedem Iterationsschritt die Differenz ausgeben lassen. Bei meiner vorigen Version wurden die Differenz kleiner und das Verfahren hat konvergiert. Jetzt aber konvergiert es zwar auch, doch die Differenzwerte werden größer. In einem Bsp. ist die dif = 2, dann = 2,5 und zuletzt 2,8... Müssten die nicht kleiner werden??? 🤔

T
512 Beiträge seit 2006
vor 17 Jahren

Ich hoffe du setzt dif auch vor der Schleife mal zurück auf 0? 😉

e.f.q.

Aus Falschem folgt Beliebiges

merc Themenstarter:in
15 Beiträge seit 2006
vor 17 Jahren

Hehe, danke. 😁