Eigenwerte und Eigenvektoren

Methoden

Motivation

Der Effekt einer linearen Abbildung, also die Multiplikation eines Vektors \(x\) mit einer Matrix \(A\), ist im Allgemeinen unübersichtlich, d. h. der Output \(Ax\) hängt vom Input auf eine komplizierte Weise ab. Wäre der Effekt eine reine Skalarmultiplikation \(\lambda x\), so wäre der Output eine Umsaklierung (Streckung, Stauchung, Umkehrung) des Inputs \(x\) und deutlich einfacher zu handhaben.

Anwendungen:

  • Entkoppeln von Gleichungen: Dynamische Systeme, Differentialgleichungen (-> Eigenfrequenzen, etc.)
  • Nicht-lineare Optimierung: Kriterium für lokales Maximum bzw. Minimum
  • Hauptachsentransformation, Principal Component Analysis
  • PageRank-Algorithmus von Google
  • Quantenmechanik
  • etc.

Definitionen

Ein Eigenvektor (EV) einer quadratischen \(n\times n\)-Matrix \(A\) ist ein nicht-Null-Vektor \(x \in \mathbb{R}^n\), sodass \[ Ax = \lambda x \] für eine Zahl \(\lambda\) gilt. Die Zahl \(\lambda\) ist der Eigenwert (EW) von \(A\) zum Eigenvektor \(x\).

Achtung: Jedes Vielfache eines Eigenvektors ist wieder Eigenvektor zum selben Eigenwert: \(A(\alpha x) = \alpha Ax = \alpha \lambda x = \lambda (\alpha x)\).

Spektralzerlegung

Wenn man eine Basis des \(\mathbb{R}^n\) aus Eigenvektoren \(v_i \in \mathbb{R}^n, i=1,\ldots,n\) bilden kann, dann heißt \(A\) diagonalisierbar und die Basis Eigenbasis. Man bildet die Matrix \(V\), die die Basis-Eigenvektoren als Spalten hat, und die Diagonalmatrix \(D\) mit den zugehörigen Eigenwerten \(\lambda_i\) auf der Diagonalen. Dann gilt \(AV = A (v_1 | v_2 | \ldots | v_n) = (\lambda_1 v_1 | \lambda_2 v_2 | \ldots | \lambda_n v_n) = VD\), also \[ AV = VD \] und, weil \(V\) als Basisvektoren-Matrix invertierbar ist, durch Umformen \(A = V D V^{-1}\) und \(V^{-1}AV = D\). Diese Formeln ergeben sich auch aus dem Basiswechsel (Koordinatentransformation) im In- und Outputraum der linearen Abbildung \(y = Ax\). Seinen \(x = Vc\) und \(y = Vd\) mit \(c\) und \(d\) die Koordinaten von \(x\) bzw. \(y\) bzgl. der Eigenbasis \(V\). Dann erhält man durch Einsetzen in \(y = Ax\) und Multiplikation mit \(V^{-1}\): \[ \begin{aligned} Vd & = AVc = VDc \\ d & = V^{-1}AVc = Dc. \end{aligned} \] Die Abbildung, die in Standardkoordinaten die Matrix \(A\) hat, hat bzgl. den angepassten Eigenkoordinaten die Diagonalmatrix \(D\).

Das Auffinden der EW und EV einer quadratischen Matrix \(A\) und die Diagonalisierung \(A = V D V^{-1}\) wird auch Spektralzerlegung von \(A\) genannt. Weitere wichtige Matrixzerlegungen sind

Berechnung

Falls \(x\) ein Eigenvektor von \(A\) zum Eigenwert \(\lambda\) ist, dann gelten: \[ \begin{aligned} Ax & = \lambda x \\ Ax - \lambda x & = 0 \\ Ax - \lambda I x & = 0 \\ (A - \lambda I) x & = 0. \end{aligned} \] Da \(x\) laut Annahme ein nicht-Null-Vektor ist, hat das letzte Gleichungssystem neben der trivialen Lösung des Nullvektors mit \(x\) eine zweite und somit nicht eindeutige Lösung. Daher muss die Determinante der Koeffizientenmatrix Null sein: \[ \det(A - \lambda I) = 0. \] Andernfalls hätte das letzte Gleichungssystem den Nullvektor als eindeutige Lösung. Aus der Bedingung \(\det(A - \lambda I) = 0\) lassen sich alle Eigenwerte berechnen, die auch komplex sein können. Anschließend kann zu jedem Eigenwert \(\lambda\) ein Eigenvektor berechnet werden, indem eine nicht-triviale Lösung von \((A - \lambda) x = 0\) bestimmt wird. Die Bedingung \(\det(A - \lambda I) = 0\) heißt charakteristische Gleichung von \(A\).

