vb@rchiv
VB Classic
VB.NET
ADO.NET
VBA
C#

https://www.vbarchiv.net
Rubrik: Variablen/Strings · Algorithmen/Mathematik   |   VB-Versionen: VB200824.07.08
Eigenwerte und -vektoren einer diagonalsymmetrischen Matrix

Der Jacobi-Algorithmus zur Bestimmung der Eigenwerte und Eigenvektoren - programmiert im 'modernen' Retro-Stil.

Autor:   Manfred BohnBewertung:  Views:  14.046 
ohne HomepageSystem:  Win2k, WinXP, Win7, Win8, Win10, Win11kein Beispielprojekt 

Das hier vorgestellte Verfahren zur Bestimmung der Eigenwerte und Eigenvektoren einer diagonal-symmetrischen Matrix ist bereits im Jahr 1846 von dem Mathematiker Jacobi veröffentlicht worden. Es gilt als das einfachste Verfahren zur 'Diagonalisierung' einer Matrix (iterativ: orthogonale Drehungen), das aber langsamer arbeitet als die Varianten des heute gebräuchlichen QR-Verfahrens.

Die Jacobi-Rotation ist ein Klassiker unter den Algorithmen und sollte in keiner Sammlung fehlen.

Auch aus der BASIC-Perspektive betrachtet ist diese Methode von besonderem Interesse.

Die ersten MS-Basic-Versionen sind vor ca. 30 Jahren als Beipack des Betriebssystems DOS ausgeliefert worden. Für den BASIC-Interpreter, den BASIC-Code und die Programm-Variablen standen insgesamt lediglich 64 KiloByte zur Verfügung. Von entscheidender Bedeutung war beim BASIC-Programmieren die Verfügbarkeit kurzer Algorithmen, die bei der Ausführung nur wenig Speicher belegten. Die Jacobi-Diagonalisierung gehörte deshalb zu den häufig eingesetzten Matrizen-Operationen.

Als Referenz an das klassische BASIC ist der Algorithmus - soweit möglich - im Stil der damaligen Zeit programmiert worden.

Wie allgemein üblich, werden die Eigenwerte nach Abschluss des Verfahrens (fallend) sortiert und die Eigenvektoren entsprechend angeordnet (vgl. Routine: Eigenwerte_Sortieren). Die Eigenvektoren sind normiert (Euklid-Norm) und reflektiert. Falls das Verfahren scheitert, stehen die Rückgabe-Arrays auf 'Nothing'. Das Modul 'modAdjust_GWB' wird benötigt, damit der Code unter Visual Basic ausführbar ist.
(Keine Sorge. Der Quellcode arbeitet auch unter den aktuellen Basic-Versionen korrekt. Und falls Sie die zahlreichen GOTOs irritieren sollten: eine ELSE-Bedingung gab es damals noch nicht. Es war auch üblich, den Datentyp 'Single' zu favorisieren, weil dadurch eine größere Anzahl von Werten im Speicher gehalten werden konnte. GW-Basic läuft auch unter Windows XP. Das Kürzel 'GW' steht für Graphics & Windows - ungelogen )

''' <summary>
''' Bestimmung der Eigenwerte und Eigenvektoren 
''' einer diagonalsymmetrischen Matrix
''' (Jacobi-Rotation)
''' </summary>
''' <param name="ds_mat">diagonalsymmetrische Matrix</param>
''' <param name="Eigenwerte">Rückgabe: 
''' Eigenwerte-Vektor</param>
''' <param name="EigenVektoren">Rückgabe: 
''' Eigenvektoren-Matrix</param>
Public Sub Jacobi(ByVal ds_mat(,) As Double, _
  ByRef Eigenwerte() As Double, _
  ByRef EigenVektoren(,) As Double)
 
  Eigenwerte = Nothing
  EigenVektoren = Nothing
  N% = UBound(ds_mat, 1)
  If Not IsDiagSym(ds_mat) Then Exit Sub
 
  Dim dsmat(,) As Double = CType(ds_mat.Clone, Double(,))
  Dim mat(N, N) As Double
 
  ' Klassischer VB-Code (GW-Basic)
 
  110: eps# = 0.00000000000001 REM Genauigkeitsschranke
  115: REM Initialisierung
  120:
  130: For j% = 0 To N
  140:   mat(j, j) = 1
  150: Next j
  155: REM Berechnung der Norm
  160: norm# = 0
  170: For i% = 1 To N
  180:   For j = 0 To i - 1
  190:     norm = norm + 2 * dsmat(i, j) ^ 2
  200:   Next j
  210: Next i
  300:
  310: norm = SQR(norm) : schranke# = eps * norm / (N + 1)
  315: itw# = norm : abb% = 0
  318:
  320: itw = itw / (N + 1)
  325:
  330: REM Jacobi-Iteration
  340: For q% = 1 To N
  350:   For p% = 0 To q - 1
  360:     If ABS(dsmat(p, q)) <= itw Then GoTo 650
  370:     abb = 1
  380:     REM Drehungen: Berechnung
  390:     el1# = dsmat(p, p)
  392:     el2# = dsmat(p, q)
  394:     el3# = dsmat(q, q)
  400:     mn# = (el1 - el3) / 2
  410:     If mn <> 0 Then GoTo 430
  420:     ro# = -1 : GoTo 440
  430:     ro = -SGN(mn) * el2 / SQR(el2 ^ 2 + mn ^ 2)
  440:     s1# = ro / SQR(2 * (1 + SQR(1 - ro / 2)))
  450:     s2# = s1 ^ 2 : c1# = SQR(1 - s2)
  460:     c2# = c1 ^ 2 : s3# = s1 * c1
  500:     REM Rotation
  510:     For i = 0 To N
  520:       help1 = dsmat(i, p) * c1 - dsmat(i, q) * s1
  530:       dsmat(i, q) = dsmat(i, p) * s1 + dsmat(i, q) * c1
  540:       dsmat(i, p) = help1
  541:       help1 = mat(i, p) * c1 - mat(i, q) * s1
  550:       mat(i, q) = mat(i, p) * s1 + mat(i, q) * c1
  560:       mat(i, p) = help1
  570:     Next i
  580:     For i = 0 To N
  590:       dsmat(p, i) = dsmat(i, p)
  595:       dsmat(q, i) = dsmat(i, q)
  600:     Next i
  610:     dsmat(p, p) = el1 * c2 + el3 * s2 - 2 * el2 * s3
  620:     dsmat(q, q) = el1 * s2 + el3 * c2 + 2 * el2 * s3
  630:     dsmat(p, q) = (el1 - el3) * s3 + el2 * (c2 - s2)
  640:     dsmat(q, p) = dsmat(p, q)
  650:   Next p
  660: Next q
  700: If abb <> 1 Then GoTo 720
  710: abb = 0 : GoTo 340
  720: If itw > schranke Then GoTo 320
 
  ' Rückgabe aufbereiten (aktuelles VB)
  ReDim Eigenwerte(N)
  For i = 0 To N
    Eigenwerte(i) = dsmat(i, i)
    CheckArg(Eigenwerte(i))
  Next i
 
  EigenVektoren = CType(mat.Clone, Double(,))
  Eigenwerte_Sortieren(Eigenwerte, EigenVektoren)
 
