A*x = b für 4*4 Matrix instabil

[thesa]

Hallo zusammen!
Für ein Mathematikprojekt in unserer Schule habe ich in C++ die Gauss Elimination für eine 4*4 Matrix geschrieben, allerdings sind die Einträge der Matrix später im Gesamtprogramm teilweise sehr klein oder sogar Null.
Dadurch funktioniert das Verfahren nur begrenzt, ich müsste ständig Zeilen oder Spalten tauschen, was aufwendig ist. Unser Lehrer hat erwähnt, dass es unter netlib.com eventuell bereits fertige Löser gibt (Stichwort pivoting??), allerdings haben wir bisher dort nichts vernünftiges gefunden.

Nun meine Frage:
Kennt jemand von euch eventuell einen Link etc. zu einem wirklich stabilen Lösungsverfahren für 4*4 Matrizen? Wäre um Anregungen sehr dankbar, bin recht verzweifelt.

Vielen Dank,
Liebe grüße,
Thesa

[Maxi]

es gibt doch formeln, wie man die inverse einer 4*4 bildet. Is zwar ziemlich viel rechnerei, aber da sihst du dann, obs lösbar is. Wenn die Determinante sehr klein wird, also fast 0, dann ists wahrscheinlich nicht eindeutig lösbar, ansonsten schon glaub ich... Such ma nach Inverse oder so, dann kannst du das machen:

geg: A*x=b
In=inv(A)

In*A*x=Inb
x=Inb

so kommst du dann auf den lösungsvektor

[Jester]

Das dürfte nicht viel Helfen. Die Lösung ist eindeutig, sobald die Determinante ungleich 0 ist. Das Problem ist, daß der einfache Gaußalgo nicht stabil ist. In ungünstigen Situationen kann eine kleine Änderung der Zahlen (also ein kleiner Rechenfehler) zu einer großen Änderung des Ergebnisses kommen. Die Matrixinversion ist übrigens aus Effizienzgründen nicht so gut geeignet. Das direkte Lösen des LGS ist deutlich schneller.

Evtl. ist hier http://library.lanl.gov/numerical/bookcpdf.html in Kapitel 2 etwas zu finden.

MfG Jester

[thesa]

Vielen Dank für die schnelle Antwort, sehr lieb von euch!
Aber wo finde ich ein Programm, dass diese Inverse der 4x4 Matrix berechnet?
Gauss Verfahren hab ich soweit prinzipiell verstanden, aber Inverse habe ich noch nie gehört...bin noch Schülerin.

[Turing]

Und inverse Matrix berechnen soll stabil sein?

NEIN!

Um genau zu sein hast Du damit einen fast immer instabilen, schlecht konditionierten, langsamen und aufwendigen Algorithmus.

Gaußelimination mit Spaltenpivotisierung sollte eigentlich ausreichen und genau so wird es in der Praxis auch gemacht.

Wenn Deine Matrix fast singulär (oder schlecht konditionniert) ist, dann wird jeder Löser instabile Ergebnisse produzieren. Mit Gaußelimination + Pivotisieruing machst Du die Ergebnisse wenigstens nicht noch instabiler als sie durch Deine Matrix ohnehin schon sind.

[Maxi]

warum soll die inversion instabil sein? was is überhaupt mit instabil gemeint? Entweder das LGS is lösbar, dann kann die inverse bilden oder es is nich. Dass das ewig dauert, hab ich ja schon gesagt...

[Turing]

was is überhaupt mit instabil gemeint?

Soviel wie: Wenn ich ein bisschen an meinen Eingangsdaten wackel, dann wackelt die Lösung ganz gewaltig.

Wo wird denn gewackelt? - Na das passiert andauernd. Beispielsweise kann man den dezimalen Wert 0.1 nicht als Fließkommazahl darstellen. Oder den Bruch 1/3 kann man nicht als Dezimalzahl darstellen.

warum soll die inversion instabil sein?