Beispiel: \(A = \begin{pmatrix} -2 & -5\\ 1 & 4 \end{pmatrix}\), \(\det(A - \lambda I) = \det\begin{pmatrix} -2 - \lambda & -5\\ 1 & 4 - \lambda \end{pmatrix} = (-2 - \lambda)(4 - \lambda) - (-5) = \lambda^2 - 2\lambda - 3 = 0\) hat die Lösungen \(\lambda_1 = -1\) und \(\lambda_2 = 3\). Das sind die Eigenwerte von \(A\).

  • Eigenwert \(\lambda_1 = -1\): \((A - \lambda_1 I) x = \begin{pmatrix} -1 & -5\\ 1 & 5 \end{pmatrix} x = 0\) hat unendlich viele Lösungen, z. B. \(v_1 = \begin{pmatrix} 5\\ -1 \end{pmatrix}\). \(v_1\) ist ein Eigenvektor von \(A\) zum Eigenwert \(\lambda_1 = -1\).
  • Eigenwert \(\lambda_2 = 3\): \((A - \lambda_2 I) x = \begin{pmatrix} -5 & -5\\ 1 & 1 \end{pmatrix} x = 0\) hat unendlich viele Lösungen, z. B. \(v_2 = \begin{pmatrix} 1\\ -1 \end{pmatrix}\). \(v_2\) ist ein Eigenvektor von \(A\) zum Eigenwert \(\lambda_2 = 3\).

EW und EV von speziellen Matrizen:

  • Für symmetrische Matrizen \(A\) gilt, dass die Eigenwerte reell sind und die Eigenvektoren zu verschiedenen Eigenwerten orthogonal sind. Aus den Eigenvektoren kann man daher eine Orthonormalbasis bilden. Die Matrix \(V\) der Orthonormalbasisvektoren ist dann orthogonal, d. h. \(V^{-1} = V^T\).
  • Die Eigenwerte einer Diagonalmatrix oder allgemeiner einer Dreiecksmatrix sind die Diagonaleinträge. Die Eigenvektoren sind die Standardbasisvektoren.

Definitheit

Eine quadratische \(n\times n\)-Matrix \(A\) heißt

  • positiv definit, falls \(x^{T}Ax > 0\)
  • positiv semidefinit, falls \(x^{T}Ax \geq 0\)
  • negativ definit, falls \(x^{T}Ax < 0\)
  • negativ semidefinit, falls \(x^{T}Ax \leq 0\)

für alle \(x \neq 0\) in \(\mathbb{R}^n\).

Eine quadratische, symmetrische Matrix ist genau dann

  • positiv definit, wenn alle Eigenwerte größer als null sind
  • positiv semidefinit, wenn alle Eigenwerte größer oder gleich null sind
  • negativ definit, wenn alle Eigenwerte kleiner als null sind
  • negativ semidefinit, wenn alle Eigenwerte kleiner oder gleich null sind und
  • indefinit, wenn positive und negative Eigenwerte existieren.

Die Definitheit der Hesse-Matrix spielt bei der Untersuchung von kritischen Stellen einer Funktion \(f\colon \mathbb{R} ^{n}\to \mathbb{R}\), also der Extremwertberechnung, eine entscheidende Rolle.

Python

Die zentrale Python-Funktion ist numpy.linalg.eig. Sie liefert die Eigenwerte und Eigenvektoren einer Matrix. Die Eigenwerte werden in einem Vektor und die Eigenvektoren als Matrix mit den zugehörigen Eigenvektoren als Spalten zurückgegeben. Die Eigenvektoren sind auf die Länge 1 normiert, d. h. sie sind Einheitsvektoren.

Für die Beispiele und Aufgaben verwenden wir folgende Python-Bibliotheken, siehe Kapitel Python Tutorial:

Code 
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

Beispiele

Berechnung von EW und EV

Für \(A = \begin{pmatrix} 0.8 & 0.3 \\ 0.2 & 0.7 \end{pmatrix}\) liefert die Bedingung \(\det(A - \lambda I) = 0\) die quadratische Gleichung \(\lambda^2 - 1.5\lambda + 0.5=0\). Deren Lösungen \(\lambda_1=1\) und \(\lambda_2=\frac{1}{2}\) sind die Eigenwerte von \(A\).

  • Die Eigenvektoren zum Eigenwert \(\lambda_1\) sind die Lösungen von \((A - \lambda_1 I)x = 0\), in unserem Beispiel ist das \(\begin{pmatrix} -0.2 & 0.3 \\ 0.2 & -0.3 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\). Die Eigenvektoren sind alle Vielfachen von z. B. \(\begin{pmatrix} 3 \\ 2 \end{pmatrix}\).
  • Die Eigenvektoren zum Eigenwert \(\lambda_2\) sind die Lösungen von \((A - \lambda_2 I)x = 0\), in unserem Beispiel ist das \(\begin{pmatrix} 0.3 & 0.3 \\ 0.2 & 0.2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}\). Die Eigenvektoren sind alle Vielfachen von z. B. \(\begin{pmatrix} 1 \\ -1 \end{pmatrix}\).
