Gauß-Verfahren mit Pivoting

Im letzten Post habe ich das Gaußverfahren zur Lösung von linearen Gleichungssysteme vorgestellt. Mittlerweile habe ich den Algorithmus um Pivoting, also dem Vertauschen von Zeilen und Spalten, ergänzt. Dadurch verhindert man eine mögliche Division durch Null und der Algorithmus läuft deutlich stabiler.


Die Lösung ist noch nicht vollständig ausgereift und getestet. Vor allem ist zu beachten, dass ich grundsätzlich zweidimensionale Datenfelder verwende (dazu später mehr), ein Vektor muss also nicht mit Vektor(10) as Double sondern mit Vektor(10,1) as Double initialisiert werden. Die Übergabe von eindimensionalen Datenfeldern führt zum Abbruch der Funktion.

Public Function GaussEliminationWithPivoting(Matrix() As Double, _
                                    Vektor() As DoubleAs Double()
    Dim L() As Double    'lower diagonal
    Dim U() As Double    'upper diagonal
    Dim y() As Double
    Dim x() As Double
    Dim Row() As Long, Column() As Long
    Dim k As Long, i As Long, j As Long
    Dim Lij As Double
    Dim temp As Double
    Dim dblMax As Double, rowMax As Long, columnMax As Long

    'only valid for symmetric matrix
    If UBound(Matrix, 1) <> UBound(Matrix, 2) Then Exit Function

    'initialize
    ReDim L(UBound(Matrix, 1), UBound(Matrix, 1))
    ReDim U(UBound(Matrix, 1), UBound(Matrix, 1))
    ReDim Row(UBound(Matrix, 1))
    ReDim Column(UBound(Matrix, 1))

    'GaussElimination method:
    'A  x  = b
    'L U x = b
    'L  y  = b   (Lij = aij/ajj, current Matrix)
    'U  x  = y

    'initialize index-vectors, upper- & lower-diagonal matrix
    For i = 1 To UBound(Matrix, 1)
        'Initialize Index-Vectors
        Row(i) = i
        Column(i) = i
        For j = 1 To UBound(Matrix, 2)
            'initialize matrix U
            U(i, j) = Matrix(i, j)
            'initialize matrix L
            If i = j Then
                L(i, j) = 1
            Else
                L(i, j) = 0
            End If
        Next
    Next

    'get upper- & lower-diagonal matrix
    For k = 1 To UBound(Matrix, 1) - 1    'k: loop over all rows
        'Search for biggest entry in current row and column
        dblMax = Abs(U(Row(k), Column(k)))
        rowMax = k
        columnMax = k
        For i = k + 1 To UBound(Matrix, 1)
            If Abs(U(Row(i), Column(k))) > dblMax Then
                dblMax = Abs(U(Row(i), Column(k)))
                rowMax = i
                columnMax = k
            End If
            If Abs(U(Row(k), Column(i))) > dblMax Then
                dblMax = U(Row(k), Column(i))
                rowMax = k
                columnMax = i
            End If
        Next
        'swap rows and columns
        temp = Row(k)
        Row(k) = rowMax
        Row(rowMax) = temp
        temp = Column(k)
        Column(k) = columnMax
        Column(columnMax) = temp
        'get upper- & lower-diagonal matrix
        For i = k + 1 To UBound(Matrix, 1)   'i: row index
            'U(i,k) = most left non zero entry in line i
            Lij = U(Row(i), Column(k)) / U(Row(k), Column(k))
            For j = 1 To UBound(Matrix, 2)   '
                U(Row(i), Column(j)) = _
                            U(Row(i), Column(j)) - Lij * U(Row(k), Column(j))
            Next
            L(Row(i), Column(k)) = Lij
        Next
    Next
    'Gauss-Elimination finished

    'get y
    'L  y  = b
    '1 0 0       y   b
    'L 1 0 times y = b
    'L L 1       y   b
    ReDim y(UBound(Vektor, 1), 1)
    For i = 1 To UBound(y, 1)
        temp = 0
        'y(i, 1) = (Vektor(i, 1)
        For j = 1 To i - 1
            temp = temp + y(Row(j), 1) * L(Row(i), Column(j))
        Next
        y(Row(i), 1) = (Vektor(Row(i), 1) - temp)
    Next

    'get x
    'U  x  = y
    'u u u       x   y
    '0 u u times x = y
    '0 0 u       x   y
    ReDim x(UBound(Vektor, 1), 1)
    For i = UBound(x, 1) To 1 Step -1
        temp = 0
        For j = UBound(Vektor, 1) To i + 1 Step -1
            temp = temp + x(Row(j), 1) * U(Row(i), Column(j))
        Next
        x(Row(i), 1) = (y(Row(i), 1) - temp) / U(Row(i), Column(i))
    Next

    'return Result
    For i = 1 To UBound(Vektor, 1)
        'y is temporary vector
        y(i, 1) = x(Row(i), 1)
    Next
    GaussEliminationWithPivoting = y
End Function

Code eingefügt mit Syntaxhighlighter 4.15
Your rating: Keine Average: 3.8 (4 votes)