Ing Mathematik: Numerisches Lösen von Eigenwertproblemen

Zurück zu Numerisches Lösen nichtlinearer Gleichungen |

Vor zu Numerisches Lösen von Ausgleichsproblemen

Insbesondere bei Schwingungen sind Eigenwerte und -vektoren von Bedeutung (siehe z.B. Schwingung oder Hagedorn, Hochlenert: Technische Schwingungslehre. Harri Deutsch, 2012, ISBN 978-3-8171-1890-8). Aber auch die Hauptspannungen in der Festigkeitslehre können als Eigenwertproblem behandelt werden (siehe z.B. Mechanische_Spannung#Hauptspannung_und_Hauptspannungsrichtung und Hagedorn, Wallaschek: Technische Mechanik Bd.2, Festigkeitslehre. Europa, 5.Aufl., 2015, ISBN 978-3-8085-5691-7).

Eigenwerte und Eigenvektoren - Allgemeines

Es sei

$Ax=\lambda x$ mit $A\in \mathbb {R} ^{n\times n},x\in \mathbb {R} ^{n},\lambda \in \mathbb {R}$ .

$Ax-\lambda Ix=0$ mit $I$ als nxn-Einheitsmatrix.

$(A-\lambda I)x=0$

Für nicht triviale Lösungen ( $x\neq 0$ ) muss gelten:

$\det(A-\lambda I)=0$ (charakteristische Gleichung oder Säkulargleichung).

$\lambda ^{n}+\alpha _{n-1}\cdot \lambda ^{n-1}+\dotsb +\alpha _{1}\cdot \lambda +\alpha _{0}=0$ . Die Lösungen dieses Polynoms heißen Eigenwerte.

$\alpha _{0}=\det A;\quad \alpha _{n-1}={\text{Spur}}A;\quad \alpha _{n}=1$

Ein Polynom vom Grad n hat maximal n Nullstellen. D.h. es gibt auch höchstens n Eigenwerte. Die zu $\lambda _{i}$ gehörenden Vektoren $x_{i}\neq 0$ heißen Eigenvektoren.

Siehe auch Eigenwerte und Eigenvektoren.

Ein Beispiel mit Python und SciPy:

import numpy as np
import scipy.linalg as la

A = np.array([[ 24,  54, -38,  -8],
              [-11, -27,  20,  -2],
              [  0,  -2,   3,  -6],
              [  6,  14,  -10, -1]])

w, v = la.eig(A)

# Eigenwerte
print(w)
# Eigenvektoren
print(v)

Ausgabe:

[-3.+0.j  2.+0.j  1.+0.j -1.+0.j]
[[-9.05821627e-01  4.08248290e-01 -9.25820100e-01  9.17662935e-01]
 [ 3.39683110e-01  4.08248290e-01  1.54303350e-01 -2.29415734e-01]
 [-1.13227703e-01  8.16496581e-01 -3.08606700e-01  2.29415734e-01]
 [-2.26455407e-01  1.93235474e-14 -1.54303350e-01  2.29415734e-01]]

Siehe auch [1]

Spektrum einer Matrix

Die Menge der Eigenwerte nennt man das Spektrum $\sigma (A)$ . Der Spektralradius ist $\rho =max(|\lambda _{i}|)$ .

Vielfachheiten

Algebraische Vielfachheit: $\kappa _{i}$ (k-facher Eigenwert von $A$ )

{\text{Spur}}A=\kappa _{1}\lambda _{1}+\dots \kappa _{r}\gamma _{r}

\det A=\lambda _{1}^{\kappa _{1}}\dots \lambda _{r}^{\gamma _{r}}

Geometrische Vielfachheit: $\gamma _{i}={\text{dim Kern}}(A-\lambda _{i}I)=n-{\text{Rang}}(A-\lambda _{i}I)$

Mehrfache Nullstellen werden entsprechend ihrer Vielfachheit gezählt.

Symmetrische Matrizen

Für reelle, symmetrische Matrizen gilt:

Eigenwerte sind reell
Eigenvektoren stehen senkrecht aufeinander
geometrische und algebraische Vielfachheit stimmen bei jedem Eigenwert überein.

Seien

p ... Anzahl der positiven Eigenwerte der symmetrischen Matrix $S$
q ... Anzahl der negativen Eigenwerte der symmetrischen Matrix $S$
d ... Anzahl der Eigenwerte die Null sind.

Dann sind

(p, q) ... Signatur von $S$
p-q ... Trägheitsindex von $S$
d ... Defekt von $S$ .

Eine symmetrische Matrix ist positiv definit, wenn alle ihre Eigenwerte > 0 sind.

Siehe auch Defekt (Mathematik).

Diagonalisierbarkeit

Eine Matrix $A$ lässt sich diagonalisieren, wenn sie n linear unabhängige Eigenvektoren $x_{1},\dots ,x_{n}$ aufweist.

$C=(x_{1},\dots ,x_{n})$

$C^{-1}AC={\text{diag}}(\lambda _{1},\dots ,\lambda _{n})$

Jordansche Normalform

Nicht jede Matrix ist diagonalisierbar. Aber jede quadratisch komplexe Matrix $A$ lässt sich auf die jordansche Normalform $J$ bringen.

$T^{-1}AT=J={\begin{pmatrix}J_{1}&&0\\&\ddots &\\0&&J_{k}\end{pmatrix}}$

mit $T$ regulär komplex, $J_{j}={\begin{pmatrix}\lambda _{j}&1&&&0\\&\lambda _{j}&1&&\\&&\ddots {}&\ddots {}\\&&&\lambda _{j}&1\\0&&&&\lambda _{j}\end{pmatrix}}\in \mathbb {C} ^{s_{j}\times s_{j}}$

Wie man die Transformation auf die jordansche Normalform durchführt (Berechnung von $J,T$ aus gegebenen $A$ ), siehe z.B. Burg, Haf, Wille, Meister: Seite 237ff. Dies ist einigermaßen aufwendig. Umgekehrt ist das natürlich sehr einfach.

Shift - Verschieben von Eigenwerten

Sind $\lambda _{i}$ die Eigenwerte der nxn-Matrix $A$ , so hat die Matrix $A_{\epsilon }=A+\epsilon I$ die Eigenwerte $\mu _{i}=\lambda _{i}+\epsilon$ . $\lambda _{i}$ und die $\mu _{i}$ haben gleiche algebraische Vielfachheit.

Hessenbergmatrix

Eine (obere) Hessenbergmatrix ist eine quadratische Matrix $H\in \mathbb {C} ^{n\times n}$ , deren Einträge unterhalb der ersten Nebendiagonalen gleich Null sind, also $h_{ij}=0$ für alle $i>j+1$ .

$H={\begin{pmatrix}h_{11}&h_{12}&h_{13}&\cdots &h_{1n}\\h_{21}&h_{22}&h_{23}&\cdots &h_{2n}\\0&h_{32}&h_{33}&\cdots &h_{3n}\\\vdots &\ddots &\ddots &\ddots &\vdots \\0&\cdots &0&h_{nn-1}&h_{nn}\end{pmatrix}}$

Siehe auch Hessenbergmatrix.

Krylov-Verfahren

Das Krylov-Verfahren eignet sich nur für kleine Matrizen ( $n\leq 8$ ) (Quelle: Burg, Haf, Wille, Meister: Seite 249). Die Matrix $A$ habe nur unterschiedliche Eigenwerte. Das Verfahren ist ziemlich instabil. Es wird deshalb in der Praxis nur selten eingesetzt.

Gegeben sei eine 4x4-Matrix $A$ . Startvektor sei $z_{0}=[1,1,1,1]$ . Funktioniert das nicht, fährt man mit $z_{0}=[-1,1,1,1]$ fort usw. Funktioniert das alles nicht, so lässt sich das Problem mit dem Krylov-Verfahren nicht lösen.

Python-Code (rudimentär):

import numpy as np
import scipy.linalg as la
import scipy.optimize as opt
import matplotlib.pyplot as plt

# 4x4-Matrix
A = np.array([[3, 1, -10, 5],
              [1, -7, 1, 2],
              [-10, 1, 3, -6],
              [5, 2,  -6, -1]])

n = A.shape[0]
z = np.zeros((n,n))
zn = np.zeros(n)

# Koeffizientenmatrix für lin. GS erstellen
z[:,0] = [ 1, 1, 1, 1]

for k in np.arange(1, n):
    z[:,k] = A @ z[:,k-1]

zn = A @ z[:,3]

# Löse lin. GS
x = la.solve(z, -(-1)**4 * zn)

# Polynom erstellen
p = np.poly1d([1, x[3], x[2], x[1], x[0]])

# Nullwerte ermitteln
eigenwert = np.zeros(n)

for k in np.arange(0, n):
    if k == n-1:
        i = 15
    else:
        i = -10 + k*4
    eigenwert[k] = opt.newton(p, i) 

# Konsolenausgabe
for k in np.arange(0, n):
    print(eigenwert[k])

# Grafik
x = np.arange(-12, 20, .1)
y = p(x)
x1 = np.arange(-10, 0, .1)
y1 = p(x1)

plt.plot(x, y)
for k in np.arange(0, n):
    plt.plot(eigenwert[k], 0, "o")

plt.axes([0.3, 0.5, 0.25, 0.25])
plt.plot(x1, y1)
for k in np.arange(0, n-1):
    plt.plot(eigenwert[k], 0, "o") 

plt.grid()
plt.show()

Konsolenausgabe:

-9.125090927816691
-5.8570090791499965
-3.5129666875859895
16.49506669455273

Grafik:

Das charakteristische Polynom für dieses Beispiel lautet: $\chi (x)=x^{4}+2x^{3}-199x^{2}-1562x-3097=0$

Jacobi-Verfahren

Dieses Verfahren eignet sich für kleine reelle symmetrische Matrizen. Es arbeitet mit Givens-Rotationen:

$A^{(0)}=A$

$R_{k}={\begin{bmatrix}1&\cdots &0&\cdots &0&\cdots &0\\\vdots &\ddots &\vdots &&\vdots &&\vdots \\0&\cdots &\cos \varphi &\cdots &\sin \varphi &\cdots &0\\\vdots &&\vdots &\ddots &\vdots &&\vdots \\0&\cdots &-\sin \varphi &\cdots &\cos \varphi &\cdots &0\\\vdots &&\vdots &&\vdots &\ddots &\vdots \\0&\cdots &0&\cdots &0&\cdots &1\end{bmatrix}},$

Damit folgt wieder die symmetrische Matrix $A^{(1)}=R_{k}^{T}AR_{k}$ . Führt man dies immer weiter durch, so ergibt sich eine Diagonalmatrix.

$A^{(k+1)}=R_{k}^{T}A^{(k)}R_{k}=\underbrace {R_{k}^{T}R_{k-1}^{T}\ldots R_{0}^{T}} _{T_{k}^{T}}A^{(0)}\underbrace {R_{0}\ldots R_{k-1}R_{k}} _{T_{k}}={\text{diag}}(\lambda _{1},\dots ,\lambda _{n})$

Siehe auch Jacobi-Verfahren (Eigenwerte), Burg, Haf, Wille, Meister: Seite 251ff, Hanke-Bourgeois: Seite 238ff oder [2].

Potenzmethode nach von Mises

Sie wird auch Vektoriteration oder Von-Mises-Iteration genannt und berechnet den betragsgrößten Eigenwert. Sie eignet sich besonders für dünnbesetzte Matrizen. Das Verfahren lässt sich den Krylow-Unterraum-Verfahren zuordnen.

Python-Code:

import numpy as np

A = np.array([[3, 1, -10, 5],
              [1, -7, 1, 2],
              [-10, 1, 3, -6],
              [5, 2,  -6, -1]]) 

r = [1,1,1,1]

for k in np.arange(0,50):
    r = A@r
    lam = np.max(np.absolute(r))
    r = r/lam

print(lam)

Ausgabe:

16.495066694547337

Siehe auch Potenzmethode, Burg, Haf, Wille, Meister: Seite 254ff, Hanke-Bourgeois: Seite 218ff

Deflation

Siehe vorerst Deflation (Mathematik) und Burg, Haf, Wille, Meister: Seite 257f.

Wielandt-Iteration

Helmut Wielandt (deutscher Mathematiker, 1910-2001)

Siehe vorerst Inverse Iteration und Burg, Haf, Wille, Meister: Seite 258f.

Lanczos-Verfahren

Siehe vorerst Lanczos-Verfahren und Hanke-Bourgeois: Seite 259ff

Weitere Verfahren

QR-Verfahren basierend auf der Hessenberg-Matrix (siehe auch QR-Algorithmus, und Hanke-Bourgeois: Seite 227ff)
Hyman-Verfahren

Eigenwerte und -vektoren mit SciPy

Mit SciPy und Python lassen sich die Eigenwerte und -vektoren sehr einfach berechnen. Dazu gibt es eine Reihe von SciPy-Funktionen. Nachfolgend wird die Verwendung einiger dieser Funktionen gezeigt. Siehe dazu auch die SciPy-Dokumentation [3].

Python-Code (gegeben sei eine symmetrische Matrix):

import numpy as np
from scipy import linalg

a = np.array([[3, 1, -10, 5],
              [1, -7, 1, 2],
              [-10, 1, 3, -6],
              [5, 2,  -6, -1]])

ev1 = linalg.eigvals(a)   # Eigenwerte für allgemeine Matrizen
ev2 = linalg.eigvalsh(a)  # Eigenwerte für symmetrische/hermitesche Matrizen
w, vl = linalg.eig(a)     # Eigenwerte und Eigenvektoren

print(ev1)
print("-"*75)
print(ev2)
print("-"*75)
print(w)
print("-"*75)
print(vl)

Ausgabe:

[16.49506669+0.j -9.12509093+0.j -5.85700908+0.j -3.51296669+0.j]
---------------------------------------------------------------------------
[-9.12509093 -5.85700908 -3.51296669 16.49506669]
---------------------------------------------------------------------------
[16.49506669+0.j -9.12509093+0.j -5.85700908+0.j -3.51296669+0.j]
---------------------------------------------------------------------------
[ 0.63731292  0.33339446  0.64983438 -0.24575529]
[ 0.03424931 -0.71822344  0.51243339  0.4694615 ]
[-0.65183109  0.5045163   0.48861626  0.28606588]
[ 0.40960403  0.34418968 -0.27637196  0.79835866]]

Gedruckte Werke (auszugsweise)

Burg, Haf, Wille, Meister: Höhere Mathematik für Ingenieure, Band II: Lineare Algebra. 7. Auflage, Springer Vieweg, 2012, ISBN 978-3-8348-1853-9
Hanke-Bourgeois: Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens. 3. Aufl., Vieweg+Teubner, 2009, ISBN 978-3-8348-0708-3

Zurück zu Numerisches Lösen nichtlinearer Gleichungen |

Hoch zu Gesamtinhaltsverzeichnis |

Vor zu Numerisches Lösen von Ausgleichsproblemen