Code 
import numpy as np

A = np.array([[0.8, 0.3],
              [0.2, 0.7]])
L, V = np.linalg.eig(A)
print(f"eigenvalues: {L}")
print(f"matrix with eigenvectors in columns:\n{V}")
eigenvalues: [1.  0.5]
matrix with eigenvectors in columns:
[[ 0.83205029 -0.70710678]
 [ 0.5547002   0.70710678]]

Bevölkerungswachstum

Die Population einer Region zerlegt sich in die Altersgruppen

  • jünger als 20 Jahre
  • 20 bis 39 Jahre und
  • über 40 Jahre.

Dynamisches System: Anfangs ist der Generationenstand durch den Vektor \(x = \begin{pmatrix} x_1\\ x_2 \\ x_3 \end{pmatrix}\) gegeben. Eine Generation später ist die Population durch zwei Effekte verändert:

  • Reproduktion: \(x_1^{\text{neu}} = F_1 x_1 + F_2 x_2 + F_3 x_3\) mit den Fertilitätsraten \(F_k\).
  • Überlebensrate: \(x_2^{\text{neu}} = P_1 x_1\) and \(x_3^{\text{neu}} = P_2 x_2\) mit den Überlebenswahrscheinlichkeiten \(P_k\).

Die Leslie Matrix \(A\) liefert \(x^{\text{neu}}\) für \(x\) via \(x^{\text{neu}}= Ax\). Betrachten Sie folgendes Beispiel: \[ x^{\text{neu}} = \begin{pmatrix} x_1\\ x_2 \\ x_3 \end{pmatrix}^{\text{neu}} = \begin{pmatrix} F_1 & F_2 & F_3 \\ P_1 & 0 & 0 \\ 0 & P_2 & 0 \end{pmatrix}\begin{pmatrix} x_1\\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} .04 & 1.1 & .01 \\ .98 & 0 & 0 \\ 0 & .92 & 0 \end{pmatrix}\begin{pmatrix} x_1\\ x_2 \\ x_3 \end{pmatrix} \]

Aufgaben:

  1. Der Generationenstand sei zu Beginn durch \(x = \begin{pmatrix} .6\\ 1 \\ .8 \end{pmatrix}\) [Mio. Menschen] gegeben. Benutzen Sie Python, um die absoluten und relativen Populationsgrößen aller drei Gruppen für die nächsten 30 Generationen zu berechnen und zu plotten.
  2. Berechnen Sie die Eigenwerte und Eigenvektoren der Leslie Matrix mithilfe des Python Befehls eig, und interpretieren Sie das Resultat. Welcher Eigenwert ist weshalb für die Bevölkerungsentwicklung entscheidend? Wie kann ein zugehöriger Eigenvektor zu diesem Eigenwert zur Prognose der sich abzeichnenden relativen Populationsgrößen verwendet werden?

Lösung:

Code 
A = np.array([[0.04, 1.1 , 0.01],
              [0.98, 0   , 0   ],
              [0   , 0.92, 0   ]])
print(A)
(L, V) = np.linalg.eig(A)

for k in range(len(L)):
  print(f"Der Eigenwert {L[k]:.2f} hat den Eigenvektor\n {V[:,[k]]}")

x0 = np.array([[0.6, 1, 0.8]]).T
# x0 = -V[:,[0]]
# x0 = V[:,[1]]
# x0 = V[:,[2]]

N = 50
X = np.zeros((3, N))
X[:,[0]] = x0

for k in range(1, N):
    X[:,[k]] = A@X[:,[k-1]]

plt.figure(figsize=(6, 4))
plt.title('Absoluter Generationenstand')
plt.plot(X[0,:], 'g', label='Alter < 20')
plt.plot(X[1,:], 'r', label='Alter 20 bis 39')
plt.plot(X[2,:], 'k', label='Alter > 40')
plt.xlabel('Generation')
plt.ylabel('Generationenstand absolut')
plt.legend()
plt.grid(True)
[[0.04 1.1  0.01]
 [0.98 0.   0.  ]
 [0.   0.92 0.  ]]
Der Eigenwert 1.06 hat den Eigenvektor
 [[-0.63392495]
 [-0.58468167]
 [-0.50624747]]
