XProfan

Deutsch

p.specht

Der deutsche Mathematiker Jordan erweiterte den Gauß'schen Eliminationsalgorithmus um die Erzeugung von Nullen nicht nur unter, sondern auch oberhalb der Hauptdiagonale. 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. Diesbezüglich besteht also noch Verbesserungspotential.

Kompilieren Markieren Separieren

WindowTitle "Gauss-Jordan Matrix-Inversion"
' Quelle: https://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