Algorithmen zur Lösung von linearen Gleichungssystemen

    • VB.NET

    Es gibt 1 Antwort in diesem Thema. Der letzte Beitrag () ist von ErfinderDesRades.

      Algorithmen zur Lösung von linearen Gleichungssystemen

      Hallo,

      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)

      Spoiler anzeigen

      VB.NET-Quellcode

      1. Public Shared Function SolveGauss(MatA As Double(,), MatB As Double(,)) As Double(,)
      2. Dim A As Double(,)
      3. Dim B As Double(,)
      4. Dim SUM1 As Double
      5. Dim i As Integer, ii As Integer, j As Integer, k As Integer, ll As Integer, Rows As Integer, Cols As Integer
      6. '#Region "LU Decompose"
      7. A = DirectCast(MatA.Clone(), Double(,))
      8. B = DirectCast(MatB.Clone(), Double(,))
      9. Rows = A.GetUpperBound(0)
      10. Cols = A.GetUpperBound(1)
      11. Dim IMAX As Integer = 0, N As Integer = Rows
      12. Dim AAMAX As Double, Sum2 As Double, Dum As Double
      13. Dim INDX As Integer() = New Integer(N) {}
      14. Dim VV As Double() = New Double(N * 10 - 1) {}
      15. Dim D As Double = 1.0
      16. For i = 0 To N
      17. AAMAX = 0.0
      18. For j = 0 To N
      19. If Math.Abs(A(i, j)) > AAMAX Then
      20. AAMAX = Math.Abs(A(i, j))
      21. End If
      22. Next
      23. If AAMAX = 0.0 Then
      24. Return Nothing
      25. End If
      26. VV(i) = 1.0 / AAMAX
      27. Next
      28. For j = 0 To N
      29. If j > 0 Then
      30. For i = 0 To (j - 1)
      31. Sum2 = A(i, j)
      32. If i > 0 Then
      33. For k = 0 To (i - 1)
      34. Sum2 = Sum2 - A(i, k) * A(k, j)
      35. Next
      36. A(i, j) = Sum2
      37. End If
      38. Next
      39. End If
      40. AAMAX = 0.0
      41. For i = j To N
      42. Sum2 = A(i, j)
      43. If j > 0 Then
      44. For k = 0 To (j - 1)
      45. Sum2 = Sum2 - A(i, k) * A(k, j)
      46. Next
      47. A(i, j) = Sum2
      48. End If
      49. Dum = VV(i) * Math.Abs(Sum2)
      50. If Dum >= AAMAX Then
      51. IMAX = i
      52. AAMAX = Dum
      53. End If
      54. Next
      55. If j <> IMAX Then
      56. For k = 0 To N
      57. Dum = A(IMAX, k)
      58. A(IMAX, k) = A(j, k)
      59. A(j, k) = Dum
      60. Next
      61. D = -D
      62. VV(IMAX) = VV(j)
      63. End If
      64. INDX(j) = IMAX
      65. If j <> N Then
      66. If Math.Abs(A(j, j)) < 10 ^ -10 Then 'GLS singulär
      67. Return Nothing
      68. 'A(j, j) = TINY
      69. End If
      70. Dum = 1.0 / A(j, j)
      71. For i = j + 1 To N
      72. A(i, j) = A(i, j) * Dum
      73. Next
      74. End If
      75. Next
      76. If Math.Abs(A(N, N)) < 10 ^ -10 Then 'GLS singulär
      77. Return Nothing
      78. 'A(N, N) = TINY
      79. End If
      80. '#End Region
      81. ii = -1
      82. For i = 0 To N
      83. ll = INDX(i)
      84. SUM1 = B(ll, 0)
      85. B(ll, 0) = B(i, 0)
      86. If ii <> -1 Then
      87. For j = ii To i - 1
      88. SUM1 = SUM1 - A(i, j) * B(j, 0)
      89. Next
      90. ElseIf SUM1 <> 0 Then
      91. ii = i
      92. End If
      93. B(i, 0) = SUM1
      94. Next
      95. For i = N To 0 Step -1
      96. SUM1 = B(i, 0)
      97. If i < N Then
      98. For j = i + 1 To N
      99. SUM1 = SUM1 - A(i, j) * B(j, 0)
      100. Next
      101. End If
      102. B(i, 0) = SUM1 / A(i, i)
      103. Next
      104. Return B
      105. End Function



      Cholesky-Algorithmus (Voraussetzung: symmetrische und positiv definite Matrizen)

      Original von Colin Caprani (colincaprani.com/programming/visual-basic/)

      Spoiler anzeigen

      VB.NET-Quellcode

      1. Public Shared Function solveCholesky(MatA(,) As Double, MatB As Double(,)) As Double(,)
      2. Dim i, j, k, n As Integer
      3. Dim sum As Double
      4. n = MatA.GetUpperBound(0)
      5. Dim p(n) As Double
      6. Dim A(,) As Double = MatA.Clone
      7. For i = 0 To n
      8. For j = i To n
      9. sum = A(i, j)
      10. For k = i - 1 To 0 Step -1
      11. sum = sum - A(i, k) * A(j, k)
      12. Next k
      13. If i = j Then
      14. If Math.Abs(sum) < 10 ^ -10 Then 'GLS singulär
      15. Return Nothing
      16. Else
      17. p(i) = sum ^ 0.5
      18. End If
      19. Else
      20. A(j, i) = sum / p(i)
      21. End If
      22. Next j
      23. Next i
      24. Dim tmpCholesky(,) As Double
      25. ReDim tmpCholesky(n, 0)
      26. 'Solve Ly = b for y
      27. For i = 0 To n
      28. sum = MatB(i, 0)
      29. For k = i To 0 Step -1
      30. sum = sum - A(i, k) * tmpCholesky(k, 0)
      31. Next k
      32. tmpCholesky(i, 0) = sum / p(i)
      33. Next i
      34. 'Solve LT x = y for x
      35. For i = n To 0 Step -1
      36. sum = tmpCholesky(i, 0)
      37. For k = i + 1 To n
      38. sum = sum - A(k, i) * tmpCholesky(k, 0)
      39. Next k
      40. tmpCholesky(i, 0) = sum / p(i)
      41. Next i
      42. Return tmpCholesky
      43. 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)

      Spoiler anzeigen

      VB.NET-Quellcode

      1. Public Shared Function solveGauCho(MatA(,) As Double, MatB(,) As Double) As Double(,)
      2. Dim i, j, a As Integer
      3. If IsSingular(MatA) = True Then Return Nothing
      4. Dim n As Integer = MatA.GetUpperBound(0)
      5. Dim m As Integer = Bandbreite(MatA)
      6. Dim f(n, m) As Double
      7. m += 1
      8. For i = 0 To n
      9. If n - i < m - 1 Then a = n - i Else a = m - 1
      10. For j = 0 To a
      11. f(i, j) = MatA(i, i + j)
      12. Next j
      13. Next i
      14. Dim b, mr As Integer
      15. mr = m - 1
      16. For i = 0 To n
      17. If i - mr > 0 Then a = i - mr Else a = 0
      18. For j = a To i - 1
      19. f(i, 0) = f(i, 0) - f(j, 0) * f(j, i - j) ^ 2
      20. Next j
      21. If i + mr < n Then a = i + mr Else a = n
      22. For c = i + 1 To a
      23. If c - mr > 0 Then b = c - mr Else b = 0
      24. For j = b To i - 1
      25. f(i, c - i) = f(i, c - i) - f(j, 0) * f(j, i - j) * f(j, c - j)
      26. Next j
      27. f(i, c - i) = f(i, c - i) / f(i, 0)
      28. Next c
      29. Next i
      30. Dim x As Double(,) = MatB.Clone
      31. For i = 0 To n 'Vorwärtseinsetzen
      32. If i - mr > 0 Then a = i - mr Else a = 0
      33. For j = a To i - 1
      34. x(i, 0) = x(i, 0) - f(j, i - j) * x(j, 0)
      35. Next j
      36. Next i
      37. For i = n To 0 Step (-1) 'Rückwärtseinsetzen
      38. x(i, 0) = x(i, 0) / f(i, 0)
      39. If i + mr < n Then a = i + mr Else a = n
      40. For j = i + 1 To a
      41. x(i, 0) = x(i, 0) - f(i, j - i) * x(j, 0)
      42. Next j
      43. Next i
      44. Return x
      45. End Function
      46. Public Shared Function IsSingular(mat(,) As Double) As Boolean
      47. Dim n, i, j, c As Integer
      48. n = mat.GetUpperBound(0)
      49. Dim z(n, n) As Double
      50. Dim d(n, 0) As Double
      51. For i = 0 To n
      52. z(i, i) = 1
      53. d(i, 0) = mat(i, i)
      54. For j = 0 To i - 1
      55. d(i, 0) = d(i, 0) - d(j, 0) * z(j, i) ^ 2
      56. Next j
      57. For c = i + 1 To n
      58. z(i, c) = mat(i, c)
      59. For j = 0 To i - 1
      60. z(i, c) = z(i, c) - d(j, 0) * z(j, i) * z(j, c)
      61. Next j
      62. z(i, c) = z(i, c) / d(i, 0)
      63. Next c
      64. Next i
      65. For i = 0 To d.GetUpperBound(0)
      66. If Math.Abs(d(i, 0)) < 10 ^ -10 Then Return True
      67. Next
      68. Return False
      69. End Function
      70. Public Shared Function Bandbreite(mat(,) As Double) As Integer
      71. Dim n As Integer = mat.GetUpperBound(0)
      72. Dim counter As Integer = 0
      73. Dim allesNull As Boolean = True
      74. For i = n To 1 Step -1
      75. allesNull = True
      76. For k = n - i To 0 Step -1
      77. If mat(k, n - counter) <> 0 Then allesNull = False : Exit For
      78. counter += 1
      79. Next
      80. If allesNull = False Then Return i
      81. counter = 0
      82. Next
      83. Return n
      84. 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:

      Spoiler anzeigen

      VB.NET-Quellcode

      1. Dim mat1 As Double(,) = _
      2. {{21000, 0, -21000, 0}, _
      3. {0, 8, 0, 4}, _
      4. {-21000, 0, 21000, 0}, _
      5. {0, 4, 0, 8}}
      6. Dim s1(,) As Double = {{0}, {2 / 3}, {0}, {-2 / 3}}
      7. Dim mat2 As Double(,) = {{8.4, 0, 4.2}, _
      8. {0, 21000, 0}, _
      9. {4.2, 0, 8.4}}
      10. Dim s2(,) As Double = {{2 / 3}, {0}, {-2 / 3}}
      11. Dim Ergebnis(,) As Double
      12. Ergebnis = SolveGauss(mat1, s1)
      13. Debug.WriteLine("Gauss:")
      14. If Ergebnis IsNot Nothing Then
      15. For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
      16. Else
      17. Debug.WriteLine("GLS singulär")
      18. End If
      19. Debug.WriteLine("")
      20. Ergebnis = solveCholesky(mat1, s1)
      21. Debug.WriteLine("Cholesky:")
      22. If Ergebnis IsNot Nothing Then
      23. For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
      24. Else
      25. Debug.WriteLine("GLS singulär")
      26. End If
      27. Debug.WriteLine("")
      28. Ergebnis = solveGauCho(mat1, s1)
      29. Debug.WriteLine("GauCho:")
      30. If Ergebnis IsNot Nothing Then
      31. For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
      32. Else
      33. Debug.WriteLine("GLS singulär")
      34. End If
      35. Debug.WriteLine("")
      36. Debug.WriteLine("")
      37. Ergebnis = SolveGauss(mat2, s2)
      38. Debug.WriteLine("Gauss:")
      39. If Ergebnis IsNot Nothing Then
      40. For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
      41. Else
      42. Debug.WriteLine("GLS singulär")
      43. End If
      44. Debug.WriteLine("")
      45. Ergebnis = solveCholesky(mat2, s2)
      46. Debug.WriteLine("Cholesky:")
      47. If Ergebnis IsNot Nothing Then
      48. For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
      49. Else
      50. Debug.WriteLine("GLS singulär")
      51. End If
      52. Debug.WriteLine("")
      53. Ergebnis = solveGauCho(mat2, s2)
      54. Debug.WriteLine("GauCho:")
      55. If Ergebnis IsNot Nothing Then
      56. For i = 0 To 2 : Debug.WriteLine("x" & (i + 1) & "=" & Ergebnis(i, 0)) : Next
      57. Else
      58. Debug.WriteLine("GLS singulär")
      59. End If


      Ausgabe:
      Spoiler anzeigen

      Quellcode

      1. Gauss:
      2. GLS singulär
      3. Cholesky:
      4. GLS singulär
      5. GauCho:
      6. GLS singulär
      7. Gauss:
      8. x1=0,158730158730159
      9. x2=0
      10. x3=-0,158730158730159
      11. Cholesky:
      12. x1=0,158730158730159
      13. x2=0
      14. x3=-0,158730158730159
      15. GauCho:
      16. x1=0,158730158730159
      17. x2=0
      18. x3=-0,158730158730159


      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“ ()

      zum einen täte ich empfehlen, immer mit Option Strict On zu programmieren.

      Dann frage ich mich, warum die 2. Matrix zweidimensional ist. Ich verstehe nix von linearen Gleichungssystemen, mir fällt nur auf, dass s zwar als 2d-Array deklariert ist, aber überhaupt keine 2. Dimension enthält.
      Ist das der mathematisch Usus?