Der Eigenwert -1.01 hat den Eigenvektor
 [[ 0.6083396 ]
 [-0.58784241]
 [ 0.53325813]]
Der Eigenwert -0.01 hat den Eigenvektor
 [[ 7.76398245e-05]
 [-9.09394695e-03]
 [ 9.99958646e-01]]

Jeder Anfangszustand \(x\) kann in der Eigenbasis darstellt werden: \(x = V y\), wobei \(y\) die Koordinaten von \(x\) in der Eigenbasis sind. Im Detail bedeutet das, dass \[ x = \sum_{i=1}^3 y_i v_i \] mit den Eigenvektoren \(v_i\) und den Koordinaten \(y_i\). Um den Zustand \(n\) Generationen später zu berechnen, wendet man die Leslie-Matrix \(A\) \(n\)-mal auf \(x\) an: \[ x^{(n)} = A^n x = \sum_{i=1}^3 y_i A^n v_i = \sum_{i=1}^3 y_i \lambda_i^n v_i. \] Der Zustand \(x^{(n)}\) ist also eine Linearkombination der Eigenvektoren \(v_i\) mit den Faktoren \(\lambda_i^n y_i\). Der Eigenvektor \(v_i\) mit dem größten Eigenwert \(\lambda_i\) bestimmt das langfristige Verhalten der Population.

Umskalieren des Eigenvektors mit größten Eigenwert:

Code 
E = V/np.sum(V, axis=0)
e1 = E[:,[0]]*100  # rescale the first eigenvector to percentage
print(f"rescaled the first eigenvector:\n{e1}")

plt.figure(figsize=(6, 4))
plt.title('Relative Generationenstände')
plt.plot(X[0,:]/np.sum(X, axis = 0)*100, 'g', label='Alter < 20')
plt.plot(X[1,:]/np.sum(X, axis = 0)*100, 'r', label='Alter 20 bis 39')
plt.plot(X[2,:]/np.sum(X, axis = 0)*100, 'k', label='Alter > 40')
plt.hlines(e1[0,0], 0, N, color='g')
plt.hlines(e1[1,0], 0, N, color='r')
plt.hlines(e1[2,0], 0, N, color='k')
plt.legend()
plt.xlabel('Generation')
plt.ylabel('Generationenstand in Prozent')
plt.grid(True)
rescaled the first eigenvector:
[[36.75238138]
 [33.89745683]
 [29.3501618 ]]

Entkoppeln von Differentialgleichungen

Wir betrachten das folgende lineare, homogene Differentialgleichungssystem erster Ordnung: \[ \begin{aligned} \dot{x}_1 (t) & = x_1 (t) + x_2 (t) \\ \dot{x}_2 (t) & = 4x_1 (t) + x_2 (t). \end{aligned} \] Man sagt, das System ist gekoppelt, weil \(\dot{x}_1 (t)\) auch von \(x_2 (t)\) und \(\dot{x}_2 (t)\) auch von \(x_1 (t)\) abhängt. Gesucht sind Lösungsfunktionen \(x_1(t)\) und \(x_2(t)\), die wir zu einem Vektor \(x(t) = \begin{pmatrix} x_1(t)\\ x_2(t) \end{pmatrix}\) zusammenfassen. Dadurch können wir das Differentialgleichungssystem in der kompakten Form \(\dot{x}(t) = Ax(t)\) schreiben, wobei \(A = \begin{pmatrix} 1 & 1\\ 4 & 1 \end{pmatrix}\). Mit Hilfe der Eigenwerte und Eigenvektoren von \(A\) können wir das Differentialgleichungssystem entkoppeln:

Code 
A = np.array([[1, 1],
              [4, 1]])
L, V = np.linalg.eig(A)
print(f"eigenvalues: {L}")
print(f"matrix with eigenvectors in columns:\n{V}")
print(f"rescaled eigenvectors:")
v_1 = V[:,[0]]/V[0, 0]
v_2 = V[:,[1]]/V[0, 1]
print(" v_1 = \n", v_1)
print(" v_2 = \n", v_2)
eigenvalues: [ 3. -1.]
matrix with eigenvectors in columns:
[[ 0.4472136  -0.4472136 ]
 [ 0.89442719  0.89442719]]
rescaled eigenvectors:
 v_1 = 
 [[1.]
 [2.]]
 v_2 = 
 [[ 1.]
 [-2.]]