1000:
End Sub
Module modAdjust_GWB
  ' Anpassung der benötigten Routinen im Modul Math
  ' an die klassischen VB-Aufrufe
  Public Function CheckArg(ByVal arg As Single)
    ' Früher war das nicht nötig
    If Single.IsNaN(arg) Then
      Throw New ArithmeticException
    End If
    If Single.IsNegativeInfinity(arg) Or _
      Single.IsPositiveInfinity(arg) Then
      Throw New ArgumentOutOfRangeException
    End If
    Return True
  End Function
  Public Function SQR(ByVal arg As Single) As Single
    CheckArg(arg)
    Return Math.Sqrt(arg)
  End Function
  Public Function ABS(ByVal arg As Single) As Single
    CheckArg(arg)
    Return Math.Abs(arg)
  End Function
  Public Function SGN(ByVal arg As Single) As Single
    CheckArg(arg)
    Return Math.Sign(arg)
  End Function
  Public Sub Eigenwerte_Sortieren( _
    ByRef Eigenwerte() As Double, _
    ByRef EigenVektoren(,) As Double)
 
    ' Eigenwerte fallend sortieren
    ' Eigenvektoren entsprechend umsortieren 
 
    N = UBound(Eigenwerte)
 
    For i As Integer = 0 To N - 1
      For k As Integer = i + 1 To N
        If Eigenwerte(i) < Eigenwerte(k) Then
          Swap(Eigenwerte(i), Eigenwerte(k))
          For l% = 0 To N
            Swap(EigenVektoren(l, i), EigenVektoren(l, k))
          Next l
        End If
      Next k
    Next i
 
    For i% = 0 To N
      minus# = 0
      For l% = 0 To N
        If EigenVektoren(l, i) < 0.0 Then
          minus += 1
        End If
      Next l
      If minus > (N + 1) \ 2 Then
        ' Eigenvektor zu Eigenwert 'i' reflektieren
        For l% = 0 To N
          EigenVektoren(l, i) *= -1
        Next l
      End If
    Next i
  End Sub
  Private Sub Swap(ByRef a As Double, ByRef b As Double)
    Dim c As Double = a
    a = b : b = c
  End Sub
  Public Function IsDiagSym(ByVal mat(,) As Double) As Boolean
    ' Ist die Matrix diagonalsymmetrisch?
    Dim eps As Double = 0.00001
    N = UBound(mat, 1)
    If UBound(mat, 2) <> N Then Return False
 
    For i As Integer = 0 To N - 1
      For k As Integer = i + 1 To N
        If Math.Abs(mat(i, k) - mat(k, i)) > eps Then Return False
      Next k
    Next i
    Return True
  End Function
End Module



Anzeige

Kauftipp Unser Dauerbrenner!Diesen und auch alle anderen Tipps & Tricks finden Sie auch auf unserer aktuellen vb@rchiv  Vol.6

Ein absolutes Muss - Geballtes Wissen aus mehr als 8 Jahren vb@rchiv!
- nahezu alle Tipps & Tricks und Workshops mit Beispielprojekten
- Symbol-Galerie mit mehr als 3.200 Icons im modernen Look
Weitere Infos - 4 Entwickler-Vollversionen (u.a. sevFTP für .NET), Online-Update-Funktion u.v.m.
 
 
Copyright ©2000-2024 vb@rchiv Dieter OtterAlle Rechte vorbehalten.


Microsoft, Windows und Visual Basic sind entweder eingetragene Marken oder Marken der Microsoft Corporation in den USA und/oder anderen Ländern. Weitere auf dieser Homepage aufgeführten Produkt- und Firmennamen können geschützte Marken ihrer jeweiligen Inhaber sein.