Das solltest Du besser einen Numeriker oder ein tolles Buch fragen, ich versuch mich mal.
Ein Maß für diese o.g. "Wackeleigenschaft" ist die Kondition einer Matrix. Ist sie niedrig - ist das gut, ist sie hoch - ist das schlecht. Invertierung dreht die Kondition (euklidisch) um (1/cond(A)), d.h. entweder hast am Anfang oder am Ende eine schlecht konditionierte Matrix und damit eine große Empfindlichkeit. Gauß + Pivotisierung versucht die Kondition der Ausgangsmatrix bei der Lösung nicht wesentlich zu verschlechtern. Bei der Invertierung hast Du immer eine wesentliche Verschlechterung.

Ganz davon abgesehen, dass Invertierung im Grunde genommen auch wieder die Lösung von n linearen Gleichungssystemen bedeutet...

[wolfgke]

Turing schrieb:

Mit Gaußelimination + Pivotisieruing machst Du die Ergebnisse wenigstens nicht noch instabiler als sie durch Deine Matrix ohnehin schon sind.

Der Korrektheit halber:

durch Gauß+Spaltenpivotisierung verschlechtert sich die Kondiditionszahl nicht schlimmer als um den Faktor n^2 (habe extra nochmal in meine Numerik-Aufzeichnungen geschaut; ich hoffe ihr wisst das zu schätzen ).

Wenn man es wirklich mit schlecht konditionierten Matrizen zu tun hat, sollte man Householder-Verfahren (verschlechtert die Konditionszahl theoretisch nicht und praktisch kaum) => benötigt aber ca. doppelt so viel Rechenzeit wie Gauß, verwenden.

Wenn es ganz grausam mit der Konditionierung aussieht, dann gibt es noch die Singulärwertzerlegung (SVD: Singular Value Decomposition). Die ist aber extrem rechenintensiv und gilt als schwer programmierbar (wir haben nur sehr wenig aus Zeitgründen über diese gelernt).

In der Praxis sollte es für 97% der Fälle Gauß+Spaltenpivotisierung (wenn diese nicht ausreicht: Totalpivotisierung) ausreichen. Für den Rest Householderverfahren. Und wenn es das auch nicht tut, kann man mal einen Blick auf "härteres" Zeug wie SVD werfen.

Ergänzung: Die Konditionszahl einer Matrix A bezüglich einer Norm $|| . ||_p$ ist definiert als

$|| A ||\_p * || A^{-1}||\_p$

Der Beweis, dass die Konditionszahl bezüglich der $L_2$-Norm sich nicht durch das Householderverfahren verschlechtert, sei Übungsaufgabe (ich erkläre es dir jedoch, falls du es brauchst).

[life]

Wenn ungenauigkeiten der Darstellung das Problem sind, dann nimm doch einfach nen Datentyp der für deinen Körper verlustfrei arbeitet. Z.B. eine Bruchklasse, wenn du im Körper der rationalen Zahlen arbeitest..

oder hab ich was nicht verstanden? :>

[PeterTheMaster]

klar. oder rechne symbolisch. nee, man will ja auch mal fertig werden.

[wolfgke]

Die Frage ist: willst du eine numerische oder algebraische Lösung?

Wenn es eine algebraische sein soll: dazu hat sich life geäußert. Ist halt in der Tat erheblich langsamer (für 4x4-Matrix aber absolut kein Thema). Man muss sich halt eine entsprechende Bibliothek, wie GMP oder vergleichbares besorgen. Von Hand derartige Bibliotheken zu programmieren ist nicht lustig (wegen der eingeschränkten Größe von den Standard-C/C++-Typen muss man eine eigene Hochgenauigkeitsarithmetik von Hand implementieren).

Zu der numerischen Variante habe ich alles gesagt. Die Version wie ich sie beschrieben habe, sollte es für bis zu einigen 10000 Variablen tun, darüber hinaus muss man gänzlich andere Verfahren (wie Gauß-Seidel, Jacobi, SOR, CG etc. benutzen, die gänzlich anders funktionieren). Damit sollten ein paar Millionen Variablen heutzutage problemlos machbar sein (Voraussetzung: A ist dünn besetzt).