Zur Entkopplung führen wir die Koordinatentransformation \(x = Vy\) durch, wobei die Matrix \(V\) die Eigenvektoren \(v_1\) und \(v_2\) als Spalten hat. Daher gilt \(AV = VD\), wobei \(D\) die Diagonalmatrix mit den Eigenwerten \(\lambda_1 = 3\) und \(\lambda_2 = -1\) auf der Diagonalen ist. Aus \(x = Vy\) folgt \(\dot{x} = V\dot{y}\), und \(\dot{x} = Ax\) wird durch Einsetzen zu \(V\dot{y} = AVy = VDy\), also \(\dot{y} = D y\), weil \(V\) invertierbar ist. Die letzte Gleichung \(\dot{y} = D y\) ist ein entkoppeltes System, d. h. die Gleichungen für \(y_1\) und \(y_2\) hängen nicht mehr voneinander ab und können separat gelöst werden: \[ \begin{aligned} \dot{y}_1 (t) & = 3y_1(t) \\ \dot{y}_2 (t) & = -y_2(t). \end{aligned} \] Die allgemeinen Lösungen sind \(y_1(t) = y_1(0) e^{3t}\) und \(y_2(t) = y_2 (0) e^{-t}\). In Vektorform: \(y(t) = \begin{pmatrix} y_1(t)\\ y_2(t) \end{pmatrix} = \begin{pmatrix} y_1(0) e^{3t} \\ y_2(0) e^{-t} \end{pmatrix}\). Die allgemeine Lösung für \(x(t)\) lautet nun durch Rücktransformation \[ x(t) = V y(t) = y_1(0) e^{3t} v_1 + y_2(0) e^{-t} v_2 = y_1(0) e^{3t} \begin{pmatrix} 1\\ 2 \end{pmatrix} + y_2(0) e^{-t} \begin{pmatrix} 1\\ -2 \end{pmatrix}. \] Wir sehen, dass der erste Summand exponentiell wächst, während der zweite exponentiell abfällt.

Wir lösen nun das Anfangswertproblem mit Anfangswert \(x(0) = \begin{pmatrix} 1 \\ -1 \end{pmatrix}\), indem wir die Anfangswerte \(y_1(0)\) und \(y_2(0)\) bestimmen und diese in die allgemeine Lösung einsetzen: \[ \begin{aligned} x(0) & = V y(0) \\ y(0) & = V^{-1} x(0) \end{aligned} \]

Code 
x0 = np.array([1, -1])
y0 = np.linalg.solve(V, x0)
print(f"y0 = {y0}")

t = np.linspace(0, 1, 100)
y = np.array([y0[0]*np.exp(3*t),
              y0[1]*np.exp(-t)])
x = V@y

plt.figure(figsize=(5, 3))
plt.plot(t, x[0,:], label='$x_1(t)$')
plt.plot(t, x[1,:], label='$x_2(t)$')
plt.xlabel('t')
plt.legend(loc='best')
plt.grid(True)
y0 = [ 0.55901699 -1.67705098]

Matrixexponential

Ein lineares Differentialgleichungssystem \(\dot{x} = Ax\) mit Anfangsbedingung \(x(0) = x_0\) kann auch direkt gelöst werden. Die Lösung lautet \(x(t) = e^{At}x_0\), wobei \(e^{At}\) ein Matrixexponential ist. Es ist analog zur Potenzreihe der Exponentialfunktion \[ e^x = \sum_{n=0}^\infty \frac{1}{n!} x^n \] für Zahlen \(x\) durch die Matrix-Potenzreihe \[ e^{At} = \sum_{n=0}^\infty \frac{1}{n!} A^n \] definiert. Zur effizienten Berechnung des Matrixexponentials verwendet man z. B. die Spektralzerlegung \(A = V D V^{-1}\): \[ \begin{aligned} e^{At} & = \sum_{n=0}^\infty \frac{1}{n!} A^n \\ & = \sum_{n=0}^\infty \frac{1}{n!} (V D V^{-1})^n \\ & = \sum_{n=0}^\infty \frac{1}{n!} V D^n V^{-1} \\ & = V \left(\sum_{n=0}^\infty \frac{1}{n!} D^n\right) V^{-1} \\ & = V e^{Dt} V^{-1} \end{aligned} \] Da die Matrix \(D\) eine Diagonalmatrix von Eigenwerten \(\lambda_i\) ist, ist \(e^{Dt}\) einfach zu berechnen: \[ e^{Dt} = \begin{pmatrix} e^{\lambda_1 t} & 0 & \ldots & 0\\ 0 & e^{\lambda_2 t} & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & e^{\lambda_n t} \end{pmatrix}. \] Mit dieser Alternative können wir das Anfangswertproblem des vorherigen Beispiels lösen. Dabei verwenden wir den SciPy-Befehl scipy.linalg.expm, der das Matrixexponential schon implementiert hat:

Code 
import scipy as sp

