Hallo,
momentan habe ich mich mit Algorithmen zur Lösung von linearen Gleichungssystemen auseinandersetzen müssen. Dabei habe ich folgende genauer betrachtet:
Da man im Netz nicht so einfach Beispiele dazu findet, wollte ich meine Lösungen hier mal vorstellen.
Vorweg:
Das Gleichungssystem wird per Matrizen dargestellt
[mathelounge.de/]
x und b sind eigentlich Vektoren, ich gebe diese dennoch als (mehrdimensionale) Matrizen aus. Das kann man natürlich ändern. Für mich war es so eben sinnvoller.
Wenn das Gleichungssystem singulär ist, also nicht eindeutig lösbar, dann liefern die Funktionen als Ergebnis eine Matrix mit dem Wert Nothing.
Da ich kein Mathematiker bin und von der Theorie wenig verstehe, habe ich bestehende Quellcodes einfach umprogrammiert und modifiziert.
Es fehlen beispielsweise Abfragen, ob Matrix A quadratisch ist oder generell, ob die Eingangsmatrix überhaupt mit dem entsprechenden Verfahren berechenbar ist o.Ä..
Gauss-Algorithmus (Voraussetzung: symmetrische oder unsymmetrische Matrizen)
Original von Anas Abidis Matrix-Klasse (hier abrufbar: github.com/bcourter/SpaceClaim…r/AddInLibrary/NMatrix.cs)
Spoiler anzeigen
Cholesky-Algorithmus (Voraussetzung: symmetrische und positiv definite Matrizen)
Original von Colin Caprani (colincaprani.com/programming/visual-basic/)
Spoiler anzeigen
GauCho-Algorithmus (Voraussetzung: symmetrische Matrizen)
(eine Mischung auch dem Gauss- und dem Cholesky-Algorithmus, welcher auch eine Bandstruktur von Matrizen ausnutzt um die Berechnung zu beschleunigen, was z.B. bei Gleichungssystemen in der Mechanik oft nützlich ist)
Original von Prof. Kindmann und Prof. Laumann, (ruhr-uni-bochum.de/stahlbau/mam/rubstahl-bericht_1-2004.pdf)
Spoiler anzeigen
Bei meinem Projekt habe ich mal gemessen bei einer schwach besetzten, symmetrischen und positiv definiten Matrix mit ca. 1800 Unbekannten:
Gauss: 35 Sekunden
Cholesky: 15 Sekunden
GauCho: 0,04 Sekunden
Beispiel:
Spoiler anzeigen
Ausgabe:
Korrekturen und Optimierungen sind gerne gesehen
Grüße
momentan habe ich mich mit Algorithmen zur Lösung von linearen Gleichungssystemen auseinandersetzen müssen. Dabei habe ich folgende genauer betrachtet:
- Gauss
- Cholesky
- GauCho
Da man im Netz nicht so einfach Beispiele dazu findet, wollte ich meine Lösungen hier mal vorstellen.
Vorweg:
Das Gleichungssystem wird per Matrizen dargestellt
[mathelounge.de/]
x und b sind eigentlich Vektoren, ich gebe diese dennoch als (mehrdimensionale) Matrizen aus. Das kann man natürlich ändern. Für mich war es so eben sinnvoller.
Wenn das Gleichungssystem singulär ist, also nicht eindeutig lösbar, dann liefern die Funktionen als Ergebnis eine Matrix mit dem Wert Nothing.
Da ich kein Mathematiker bin und von der Theorie wenig verstehe, habe ich bestehende Quellcodes einfach umprogrammiert und modifiziert.
Es fehlen beispielsweise Abfragen, ob Matrix A quadratisch ist oder generell, ob die Eingangsmatrix überhaupt mit dem entsprechenden Verfahren berechenbar ist o.Ä..
Gauss-Algorithmus (Voraussetzung: symmetrische oder unsymmetrische Matrizen)
Original von Anas Abidis Matrix-Klasse (hier abrufbar: github.com/bcourter/SpaceClaim…r/AddInLibrary/NMatrix.cs)
VB.NET-Quellcode
- Public Shared Function SolveGauss(MatA As Double(,), MatB As Double(,)) As Double(,)
- Dim A As Double(,)
- Dim B As Double(,)
- Dim SUM1 As Double
- Dim i As Integer, ii As Integer, j As Integer, k As Integer, ll As Integer, Rows As Integer, Cols As Integer
- '#Region "LU Decompose"
- A = DirectCast(MatA.Clone(), Double(,))
- B = DirectCast(MatB.Clone(), Double(,))
- Rows = A.GetUpperBound(0)
- Cols = A.GetUpperBound(1)
- Dim IMAX As Integer = 0, N As Integer = Rows
- Dim AAMAX As Double, Sum2 As Double, Dum As Double
- Dim INDX As Integer() = New Integer(N) {}
- Dim VV As Double() = New Double(N * 10 - 1) {}
- Dim D As Double = 1.0
- For i = 0 To N
- AAMAX = 0.0
- For j = 0 To N
- If Math.Abs(A(i, j)) > AAMAX Then
- AAMAX = Math.Abs(A(i, j))
- End If
- Next
- If AAMAX = 0.0 Then
- Return Nothing
- End If
- VV(i) = 1.0 / AAMAX
- Next
- For j = 0 To N
- If j > 0 Then
- For i = 0 To (j - 1)
- Sum2 = A(i, j)
- If i > 0 Then
- For k = 0 To (i - 1)
- Sum2 = Sum2 - A(i, k) * A(k, j)
- Next
- A(i, j) = Sum2
- End If
- Next
- End If
- AAMAX = 0.0
- For i = j To N
- Sum2 = A(i, j)
- If j > 0 Then
- For k = 0 To (j - 1)
- Sum2 = Sum2 - A(i, k) * A(k, j)
- Next
- A(i, j) = Sum2
- End If
- Dum = VV(i) * Math.Abs(Sum2)
- If Dum >= AAMAX Then
- IMAX = i
- AAMAX = Dum
- End If
- Next
- If j <> IMAX Then
- For k = 0 To N
- Dum = A(IMAX, k)
- A(IMAX, k) = A(j, k)
- A(j, k) = Dum
- Next
- D = -D
- VV(IMAX) = VV(j)
- End If
- INDX(j) = IMAX
- If j <> N Then
- If Math.Abs(A(j, j)) < 10 ^ -10 Then 'GLS singulär
- Return Nothing
- 'A(j, j) = TINY
- End If
- Dum = 1.0 / A(j, j)
- For i = j + 1 To N
- A(i, j) = A(i, j) * Dum
- Next
- End If
- Next
- If Math.Abs(A(N, N)) < 10 ^ -10 Then 'GLS singulär
- Return Nothing
- 'A(N, N) = TINY
- End If
- '#End Region
- ii = -1
- For i = 0 To N
- ll = INDX(i)
- SUM1 = B(ll, 0)
- B(ll, 0) = B(i, 0)
- If ii <> -1 Then
- For j = ii To i - 1
- SUM1 = SUM1 - A(i, j) * B(j, 0)
- Next
- ElseIf SUM1 <> 0 Then
- ii = i
- End If
- B(i, 0) = SUM1
- Next
- For i = N To 0 Step -1
- SUM1 = B(i, 0)
- If i < N Then
- For j = i + 1 To N
- SUM1 = SUM1 - A(i, j) * B(j, 0)
- Next
- End If
- B(i, 0) = SUM1 / A(i, i)
- Next
- Return B
- End Function
Cholesky-Algorithmus (Voraussetzung: symmetrische und positiv definite Matrizen)
Original von Colin Caprani (colincaprani.com/programming/visual-basic/)
VB.NET-Quellcode
- Public Shared Function solveCholesky(MatA(,) As Double, MatB As Double(,)) As Double(,)
- Dim i, j, k, n As Integer
- Dim sum As Double
- n = MatA.GetUpperBound(0)
- Dim p(n) As Double
- Dim A(,) As Double = MatA.Clone
- For i = 0 To n
- For j = i To n
- sum = A(i, j)
- For k = i - 1 To 0 Step -1
- sum = sum - A(i, k) * A(j, k)
- Next k
- If i = j Then
- If Math.Abs(sum) < 10 ^ -10 Then 'GLS singulär
- Return Nothing
- Else
- p(i) = sum ^ 0.5
- End If
- Else
- A(j, i) = sum / p(i)
- End If
- Next j
- Next i
- Dim tmpCholesky(,) As Double
- ReDim tmpCholesky(n, 0)
- 'Solve Ly = b for y
- For i = 0 To n
- sum = MatB(i, 0)
- For k = i To 0 Step -1
- sum = sum - A(i, k) * tmpCholesky(k, 0)
- Next k
- tmpCholesky(i, 0) = sum / p(i)
- Next i
- 'Solve LT x = y for x
- For i = n To 0 Step -1
- sum = tmpCholesky(i, 0)
- For k = i + 1 To n
- sum = sum - A(k, i) * tmpCholesky(k, 0)
- Next k
- tmpCholesky(i, 0) = sum / p(i)
- Next i
- Return tmpCholesky
- End Function
GauCho-Algorithmus (Voraussetzung: symmetrische Matrizen)
(eine Mischung auch dem Gauss- und dem Cholesky-Algorithmus, welcher auch eine Bandstruktur von Matrizen ausnutzt um die Berechnung zu beschleunigen, was z.B. bei Gleichungssystemen in der Mechanik oft nützlich ist)
Original von Prof. Kindmann und Prof. Laumann, (ruhr-uni-bochum.de/stahlbau/mam/rubstahl-bericht_1-2004.pdf)
VB.NET-Quellcode
- Public Shared Function solveGauCho(MatA(,) As Double, MatB(,) As Double) As Double(,)
- Dim i, j, a As Integer
- If IsSingular(MatA) = True Then Return Nothing
- Dim n As Integer = MatA.GetUpperBound(0)
- Dim m As Integer = Bandbreite(MatA)
- Dim f(n, m) As Double
- m += 1
- For i = 0 To n
- If n - i < m - 1 Then a = n - i Else a = m - 1
- For j = 0 To a
- f(i, j) = MatA(i, i + j)
- Next j
- Next i
- Dim b, mr As Integer
- mr = m - 1
- For i = 0 To n
- If i - mr > 0 Then a = i - mr Else a = 0
- For j = a To i - 1
- f(i, 0) = f(i, 0) - f(j, 0) * f(j, i - j) ^ 2
- Next j
- If i + mr < n Then a = i + mr Else a = n
- For c = i + 1 To a
- If c - mr > 0 Then b = c - mr Else b = 0
- For j = b To i - 1
- f(i, c - i) = f(i, c - i) - f(j, 0) * f(j, i - j) * f(j, c - j)
- Next j
- f(i, c - i) = f(i, c - i) / f(i, 0)
- Next c
- Next i
- Dim x As Double(,) = MatB.Clone
- For i = 0 To n 'Vorwärtseinsetzen
- If i - mr > 0 Then a = i - mr Else a = 0
- For j = a To i - 1
- x(i, 0) = x(i, 0) - f(j, i - j) * x(j, 0)
- Next j
- Next i
- For i = n To 0 Step (-1) 'Rückwärtseinsetzen
- x(i, 0) = x(i, 0) / f(i, 0)
- If i + mr < n Then a = i + mr Else a = n
- For j = i + 1 To a
- x(i, 0) = x(i, 0) - f(i, j - i) * x(j, 0)
- Next j
- Next i
- Return x
- End Function
- Public Shared Function IsSingular(mat(,) As Double) As Boolean
- Dim n, i, j, c As Integer
- n = mat.GetUpperBound(0)
- Dim z(n, n) As Double
- Dim d(n, 0) As Double
- For i = 0 To n
- z(i, i) = 1
- d(i, 0) = mat(i, i)
- For j = 0 To i - 1
- d(i, 0) = d(i, 0) - d(j, 0) * z(j, i) ^ 2
- Next j
- For c = i + 1 To n
- z(i, c) = mat(i, c)
- For j = 0 To i - 1
- z(i, c) = z(i, c) - d(j, 0) * z(j, i) * z(j, c)
- Next j
- z(i, c) = z(i, c) / d(i, 0)
- Next c
- Next i
- For i = 0 To d.GetUpperBound(0)
- If Math.Abs(d(i, 0)) < 10 ^ -10 Then Return True
- Next
- Return False
- End Function
- Public Shared Function Bandbreite(mat(,) As Double) As Integer
- Dim n As Integer = mat.GetUpperBound(0)
- Dim counter As Integer = 0
- Dim allesNull As Boolean = True
- For i = n To 1 Step -1
- allesNull = True
- For k = n - i To 0 Step -1
- If mat(k, n - counter) <> 0 Then allesNull = False : Exit For
- counter += 1
- Next
- If allesNull = False Then Return i
- counter = 0
- Next
- Return n
- End Function
Bei meinem Projekt habe ich mal gemessen bei einer schwach besetzten, symmetrischen und positiv definiten Matrix mit ca. 1800 Unbekannten:
Gauss: 35 Sekunden
Cholesky: 15 Sekunden
GauCho: 0,04 Sekunden
Beispiel:
VB.NET-Quellcode
- Dim mat1 As Double(,) = _
- {{21000, 0, -21000, 0}, _
- {0, 8, 0, 4}, _
- {-21000, 0, 21000, 0}, _
- {0, 4, 0, 8}}
- Dim s1(,) As Double = {{0}, {2 / 3}, {0}, {-2 / 3}}
- Dim mat2 As Double(,) = {{8.4, 0, 4.2}, _
- {0, 21000, 0}, _
- {4.2, 0, 8.4}}
- Dim s2(,) As Double = {{2 / 3}, {0}, {-2 / 3}}
- Dim Ergebnis(,) As Double
- Ergebnis = SolveGauss(mat1, s1)
- Debug.WriteLine("Gauss:")
- If Ergebnis IsNot Nothing Then
- For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
- Else
- Debug.WriteLine("GLS singulär")
- End If
- Debug.WriteLine("")
- Ergebnis = solveCholesky(mat1, s1)
- Debug.WriteLine("Cholesky:")
- If Ergebnis IsNot Nothing Then
- For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
- Else
- Debug.WriteLine("GLS singulär")
- End If
- Debug.WriteLine("")
- Ergebnis = solveGauCho(mat1, s1)
- Debug.WriteLine("GauCho:")
- If Ergebnis IsNot Nothing Then
- For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
- Else
- Debug.WriteLine("GLS singulär")
- End If
- Debug.WriteLine("")
- Debug.WriteLine("")
- Ergebnis = SolveGauss(mat2, s2)
- Debug.WriteLine("Gauss:")
- If Ergebnis IsNot Nothing Then
- For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
- Else
- Debug.WriteLine("GLS singulär")
- End If
- Debug.WriteLine("")
- Ergebnis = solveCholesky(mat2, s2)
- Debug.WriteLine("Cholesky:")
- If Ergebnis IsNot Nothing Then
- For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
- Else
- Debug.WriteLine("GLS singulär")
- End If
- Debug.WriteLine("")
- Ergebnis = solveGauCho(mat2, s2)
- Debug.WriteLine("GauCho:")
- If Ergebnis IsNot Nothing Then
- For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
- Else
- Debug.WriteLine("GLS singulär")
- End If
Ausgabe:
Korrekturen und Optimierungen sind gerne gesehen
Grüße
Für ein Mindestmaß an Rechtschreibung, Interpunktion und Majuskeln!
Dieser Beitrag wurde bereits 1 mal editiert, zuletzt von „bla“ ()