Aufgaben
Lösung zu Aufgabe 10.12
Um das $\textrm{CO}_2$-Molekül zu simulieren, muss man zunächst die Geometrie festlegen. Da wir keine Stützpunkte haben, erzeugen wir für die Stützpunkte eine leere Liste. Auf die absoluten Abstände der Atome voneinander kommt es nicht an, sodass wir hier willkürlich eine Einheit wählen. Anstelle der Steifigkeiten der Stäbe müssen wir nun Federkonstanten angeben. Dementsprechend muss man die Zeile für die Berechnung der Matrix $\hat A$ anpassen, da man hier nun nicht mehr durch die Länge des Stabes dividieren darf. Für die Ausgabe der Frequenz eignet sich die Einheit THz (1 THz = $10^{12}\,\textrm{Hz}$).
Man erkennt in der Animation 6 Eigenmoden von denen die ersten 4 eine Eigenfrequenz von Null haben. Diese haben die folgenden Bedeutungen:
- Die Mode $f_1$ ist eine Starrkörperverschiebung entlang der x-Achse.
- Die Summe der Moden $f_2$, $f_3$ und $f_4$ ist eine Starrkörperverschiebung entlang der y-Achse.
- Die Differenz der Moden $f_2$ und $f_4$ ist eine Starrkörperdrehung des Moleküls in der Bildebene.
- Die Mode $f_3$ beschreibt eine Biegemode des Moleküls. Diese hat hier die Frequenz Null, da die Federn in unserem Modell keinerlei Biegesteifigkeit haben.
Es ist interessant zu bemerken, dass das Verhältnis der beiden Schwingungsfrequenzen $f_6$ zu $f_5$ in dieser Berechnung ungefähr $1{,}92$ ist, während experimentell Messungen einen Wert von $1{,}76$ ergeben. Die Ursache ist in einer sogenannten Fermi-Resonanz mit der Biegemode des Moleküls zu sehen, auf die im Rahmen dieses Buches aber nicht weiter eingegangen werden kann.
Schwingungen/Loesungen/eigenmoden_co2.py
"""Eigenmoden eines CO2-Moleküls in linearer Näherung."""
import numpy as np
import matplotlib as mpl
import matplotlib.animation
import matplotlib.pyplot as plt
# Lege die Positionen der Punkte fest [m].
punkte = np.array([[-1.0, 0], [0.0, 0], [1.0, 0]])
# Erzeuge eine Liste mit den Indizes der Stützpunkte.
indizes_stuetz = []
# Jeder Stab verbindet genau zwei Punkte. Wir legen dazu die
# Indizes der zugehörigen Punkte in einem Array ab.
staebe = np.array([[0, 1], [1, 2]])
# Federkonstante der Verbindungen [N/m].
federkonstanten = np.array([1.5e3, 1.5e3])
# Massen der einzelnen Punkte [kg].
knotenmassen = np.array([16.0, 12.0, 16.0]) * 1.6605e-27
# Amplitude, mit der die Eigenmoden dargestellt werden [m].
amplitude = 0.3
# Definiere die Dimension sowie die Anzahl der Punkte, Stäbe, etc.
n_punkte, n_dim = punkte.shape
n_staebe = len(staebe)
n_stuetz = len(indizes_stuetz)
n_knoten = n_punkte - n_stuetz
# Erzeuge eine Liste mit den Indizes der Knoten.
indizes_knoten = list(set(range(n_punkte)) - set(indizes_stuetz))
def ev(i_pkt, i_stb, koord=punkte):
"""Bestimme den Einheitsvektor in einem Punkt für einen Stab.
Args:
i_pkt (int):
Index des betrachteten Punktes.
i_stb (int):
Index des betrachteten Stabes.
koord (np.ndarray):
Koordinaten der Punkte (n_punkte × n_dim).
Returns:
np.ndarray: Berechneter Einheitsvektor oder der Nullvektor,
wenn der Stab den Punkt nicht enthält.
"""
stb = staebe[i_stb]
if i_pkt not in stb:
return np.zeros(n_dim)
if i_pkt == stb[0]:
vektor = koord[stb[1]] - koord[i_pkt]
else:
vektor = koord[stb[0]] - koord[i_pkt]
return vektor / np.linalg.norm(vektor)
def laenge(i_stb, koord=punkte):
"""Berechne die Länge eines Stabes.
Args:
i_stb (int):
Index des betrachteten Stabes.
koord (np.ndarray):
Koordinaten der Punkte (n_punkte × n_dim).
Returns:
float: Länge des Stabes.
"""
i1, i2 = staebe[i_stb]
return np.linalg.norm(koord[i2] - koord[i1])
# Stelle das Gleichungssystem für die Kräfte auf.
A = np.zeros((n_knoten, n_dim, n_knoten, n_dim))
for n, k in enumerate(indizes_knoten):
for m, j in enumerate(indizes_knoten):
for i in range(n_staebe):
A[n, :, m, :] -= (federkonstanten[i]
* np.outer(ev(k, i), ev(j, i)))
A = A.reshape((n_knoten * n_dim, n_knoten * n_dim))
# Erzeuge ein Array, das die Masse für jede Koordinate
# der Knotenpunkte enthält.
massen = np.repeat(knotenmassen, n_dim)
# Berechne die Matrix Lambda.
Lambda = -A / massen.reshape(-1, 1)
# Bestimme die Eigenwerte und die Eigenvektoren.
eigenwerte, eigenvektoren = np.linalg.eig(Lambda)
# Eigentlich sollten alle Eigenwerte reell sein.
if np.any(np.iscomplex(eigenwerte)):
print('Achtung: Einige Eigenwerte sind komplex.')
print('Der Imaginärteil wird ignoriert')
eigenwerte = np.real(eigenwerte)
eigenvektoren = np.real(eigenvektoren)
# Eigentlich sollte es keine negativen Eigenwerte geben.
eigenwerte[eigenwerte < 0] = 0
# Sortiere die Eigenmoden nach aufsteigender Frequenz.
indizes_sortiere_eigenwerte = np.argsort(eigenwerte)
eigenwerte = eigenwerte[indizes_sortiere_eigenwerte]
eigenvektoren = eigenvektoren[:, indizes_sortiere_eigenwerte]
# Berechne die Eigenfrequenzen.
eigenfrequenzen = np.sqrt(eigenwerte) / (2 * np.pi)
# Erzeuge eine Figure.
fig = plt.figure()
fig.set_tight_layout(True)
# Anzahl der darzustellenden Eigenmoden.
n_moden = eigenfrequenzen.size
# Erzeuge ein geeignetes n_zeilen × n_spalten - Raster.
n_zeilen = int(np.sqrt(n_moden))
n_spalten = n_moden // n_zeilen
while n_zeilen * n_spalten < n_moden:
n_spalten += 1
# Erzeuge eine Liste der animierten Grafikobjekte jeder Eigenmode.
plots = []
# Erstelle die Plots für jede Eigenmode in einer eigenen Axes.
for mode in range(n_moden):
ax = fig.add_subplot(n_zeilen, n_spalten, mode + 1)
ax.set_title(
f'$f_{{{mode+1}}}$={eigenfrequenzen[mode] / 1e12:.1f} THz')
ax.set_xlim(-1.5, 1.5)
ax.set_ylim(-2, 2)
ax.set_aspect('equal')
ax.axis('off')
# Erzeuge ein Dictionary für die animierten Plot-Objekte
# dieser Mode und hänge dieses an die Liste plots an.
plot_objekte = {}
plots.append(plot_objekte)
# Lege einen Plot für die Knotenpunkte in Blau an.
plot_objekte['knoten'], = ax.plot([], [], 'bo', zorder=5)
# Lege Plots für die Stäbe an.
plot_objekte['staebe'] = []
for stab in staebe:
plot_stab, = ax.plot([], [], color='black', zorder=4)
plot_objekte['staebe'].append(plot_stab)
# Plotte die Stützpunkte in Rot.
ax.plot(punkte[indizes_stuetz, 0], punkte[indizes_stuetz, 1],
'ro', zorder=5)
# Plotte die Ausgangslage der Knotenpunkte hellblau.
ax.plot(punkte[indizes_knoten, 0], punkte[indizes_knoten, 1],
'o', color='lightblue', zorder=2)
# Plotte die Ausgangslage der Stäbe Hellgrau.
for stab in staebe:
ax.plot(punkte[stab, 0], punkte[stab, 1],
color='lightgray', zorder=1)
# Zeitachse, die 60 Punkte im Bereich von 0 ... 2 pi enthält.
t = np.radians(np.arange(0, 360, 6))
def update(n):
"""Aktualisiere die Grafik zum n-ten Zeitschritt."""
for mode in range(n_moden):
# Stelle den zu dieser Mode gehörenden Eigenvektor
# als ein n_knoten × dim - Array dar.
ev = eigenvektoren[:, mode].reshape(n_knoten, n_dim)
# Berechne die aktuellen Positionen p aller Punkte.
p = punkte.copy()
p[indizes_knoten] += amplitude * np.sin(t[n]) * ev
# Aktualisiere die Positionen der Knotenpunkte.
plots[mode]['knoten'].set_data(p[indizes_knoten].T)
# Aktualisiere die Koordinaten der Stäbe.
for linie, stab in zip(plots[mode]['staebe'], staebe):
linie.set_data(p[stab, 0], p[stab, 1])
# Gib eine Liste aller geänderten Objekte zurück.
geaendert = []
for p in plots:
geaendert.append(p['knoten'])
geaendert += p['staebe']
return geaendert
# Erzeuge das Animationsobjekt und starte die Animation.
ani = mpl.animation.FuncAnimation(fig, update, frames=t.size,
interval=30, blit=True)
plt.show()