x = np.zeros((2, 100))
for k in range(100):
    x[:,k] = sp.linalg.expm(A*t[k])@x0

plt.figure(figsize=(5, 3))
plt.plot(t, x[0,:], label='$x_1(t)$')
plt.plot(t, x[1,:], label='$x_2(t)$')
plt.xlabel('t')
plt.legend(loc='best')
plt.grid(True)

Aufgaben

Verständnisfragen und kurze Aufgaben

  1. Beschreiben Sie die Rechenschritte zur Bestimmung der Eigenvektoren einer Matrix.
  2. Welche Eigenwerte hat die Matrix \(\begin{pmatrix} 1 & 3\\ 0 & -2 \end{pmatrix}\)?
  3. Welche Eigenwerte hat die Matrix \(\begin{pmatrix} 1 & -2\\ 1 & 4 \end{pmatrix}\)?

Ergebnisse:

  1. Siehe Abschnitt Berechnung
  2. \(\lambda_1 = 1, \lambda_2 = -2\)
  3. \(\lambda_1 = 2, \lambda_2 = 3\)

Eigenwerte und -vektoren

Berechnen Sie Eigenwerte und Eigenvektoren der Matrix \(\begin{pmatrix} 3 & 0.25 \\ 1 & 3 \end{pmatrix}\), und überprüfen Sie Ihr Ergebnis in Python.

Ergebnisse: Die Eigenwerte sind \(\lambda_1 = 3.5\) und \(\lambda_2 = 2.5\). Die Eigenvektoren sind z. B. \(v_1 = \begin{pmatrix} 1\\ 2 \end{pmatrix}\) und \(v_2 = \begin{pmatrix} 1\\ -2 \end{pmatrix}\).

Diskrete Dynamische Systeme, Eigenwerte und -vektoren

Jedes Jahr ziehen etwa 10 % der Bevölkerung einer Stadt in die umliegenden Vororte, und etwa 20 % der Vorstadtbevölkerung in die Stadt. Im Jahr 2016 gibt es 10 Millionen Menschen in der Stadt und 8 Millionen Menschen in den Vororten, was wir mit dem Zustandsvektor \(x_0 = \begin{pmatrix} 10 \\ 8 \end{pmatrix}\) modellieren.

  1. Bestimmen Sie die Übergangsmatrix \(M\), so dass \(Mx\) der Einwohnerverteilung ein Jahr nach \(x\) entspricht.
  2. Analysieren Sie das Langzeitverhalten des dynamischen Systems mit Hilfe der Eigenwerte und -vektoren von \(M\), d. h. wohin tendiert der Zustandsvektor für große Zeiten?

Ergebnisse: Die Übergangsmatrix ist \(M = \begin{pmatrix} 0.9 & 0.2\\ 0.1 & 0.8 \end{pmatrix}\). Der Zustand tendiert für große Zeiten gegen \(x = \begin{pmatrix} 12\\ 6 \end{pmatrix}\).

Ein Räuber-Beute-System

Tief in den Redwood-Wäldern von Kalifornien liefern die Dunkelfuß-Buschratten bis zu 80 % der Nahrung für den Fleckenkauz, den Hauptfeind der Buschratten. Wir verwenden ein lineares dynamisches System, um die Dynamik des Eulen- und Rattenbestands zu modellieren. Das Modell ist in mehrfacher Hinsicht unrealistisch, aber es kann als Ausgangspunkt für die Untersuchung komplizierterer nichtlinearer Modelle verwendet werden. Mit \(x_t = \begin{pmatrix} f_t \\ b_t \end{pmatrix}\) bezeichnen wir für den Monat \(t\) den Bestand des Fleckenkauzes sowie den Bestand an Buschratten gemessen in Tausenden. Ein Monat später sei der Bestand gegeben durch \[ \begin{aligned} f_{t+1} & = 0.5 f_t + 0.4 b_t, \\ b_{t+1} & = -0.104 f_t + 1.1 b_t. \end{aligned} \] Der Term \(0.5 f_t\) in der ersten Gleichung besagt, dass ohne Buschratten als Nahrung nur die Hälfte der Fleckenkauz jeden Monat überleben wird, während der Term \(1.1 b_t\) in der zweiten Gleichung besagt, dass ohne Fleckenkauz als Räuber die Buschrattenpopulation um 10 % pro Monat zunehmen wird. Wenn es viele Buschratten gibt, führt der Term \(0.4 b_t\) zu einem Anstieg der Fleckenkauzpopulation, während der negative Term \(-0,104 f_t\) den Tod von Buschratten durch den Fleckenkauz beschreibt. Ein Fleckenkauz frisst also pro Monat im Schnitt 104 Buschratten.

