Diese Seite bezieht sich auf die erste Auflage (2020) des Buches.   ➜   Zur aktuellen Auflage.

Aufgaben

Lösung zu Aufgabe 9.1

Um die analytische Lösung für die Audioausgabe einer gedämpften Schwingung zu verwenden, berechnen wir zunächst den Abklingkoeffizienten $\delta$ und die Eigenkreisfrequenz $\omega_0$, überprüfen ob es sich um einen Schwingfall handelt und berechnen die Kreisfrequenz \[ \omega = \sqrt{\omega_0^2 - \delta^2} \] der gedämpften Schwingung. Anschließend müssen wir die Amplitude $A$ und den Phasenwinkel $\varphi$ in der analytischen Lösung \[ x(t) = A \textrm{e}^{-\delta t}\cos(\omega t + \varphi) \] aus den Anfangsbedingungen \[ \begin{align} x(0) &= x_0 \\ \dot x(0) &= v_0 \end{align} \] bestimmen. Dazu berechnen wir $x(0)$ und $\dot x(0)$ aus der analytischen Lösung und erhalten ein Gleichungssystem: \[ \begin{align} A \cos(\varphi) &= x_0 \\ -A \delta \cos(\varphi) - A \omega \sin(\varphi) &= v_0 \end{align} \] Nach einigen Umformungen erhält man daraus die Gleichungen: \[ \begin{align} \cos(\varphi) &= \frac{x_0}{A} \\ \sin(\varphi) &= - \frac{v_0 + \delta x_0}{A\omega} \end{align} \]

Es ist ein häufig gemachter Fehler, nun die beiden Gleichungen durcheinander zu teilen und anschließend den Arkustangens zur Berechnung des Winkels $\varphi$ zu benutzen, da der Arkustangens immer nur Winkel im Bereich $-\pi/2\dots+\pi/2$ ergibt, was zu einem falschen Ergebnis bei negative Anfangsauslenkung führt. In NumPy gibt es für genau dieses Problem die Funktion np.arctan2(y, x). Diese Funktion berechnet einen Winkel $\varphi$, sodass $\cos(\varphi)=x$ und $\sin(\varphi)=y$ ist. Der mathematische Hintergrund der Funktion arctan2 ist beispielsweise auf Wikipedia beschrieben.

Die Amplitude $A$ kann man bestimmen, indem man die beiden Gleichungen für $\cos(\varphi)$ und $\sin(\varphi)$ quadriert und addiert. Wegen $\sin^2(\varphi)+\cos^2(\varphi)=1$ ergibt sich daraus direkt: \[ A = \sqrt{x_0^2 + \left(\frac{v_0 + \delta x_0}{\omega}\right)^2} \] Das vollständige Programm, das diese Gleichungen umsetzt, ist unten abgedruckt.

Download    Schwingungen/Loesungen/gedaempfte_schwingung_analytisch.py

"""Audioausgabe einer gedämpften Schwingung. """

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile
import sounddevice

# Zeitdauer der Simulation [s].
T = 3.0

# Abtastrate für die Tonwiedergabe [1/s].
rate = 44100

# Federkonstante [N/m].
D = 7643.02

# Masse [kg].
m = 1e-3

# Reibungskoeffizient [kg/s].
b = 0.005

# Anfangsauslenkung [m].
x0 = 1e-3

# Anfangsgeschwindigkeit [m/s].
v0 = 0

# Berechne die Eigenkreisfrequenz und den Abklinkoeffizienten.
omega0 = np.sqrt(D / m)
delta = b / (2 * m)

# Brich das Programm ab, wenn es sich nicht um einen
# Schwingfall handelt.
if delta >= omega0:
    print('ES handelt sich nicht um einen Schwingfall.')
    exit(0)

# Lege die Zeitachse fest.
t = np.arange(0, T, 1 / rate)

# Berechne die Kreisfrequenz der Schwingung.
omega = np.sqrt(omega0**2 - delta**2)

# Berechne die Amplitude.
A = np.sqrt(x0**2 + (x0 + delta * v0)**2 / omega**2)

# Berechne den Phasenwinkel.
phi = np.arctan2(-(v0 + delta * x0) / omega, x0)

# Werte die analytische Lösung an den Zeitpunkten aus.
x = A * np.exp(-delta * t) * np.cos(omega * t + phi)

# Skaliere das Signal so, dass es in den Wertebereich von
# ganzen 16-bit Zahlen passt (-32768 ... +32767) und wandle
# es anschließend in 16-bit-Integers um.
dat = np.int16(x / np.max(np.abs(x)) * 32767)

# Gib das Signal als Audiodatei im wav-Format aus.
scipy.io.wavfile.write('output.wav', rate, dat)

# Gib das Signal als Sound aus.
sounddevice.play(dat, rate, blocking=True)

# Erzeuge eine Figure und eine Axes und plotte den
# Zeitverlauf der Auslenkung.
fig = plt.figure()
fig.set_tight_layout(True)
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('t [s]')
ax.set_ylabel('x [mm]')
ax.plot(t, 1e3 * x)

# Erzeuge eine Ausschnittsvergößerung.
axins = ax.inset_axes([0.55, 0.67, 0.4, 0.25])
axins.set_xlabel('t [s]')
axins.set_ylabel('x [mm]')
axins.set_xlim(0.5, 0.52)
axins.set_ylim(-0.4, 0.4)
axins.plot(t, 1e3 * x)

# Zeige den Plot an.
plt.show()