Abt. Matrixinversion mittels Gauss-Jordan-Algorithmus
==================================
Zit. aus Wikipedia: Wilhelm Jordan (* 1. März 1842 in Ellwangen, Württemberg; † 17. April 1899 in Hannover) war ein bedeutender deutscher Geodät und Mathematiker. Jordan war zunächst Professor am Polytechnikum Stuttgart (1865–68) und Karlsruhe (bis 1881), wo sein erstes Taschenbuch der Praktischen Geometrie (1873) herausgegeben wurde. Zit. Ende.
Seine Erweiterung des Gauß'schen Eliminationsalgorithmus um die Erzeugung von Nullen nicht nur unter, sondern auch oberhalb der Hauptdiagonale baute er später zu einem direkten, also nicht-iterativen Verfahren der Matrixinversion aus.
Besonders elegant wird sein Verfahren, wenn zusätzlich ...
1. die Matrizen vorab auf numerische Stabilität und Invertierbarkeit geprüft werden ("Explosionsschutz"),
2. im Fall von Nullelementen die aktuelle Spalte gegen jene mit dem Pivotelement der Zeile vertauscht wird (Pivot= Angelpunkt, d.i. das Element mit dem maximalen Absolutbetrag),
3. statt unnötigem "Division durch Null"-Fehlerabbruch der Rechner-Underflow-Grenzwert errechnet und statt der Null eingesetzt wird, sodaß im Ergebnis allenfalls "Rechnerunendlich" auftaucht, was immerhin noch die Chance auf das Ablesen einer analytischen Lösungsgleichung ergibt.
4. das Ganze sogar "in Place", also ohne weiteren Speicherplatz im RAM, erfolgen könnte.
WindowTitle "Gauss-Jordan Matrix-Inversion"
' Quelle: http://www.cs.berkeley.edu/~wkahan/MathH110/gji.pdf
' Aus IBM BASIC übersetzt nach XProfan 11.2a, P. Specht 2012-04
' Ausschließlich für Demo-Zwecke, keine Gewähr!
' Verwendung samt und sonders auf Risiko des Anwenders!!!
' Enthält Überprüfungen auf exzessives Wachstum trotz Spaltenpivotierung.
' und eine Anpassung zur Vermeidung von Nullen als Pivotelemente.
Window 0,0 - %maxx,%maxy-52
Font 2:Randomize:Cls Rnd(8^8)
set("decimals",6):set("numwidth",16)
Var n&=5 ' Für Testmatrix nötige Maximale Zeilen- bzw. Spaltenzahl
Declare A![N&,N&],X![N&,N&],P![N&]
A![1,1]=1 :A![1,2]=0 :A![1,3]=0 :A![1,4]=0 :A![1,5]=0
A![2,1]=2 :A![2,2]=1 :A![2,3]=0 :A![2,4]=0 :A![2,5]=0
A![3,1]=3 :A![3,2]=2 :A![3,3]=1 :A![3,4]=0 :A![3,5]=0
A![4,1]=0 :A![4,2]=0 :A![4,3]=0 :A![4,4]=1 :A![4,5]=0
A![5,1]=0 :A![5,2]=0 :A![5,3]=0 :A![5,4]=0 :A![5,5]=1
n&=3 ' Diese Zeilen/Spaltenzahl hier soll tatsächlich verwendet werden
MatrInvs n&
Show n&
WaitInput
End
proc Show
parameters n&
declare i&,j&
Print " X = "
WhileLoop n&:i&=&Loop
WhileLoop n&:j&=&Loop
print X![i&,j&],
Endwhile
print
Endwhile
print
endproc
Proc MatrInvs :parameters n&
Declare UFL!,EPS!,G!,P!,Q!,T!,i&,j&,k&,l&,m&
' I bis N sind Ganzzahlvariablen, die anderen Doubleprecision Floats!
' Vorab wird der Rundungsfehler sowie Over- und Underflow-Wert festgelegt.
UFL! = val("5.9E-39") '...= max{underflow,1/overflow}-Grenzwerte.
G!=4 : G!=G!/3 : G!=G!-1 ' ... = 1/3 + Rundungswert auf 4/3
EPS! = ABS( ((G!+G!) - 1) + G! ) ' ... = Rundungspegel
G! = 1 'Neue Verwendung von G:
' G zeichnet nun Wachstumsrate des Pivotelementes auf!
' Kopiere Matrix A auf X und speichere das betragsgrößte Argument der Spalte.
' ACHTUNG: In immer-der-selben Spalte wird gesucht, in dem der Zeilenindex läuft!!!
WhileLoop n&:j&=&Loop
P![J&]=0
WhileLoop n&:i&=&Loop
T! = A![I&,J&]
X![I&,J&] = T!
T! = ABS(T!)
Case T! > P![J&] : P![J&] = T!
EndWhile
EndWhile
' Die P![Zeilenindex] beinhalten jew. den größten Betrag dieser Spalte
WhileLoop n&:k&=&Loop ' ... Elimination in Zeile k
Q!=0
J&=K& 'Annahme: Pivotzeile in Diagonale, daher gleich Spalte
' Suche ab & unterhalb k. Spalte das Pivot (Betragsmaximum)
WhileLoop k&,n&,1:i&=&Loop
T!=ABS(X![I&,K&])
if T!>Q!
Q!=T!
J&=I& ' J speichert offenbar die Zeile mit dem Pivot
'ACHTUNG FEHLERQUELLE: IBM-Basic if: Das 'then : : ´' bezieht sich auf alle : : !!!
endif
EndWhile
if Q!=0
Q!=EPS!*P![K&]+UFL!
X![K&,K&]=Q!
endif
if P![K&]>0
Q!=Q!/P![K&]
if Q!>G!
G!=Q!
endif
endif
Case G!<=(8*K&):Goto "OK"
PRINT "Wachstumsfaktor g = ";G!;" geht über Alarmgrenze ";8*K
PRINT "Versuchen Sie, Spalte ";k&;" von A als Spalte 1 zu setzen!"
END ' ... eventuell hilft eine andere Umordnung der Spalten von A.
OK:
P![k&]=j& ' ... speichere die gefundene Pivot-Zeile der Spalte k.
' ...das P![]-Array ist jetzt frei dafür, warum also nicht trotz Float verwenden...
Case J&=K& : GOTO "Skip" ' da vertauschen mit sich selbst sinnlos.
WhileLoop n&:L&=&Loop
Q!=X![J&,L&]
X![J&,L&]=X![K&,L&]
X![K&,L&]=Q!
EndWhile
Skip:
Q! = X![K&,K&]
X![K&,K&] = 1
WhileLoop n&:j&=&Loop
X![K&,J&] = X![K&,J&]/Q!
EndWhile
WhileLoop n&:i&=&Loop
Case I&=K&:Continue
Q!=X![I&,K&]
X![I&,K&]=0
WhileLoop n&:j&=&Loop
X![I&,J&] = X![I&,J&] - X![K&,J&] * Q!
EndWhile
EndWhile
EndWhile
WhileLoop n&-1,1,-1:k&=&Loop
' ... Rücktausch der Spalten von X
J&=P![K&]
Case J&=K&:Continue
WhileLoop n&:i&=&Loop
Q!=X![I&,K&]
X![I&,K&]=X![I&,J&]
X![I&,J&]=Q!
EndWhile
EndWhile
EndProc
Alles anzeigen
Gruss
P.S.: Die Übertragung von IBM-Basic nach XProfan11 war ziemlich fehlerträchtig, konnte durch ausgiebige Tests aber jeweils geklärt werden. Hinweise dazu im Sourcetext.
PPS.: Der Underflow-Grenzwert bezieht sich auf Systeme mit 4-Byte REAL Floats, XProfan verwendet 8-Byte Floats. Bin am Testen...