Wir beschreiben die Bestandsdynamik kann durch die Vektorgleichung \(x_{t+1} = M x_t\). Die Übergangsmatrix \(M\) lautet \[ M = \begin{pmatrix} 0.5 & 0.4 \\ -0.104 & 1.1 \end{pmatrix}. \]

  1. Bestimmen Sie in Python die Eigenwerte und Eigenvektoren von \(M\). Skalieren Sie die Eigenvektoren sinnvoll.
  2. Bestimmen Sie daraus den langfristigen Bestand der Eulen und Buschratten.
  3. Plotten Sie die Bestandsentwicklung der Eulen und Buschratten für die nächsten Monate für unterschiedliche Anfangsbedingungen in der \(f\)-\(b\)-Ebene.

Ergebnisse:

  1. Die Eigenwerte sind \(\lambda_1 = 0.58\) und \(\lambda_2 = 1.02\). Die Eigenvektoren sind z. B. \(v_1 = \begin{pmatrix} 1\\ 0.2 \end{pmatrix}\) und \(v_2 = \begin{pmatrix} 1\\ 1.3 \end{pmatrix}\).
  2. Langfristig gibt es pro Fleckenkauz 1300 Buschratten.

Eigenwerte und -vektoren, Invertierbarkeit, Eindeutigkeit der Lösung

Gegeben ist die Matrix \(A = \begin{pmatrix} 1 & 0 & 2 \\ 0 & 2 & 0 \\ 0 & 1 & 3 \end{pmatrix}\).

  1. Berechnen Sie von Hand und am Computer die Eigenwerte und Eigenvektoren von A.
  2. Ist die Matrix A invertierbar? Begründen Sie Ihre Antwort.
  3. Besitzt das Gleichungssystem \(Ax = b\) eine eindeutige Lösung? Begründen Sie Ihre Antwort.

Hinweis: \[ \det \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix} = a_{11} a_{22} a_{33} +a_{12} a_{23} a_{31} + a_{13} a_{21} a_{32} - a_{13} a_{22} a_{31} - a_{12} a_{21} a_{33} - a_{11} a_{23} a_{32}. \]

Ergebnisse:

  1. Eigenwerte von \(A\): \(\lambda_1 = 1, \lambda_2 = 2, \lambda_3 = 3\). Zugehörigen Eigenvektoren sind z. B. \(v_1 = (1, 0, 0)^T, v_2 = (1, 0, 1)^T, v_3 = (-2, 1, -1)^T\).
  2. Ja, denn die Determinante ist nicht Null.
  3. Ja, denn die Determinante ist nicht Null.

Eigenwerte und -vektoren

Bestimmen Sie von Hand und in Python alle Eigenwerte und -vektoren folgender Matrix. \[ \begin{pmatrix} 2 & 7 \\ 7 & 2 \end{pmatrix} \]

Ergebnisse: Eingenwerte von \(A\): \(\lambda_1 = -5, \lambda_2 = 9\). Zugehörige Eigenvektoren sind z. B. \(v_1 = (-1, 1)^T, v_2 = (1, 1)^T\).

Eigenwerte und -vektoren

Bestimmen Sie von Hand und in Python alle Eigenwerte und -vektoren folgender Matrix. \[ \begin{pmatrix} -2 & -5 \\ 1 & 4 \end{pmatrix} \] Quelle: Papula, Band 2, 7.2, S. 128 f.

Ergebnisse:

Eingenwerte von \(A\): \(\lambda_1 = -1, \lambda_2 = 3\). Zugehörige Eigenvektoren sind z. B. \(v_1 = (-5, 1)^T\), \(v_2 = (-1, 1)^T\).

Eigenwerte und -vektoren

Bestimmen Sie von Hand die Eigenwerte und Eigenvektoren von \(A = \begin{pmatrix} 5 & 1 \\ 4 & 2 \end{pmatrix}\).

Ergebnisse: \(\lambda_1 = 1, \lambda_2 = 6\), \(v_1 = (1, -4)^T\), \(v_2 = (1, 1)^T\).

Eigenwerte und -vektoren

Bestimmen Sie von Hand und am Computer die Eigenwerte und Eigenvektoren von \(A = \begin{pmatrix} -1 & 2 \\ 4 & 1 \end{pmatrix}\)?

Ergebnisse: \(\lambda_1 = 3, \lambda_2 = -3\), \(v_1 = (1, 2)^T\), \(v_2 = (1, -1)^T\).

Gekoppelter Radioaktiver Zerfall

