| |
|
|
p.specht
|
WindowTitle "Matrixinversion"
Font 2:randomize:cls rnd(8^8)
' Q: https://www.rhirte.de/vb/gleichsys.htm#mat
' Für XProfan adaptiert von P. Specht 2012-04
' Demoware, keine wie immer geartete Gewähr!
Var n&=12'= Zeilen, Spalten (Testmatrix-Größe)
Declare A![n&,n&],erg$
@MatrixAufruf n&
' Testmatrix mit Zufallszahlen belegen.
' Anschließend:
' Ob die Matrixinversion richtig rechnet, kann man einfach
' testen, denn die Inverse der Inversen muß wieder die
' Ausgangsmatrix ergeben. Im Beispiel wird deshalb die
' größte absolute Abweichung ausgegeben.
Proc @MatrixAufruf : Parameters n&
Declare i&,j&,Max!,erg$,B![n&,n&],k!
Whileloop n&:i&=&Loop
Whileloop n&:j&=&Loop
A![i&,j&]=(Rnd()-0.5)*1000000
B![i&,j&]=A![i&,j&]
Endwhile
Endwhile
print "TESTMATRIX:" : @Show A![],n&
MatInv A![],N&
print "INVERSE:" : @Show A![],n&
MatInv A![],N&
print "INVERSE RÜCKINVERTIERT:" : @Show A![],n&
' Fehlerbestimmung und -Ausgabe
Max! = -1
Whileloop n&:i&=&Loop
Whileloop n&:j&=&Loop
If Abs(A![i&,j&] - B![i&,j&]) > Max!
Max! = Abs(A![i&,j&] - B![i&,j&])
endif
'erg$ = erg$ + Format$("%e",A![i&,j&] - B![i&,j&]) + " "
endwhile
'erg$ = erg$ + chr$(10)+chr$(13)
endwhile
erg$ = erg$ + "\n Größter Fehler: "+format$("%e",Max!)
print "DIFFERENZ:" :print erg$
waitinput
Clear B![]
Endproc
' Eigentliche Inversion
Proc MatInv :parameters Mat![],N&
Declare Hlp1&[n&],Hlp2&[n&],Hlp3&[n&]
Declare Max!,T!,i&,j&,k&,ix&,iy&
Whileloop n&:i&=&Loop
Hlp3&[i&]=0
Endwhile
Whileloop n&:i&=&Loop
' Suche das größte Element
Max! = 0
Whileloop n&:j&=&Loop
If Hlp3&[j&]<>1
Whileloop n&:k&=&Loop
If (Hlp3&[k&]<>1) AND (Max! <= Abs(Mat![j&,k&]))
iy& = k&
ix& = j&
Max! = Abs(Mat![j&,k&])
EndIf
endwhile
EndIf
endwhile
Hlp3&[iy&] = Hlp3&[iy&] + 1
'Pivotisierung
If ix&<>iy&
Whileloop n&:j&=&Loop
t!=Mat![ix&,j&]
Mat![ix&,j&]=Mat![iy&,j&]
Mat![iy&,j&]=t!
Endwhile
EndIf
Hlp1&[i&] = ix&
Hlp2&[i&] = iy&
'Kontrolle auf mögliches Verschwinden eines Hauptachsenwertes
If Abs(Mat![iy&,iy&]) < 10^-300
Print "Abbruch, Inversion nicht möglich!"
Waitinput :End
Else
T! = Mat![iy&,iy&]
Mat![iy&,iy&] = 1
Whileloop n&:j&=&Loop
Mat![iy&,j&] = Mat![iy&,j&] / T!
EndWhile
Whileloop n&:j&=&Loop
If j&<>iy&
T! = Mat![j&,iy&]
Mat![j&,iy&] = 0
Whileloop n&:k&=&Loop
Mat![j&,k&] = Mat![j&,k&]- Mat![iy&,k&] * T!
endwhile
EndIf
endwhile
EndIf
endwhile
'Rücktausch
Whileloop n&:i&=&Loop
j& = N& + 1 - i&
If Hlp1&[j&]<>Hlp2&[j&]
ix& = Hlp1&[j&]
iy& = Hlp2&[j&]
Whileloop n&:k&=&Loop
T!=Mat![k&,ix&]
Mat![k&,ix&]=Mat![k&,iy&]
Mat![k&,iy&]=T!
endwhile
EndIf
endwhile
'Hilfsspeicher freigeben
Clear Hlp1&[],Hlp2&[],Hlp3&[]
'zur Ausgabe ...
EndProc
' Anzeigen der Matrix
Proc Show :parameters A![],n&
declare i&,j&
Whileloop n&:i&=&Loop
Whileloop n&:j&=&Loop
erg$ = erg$ + Format$("%e",A![i&,j&])+" "
Endwhile
erg$ = erg$+chr$(10)+chr$(13)
Endwhile
print erg$
waitinput 1000
erg$=""
EndProc
' Eine wichtige Anwendung der Matrizeninversion ist die
' Lösung von linearen Gleichungssystemen. Dieses
' Lösungsverfahren hat ja den immensen Vorteil, sofern man
' die invertierte Matrix kennt, besonders elegant zu sein.
' Denn wenn Ax = b das Gleichungssystem in vektorieller
' Schreibweise beschreibt, dann ist x = inv(A)*A*x = inv(A)*b
' bereits die Lösung. Deshalb:
Proc InvMat : parameters n&,a![],x![]
Declare i%,j%
MatInv a![],n&
WhileLoop n&:i&=&Loop
x!(i&)=0
Whileloop n&:j&=&Loop
x!(i&)=x!(i&)+a![i&,j&] * a![j&,n&+1]'<<< Rechte Seite d.LGS
Endwhile
Endwhile
Endproc
|
|
|
| XProfan 11Computer: Gerät, daß es in Mikrosekunden erlaubt, 50.000 Fehler zu machen, zB 'daß' statt 'das'... | 07.05.2021 ▲ |
|
|
|