Das radioaktive Isotop Jod zerfällt in das radioaktive Isotop Xenon, welches wiederum in etwas anderes (nicht wichtig) zerfällt. Die Halbwertszeiten von Jod und Xenon sind 6.7 h und 9.2 h, vgl. Halbwertszeit.

  1. Bestimmen Sie das System von DGL, das die Mengen \(J(t)\) von Jod und \(X(t)\) Xenon in Gramm zu jedem Zeitpunkt \(t\) in Stunden beschreibt.
  2. Bestimmen Sie die allgemeine Lösung des DGL-Systems, indem Sie das DGL-System in Matrix-Form schreiben und es über Eigenwerte und Eigenvektoren lösen.
  3. Bestimmen Sie die spezielle Lösung des DGL-Systems, wenn zu Beginn 100 g Jod und 0 g Xenon vorhanden sind.

Quelle: Farlow, Stanley J.: An Introduction to Differential Equations and Their Applications. p. 353, Problem 31.

Ergebnisse: Die allgemeine Lösung des DGL-Systems lautet: \[ \begin{pmatrix} J(t) \\ X(t) \end{pmatrix} = y_1(0) e^{-c_1 t} \begin{pmatrix} 1 \\ \frac{c_1}{c_2 - c_1} \end{pmatrix} + y_2 (0) e^{-c_2 t} \begin{pmatrix} 0 \\ 1 \end{pmatrix}. \] mit \(c_1 = \frac{\ln(2)}{6.7}\) und \(c_2 = \frac{\ln(2)}{9.2}\). Die spezielle Lösung des DGL-Systems hat den in Abbildung 1 dargestellten Verlauf.

Abbildung 1: Gekoppelter radioaktiver Zerfall

Romeo und Julia

Ein bisschen weit hergeholt könnte das folgende System von zwei linearen Differentialgleichungen \[ \begin{aligned} \frac{\text{d}r}{\text{d}t}(t) & = -aj(t), \quad (a>0) \\ \frac{\text{d}j}{\text{d}t}(t) & = br(t), \quad (b>0) \end{aligned} \] die Zuneigung zweier Liebender, Romeo und Julia, darstellen, wobei die Werte \(r(t)\) und \(j(t)\) ihre jeweilige Liebe für den anderen zum Zeitpunkt \(t\) messen. Man beachte, dass Romeo, je mehr Julia ihn liebt, anfängt, sie nicht mehr zu mögen, aber wenn sie das Interesse verliert, erwärmen sich seine Gefühle für sie. Sie hingegen macht die gegenteilige Erfahrung. Ihre Liebe zu ihm wächst, wenn er sie liebt, und verwandelt sich in Hass, wenn er sie hasst.

  1. Wie wird unter diesen Bedingungen die Lösung dieses Gleichungssystems qualitativ aussehen? Beschreiben Sie Ihre Meinung in Worten.

  2. Lösen Sie das DGL-System für den Fall \(a = b = 1\) und \(r(0) = j(0) = 2\) durch Entkoppeln in Python.

  3. Allgemeiner als bisher, könnte eine lineare Liebesbeziehung zwischen Romeo und Julia durch das folgende lineare DGL-System beschrieben werden: \[ \begin{aligned} \frac{\text{d}r}{\text{d}t}(t) & = a_{11}r(t) + a_{12}j(t) \\ \frac{\text{d}j}{\text{d}t}(t) & = a_{21}r(t) + a_{22}j(t) \end{aligned} \] Dabei können die “Liebe/Hass”-Parameter \(a_{ij}\) positiv oder negativ sein. Wenn z. B. \(a_{11} > 0\) und \(a_{12} > 0\) sind, dann wird Romeos Liebe zu Julia durch seine eigenen Gefühle und Julias Liebe zu ihm angespornt. In diesem Fall könnten wir Romeo als “Streberliebhaber” bezeichnen. Beschreiben Sie die folgenden anderen drei Liebhabertypen:

    1. \(a_{11} > 0\), \(a_{12} < 0\),
    2. \(a_{11} < 0\), \(a_{12} > 0\),
    3. \(a_{11} < 0\), \(a_{12} < 0\).

Ergebnisse:

  1. periodischer Wechsel zwischen Liebe und Hass
  2. Zwei phasenverschobene Schwingungen
  3. Beschreibung der Liebhabertypen:
    1. Romeo wird durch Julias Annäherungsversuche (\(a_{12} < 0\)) abgeturnt, aber immer noch von seinen eigenen Gefühlen (\(a_{11} > 0\)) angespornt.
    2. Romeo ist von Julias Liebe angeturnt (\(a_{12} > 0\)), aber er hat Angst vor seinen eigenen Gefühlen (\(a_{11} < 0\)).
    3. Romeo ist von Julias Annäherungsversuchen abgeschreckt (\(a_{12} < 0\)) und hat Angst vor der Liebe (\(a_{11} < 0\)).